From meniscus formation to accelerated H− beam: coupling of 3D-PIC and ion-optics simulations

The ITER NBI is based on negative hydrogen ions extracted from caesiated ion sources. The 3D particle-in-cell Monte Carlo collision code Orsay negative ion extraction (ONIX) models the beamlet formation of negative ions in such sources where surface production plays an important role. A coupling scheme between ONIX and the ion-optics code ion beam simulator (IBSimu) has been developed and compared to other particle simulation approaches. This extends the computational domain such that the complete grid system can be included while only marginally increasing the computational cost. The properties of the accelerated ONIX beamlet are studied and compared to standalone IBSimu calculations, which are based on a simplified plasma model. The comparison provides insight about the effect of approximations made in ion-optics codes, which were also used to design the ITER NBI grid systems. ONIX volume and surface produced negative ions have a different angular distribution in the accelerated beamlet. The ONIX volume produced particles have a similar core divergence compared to standalone IBSimu calculations, but there is more halo in the IBSimu angular distribution. In the ONIX simulations, a Debye sheath is formed between the plasma and the grid, which repels negatively charged particles. The sheath decreases the extracted current density at the edge of the aperture for volume produced ions. Contrarily, surface produced particles are directly extracted near the edge of the aperture. Particles extracted near the edge of the aperture are highly divergent at the end of the grid system, independent of their initial angle. To summarize, the presence of the plasma sheath around the apertures in the plasma grid as calculated by ONIX decreases the halo from volume produced particles compared to standalone IBSimu.

The ITER NBI is based on negative hydrogen ions extracted from caesiated ion sources. The 3D particle-in-cell Monte Carlo collision code Orsay negative ion extraction (ONIX) models the beamlet formation of negative ions in such sources where surface production plays an important role. A coupling scheme between ONIX and the ion-optics code ion beam simulator (IBSimu) has been developed and compared to other particle simulation approaches. This extends the computational domain such that the complete grid system can be included while only marginally increasing the computational cost. The properties of the accelerated ONIX beamlet are studied and compared to standalone IBSimu calculations, which are based on a simplified plasma model. The comparison provides insight about the effect of approximations made in ion-optics codes, which were also used to design the ITER NBI grid systems. ONIX volume and surface produced negative ions have a different angular distribution in the accelerated beamlet. The ONIX volume produced particles have a similar core divergence compared to standalone IBSimu calculations, but there is more halo in the IBSimu angular distribution. In the ONIX simulations, a Debye sheath is formed between the plasma and the grid, which repels negatively charged particles. The sheath decreases the extracted current density at the edge of the aperture for volume produced ions. Contrarily, surface produced particles are directly extracted near the edge of the aperture. Particles extracted near the edge of the aperture are highly divergent at the end of the grid system, independent of their initial angle. To summarize, the presence of the plasma sheath around the apertures in the plasma grid as calculated by ONIX decreases the halo from volume produced particles compared to standalone IBSimu.
Keywords: ONIX, IBSimu, neutral beam injection, ion sources, ion optics, 3D particle-in-cell, surface production (Some figures may appear in colour only in the online journal) * 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.

