Influence of positive ions on the beamlet optics for negative-ion neutral beam injectors

Neutral beam injectors are based on the neutralization of ion beams accelerated at the desired energy. In the case of the ITER heating and diagnostic neutral beams, the target heating power translates into stringent requirements on the acceptable beamlet divergence and aiming to allow the beam to reach the fusion plasma. The beamlets composing the accelerated beam are experimentally found to feature a transverse velocity distribution exhibiting two Gaussian components: the well-focused one is referred to as the core component while the rest of the beam, the halo, describes beam particles with much worse optics. The codes that simulate beam extraction and acceleration usually assume that the negative ions move towards the plasma meniscus with a laminar flow (no transverse velocity) or that the transverse velocity distribution can be modelled as a Maxwellian and that the current density is uniformly illuminating the meniscus; under such approximations, the presence of highly divergent components cannot be explained. In this work, we develop a simple test-particle tracing code with Monte Carlo collisions, named ICARO (for Ions Coming Around), to study the transport of negative ions in the extraction region and derive the spatial and velocity distribution of the negative ions at the meniscus (i.e. the plasma boundary where a beamlet is extracted). In particular, the origin of the beamlet halo and its dependence on the source parameters are discussed, highlighting as a key parameter the energy distribution of positive ions in the source plasma.


Introduction
The performance of negative-ion sources for fusion relies on the surface conversion of precursors (atoms or positive ions) * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. into negative ions, at a converter surface of low work function, close to the extraction region [1]. In the case of the ITER heating neutral beam (HNB) and diagnostic neutral beam a weakly ionized hydrogen/deuterium plasma will be ignited by an inductively coupled discharge. The design of the ITER ion source and its prototype SPIDER [2] is based on the prototype source developed at the Institut für Plasma Physik (IPP, Garching) [3]. An expansion region with a length of about 25 cm separates the RF drivers from the extraction region. In the expansion region a transverse magnetic field B filter , named Profiles of the plasma potential along the SPIDER source as measured by movable electrostatic probes (PG is located at z = 0 cm). P RF = 50 kW/driver, B filter = 4.2 mT at the PG surface, PG potential set to 29.6 V, and plasma discharges in deuterium. Continuous lines give a polynomial fit to the measurements; the black dashed line represents the separation between the driver (on the right) and the expansion region (on the left). the filter field, is generated with the aim of preventing hot electrons from reaching the extraction region and thus effectively neutralizing the negative ions. A typical example of the plasma potential profile as measured in SPIDER by means of electrostatic probes is given in figure 1. It is worth highlighting that by increasing the source filling pressure, the plasma potential in the RF driver is strongly reduced and so is the potential drop between the driver region and the extraction region close to the plasma grid (PG).
Reducing T e from 10-20 eV in the driver to about 1-2 eV in the extraction region and maximizing the extent of the latter region are fundamental requirements for reducing the destruction rate of the negative ions and increasing their probability of being extracted. Negative ions can be effectively generated on the source walls thanks to the surface conversion of either H 0 atoms or positive ions (H + /H 2 + /H 3 + ). Thanks to proximity to the extraction apertures, the most effective generation occurs on the PG surface, which constitutes the boundary between the ion source and the electrostatic accelerator. Ions are extracted by applying a potential difference between the multi-aperture PG and a second electrode, called the extraction grid (EG). The latter is used as a dump surface for the co-extracted electrons, which are deflected by means of permanent magnets embedded in the grid itself. Such magnets are usually referred to as electron deflection magnets (EDMs) or co-extracted electron suppression magnets (CESMs).
Fluid and particle-in-cell (PIC) codes are the most suitable tools to investigate the plasma of the source [4,5]. In most cases, PIC codes are run with a scaled plasma density or with a scaled vacuum permittivity due to the very long computational time that they require even when running on supercomputers. Such codes were also applied to investigate the extraction region in more detail [6][7][8].
Beam extraction and acceleration are usually investigated by means of ray-tracing codes such as SLACCAD [9], OPERA3D [10] and IBSimu [11]. Such codes solve the Poisson equation along the accelerator, taking into account, besides the boundary conditions given by the electrode potential, also the space charge of the beam particles. The beam at a steady state is modelled starting from a laminar flow of particles moving towards the PG, and the space charge is zeroed upstream of the so-called plasma meniscus, i.e. the surface separating the plasma where the charge neutrality condition applies from the accelerator. In most cases, a transverse velocity distribution can be considered to take into account the energy spread of the H − coming from the plasma. As ray-tracing codes were originally developed for positive ion beams, a Maxwellian distribution is commonly assumed at the meniscus, as it well describes thermalized positive ions at the edge of the plasma sheath. However, this might not be the case for negative ions originating at surfaces at a short distance from the meniscus in comparison to their characteristic mean free path.
In experiments, the accelerated negative ion beam is found to be composed of two distributions: one with small divergence, called the 'core', and a broader component, called the 'halo'. Due to its high divergence, the halo fraction is usually poorly transmitted to the torus and causes heat loads on the beamline components and the grids of the accelerator [12]. The large transverse velocity of the halo fraction does not result from aberration of the electrostatic lenses at high-energy electrodes, but needs to be traced back directly to the initial distribution function of the ions (i.e. their position and velocity vector when they were first subject to the acceleration from the extraction field).
From the modelling point of view, an effective way to reproduce the experimental halo in ray-tracing simulation was to generate negative ions on the downstream side of the PG [12], an expedient also having a plausible physical explanation (i.e. concomitant influx of caesium and hot atoms diffusing from the source at those surfaces). In the code, an ensemble of particles carrying a current about 5%-15% of the total is thus emitted from that area, and it gains from the electric field a substantial transverse energy, which then turns into a higher divergence. On the other hand, PIC codes have shown the direct extraction of negative ions from the edge of PG apertures, with negative ions either crossing too short a distance in the plasma before the meniscus or being directly taken from the surface by the extraction field, leading to substantial aberrations of the beam optics [13].
Another approach to the numerical description of beamlet formation was attempted by tracing the negative ions from their point of origin at the caesiated surface, in a given plasma background and electromagnetic fields, up to their destruction or-optimistically-their extraction [14]. Despite lacking self-consistency, this approach is extremely useful for developing an insight into the physics of the negative-ion extraction and its dependence on the source parameters.
In this work we follow this test particle approach to discuss the features of the H − velocity distribution at the meniscus starting from the formation of the H − ions at the PG surface. To such purpose, H − ions are emitted from the PG with different initial energies and are traced in the source and accelerator magnetic fields with Monte Carlo collisions with the source gas and plasma. The spatial distribution j(x, y) of the extracted current density, the negative-ion transverse velocity distribution and the resulting emittance at the meniscus calculated from such a test particle model are discussed in terms of their dependence on the plasma parameters. In particular, it is highlighted that a substantial high-energy tail can be present in some experimental conditions. The paper is organized as follows: section 2 describes the model in terms of computational aspects and physics assumptions. Section 3 presents the typical features of the extracted particles, while in sections 4-6 the most relevant plasma source and computational parameters are scanned. In section 7 a comparison with experiments is provided. Finally, some conclusions are drawn in section 8 and the outlook of the work is discussed in section 9.

Description of the ICARO model
Negative ions are emitted from the chamfered portion of the PG with a given energy E 0 and with a cosine angular distribution. Their collisions with H 0 , H 2 , H + and electrons are modelled with a null collision method [15]. The collisional processes modelled here can lead to either neutralization of the H − (loss) or momentum transfer. A summary of the modelled collisions is given in table 1. Differently from [6,14], Coulomb collisions (CCs) with negative ions are also considered, assuming T H− = 1 eV with no drift velocity. For simplicity, the plasma is assumed to be formed by H + only. The plasma domain extends from the PG up to 150 mm upstream of the PG, a relatively large distance from the extraction region, at which the H − is considered to be lost. The  [20] Coulomb collisions with H + (CCs) plasma parameters are assumed to be constant throughout the domain. Such an assumption is somewhat justified by the almost flat profile of T e and n plasma in the last 150 mm of the expansion region, which is found experimentally in standard operating conditions [16]. For the gas parameters, the measurements by Briefi et al [17] on an RF-driven prototype source of identical design are assumed as a reference. In particular, for the H 2 we consider a cold population at T H2,cold = 630 K and a hot population with T H2,hot = 4000−5000 K. For the H 0 , a cold population with T H0,cold = 0.19 eV and a hot population with T H0,hot = 2.5 eV are considered. The relative abundance for the two atom population is set to 50:50; while for the molecules it is set to 67:33. Such numbers were obtained in the case of the BATMAN Upgrade ion source [18] for anRF power P RF = 70 kW/driver, a source filling pressure p fill = 0.3 Pa and a filter field B filter = 3 mT in front of the PG; in our model, to derive the gas density corresponding to the measured temperatures and abundances of H 0 and H 2 , gas flow conservation in SPIDER is considered. The derivation of n H0 and n H2 is described in appendix A. The calculated values of n H0 and n H2 are found to be much lower than the ones assumed by Gutser et al [14], which were derived from spectroscopic results available at the time. The reference value for the positive ion temperature is set to be T H+ = 0.8 eV as in other works [14,16]. Such an assumption is usually motivated by the fact that H + ions are generated mostly in the RF drivers by dissociative ionization of H 2 , direct ionization of H 0 or charge exchange reactions with H 0processes with a typical birth energy for the H + ions around 1 eV [16]. Nonetheless, in negative-ion sources based on RF drivers (like BATMAN, ELISE [21] or SPIDER) such positive ions can be accelerated towards the PG by the electric field associated with plasma potential profiles similar to the one reported in figure 1; collisions with the background gas should instead slow down such ions in the expansion region [16]. The Table 2. Plasma and gas parameters used for the reference case. The parameters utilized in the work by Gutser et al [14] are given for comparison.
Gutser et al [14] This work role of the assumption on the temperature of positive ions will be discussed in section 5. Differently from Gutser et al, we assume a plasma density of n H+ = 4 × 10 17 m −3 and we assume n e = 0.8n H+ and n H− = 0.2n H+ . The ratio n H− /n H+ is surely underestimated in the last 20-30 mm upstream of the PG, where in good caesiation conditions an ion-ion plasma is almost achieved [21], but it is considered a reasonable average value over the length of the domain. The electron temperature is assumed to be T e = 2 eV, implying that the rate of H − destruction via electron detachment is moderately high; it is about four times higher than at T e = 1 eV, and four times lower than at T e = 5 eV. A summary of the reference plasma and gas parameters utilized in this work is given in table 2.
The parameters utilized in [14] are also reported for comparison; it is worth stressing that, since the neutral density is much smaller in this work, the role of the gas in slowing the surfaceproduced H − is less relevant in this case. The rates for the reactions listed in table 1 are calculated, for the parameters listed in table 2, for different energies of the H − ions. The corresponding mean free paths are shown in figure 2. As one can see, the reaction with the shortest mean free path is the charge exchange reaction with the H 0 atoms. Momentum transfer collisions with H 2 (MT2) are much less likely and they decrease their effectiveness with increasing H − energy. The collisions that affect the H − momentum more are the charge exchange ones and the CCs. The effectiveness of CCs increases for low H − energy, with a typical v REL −3 = |v H--v H+ | −3 dependence [22]. As mentioned earlier, as a possible consequence of the electric field present at the exit of the driver, the presence of a drift velocity for the positive ions towards the PG is also considered in the model. As shown by Nishioka et al [6] in fact, CCs of H − with protons described by a drifting Maxwellian can be very effective in increasing the extraction probability for negative ions. Concerning the destruction of negative ions, the most relevant channels are ED and MN. It is worth highlighting that MN is sensitive to the H − energy, while ED is not. In addition, since MN in the  relevant energy range is equally (or slightly more) effective with H 2 + targets, and it is also present with H 3 + targets, the assumption on the ion fractions (100% H + ) does not play a significant role.
The most relevant physics and computational aspects to be dealt with in the present work are described in the following, with one subsection dedicated to each of these aspects.

The plasma domain
A sketch of the plasma domain in the XZ plane is given in figure 3. The aperture position and the geometry of the apertures are those of the SPIDER source. The PG of SPIDER features 1280 apertures, divided into 16 beamlet groups with 5 × 16 apertures each. Negative ions are emitted from the surfaces around the central aperture of the beamlet group, whose centre coordinate is set to X = Y = 0 mm. The domain length in the Z direction is 150 mm; as discussed previously, we applied uniform plasma parameters T e , n e over this length. Also the plasma potential V plasma is assumed uniform over such an extent, a condition that is verified experimentally in the case of relatively high electric biasing of the PG with respect to the walls of the plasma source. The role of electric fields within the plasma could also be included as an additional input parameter to the test particle model; this will be investigated in future work. In the absence of electric fields, only the magnetic field and the collisional processes influence the negativeion trajectories. In the vicinity of the walls, it is not necessary to include the electric fields in the Debye sheath because the sheath thickness is negligible; the sheath voltage is accounted for in the initial velocity of the ion, as discussed in section 2.2. The extent in the XY plane is limited only at the PG level, so that particles reaching the grid surface outside the beamlet group are stopped and counted as 'Side' losses. To mimic a source bigger than one beamlet group, the domain size is instead infinite in the XY plane from the PG surface up to 150 mm upstream.