Introduction
The ITER neutral beam injection system is based on H − /D − ions [1]. The negative ions are generated in an RF-driven plasma source, mostly on the plasma facing surfaces by conversion of atoms and positive ions [2]. The ions are extracted from the quasi-neutral plasma through apertures in the plasma grid (PG) by the extraction potential of approximately 10 kV, applied between the PG and extraction grid (EG).
After extraction, the ions are accelerated to the full energy in five acceleration steps, which generates a low divergence beamlet. The ITER beamlines are designed assuming 85% of the beam power at a 3-7 mrad core divergence, and 15% at a 15-30 mrad halo divergence [1,3]. The e-folding divergence convention is used throughout this paper.
The test facility extraction from a large ion source experiment (ELISE) has a half-size ITER ion source, and aims to demonstrate the ITER-required current density and coextracted electron to ion ratio. The ELISE grid system has one acceleration stage, between the EG and the grounded grid (GG). The total beam energy is limited to 60 keV, which puts strict limits on the perpendicular energy distribution of the accelerated beamlet to satisfy the 1 MeV ITER divergence requirements [4]. The ELISE experiments are supported by modelling, in which beam formation is studied using several methods.
Particle in cell (PIC) codes take a physical approach where the system under study, typically a part of the plasma and the extraction gap, is evolved over time until steady state is reached. Complicated phenomena such as the evolution of the 3D plasma potential, and sheath formation in the presence of a magnetic field and surface production of negative ions, follow self-consistently from the model. The drawback of this modelling approach is the high computational cost, due to the requirement to resolve the plasma frequency and Debye length, which also puts a high requirement on the number of particles to simulate. Due to the small computational domain, beamlet properties after the full grid system including acceleration stages cannot be directly calculated.
Gun-type ion-optics codes only track the negative ions towards the extraction aperture in the source plasma. The compensating positive ions are approximated by an analytical charge density, and the collisional behaviour of the particles in the plasma is mimicked by locally suppressing the magnetic field. Although this leaves out many physical processes, such as for example plasma generation and sheath formation, it is a common approximation which was successfully used in the design of many negative ion extraction systems [5][6][7][8].
In this paper, the PIC computational domain is extended by coupling with an ion-optics code, self-consistently including the effects of space charge on the potential. In this way, the accelerated beamlet properties can be calculated for a PIC beamlet, and directly compared to gun-type simulations to assess the impact of the approximations made in those simulations. The potential map from the PIC extraction domain and ion-optics domain are combined to track diagnostic particles a posteriori. Since full trajectories are followed, this allows linking upstream to downstream particle properties so that e.g. the impact of the initial velocity distribution on the divergence can be studied.