Emission from the PG
The PG is a thick copper grid whose apertures are chamfered on both sides [23]. The location where the aperture radius is minimum is usually referred to as the PG knife. The pitch between the apertures is h x = 20 mm in the horizontal direction and h y = 22 mm in the vertical direction. The chamfer angle is 40 • on the upstream side of the grid and 45 • on the downstream side, and the radius of the PG aperture at the PG knife is 7 mm. As a consequence, two neighbouring apertures overlap in both the horizontal and the vertical. The area of the chamfered portion around one aperture in the bulk of the beamlet group is A chamfer = 389 mm 2 , while the aperture area is 154 mm 2 and the area of the rhomboidal flat surfaces on the upstream side of the PG is about 36 mm 2 . The study focused on the extraction probability of negative ions emitted from the chamfered portion, as it has the largest area and it surrounds the extraction aperture. The bias plate is also not modelled, as it is present only at the boundaries of the domain. Particles are emitted from given points uniformly spaced on the PG chamfer surface.
According to Seidl et al [24], it can be assumed that the H 0 impinging on the PG bounce back with a certain fraction of their initial energy E in (E out = R E E in , with R E = 0.7) and with a cosine distribution. If a Maxwellian distribution is assumed for the impinging H 0 , then a half-Maxwellian with T H− = R E T H0 can be used to approximate the energy distribution of the emitted particles. Following a test-particle Monte Carlo approach, particles at fixed emission energies will be studied first (sections 3-6) in order to investigate the energy dependence of the results. Such results for the monochromatic emission are then integrated a posteriori over the negative ion distribution at the PG (derived in appendix C), as discussed in section 7.
A certain positive potential difference V sheath − V PG is assumed between the PG and the plasma bulk, as is common practice in experiments. Such a potential difference increases the particle velocity in the direction perpendicular to the PG, thus modifying the distribution of the H − outflowing from the PG towards the plasma. In experiments V sheath − V PG is usually kept around 1-2 V to optimize the source performance. In the following, when not stated otherwise, V sheath − V PG = 1.4 V will be assumed.
In particular, with xyz being a frame of reference with x and y parallel to the PG at the emission point and z the normal to the PG surface at that point, the particle energy E 0 is first divided in the three directions imposing a cosine distribution. The velocity of the test particle in the z direction is then increased from v z0 to v z according to equation (1), in which e is the quantum of charge. Finally the velocity vector is moved to the XYZ frame of reference, in which Z is along the beamlet axis and XY is the plane perpendicular to the beam.

Particle motion
Negative ions move in the magnetic field given by the combination of the CESM magnets and the filter field. The time step is adapted at any iteration to keep the displacement at each step below ds = 0.5 mm. Such a choice allows one to properly describe collisions (the shorter mean free path is about two orders of magnitude larger) and the motion in the magnetic field. For comparison, the Larmor radius of an H − ion with 1 eV energy in a 17 mT magnetic field (the field due to the CESM at the PG knife) is about 8.2 mm. At each time step Monte Carlo collisions are implemented (see section 2.4) and particles might eventually be neutralized or exchange momentum with gas and plasma particles. Particles exiting from the upstream side of the domain are counted as losses. Particles aiming towards the PG can be lost on its surface only if they have a perpendicular kinetic energy larger than e(V sheath − V PG ), otherwise they are simply reflected back with an opposite velocity in the direction perpendicular to the PG. Particles are counted as extracted if they overcome the surface of the plasma meniscus; such a surface depends on the negative ion current density and the applied extraction voltage via the beamlet perveance, and therefore it is obtained from ray-tracing simulations (see section 2.5). If a particle reaches the PG outside the projection of the beamlet group it is also counted as a loss.

Collisions
Collisions are treated with a Monte Carlo method. At each time step, the probability for a collision to occur P coll is calculated from equation (2) where λ TOT = (1/λ CX + 1/λ MT2 + 1/λ LOSS ) −1 is the total mean free path, V is the H − ion velocity and dt is the time step. λ CX , λ MT2 and λ LOSS give the mean free paths for CX, MT2 and LOSS reactions. Without loss of generality, the ith mean free path is calculated from the collision frequency ν i as with v H − being the negative-ion velocity. A random number r 1 is generated and compared to P coll . If r 1 < P coll a collision occurs and to select which of the three events has to be considered another random number r 2 is generated and compared with the cumulated probability of the three events, each one calculated as Each of the three terms is obtained by summing over its subprocesses. For CX collisions, for example, collisions with the two H 0 populations are treated as collisions with different particles and the choice of which Maxwellian the H 0 has to be sampled from is based again on comparing random numbers with the ratio of mean free paths. The same is done for MT2 collisions with H 2 molecules. MT collisions with the gas (i.e. MT2 and CX) are treated with the small-angle approximation as in [15]. If a CX collision occurs, the scattering of the two particles is performed and then the H − is substituted by the scattered H 0 . CCs are treated separately, following the algorithm reported in [25]. (2)

Plasma meniscus and magnetic fields
The meniscus surface can be defined by the user and H − ions are counted as extracted if during their motion they overcome such a surface. A convenient choice is to import the solution of ray-tracing codes such as OPERA3D or IBSimu. In this work, the meniscus surface calculated by OPERA3D for a case of matching perveance has been adopted. The adopted extracted current density and extraction voltage are j ext = 250 A m −2 and U ex = 7.9 kV, respectively, corresponding to a perveance Π = 5.5 × 10 −8 A V −3/2 . The magnetic field of magnets embedded in the EG is also calculated by OPERA3D; the filter field is instead set to B = (2, 0, 0) mT anywhere. A filter field of about 2 mT is in fact found to be the optimum for operation in hydrogen at both ELISE [21] and SPIDER [26]. Since B x > 0 upstream of the PG, the line integral of the filter field along the motion of the extracted H − ions is such that the Lorentz force deflects the particles upwards.

Distribution of positive ions
Since the cross section for mutual neutralization of the H − and H 2 + /H 3 + is not known, it is assumed that all of the positive ions are H + . Their velocity is sampled from a Maxwellian distribution with T H+ = 0.8 eV. An average drift velocity v drift towards the PG is also considered. The value of this drift velocity is very important in determining the capability of CCs to exchange momentum with the H − and thus to cause the reversal of their motion and increase their extraction probability. Since ν CC ∼ |v H--v H+ | −3 = v REL −3 is the collision frequency for CCs, the dependence of the extraction probability on v drift is expected to exhibit a maximum. This expectation will be verified in section 5. In the following, when not stated otherwise, the drift velocity is set to v drift = 20 km s −1 , corresponding to an energy of 2 eV.

Virtual cathode
Due to the large negative-ion current density generated by a well caesiated surface a so-called virtual cathode might arise [27]. In practice, the space charge of the emitted negative ions could be so large as to induce a potential well preventing other negative ions from being extracted. The depth of the virtual cathode depends on several parameters, such as the negativeion emitted current density, the density of positive ions, and the potential difference V sheath − V PG [28]. The onset of a virtual cathode defines the limit of current density that can reach the sheath edge j em,s and enter the plasma. In our approach, the extraction probability p extr is calculated for particles starting at the sheath edge. The transmission of the negative ions through the sheath-i.e. the ratio p 0s = j em,s /j em,0 , which depends on the local plasma parameters-can be estimated with the model presented in [28]. The reader may note that, in the presence of a virtual cathode, H − with higher perpendicular energies can be preferred, with consequences also for the extraction probability. A more accurate description of the effect of a virtual cathode on the extraction probability is left for a future work. If the depth of the virtual cathode is low, however, we can assume that, to first order, this effect will be small. Therefore, the overall extraction probability can be calculated by multiplying the two factors. A brief discussion of the virtual cathode is given in appendix B.

Main features of the extracted particle distribution
The fraction of extracted particles and their spatial and velocity distributions are expected to be strongly influenced by the parameters of neutrals and plasma particles and by the boundary conditions. Also the Ions Coming Around (ICARO) was found to be very sensitive to the input parameters. While most such dependences will be discussed in the following sections, we highlight first some general features. To such purpose, we consider a fixed H − emission energy E 0 = 3 eV. For the parameters reported in table 2, we obtain a fraction of neutralized H − of 38.8%; 25.8% of the test particles impinge on the PG, 9.7% are lost in the plasma and 15.6% of the emitted particles are extracted.
The corresponding density map of the extracted particles in the XY plane is given in figure 4. By plotting the map of the extracted particles on a logarithmic scale it is possible to highlight the role of the filter field (see figure 4(a)): negative ions are transported preferentially towards apertures located above the emission aperture. Most of the extracted particles are extracted from the original aperture and with an asymmetric distribution (see figure 4(b)), resembling the one calculated in [14]. Such a characteristic is due to the magnetic field of the magnets embedded in the EG, which bends the negative ions on one side of the aperture towards the aperture itself while pushing those on the other side towards the bulk of the plasma, as exemplified in figure 3. Due to this effect a leftright asymmetry is found in the extracted particle density also on the beamlet group scale. Negative ions emitted from the aperture above or below the one we consider would show a reversed pattern in the horizontal direction, which changes the sign of B y from even to odd aperture rows (indicated as B CESM in figure 3). The presence of ions aiming directly towards the meniscus yields a higher current density in the external annulus of the aperture. This effect is magnified in the figure since the crossing angle of the particles with respect to the meniscus was not considered.

Energy dependence
The emission energy E 0 of the particles plays a key role in determining their extraction probability and their energy at the meniscus since it affects the reaction rates and the Larmor radius for gyration in the magnetic field. To provide an example, the probabilities of extraction, neutralization and losses (on the PG, outside the PG or in the plasma upstream of the simulated domain) is given in figure 5 for different emission energies E 0 .
It can be noted that, with increasing E 0 , the fractions of extracted and neutralized particles decrease while more particles are lost in the plasma and on the PG. The reason for this behaviour lies in the increase in the mean free path for collisions with the gas and the plasma. More and more H − thus return to a region where electrons have higher density and temperature and where they will be neutralized. The fraction of ions that reach the PG outside the beamlet group region (dubbed 'Side' in the figure) is of the order of 5%-10%. The velocity distributions at the meniscus generated from an ensemble of H − with E 0 = 5 eV and E 0 = 1 eV are compared in figure 6. In the former case two peaks around the emission energy can be clearly identified and correspond to the particles that reach the meniscus without any collisions with the background gas. Such peaks are non-monochromatic due to the distributed CCs. The velocity distribution of the extracted particles was separated into two components: those extracted from the aperture around which they were originally emitted from the chamfered surfaces (dubbed 'APERTURE' in the figure) and those that travel inside the plasma before being extracted from a different aperture than the starting one (dubbed 'REST' in the figure). This latter population is obtained as the summation of particles extracted from all apertures, and therefore they well represent the case of one aperture in the middle of the beamlet group, surrounded on all sides by a pair of apertures at least. As one can see, for large E 0 (figures 6(a) and (c)), the vast majority of the extracted particles exit through the emission aperture and with a velocity distribution featuring two large velocity peaks; in practice, negative ions have to be emitted towards the meniscus in order to contribute to the extracted current. For low H − energy instead, H − ions interact more with the gas and the plasma. As a consequence, they are less likely to be lost on the upstream side of the domain; such particles experience collisions with the background gas and tend to thermalize towards lower energy due to the interaction with H 0 , H 2 and H + , also inverting their velocity back toward the PG and its extraction apertures (as indicated by the 'REST' population in figures 6(b) and (d)). The fraction of H − that is extracted at the same aperture around which they were generated decreases from 63.2% to 30.9% when reducing E 0 from 5 eV to 1 eV ('APERTURE' population in figure 6(b)).
For every energy E 0 , the distribution of the particles being extracted from other apertures (REST) is more symmetric and generally narrower. The velocity distribution in the horizontal direction (V x ) is less symmetric than the one in the vertical direction (V y ) due to the vertical magnetic field generated by the CESM magnets.
By integrating the particle velocity distribution, the integrated extraction probability p extr can be calculated for every emission energy E 0 . An interesting feature is that, by increasing the H − energy E 0 , the extraction probability saturates to a minimum non-zero value, which is not negligible with respect to the extraction probability at low E 0 . Such a featuresomehow related to the direct interception of H − by the meniscus immediately after their extraction-is very important because it means that the fraction of high-energy H − ions at the meniscus is comparable to the fraction of high-energy H − ions emitted from the PG (i.e. negative ions with 1 eV have less than twice the probability of extraction of a negative ion emitted with 10 eV). The dependence of the calculated weighted extraction probability on the emission energy is compared to the one calculated by Gutser et al [14] in figure 7. The trend is very similar, while minor discrepancies can be attributed to the different grid geometry and to the slightly different gas and plasma parameters.
As mentioned in section 2.1, the electric field within the plasma is assumed to be zero. In reality, the extraction probability with respect to the emission energy might be influenced by the PG bias and by the presence of the bias plate, via the development of electric fields E z perpendicular to the PG in the plasma, favouring either slow or fast ions at emission depending on the sign of E z . It is true that even though the test-particle method cannot provide a selfconsistent electric field in the plasma, it can anyway consider given electric fields; nevertheless, in the experiment E z changes depending on bias voltages, and it can be minimized by a suitable choice of bias voltages (e.g. determining the variation of spatial potential below the characteristic H − energies) but also maximized or reversed. However, as discussed earlier, studying the influence of E z is not the purpose of this work.  [14] for different filter field strengths.