Computational methods
The particle-in-cell Monte Carlo collision (PIC-MCC) [9] code Orsay negative ion extraction (ONIX) [10] was developed to simulate negative ion extraction in sources which are dominated by surface production of negative ions, which is the case for caesiated H − /D − sources [2,11,12]. The ONIX code self-consistently simulates particle transport in a hydrogen or deuterium plasma, including dissociation, the effect of collisions, and the influence of an external magnetic field. The collisional processes included in ONIX are described in [13]. Figure 1 shows the ONIX simulation domain for the ELISE grid system [4], the domain covers one central aperture of the PG and ends at the EG with periodic boundary conditions in the directions perpendicular to the beam axis. A filter field B FF is generated by a current flowing through the PG. Due to the geometry of the extraction aperture, the current density distribution in the PG and the resulting B FF in the extraction region are highly non-uniform. In order to take these non-uniformities into account, the B FF is simulated by FEM using Ansys R Academic Research Mechanical [14], with a spatial resolution of ∼0.5 mm in the ONIX domain. The B FF distribution is included in the ONIX simulations for a PG-current value of 2.5 kA. The B FF , which in the extraction region has a strength of a few mT, magnetizes the electrons and causes a reduction in both the electron temperature and density in front of the PG [15]. The B FF also causes a small vertical deflection of the negative ions.
Permanent magnets embedded in the EG generate a deflection field with a maximum on-axis field strength of ≈60 mT. The B DF extends into the plasma and influences the trajectories of both electrons and ions. The co-extracted electrons are deflected out of the beamlet in the gap between the PG and EG; the extracted H − ions leave the grid system with a horizontal deflection. The B DF is simulated using Ansys R and added to the B FF in the simulations. An overview of the magnetic field topology and the impact on the extraction physics is available in reference [16]. In ELISE, in order to reduce the amount of co-extracted electrons, there is a bias plate (BP) installed upstream of the PG. The BP is not included in the computational domains since a central aperture of a beamlet group is simulated, which is far from the BP.
In ONIX, charged particles (electrons, negative and positive ions) are represented by macro-particles. Their charge is interpolated onto a regular mesh. The upper limit of the mesh size is dictated by the electron Debye length, which needs to be resolved to avoid nonphysical numerical heating (increase of kinetic energy of the charged particles until the effective Debye length is of the same order as the grid size). This practically restricts the simulation domain to one aperture, thus the typical large inhomogeneities in plasma density and extracted current in large ion sources are neglected. Particle Schematic horizontal section of the ELISE ion source and the ONIX/IBSimu simulation domains, marked in red and blue, respectively. In the ELISE source, the plasma generated in the drivers expands onto the grid system. The ONIX domain contains the extraction region of a single aperture, whereas the IBSimu domain includes the extraction and acceleration stages. The deflection magnets inside the EG are not visible in this horizontal section, but the magnetic field direction is indicated. drifts due to the magnetic field on the scale of the ion source are also not well reproduced in a one-aperture simulation domain. In the Poisson-equation solver, Dirichlet conditions are set for the PG and the upstream and downstream boundaries, the other boundaries are periodic. The downstream boundary of the ONIX domain contains the EG aperture, and thus the local potential depends on the space charge of the beamlet itself. The vacuum solution of the potential is used as an initial guess, the impact of this approximation and possible improvements are explored in section 3.1. The Boris-Buneman Leapfrog method is used to step the particles [9]. Complex kinetic interactions between plasma particles (mutual neutralization, associative detachment, non-associative detachment, charge exchange, electron dissociative attachment, electron detachment and electron elastic collisions) are included and treated using MC method. A detailed description of ONIX can be found in [10,[17][18][19] and validation in [10,20].
Protium, further referred to as hydrogen, is used as isotope in this paper in order to connect to previous studies [10,16,20]. The H − ions are produced in the plasma (volume production) and on caesiated surfaces (surface production) [2,11,12].
Surface produced negative ion (H − SP ) are injected on the plasma facing side of the PG. There are several extraction paths for surface produced ions, such as direct and indirect extraction, as schematically illustrated in figure 2.
Electrons and positive ions are injected at the upstream boundary; they represent the plasma source term for the simulation.
Experimental investigations have shown that the negative ion density increases by an order of magnitude during caesiation, due to emission of H − SP ions from the PG [21]. This leads to a formation of an electronegative plasma, i.e. with n H − n e , in front of the PG. In order to reproduce this in H − that exist as a consequence of the caesiation. From which surfaces of the source the latter are created is not completely understood, and for the ONIX calculations it is simply assumed that all the H − entering the simulation domain from bulk plasma have the same temperature and spatial distribution. That significantly simplifies the calculations and is justified as ≈90% of those H − are of the same origin, i.e. they are 'new' H − .
The ONIX input parameters are mainly based on measurements in the ELISE and BATMAN Upgrade ion sources. The electron density n e is 6 × 10 16 m −3 and the electron temperature T e is 2.0 eV [22]. The density of negative ions is highly dependent on the caesiation of the source, as this lowers the work function of the walls, allowing for increased negative ion emission by surface production [23]. The negative ion density 1 cm from the PG can increase by a factor 10 during the caesiation [21]. Here, the negative ion density close to the upstream boundary was set equal to the electron density; a 0.8 eV H − VP temperature was used, as in previously published ONIX calculations [16,20].
The influx of the particles from the upstream boundary is regulated to match the user requested particle densities. Each injected particle species is injected with a Maxwellian distribution. The emission rate of H − SP is set to 200 A m −2 , which is based on work function measurements of 2.1 eV [24]. A 0.6 eV H − SP temperature was used, same as in previous ONIX publications [16,20].
Modelling the plasma from first principles is computationally expensive. Gun-type codes, such as the ion-optics code ion beam simulator (IBSimu), make more approximations but are many orders of magnitude faster [25,26]. In a standalone IBSimu simulation, the full grid system is contained, starting 2 mm before the upstream face of the PG, and continuing 10 mm after the downstream side of the GG. IBSimu includes the modification of the electric field by the beam itself, through the space charge of the previous iteration which is used in the solution of the Poisson equation. IBSimu only tracks the beam particles explicitly, which leads to an artificially high space charge in the plasma region. To avoid a non-physical solution to the Poisson equation, there are compensating charges such that the net space charge in the plasma is zero. These compensating charges are not explicitly tracked particles but modeled by an analytical charge density function corresponding to thermal positive ions of 0.8 eV [26]. Collisions are neglected throughout the grid system. In the plasma, the particles are collisional and obey the magnetic field. It is difficult to model the particle dynamics in a gun-type code and therefore the magnetic field is suppressed inside the plasma region so that the particles trajectories are straight and the aperture is evenly illuminated. Drifts of the negative ions in the plasma, or an uneven illumination of the aperture are possible reasons for differences between simulations that neglect such effects and measurements [27]. Because the physics of plasma creation and wall losses is not included, the plasma potential and the sheath at the PG do not follow from the model: the potential is 0 V in the plasma region. The beam is initialized on the upstream boundary of the domain, so that all the injected particles are effectively from the volume. The current density is an input parameter. To enable the particles to cross the field free region upstream of the PG, they are provided with a starting energy of 3 eV in the beam direction. The particles are created with a perpendicular temperature T ⊥ of 1 eV, which provides an initial divergence of T ⊥ /U tot . The perpendicular temperature is chosen above the normally assumed 0.8 eV, so that the divergence estimates are conservative [28].
The properties of the accelerated ONIX beamlet are calculated by coupling IBSimu to the ONIX domain. In a coupled ONIX-IBSimu simulation, the IBSimu domain starts at the axial plane where the ONIX domain ends, i.e. at the upstream EG surface. The axial positions are given with respect to the boundary plane between the IBSimu and ONIX simulation domains, see figure 1. Particles exiting the downstream boundary of the ONIX domain are further tracked in IBSimu. The particles of 100 ONIX timesteps are used to improve the statistics; the particle weights are adjusted to keep the same current density.
Since the ONIX simulation only stores particle properties at the current timestep, the history of the particles is lost, and it is not possible to link specific particle properties at the exit of the grid system to upstream phenomena. To enable this type of analysis, diagnostic particles are tracked a posteriori through the combined ONIX-IBSimu fields. In the a posteriori analysis, the H − VP ions are injected at the upstream boundary, 15 mm upstream of the PG, with a temperature of 0.8 eV and the H − SP ions are produced uniformly on the conical surface of the PG with a cosine-distribution and an initial temperature of 0.6 eV. These parameters are the same as in the ONIX simulation. The number of tracked particles is on the order of millions to obtain data with a high resolution. The fields are interpolated on a 0.1 mm mesh so that the sheath is resolved.
There is a maximum extractable current density j SC due to space charge limited emission [29]: where 0 is the vacuum permittivity, q is the charge number, e is the elementary charge, m the mass of the extracted species, U ext is the applied extraction voltage, d PG-EG is the distance between the opposing faces of the PG and EG. The normalized perveance is the extracted current density j ext as fraction of the space charge limited emission j SC . The equipotential surface which separates the quasi-neutral plasma from the extracted negative ion beam is called the meniscus. Neither the meniscus nor the electrostatic equipotential surfaces downstream of the meniscus are flat, and the ions leaving the meniscus are convergent in the examples studied, which has a large impact on the beam divergence. Because of the large potential gradient in the extraction gap, the choice of the meniscus potential has little impact on the spatial position when chosen in a reasonable range: at 10 kV extraction potential there is approximately a 70 V increase per grid cell. To avoid incorrectly assigning the meniscus position based on statistical fluctuations of the potential, the 30 V surface is taken as meniscus.
In the context of NBI systems, the width of the beamlet's angular distribution is of particular interest, since it largely determines the transmission of the beam through the beamline. A single Gaussian divergence θ stat. is calculated from the particle angles as: where y is perpendicular to the axial direction. The divergence can be calculated in the horizontal and vertical direction, but the distribution should be centred before using equation (2). When the local v is used, the divergence decreases throughout the grid system as the beamlet is accelerated. Here the final accelerated velocity is used for v , so that the divergence evolution is not obscured by the increase in axial velocity. For more insight into the shape of the angular distribution, a double Gaussian is fitted: a core Gaussian with a smaller divergence θ core , and a halo Gaussian with a larger divergence θ halo . Adding the two components geometrically, weighted by their respective fractions f core and f halo , recovers the single Gaussian divergence θ stat. .

ONIX-IBSimu coupling
IBSimu is coupled with ONIX to extend the computational domain. At the boundary plane, which is the upstream surface of the EG, the particles exit the ONIX domain and are further tracked in IBSimu, self-consistently treating the space charge of the beam with a mean-field approach. The ONIX simulation was initially run until steady-state with the vacuum boundary condition φ vacuum x=0 on the downstream edge. φ vacuum x=0 was obtained by solving the Laplace equation for the full grid system with the IBSimu field solver. When using the vacuum solution of the Poisson equation at the interface between the two domains, there is an unphysical discontinuity in the potential gradient due to the lack of space charge on the boundary, as shown in figure 3(a) for U ext = 10 kV. To obtain the potential at the EG plane including the effect of space charge, the ONIX domain can be extended towards the downstream side, or the IBSimu domain can be extended towards the upstream side. The second option was chosen since it is computationally cheaper. The coupling plane was moved 6 mm upstream so that there is overlap between the IBSimu and ONIX domain to determine the potential on the boundary. 6 mm is chosen since it is sufficiently upstream such that the effect of the vacuum boundary condition is negligible, while still being outside the plasma. Particle properties and the potential on this plane were then exported from ONIX and used at the upstream boundary in IBSimu. Lastly, the non-vacuum potential including the effect of space charge at the EG-plane φ non−vacuum x=0 was exported from IBSimu and used as an updated boundary condition in ONIX.
The resulting on-axis potentials for the coupled ONIX-IBSimu simulations are shown in figure 3(a). The contour lines in the horizontal plane, shown in figure 3(b), illustrate the effect of space charge on the potential, and demonstrate that 6 mm is sufficiently far upstream. The penetration and shape of the meniscus is virtually unchanged. The vacuum boundary condition switches off the space charge on the boundary, which is unphysical, and leads to flatter potential surfaces in the extraction gap. Although the on-axis potential at the EG surface changes almost 1 kV when including the effect of space charge on the boundary instead of the vacuum solution, an additional iteration between the two models (ONIX and IBSimu) yielded a change of only a few volts, indicating that no additional iterations with an updated non-vacuum field are necessary.

ONIX-IBSimu beamlet properties
The extended computational domain allows linking 3D-PIC simulated extraction scenarios to beamlet properties. Two coupled ONIX-IBSimu simulations are presented. In addition standalone IBSimu simulations with matched extracted current are shown. The first case, which is presented in section 3.1, is simulated with U ext = 10 kV and an acceleration voltage U acc = 55 kV, close to the ELISE design operating potentials determined with Kobra3 [30], a gun-type ion optics  code, but below the optimum perveance. The second case has the potentials scaled by half so that U ext = 5 kV and U acc = 27.5 kV to investigate the impact of the changing extraction field. Figure 4 shows a horizontal slice of the negative ion density throughout both simulations domains. The PG, EG, and GG are indicated in gray. The equipotential surface separating the quasi-neutral plasma and the ion beam, the meniscus, is drawn in blue. The extraction paths of H − SP are distinct from the paths taken by the extracted H − VP ions. Many of the extracted H − SP ions are directly extracted near the knife edge of the PG without entering the plasma, as illustrated by trajectory (a) in figure 2. This leads to a different initial phase-space distribution compared to the particles extracted from the plasma volume, and thus a different spatial and angular distribution at the end of the grid system. For the H − SP ions produced upstream of the meniscus to be extracted, their velocity needs to be reversed, this can for example happen due to electrostatic effects, due to bending of the trajectory by the magnetic field [30], or due to collisions with neutrals [31].
The ITER heating NBI extracted current density requirements for H − and D − are 329 and 286 A m −2 , respectively [32]. In the simulation, an extraction potential of 10 kV leads to an extracted H − current density j ex of 149 A m −2 , well below the ITER target. The current density translates to a normalized perveance of 0.10, far below the design value of approximately 0.26 for the divergence optimum of this grid system [30]. For U ext to 5 kV, j ex = 88 A m −2 , which corresponds to a normalized perveance of 0.16 which is closer to the optimum. Due to the shallower penetration of the meniscus in the U ext = 5 kV case, a smaller area of the H − SP emitting PG is located downstream of the meniscus. This leads to a decrease in direct extraction of H − SP in the U ext = 5 kV case (26 A m −2 ) compared to the U ext = 10 kV case (57 A m −2 ). A decrease of H − SP closer to the perveance optimum was also observed  in [33]. The co-extracted electron current was 105 A m −2 and 90 A m −2 for the U ext = 10 kV and U ext = 5 kV, respectively.
In the U ext = 10 kV case, the meniscus curvature corresponds to a large electric field component perpendicular to the beamlet axis, which strongly focuses the beamlet.  depend on the shape of the meniscus, where on the meniscus particles are extracted, and the particle properties at this location. For H − SP ions in particular, these three variables are not independent. The relation between particle properties at the meniscus and in the accelerated beamlet is presented in section 3.5. Figure 5 shows the horizontal divergence, calculated with equation (2), throughout the grid system. H − VP and H − SP ions are shown for the two coupled ONIX-IBSimu simulations, along with standalone IBSimu simulations with matched extracted current. The value and slope of the divergence match well on the boundary between the coupled simulation domains.
In the standalone IBSimu simulation (particle density not shown), only volume produced particles are considered. Upstream of the meniscus, in the plasma region, the ONIX beamlet is more divergent than the standalone IBSimu beamlet, which simply has the T ⊥ /U tot contribution to the divergence. The surface produced particles are more divergent than the volume produced particles, which are qualitatively similar to a standalone IBSimu simulation. The divergence in the ONIX simulations follows from a self-consistent evaluation of the plasma region, whereas in IBSimu the divergence depends on the user-selected perpendicular temperature T ⊥ .
The divergence evolution of the H − SP is different from the H − VP throughout the grid system. The H − VP converges when it is accelerated, which leads to an increased divergence in the extraction and acceleration gap since equation (2) does not differentiate between a converging and diverging distribution. The initially converging H − SP changes to a diverging beamlet after the trajectories cross the middle of the aperture. In the 10 kV case the diverging beamlet and the focusing acceleration field briefly cancel each other leading to a divergence minimum in the acceleration gap. The divergence of the IBSimu beamlet is higher than that of the H − VP ions, but follows the same evolution along the beam axis. The divergence is lower in the 5 kV simulation because it is closer to the perveance optimum.
The divergence determination (equation (2)) assumes that the angular distribution is a single Gaussian; large particle angles heavily bias the divergence value. Figure 6 shows the vertical angular distribution and the divergence values at the end of the grid system at an extraction potential of 10 kV for the volume produced ONIX particles, and for the standalone IBSimu simulation. The horizontal angular distribution has a divergence within 1% of the vertical divergence. The ONIX-IBSimu H − VP and standalone IBSimu angular distribution are not strictly single Gaussians; a double Gaussian is a better fit to the distribution. The volume produced particles have a similar core divergence θ core , which represents a fraction f core of the distribution. The halo divergence θ halo is higher in the IBSimu standalone simulation than for the coupled ONIX-IBSimu H − VP ions. The higher single Gaussian divergence θ stat. for the IBSimu standalone calculation compared to the H − VP is mostly caused by the different amount and width of the halo component.

Code-to-code benchmark
A code-code benchmark was performed between the ONIX-IBSimu and ONIX-ONAC (Orsay negative ion acceleration) results, tracking particles from the ONIX simulation over the same domain, from the upstream side of the EG until 10 mm past the GG with U ext = 10 kV and U acc = 55 kV. The same particle properties as in the ONIX-IBSimu coupling was used, with j ex = 144 A m −2 . ONAC is a 3D PIC-MCC code dedicated to the modelling of negative ion acceleration in ITER-like accelerators [34]. The ONAC algorithms are similar to the ONIX algorithms but both codes have been developed independently. ONAC results presented here neglect the effect of collisions (such as stripping, double stripping etc) on the beam transport, consistent with IBSimu configuration. This comparison is done in order to compare two different models for beam transport, and to show that the results produced by the ONIX-IBSimu coupling scheme are model independent and can be reproduced by coupling with a beam model using 3D-PIC. Figure 7 shows the angular distribution of H − SP at the exit plane for ONIX-IBSimu coupling and ONIX-ONAC coupling. The angular distribution is similar for both cases, both for H − SP ions and H − VP ions (not shown). The surface produced negative ions do not have a Gaussian distribution. A separation between the two beam components, beam core and beam halo, is observed [35]. Although the simulations are run at ITER relevant voltages, the results shown in figure 7 are not in the perveance optimum, the high current density in the beam halo would therefore be significantly reduced when operating in optimum divergence. Finally, the potential maps from the ONIX and ONIX-IBSimu simulations were extracted and stitched together for the a posteriori analysis. Negative ions from the volume and the surface of the PG were tracked until the end of the grid system, 10 mm past the GG. The angular distribution at the end of the grid system generated by the a posteriori analysis was compared with the coupled ONIX-IBSimu and ONIX-ONAC runs, see figure 7. The angular distribution of the a posteriori run agrees well with the angular distribution of the coupled ONIX-IBSimu and ONIX-ONAC runs, showing that equivalent results can be derived using the a posteriori approach. This method is used for an in-depth analysis of the beamlet and meniscus properties in sections 3.4 and 3.5.

Extraction probability and current density profiles
The a posteriori analysis is used to study the extraction probability and the resulting current density profile of the negative ions. The potential map from the U ext = 10 kV, U acc = 55 kV case is used for the analysis. Approximately 95% of the H − SP is extracted in the aperture where they are generated; the average flight distance from injection to crossing the meniscus was 4 mm, indicating that most of the extracted H − SP ions are produced close to the meniscus and directly extracted, see figure 2. The extraction aperture of the H − VP depends on the injection location of the particles: moving the injection location further upstream increases the fraction extracted by neighbouring apertures. The H − VP ions were injected uniformly on the upstream horizontal-vertical boundary, 71% of the particles were extracted through the same aperture as the domain in which they were injected, their average flight distance until crossing the meniscus is 26 mm.
A local minimum of the potential, commonly called a virtual cathode, of approximately 1 V depth is formed close to the PG due to the emission of H − SP ions from the PG [16,36]. The initial location, energy and the depth of the virtual cathode has a large impact on the extraction probability of the H − SP ions. The H − SP ions emitted from the PG are either extracted, reflected back towards the PG, lost by a destructive collision or lost by exiting the upstream side of the simulation domain. The negative ion density, shown in figure 4, suggests that a large fraction of the extracted H − SP ions are produced near the PG knife and directly extracted. This is confirmed by the a posteriori analysis: although the overall extraction probability for H − SP emitted from the PG surface is ∼20%, nearly all H − SP which are emitted near the PG knife, where no virtual cathode is present, are directly extracted. The remaining 80% are mainly directly reflected back towards the PG by the virtual cathode, or exit the simulation domain on the upstream boundary towards the bulk plasma.
Surface produced particles that start further from the PG knife edge need to overcome the virtual cathode, resulting in a lower extraction probability at lower energies. The effect of the virtual cathode on the extraction probability can be seen in figure 8, which shows the probabilities for various outcomes for test particles of H − SP ions as a function of their initial energy E 0 . The depth of the virtual cathode is ∼1 V, thus H − SP ions with initial energy lower than ∼1 eV are immediately reflected back towards the PG. For an initial energy below approximately 1 eV, the majority of the H − SP ions return back to the PG (see the red line in figure 8). The depth of the virtual cathode is not uniform over the PG, so the probability of extraction depends both on the initial location and initial energy of the H − SP ions. In [30], a tendency of reduced extraction probability for higher H − SP energy was shown. However, the effects of the virtual cathode and electric field in the plasma were neglected. H − SP which enter the plasma are accelerated by the plasma sheath (≈2.5T e = 5 V) [37] are therefore rarely extracted, instead they exit the simulation domain on the upstream boundary. 5% of the extracted H − SP are extracted after passing through the periodic boundaries on the horizontal and vertical sides of the simulation domain. Figure 9 shows the average current density as a function of the horizontal position on the meniscus for H − SP and H − VP ions, and the meniscus from the IBSimu standalone simulation. The averaged current density is calculated using the radius of the meniscus, which is 7.5 mm (0.5 mm wider than the radius of the PG aperture). Due to the sheath formation near the PG, the H − VP current density decreases near the edge of the meniscus. In the IBSimu standalone simulations the sheath is neglected, and thus the current density over the meniscus is uniform. For the H − SP ions, a clear peak can be observed for horizontal position ±7 mm, which corresponds to the directly extracted particles close to the aperture edge, the majority of H − SP is extracted from this region. Figure 10 shows a side and front view of the meniscus coloured by the summed current density of H − SP and H − VP ions, and a standalone IBSimu simulation. The current densities in the side-view are high since the plots are line-integrated. For both ONIX and standalone IBSimu, the meniscus is penetrating 3.6 mm into the plasma from the PG knife edge. There is a difference in the distribution of current through the meniscus for the three cases. The directly extracted H − SP ions can be seen in figure 10 as a ring of increased current density close to the PG knife edge. In the case for H − VP ions, the effect of the plasma sheath towards the PG can be seen, i.e. the H − VP ions do not reach the PG. Since the sheath is not included in the IBSimu standalone run, the negative ions are able to cross the meniscus at larger radii close to the PG. The current density distribution over the meniscus influences the properties of the accelerated beamlet, as will be shown in the next section.

From meniscus to accelerated beamlet
The a posteriori analysis is used to link the particle properties at the meniscus to particle properties in the accelerated beamlet. The position where the particles cross the meniscus largely determines the particle angle with respect to the beamlet normal (final angle) at the end of the grid system, as shown in figure 11 for (a) H − VP and (b) H − SP . The directly extracted ions are extracted through the outer edge of the meniscus, resulting in a final angle of >100 mrad, corresponding to the halo in figure 4. The H − SP ions that are extracted through the center region of the meniscus (either by transport to a neighbouring aperture or extraction through the center of the aperture in which they were created) have much smaller final angle. As shown in figure 9, most H − VP ions are extracted through the center of the meniscus, which results in a lower angle at the end of the grid system. Alternatively, one can look at the final angle in terms of the distance the particles travel before extraction. H − SP ions that cross the meniscus after traveling only a few mm from their creation position (direct extraction, see figure 2) have a much higher final angle than H − VP ions. However, the final angle for same-aperture indirect extraction is significantly lower. For the H − SP ions extracted via neighbouring aperture extraction, the divergence matches those of H − VP ions. The position where the particles cross the meniscus strongly impacts the angle between the particles' velocity and the beam axis at the end of the grid system. To check whether this is purely a spatial effect, or linked to the initial particle properties, figure 12 shows the final angle as function of the radial crossing point on the meniscus and the perpendicular energy on the meniscus for H − VP and H − SP ions. The H − VP and H − SP plots are qualitatively similar, which is expected since the particles go through the same electrostatic potential. Particles with the highest final angle are extracted at the largest radial meniscus position, almost independent of their initial perpendicular energy. Particles going through the central region of the meniscus tend to have a low angle at the end of the grid system, unless their perpendicular energy is above a critical value. For H − SP ions, there are more particles with a higher perpendicular energy component in the central part of the meniscus. These particles are accelerated by the PG sheath and travel almost perpendicular to the meniscus until extraction; for H − VP , such an extraction path is highly unlikely. Repeating this analysis for the standalone IBSimu simulation gives the same qualitative picture. Since the sheath is not included in the IBSimu standalone simulations, negative ions are able to cross the meniscus close to the PG. This generates divergent particles at the end of the grid system, which can be observed as a halo in the angular distribution as shown in figure 6. These results show that the Debye sheath that forms in front of the PG has an impact on the ion-optics, and cannot safely be neglected.

Conclusions
Negative ion beamlet formation, extraction and acceleration was studied with the 3D-PIC MCC code ONIX. The ONIX code was coupled with the ion-optics code IBSimu to extend the computational domain to contain an adequate portion of the plasma in the extraction region and the full grid system, a prerequisite for calculating the properties of the accelerated beamlet. With an a posteriori approach, trajectory data was obtained, which enables connecting particle properties in the accelerated beamlet with upstream phenomena.
The ONIX input parameters are informed by experiments, but the extracted current densities are below the ITER target, resulting in a highly divergent beamlet in the simulations. ITER-HNB grid design calculations were done with codes that do not have a physical plasma model and surface production of negative ions. The current results provide insights into the effects of these approximations in low perveance conditions. Particles extracted through the edge of the aperture end up highly divergent in the accelerated beamlet, independent of the perpendicular energy when crossing the meniscus; the perpendicular energy does play a role for the particles extracted through the center of the aperture. The H − SP ions are highly divergent at the end of the grid system, because they are mostly directly extracted near the edge of the aperture. The angular distribution of the H − SP ions is not Gaussian, instead a core component and two highly divergent peaks are observed. The H − VP ions are extracted more evenly over the aperture, but the repelling Debye sheath around the PG decreases the current density at the edges; IBSimu simulations that neglect sheath physics have a flat current density profile. The different extracted current density profiles have an impact on the accelerated beamlet, which in both cases is well described by a double Gaussian. Although the core divergence matches well between ONIX-IBSimu and standalone IBSimu simulations, there is more halo in the standalone IBSimu simulations since more particles are extracted through the edge.
The presented methods enable the study of accelerated beamlets with the 3D-PIC MCC codes. It was shown that the surface produced particles should be included for an accurate estimation on the angular distribution of the accelerated beamlet. It was shown that the Debye sheath in front of the PG impacts the ion-optics and thus cannot be safely neglected. The Debye sheath influences the current density profile at the edge of the aperture, which influences the amount and divergence of the beam halo. Future work includes taking steps towards investigation of the effect of the bias potential on the sheath formation and beam divergence by changing the injection method of the particles, studying the influence of the PG shape on the beamlet divergence, and comparison with experimental results.