The role of positive ions
As we discussed earlier, positive ions are important both for the neutralization process and for momentum transfer via CCs. Since CCs have mean free paths that increase sharply with the H − energy, the positive ions have far more effect on the extraction probability of low-energy H − ; also the background of plasma ions exhibiting relatively high drift velocity can be less effective in this regard. To highlight the role of CCs, they were turned off in dedicated runs, which are compared to the cases with CCs. CCs help in extracting more negative ions (see figure 8(a)) and we can conclude that the most relevant contribution to the improvement of the velocity distribution of the extracted negative ions is given by CCs; an example for E 0 = 1 eV is given in figure 8(b). These results are in agreement with previous studies by Nishioka et al [6], who, however, used a particle-in-cell model to investigate those processes; the fact that test-particle results of ICARO are coherent with PIC results is especially important, as it is a confirmation that ICARO approximations hold and that ICARO is capable of grasping the main physical processes determining the ion extraction (despite the much lower computational cost in comparison to the PIC model).
Since CCs are found to be so relevant in affecting the extracted beamlet, it is important to discuss the features of the velocity distribution of positive ions in the expansion region. Their energy distribution depends on their temperature in the driver region, where most ionization events occur, on the potential profile along the source and on the density and velocity of the various species constituting a target for their collision along their path to the extraction region. The plasma potential profile in an RF-driven ion source was thoroughly characterized by McNeely et al [29] as a function of the source parameters, and it was shown that the potential difference between the RF driver and the expansion region, which provides lot of kinetic energy to  [16].
Experimentally, a positive ion flow with v drift ∼ 10 km s −1 was measured by Bandyopadhyay et al [30] in the expansion region of BATMAN for a filling pressure p fill = 0.7 Pa.
To have a practical characterization of the EDF of H + ions, a Maxwellian distribution with temperature T H+ is considered in the plane perpendicular to the PG, while a drifting Maxwellian with the same temperature and a drift velocity v drift is considered in the direction perpendicular to the grid. The EDF is thus described by three parameters: the density of positive ions n plasma , their temperature T H+ and the drift velocity v drift . For the sake of simplicity, the latter is assumed to be directed towards the PG, neglecting the velocity component due to plasma vertical drift, which is expected to have a lower magnitude than the one from the drivers towards the PG. For T H+ a range between half and double the reference value reported in table 2 is considered. Such a range ([0.4, 1.6] eV) is compatible with the one suggested by McNeely and Schiesko [31].
In this section, we discuss the effect of such parameters on the extraction probability and on the features of the extracted ions. The dependence of p extr on v drift for n plasma = 4 × 10 17 m −3 and different temperatures T H+ is shown in figure 9.
As expected, v drift > 0 improves the extraction probability and it is found that any v drift below 30 km s −1 has a beneficial effect on p extr . If v drift is increased too much, however, the extraction probability is reduced again; there is thus an optimum v drift , which is found to increase as T H+ increases. For example, for T H+ = 0.8 eV the optimum is located around v drift = 10 km s −1 while for T H+ = 1.6 eV the optimum is found around 20 km s −1 .
It is worth highlighting that the lower T H+ , the higher the extraction probability p extr that can be obtained. Such behaviour can be explained in terms of the |v H--v H+ | −3 dependence of the CC frequency ν CC . The lower T H+ in fact, the higher the collision frequency. When T H+ is large, collisions occur seldom and have little effect on the reversal of the H − motion. In addition, since such collisions are less frequent, they can affect the motion of the H − only if the momentum of the H + towards the PG is large. This should explain why the optimum is found for larger v drift .
To scan the plasma density, a simplistic assumption is made: we assume the same gas density, temperature and degree of dissociation while n plasma is simply increased with a fixed n e /n plasma . The dependence of p extr on E 0 for different n plasma is given in figure 9(b). As expected, the extraction probability increases mostly for low E 0 . Particle velocity distributions at the meniscus for various v drift and n plasma are compared in figure 10. The increased extraction probability is found to be due to particles with low transverse velocity, which has a beneficial effect on the optics of the extracted particles. In practice, a denser plasma can 'select' particles with lower transverse velocity for extraction.

Pressure dependence
A typical feature that is observed experimentally is the dependence of the measured beamlet divergence on the source filling pressure even though perveance is kept fixed [32]. The filling pressure affects the plasma discharge in many ways, both in the drivers and in the expansion region. We will use ICARO to discuss the differences between two pressure levels, trying to include the change of parameters in the boundary conditions and study the influence on the extraction probability. In experiments, p fill acts on several parameters, such as the gas density, the plasma density, the electron temperature and the energy of positive ions, due to the reduced potential difference V driver − V expansion and increased collision frequency for collisions in the expansion chamber at large pressure. A possible effect one might think of is that the higher the filling pressure, the more H − ions will slow down by colliding with the gas before being extracted. The cases selected for our investigation, listed in table 3, are based on known trends of the plasma parameters as a function of the filling pressure. In particular, as a first step, we tried different combinations of parameters while keeping a constant n plasma . Namely, the densities of H 0 and H 2 are increased while keeping constant their temperatures and the degree of dissociation. The measurements by Fantz et al [33] suggest that this assumption, at least concerning the temperatures and relative abundances of neutrals, should be reasonable. The standard case simulated up to now (Case 0) was compared with a case with doubled gas density and reduced T e (Case 1), in which T e was set to T e = 1 eV to reproduce the typical measurements by the electrostatic probes for p fill = 0.6 Pa [34]. The variation of the extraction probability by such a combination is found to be quite limited. This outcome confirms that CX collisions are much less effective than CCs in affecting the extraction probability and that the particle velocity distribution at the meniscus could be affected by pressure mostly via a variation of the velocity distribution of positive ions.
The temperature of positive ions was thus varied around the value T H+ = 0.8 eV utilized in the previous section. A case with p fill = 0.3 Pa was also modelled assuming T H+ = 1.6 eV and keeping v drift = 20 km s −1 (Case 2). In Case 3, the drift velocity was reduced to 10 km s −1 , consistently with the measurements reported in [30] at p fill = 0.7 Pa, and T H+ was set to T H+ = 0.4 eV. As said before, the tested range of T H+ is compatible with the range [0.5, 2.2] eV proposed by McNeely and Schiesko [31] from measurements at BATMAN.
The dependence of p extr on E 0 in the four cases is given in figure 11(a). As expected, the extraction probability increases in Case 1 and Case 3, but this happens mostly for low E 0 . The particle velocity distributions at the meniscus for the four cases are compared in figure 11(b). Again, higher collision rates help in increasing the fraction of extracted particles with low transverse velocity. It is worth noticing that changing only the electron temperature and the gas density has a small effect on the velocity distribution of the extracted negative ions (Case 0 vs Case 1), while a much larger effect is obtained when considering the influence that the filling pressure might have on the parameters that are used to describe the EDF of the positive ions (Case 2 vs Case 3).

Application to experiments
In this section we apply the results of previous sections in the context of experiments. In particular, in section 7.1 we derive the extraction probability and the velocity distribution at the meniscus for the negative ions that are generated from a Maxwellian distribution of H 0 atoms impinging on the PG with T H0,hot = 2.5 eV, taking as reference the cases discussed in section 6. In section 7.2 we discuss the contribution by positive ions and in section 7.3 we compare the outcomes from ICARO to experimental data.

Surface conversion of neutrals
The energy distribution of the negative ions emitted from the PG surface after the conversion of the H 0 population with 2.5 eV is determined as described in appendix C. The average extraction probability is obtained as the weighted average of the calculated probabilities at various E 0 (e.g. figure 11(a)). The derivation of the weights is also described in appendix C. By such averaging we obtain p extr,avg = 16.0% for Case 0 and 14.4% for Case 2 of the previous section.
The resulting velocity distributions for Case 2 and Case 3 are given in figure 12. As discussed earlier, the velocity distribution in the horizontal direction (see figures 12(a) and (b)) tends to be broad and asymmetric. For this reason we discuss the width of the velocity distribution for the vertical component of the velocity V y (see figures 12(c) and (d)). The plots of V y report a Gaussian fit applied to data above 0.4, in order to highlight the non-Gaussian tails at high transverse velocities. The difference between the two cases is quite significant: transforming the 1/e width into a temperature, one gets T H− = 2.4 eV for Case 2, while T H− is reduced to 0.9 eV for Case 3. For Case 2 the fraction of particles lying outside the fitted Gaussian distribution is 7.9% while it is increased to 15.8% for Case 3. In practice, high pressure (Case 3) yields a narrower distribution, in which tails, generated by H − directly extracted from the PG surface, and therefore characterized by a large transverse velocity, are clearly distinguishable. At lower pressure (Case 2), the whole distribution is much broader, so that the tails due to direct extraction are harder to distinguish. It is interesting to discuss these findings in terms of effects on the beamlet optics. The initial transverse velocity, being linearly independent-at first-order approximation-from the transverse velocity developed by the electric field in the accelerator, can be combined with the latter as the root sum of squares. In practical terms, a 1/e divergence of 7 mrad for a beam energy of 60 keV obtained in the ideal case of cold H − would become, for the low-pressure case (Case 2), a divergence of 10 mrad.
This approach offers a value to compare with ray-tracing codes like IBsimu, in which a Maxwellian distribution can be assumed for the transverse velocity of the extracted negative ions.
The value obtained for Case 3 is quite close to the typical assumption by Den Harder et al (T t = 1 eV) when simulating the negative ion extraction and acceleration from BATMAN Upgrade or ELISE [32]. The very different value obtained for T H− in the cases of small and large p fill is expected to strongly affect the optics, thus possibly explaining the experimental evidence of beams with larger divergence when the source is operated at low pressure.
Exploiting symmetries, the particles emitted from the central aperture of the beamlet group and extracted from other apertures can be combined to describe the particles that come from other apertures and enter the accelerator through the central aperture. By doing so, it is possible to plot a current density map at the meniscus and a beamlet emittance (see figure 13). The current density is found to be peaked in the side annulus of the aperture and unbalanced in the horizontal direction (the one perpendicular to the local magnetic field due to EDM). Such results are in agreement with those of backward-tracing models applied to emittance measurement performed at QST (Japan) [35,36]. As it can be noticed, the high transverse velocity tails highlighted in figure 12 mainly correspond to convergent particles, populating the odd quadrants of the emittance plot (see figure 13(b)). These particles, which we can think to be responsible for the beamlet halo, are not accounted for in typical gun codes used in the design of electrostatic accelerators (neither in terms of current density distribution nor initial velocity), while they are usually obtained from PIC codes [8,37]. Such particles are of high relevance in the case of multigrid accelerator, as it is the case of the ITER HNB, as they can be responsible for huge heat loads on the last grids of the accelerator and on the beamline components.

Surface conversion of positive ions
In order to assess the contribution of the conversion of positive ions to the extracted current density a measurement of the energy distribution of the positive ions impinging on the PG is needed. At present, we limit ourselves to some general considerations based on the outcomes from ICARO presented in this work.
The H + ions that are accelerated from the driver towards the PG usually have only a few eV of energy when they reach the PG since they have a high rate for slowing-down collisions [16]. The molecular positive ions (H 2 + and H 3 + ), instead, are expected to reach the PG with much larger energy as they come directly from the driver region or from the first part of the expansion region and have a much lower collisionality [38]. Such ions could thus lead to the emission of H − with a large birth energy E 0 . Measurements with a retarding field energy analyser in SPIDER confirm the presence of a substantial flux of high-energy positive ions for increasing RF power [39]. More details about the yield of negative ions from positive ions are given in appendix C.
To provide some reference numbers, in experimental campaigns at SPIDER, the flux of positive ions towards the PG was estimated to be about Γ + = 2.7 × 10 21 m −2 s −1 with P RF = 50 kW/driver and p fill = 0.35 Pa. For P RF = 70 kW/driver and p fill = 0.3 Pa we can assume the flux to scale linearly with power and we thus get Γ + = 3.8 × 10 21 m −2 s −1 . The fraction of such flux due to molecular ions can be expected to be about Γ H2+,H3+ = 0.4Γ + . Such a number is assumed considering that the typical fraction of molecular ions accelerated from an RF-driven source of positive ions is about 30% [40], but such a source can operate at much larger pressure and thus have a larger H + fraction thanks to the higher rates for dissociative recombination of molecular ions at low electron temperature. If the average yield of negative ions for the molecular ions is assumed to be Y avg,+ = 0.35 (corresponding to an energy of about 10 eV), an emitted H − current density of about j H2+,H3+ H− ∼ 60 A m −2 can be estimated. Molecular positive ions can thus generate a relevant number of negative ions with an energy of several eV. From the calculations presented in this work, it is expected that such ions can be extracted with a non-negligible probability with respect to the low-energy ones and that this extraction occurs mostly through the original aperture with a convergent velocity distribution. Due to their large emission energy, such H − are unlikely to be stopped by the virtual cathode and will be present in the beam whatever the virtual cathode depth. An extraction probability of 7%-10 % might be applicable to these ions, which is about the saturation value we found for p extr at high energy: we find that they contribute about 11-16 A m −2 to the extracted current, which is about 3%-5% of the target current density for the ITER HNB.
Concerning the H + , if a drift velocity of 20 km s −1 (as in Case 0) and V sheath − V PG = 2 V are assumed, the energy of the H + impinging on the PG is about 4 eV, resulting in a conversion yield about 0.16 (see appendix C). Assuming Γ H+ = 0.6Γ + , the negative-ion current density emitted by conversion of H + is estimated to be about j H+ H− ∼ 60 A m −2 . The average energy of such H − is expected to be R E × 4 eV = 2.8 eV. The corresponding p extr from ICARO in Case 0 is about 15% (see figure 11(a)), yielding an extracted current density of about 20 A m −2 .

Comparison with experiments
The most suitable measurements to compare the results from ICARO with are those of the IPP facilities BUG or ELISE, which feature the same PG geometry as SPIDER and the same PG-EG gap. The comparison should be done for the operational parameters for which the gas density and temperature reported in table 2 were measured, that is for P RF = 70 kW/driver and p fill = 0.3 Pa. The maximum current density that was measured at ELISE with such parameters was about j ELISE ∼ 300 A m −2 [41]. Such high current density was obtained with a large extraction voltage (U ex > 9.5 kV) so that the beam was in under-perveance conditions. While higher current densities will be achieved with higher available RF power, it may be argued that, at the present stage, ELISE could have reached an extracted current density about j ext ∼ 250 A m −2 in perveance match conditions for source parameters around those for which we used the gas parameters provided in [17]. With the parameters of table 2 and assuming a work function W f = 1.5 eV, the formulae of appendix C can be used to estimate a total emitted current density j em,0 = 610 A m −2 . In figure 14 we show the dependence of the average extraction probability p extr,avg and of the transmission probability through the sheath p 0,s (see appendix B) on the sheath voltage for Case 0. Note that p 0,s /5 is plotted instead of p 0,s to increase the readability of the figure.
The product P of these two variables is maximum for V sheath − V PG ∼ 1.5-2 V and its value is about 9%, which yields an extracted current density of about j ext,ICARO = P × j em,0 × A chamfer /A aperture ∼ 140 A m −2 . As mentioned in section 2.7, when p 0,s is very far from 1, the virtual cathode depth is larger and p extr,avg is overestimated in our calculations. The effect of this overestimation should anyhow be quite small in the range where P is maximum, p 0,s being around 0.7. If we also consider that about 25-35 A m −2 of the total extracted current originates in the plasma volume (this being the current density achievable in operation without caesium) and we sum the contribution to the extracted current density from the conversion of positive ions, which we estimated to be about 35 A m −2 in section 7.2, we obtain j ext,tot = 200-210 A m −2 , which is 80%-84% of the experimental extracted current density. The agreement is quite good and the gap between experiments and these calculations could be closed by slightly tuning the plasma parameters.
The most relevant qualitative modification of ICARO to be implemented in a future work is the implementation of profiles of n H+ , n H− and n e along the plasma source. The presence of an ion-ion plasma close to the PG is in fact expected to significantly increase the extraction probability due to the reduced electron density.

Summary
In this work, a particle tracing code with Monte Carlo collisions, named ICARO, was implemented similarly to what was done by Gutser et al [14] in order to describe the transport of the negative ions produced on the upstream side of the PG in the extraction region up to their eventual extraction. As demonstrated by the case studies reported in this paper, despite its relative simplicity the test-particle approach is capable of modelling the key processes influencing the extraction of negative ions, yielding results similar to fully consistent simulation codes, which are normally extremely time-consuming when three-dimensional geometry must be treated. The effect of most source parameters can be easily studied with this model, yielding important indications for the operation of the sources.
In a typical source of negative ions for neutral beam injectors (NBIs), most of the negative ions are generated by surface conversion of thermal atoms. In the case of RF sources, a hot population with T H0 = 2.5 eV is identified in a wide range of source parameters [33]. For this population, we can assume an average yield of negative ions of about 0.2-0.3 [24] and an energy reflection coefficient of about R E = 0.7. As a result, negative ions are emitted with an energy distribution roughly corresponding to T H− = 1.75 eV. The average extraction probability calculated by ICARO, around 10% as discussed in section 7, was found to be consistent with the extracted current density from experiments starting from the plasma and gas parameters.
In this work, we have studied the conditions for which the high-energy H − ions have a non-negligible probability of being extracted. Such particles usually have a long mean free path for collisions with the gas and the plasma particles, so that they are unlikely to be slowed down or pushed towards the PG by collisions. Due to their large energy, also their Larmor radius is large and the magnetic field is thus quite ineffective in ensuring the reversal of their motion. Such particles can therefore be extracted only if they are emitted with a velocity towards the meniscus. They are usually extracted with a convergent distribution, possibly leading to the formation of a socalled beam halo.
The low-energy H − ions instead have shorter mean free path and can be more easily extracted also from apertures other than the emission one. Clearly, precursors (either atoms or positive ions) impacting the converter surface with a low energy would be preferred for this reason, but the conversion yield is low at low impact energy. As discussed in section 7.2, for the same reason a reduced sheath voltage would help to increase the extraction probability, but a local optimum around 1.5-2 V is found when the formation of a virtual cathode is also considered.
CCs with positive ions were found to be the most relevant in affecting the extraction probability and the velocity distribution of the H − at the meniscus. This feature is very peculiar since the collision frequency of CCs has a sharp dependence on the relative velocity between the test negative ion and the positive ions. As a consequence, the extraction of low-energy H − can be favoured. Increasing the density of positive ions could thus be a valuable option to increase the fraction of low-energy H − within the extracted beamlets, thus yielding a narrower transverse velocity distribution for the H − at the meniscus (i.e. better beamlet optics) and also increasing the allowable negative-ion current. A higher density of positive ions also prevents the onset of the virtual cathode before the converter surface. The higher efficacy is anyhow obtained acting on the energy distribution of the positive ions themselves. In particular, a relatively slow motion towards the PG (v drift ∼ 5-10 km s −1 ) with a low T H+ was found to be the best combination to limit the most the fraction of tails within the velocity distribution of negative ions at the meniscus.
From the point of view of improving the beamlet optics by means of CCs it appears thus more effective to reduce the energy of positive ions than to simply increase their density. This is found by maintaining similar parameters for the H − emission from the converter surface (which might not be the case when changing the plasma density). Such a statement applies also to molecular ions (H 2 + , H 3 + ) from the drivers and the expansion region, whose energy can reach tens of eV. Such particles in fact, besides being of little help in reversing the motion of the H − emitted from the PG, also tend to generate energetic negative ions, which sustain the highly divergent component of the beam. In experiments, the energy distribution of positive ions is mostly influenced by the potential profile along the source and by the thickness of the gas target that opposes the positive ions in their motion from the plasma driver to the PG. This explains why the beamlet divergence is found to be sensitive to the source filling pressure, which acts on both these aspects. As shown in figure 12, according to ICARO, the fraction of cold ions at the extraction is much improved in the high-pressure case, and so also is the quality of the beam. It must be noted that this result assumed identical plasma density; in RF sources of negative ions, when increasing the pressure, the diffusion of positive ions from RF drivers to the PG through the transverse filter field is reduced, but this is compensated by the increase in plasma density inside the drivers.
As the reader would note, the potential profile is also influenced by the PG bias: a high bias would flatten the potential profile, but remove the contribution of the electric field to the plasma diffusion towards the extraction region and therefore the density of positive ions at the PG; so a compromise is necessary. Moreover, the potential profile cannot be completely flattened by increasing the PG bias [39].
The bulk of the velocity distribution of extracted negative ions resembles a Mawellian distribution with a transverse temperature about T t = 1-2 eV, but a fraction of the order of 10%-20% of the extracted particles have a higher transverse velocity, which is not taken into account in standard ray-tracing codes. Such particles are mostly located in the convergent quadrant of the emittance plot as discussed in section 7. The results of our study, summarized here, contribute to the interpretation of the overall picture that is already found experimentally, when comparing the design and operation of large RF sources to large filament-arc sources. The discussed arguments could explain why the beamlet divergence achieved from RF sources is larger than that from arc sources, where the temperature of negative ions in the extraction region was measured to be about T H− = 0.15 eV [42]. After all, arc sources typically have a substantially higher degree of dissociation [43] (i.e. a larger H 0 target thickness in the expansion region), a larger plasma density [43] and a lower temperature of ositive ions (T H+ = 0.25-0.5 eV) [44,45] than RF sources, likely due to the much smaller plasma potential difference between the driver and the extraction region.
In this work, we highlighted that the influence that positive ions have on the beam optics is given by several competing effects. The most relevant are reported in table 4.

Outlook
The first promising results obtained with ICARO suggest that such a tool can be helpful to understand the transport of negative ions in the extraction region and to provide input to experiments. We have shown that low-energy positive ions have to be favoured since they can effectively help the extraction of negative ions with low transverse velocities, ensuring a lower beamlet divergence and a smaller fraction of highly divergent particles (beamlet halo). Since the energy and the flux of positive ions impinging on the PG depend significantly on the potential difference V driver − V exp between the driver and the expansion region, and thus on the electron temperature in the RF drivers and on the applied filter field opposing the plasma diffusion towards the PG, operation at high pressure and low power usually yields a better beamlet optics. In the case of the multi-grid electrostatic accelerator of the ITER NBIs (and related facilities), however, p fill is fixed at 0.3 Pa in order to limit the stripping losses, and since a current density of the order of 330-350 A m −2 has to be extracted at perveance match, RF power larger than 70 kW/driver will have to be used [41]. The effect of using higher RF power on V driver − V exp will be to increase such a potential difference. This could lead to a further increased halo fraction.
A possible approach to the mitigation of tails when operating at low pressure could thus be to improve the plasma confinement in the drivers (to reduce the plasma potential in the driver and thus the acceleration of the positive ions) or to have a longer expansion chamber (so as to increase the thickness of the gas target at constant pressure, even though this might also influence the plasma diffusion, yielding lower plasma densities at the PG). In order to improve the plasma confinement in the RF drivers, experimental campaigns with cusp magnets around the RF drivers will be performed in SPIDER.
Concerning the development of ICARO, future steps will be focused on the implementation of the virtual cathode, on the possibility of having profiles for the plasma parameters along the expansion chamber and on the coupling of ICARO with some ray-tracing codes, so as to have more self-consistent particle tracing. and atoms as in equations (3) and (4) with T H0,1 = 0.19 eV, T H0,2 = 2.5 eV, T H2,1 ∼ 630 K and T H2,1 ∼ 4500 K. The degree of dissociation was estimated to be n H0 /n H2 = 0.3. n H2 = n H2 ( n H2,1 (T H2,1 ) n H2 + n H2,2 (T H2,2 ) n H2 n H0 = n H0 ( n H0,1 (T H0,1 ) n H0 + n H0,2 (T H0,2 ) n H0 Following [17,33], the fractions A 1 and A 2 were set to 0.67 and 0.33, respectively. The atomic fractions B 1 and B 2 were both set to 0.5. In the case of SPIDER, the gas conductance from the beam source to the vessel was estimated to be C BS = 30 m 3 s −1 . The corresponding probability for a particle crossing the PG to be transmitted through the accelerator is found to be F = 0.20. For p fill = 0.3 Pa and T H2 = 293 K a density n H2,RT = 7.4 × 10 19 m −3 is estimated. Setting a 1 = n H2,1 /n H2,RT , a 2 = n H2,2 /n H2,RT , b 1 = n H0,1 /n H2,RT and b 2 = n H0,2 /n H2,RT , we can write equation (5), where v ij are the average thermal speeds of the various species. Imposing the degree of dissociation to be 0.3 turns yields equation (6). Solving these two equations for the selected values of A 1 , A 2 , B 1 and B 2 yields the densities of the various neutrals in the case of p fill = 0.3 Pa. Such values are reported in section 2. Larger pressures were simulated by simply scaling linearly the neutral density with pressure.

Appendix B. The sheath voltage
In this appendix we follow the approach by McAdams et al [28] to derive the maximum negative-ion current density that can reach the plasma sheath in two conditions applicable to our investigation. The limit of the negative-ion current density that can reach the plasma sheath j em,s is roughly given by the onset of a virtual cathode, with only a mild increase in j em,s when increasing further the current density emitted at the converter surface j em,0 . Figure 15 shows, for varying sheath voltage V sheath − V PG , the emitted current density at which a virtual cathode appears, for relevant plasma properties (named Case 2 and Case 3 in table 3, section 6). If we consider an emission from the wall j em,0 = 610 A m −2 , at a sheath voltage V s − V PG = 1.4 V, we obtain p 0s = 0.87 and p 0s = 0.58 for Case 3 and Case 2, respectively; already at V sheath − V PG = 2 V, almost full transmission would occur in Case 3, while one quarter of the generated H − are reflected back to the emitter surface in Case 2 (p 0s = 0.74). According to recent measurements in SPIDER [46], the sheath voltage V sheath − V PG was measured to be about 2 V in a caesium condition providing intermediate beam current densities (up to 160 A m −2 ), and quite independent of the caesium injection rate (which is commonly considered to help increase the surface conversion yield in the presence of low-pressure background oxidants and impurities). To ensure p 0s = 1 for Case 2 a potential difference V sheath − V PG = 4 V is required. Note that in such a case also positive ions, after being accelerated in the sheath, would give a relevant contribution to the generation of negative ions; according to [24] the yield of negative ions for an H + ion with 4 eV could be about Y = 0.16 (see appendix C). Such a contribution is not considered in this work.

Appendix C. Weights for the energy distribution
The yield for the conversion of a hydrogen atom into a negative ion at the PG surface can be calculated following Seidl et al [24]. Assuming a perfectly caesiated surface, the yield of negative ions for hydrogen can be calculated from equation (7) with R N η 0 = 0.42, E th /R E = 1.05 and E in the energy of the incident particle: For positive ions, the same formula can be used, but with R N η 0 = 0.30 and E th /R E = 2.0. For molecular ions, it is assumed that they dissociate before being converted into H − , so that an H 2 + with 10 eV is considered like 2 H + ions with E in = 5 eV. The calculated yield of negative ions for atoms and for the three positive ion species is given in figure 16.
Assuming a Maxwellian distribution with T H0 for the impinging atoms, the distribution of the emitted negative ions can be calculated by weighting over the yield. To investigate this aspect, a sample of 200 000 particles with a Maxwellian distribution of temperature T H0 was considered. The particle energy was reduced to E out = R E E in to consider the impact on the wall, and they were re-emitted with a cosine distribution. The yield for each particle was calculated from equation (7) and an acceptance-rejection method was applied by comparing random numbers uniformly distributed between 0 and 1 with the calculated yield for each particle. In such a way, it was possible to weight the distribution over the yield. The energy of the particles that generated a negative ion was reduced to E = E out − E th . In section 7, the monochromatic solutions of sections 3-6 were weighted over the distribution calculated as described in this appendix for T H0 = 2.5 eV using bins with a width of 1 eV.