3D-PIC modelling of a low temperature plasma sheath with wall emission of negative particles and its application to NBI sources

The 3D particle-in-cell (PIC) code ONIX (orsay negative ion extraction) is a versatile tool for simulating the formation and extraction of negative hydrogen ions and co-extracted electrons in caesiated negative ion sources. Mandatory for modelling these processes is a self-consistent description of the plasma-wall interface, i.e. the plasma sheath, in the presence of surface emitting charged particles. For simplified test cases, the description of the plasma sheath by ONIX is critically validated. Additionally, a set of numerical parameters for reducing the computational cost while still accurately simulating the sheath are presented. Finally, the updated version of ONIX is applied to the plasma of a negative ion source test facility for the ITER NBI system. Results of basic investigations on the influence of the magnetic fields and the meniscus shape on the extraction of negative ions, co-extraction of electrons and the beamlet quality are presented.


Introduction
Simulating plasma phenomena is a necessary approach to understand physical processes when measurements at the experiment are not possible or analytical models are not an option due to the complexity of the system. A particular type for self-consistent simulations of plasma phenomena is the particle-in-cell (PIC) model [1]. Due to their self-consistent character, these codes can be applied, for example, to simulations of the electrostatic plasma sheath.
The plasma sheath is one of the most fundamental phenomena in plasma physics [2,3]. The plasma sheath is a layer of positive charge accumulation in the plasma-wall interface that is created due to the higher velocity of electrons due to their smaller mass. This charge accumulation generates a potential drop (sheath potential) that balances the fluxes of charged particles to the wall facing the plasma. This non quasineutral region of the plasma is very small (extension: several Debye lengths), but of utmost importance for the physics of the whole plasma volume, since it defines (and is defined by) the fluxes of charged particles from and to the surface, as well as plasma parameters like the electron temperature and the electron density.
The physics of the plasma sheath becomes more complex when charged particles are emitted from the surface, as in caesiated negative hydrogen ion sources used in neutral beam injection (NBI) systems [4][5][6][7] or at particle accelerators [8]. Currently ongoing is a strong research process in negative hydrogen ion sources for the heating and current drive system Journal of Physics D: Applied Physics 3D-PIC modelling of a low temperature plasma sheath with wall emission of negative particles and its application to NBI sources at the international fusion experiment ITER [7]. Required from the ITER NBI ion source is an extracted ion current density of 329 A m −2 in hydrogen and 286 A m −2 in deuterium, respectively [9]. In these sources, production of negative ions takes place predominately by conversion of impinging hydrogen atoms on the caesiated low-work function surface of the plasma facing grid containing the extraction apertures [10][11][12]. The plasma sheath can be affected by the surface produced negative ions and a virtual cathode can emerge. This virtual cathode limits the negative ion emission from the surface into the plasma volume [13].
Fully understanding the physics of the plasma sheath in these ion sources is essential since the sheath is related to numerous physical effects such as the co-extraction of electrons, which can strongly limit the source performance. Within this paper, the sheath close to the extraction system of test facilities for the ITER NBI system is studied using the electrostatic 3D PIC-Monte Carlo code (MCC) ONIX [14]. Used as input parameters are particle densities and temper atures measured at the experiments. Up until now, no extensive study of the different available particle injection/re-injection schemes has been done for a case representing negative ion extraction. The first aim of the present investigation is to validate results of the different options used in existing models. Starting from these results, an example application case using the most accurate and computational efficient option will be presented.

Modelling the plasma sheath in NBI sources
The design of the ITER NBI ion source is based on the prototype source developed at IPP Garching. The two experimental test facilities BATMAN (Bavarian test machine for negative ions) and ELISE (extraction from a large ion source experiment) are integral part of the development roadmap [15,16] defined by the European domestic agency F4E for the development of the ITER NBI. Used at BATMAN is the prototype source with 1/8 size (0.3 × 0.6 m 2 , extraction area ≈6 · 10 −3 m 2 ) of the ion source for ITER, while the ELISE ion source has half the size (0.9 × 1 m 2 , extraction area ≈0.1 m 2 ) of the ITER source [9].
This kind of source can be divided into three segments: driver, expansion chamber and extraction system [9]. In the driver the plasma (T e ≈ 10 eV, n e ≈ 10 18 m −3 ) is created by inductive RF coupling (f = 1 MHz). Then the plasma diffuses into the expansion chamber where the electron density and temperature are reduced (T e ≈ 1 eV, n e ≈ 10 17 m −3 ) by a magnetic filter field (FF) in order to reduce the loss rate of negative hydrogen ions by electron collisions [17]. The FF can be produced by permanent magnets or by a current flowing through the plasma-facing grid containing the extraction apertures [18].
The extraction system of these test facilities consists of three multi-aperture grids: the plasma grid (PG), the extraction grid (EG) and the grounded grid (GG). A positive extraction potential in the range of 5-10 kV is applied between the PG and EG in order to extract negative ions. Co-extracted electrons are deflected by a magnetic deflection field (DF), generated by permanent magnets embedded inside the EG, onto the surface of the EG. Finally, the negative ions are fully accelerated by positive electric potential difference set to the GG with respect to the EG (up to a total high voltage of ≈22 kV at BATMAN and 60 kV at ELISE).
In order to restrict the power deposition of electrons onto the EG another requirement for the source is a current of coextracted electrons smaller than the extracted negative ion current [7]. In order to minimize the destruction of negative ions in the acceleration system due to collisions with the background gas, the filling pressure must be 0.3 Pa or below. Up until now, these requirements have not been achieved simultaneously for the required pulse length (1000 s in hydrogen and one hour in deuterium).
The source performance, in particular in deuterium operation, is limited mainly by the co-extracted electron current and its increase during long pulses [9]. Understanding the plasma dynamics-especially the dynamics of the co-extracted electrons-close to the extraction apertures is crucial in order to identify measures to overcome this limitation. Measurements of plasma parameters by means of diagnostics are possible up to ~2 cm upstream of the PG. Physical effects taking place closer to the PG can only be investigated by the application of computational simulations.
Modelling the region close to the PG requires a method that properly describes the plasma sheath and negative particle emission from the surface. Moreover, models should self-consistently evaluate the extraction surface, the so-called meniscus. The shape of the meniscus is of high importance, as it defines the initial angle of the extracted and accelerated particles with respect to the beamlet axis and thus is of high relevance for the beam optics [19]: the meniscus is equivalent to the surface of an electrostatic lens. Particles entering this lens with too high an angle can cause a strong beamlet halo and thus a bad transmission of the beamlet through the extraction system. The physics of the meniscus in negative ion sources is particularly complex for analytical analysis due to the presence of surface produced negative ions and the 3D magnetic field topology. Therefore, simulations determining self-consistently the shape of the meniscus are mandatory to understand the physics of negative ion extraction, the plasma sheath and the beamlet formation.
Such calculations can be done best by electrostatic PIC Monte Carlo codes (PIC-MCC) [1]. These solve self-consistently the non-equilibrium physics and have already been applied (with some simplifications to the system) to model the plasma in such ion sources [13,14,[20][21][22][23][24][25]. Due to a very high computational cost, the domain is usually restricted to a plasma volume close to one extraction aperture. Calculations are performed in either 2D or 3D, but for a full description 3D models have to be applied. The reason is that close to the PG the FF and the DF are perpendicular to each other, breaking any spatial symmetry in the system.

Particle in cell models
Investigated within the scope of this work is the low temperature plasma close to the PG of ITER-relevant sources for negative ions (T e ~ 1-2 eV, n e ~ 10 17 m −3 ). In this energy range electrostatic PIC models are commonly used [1]. These models track macro particles which represent a spatial distribution of real particles, and use a mesh to calculate density distributions, the electrostatic potential and the electric field. The method is time iterative. For each time step Δt the following cycle, schematized in figure 1, is realized: (1) projection of the charge density from the macro particles onto the mesh; (2) the electrostatic potential is obtained by solving the Poisson's equation; (3) the electric field is calculated in the mesh; (4) the Lorentz force is calculated at the position of the particles; (5) the dynamic equation of the particles is solved for Δt; (6) collisions and chemical reactions are calculated by the Monte-Carlo method [26,27].
The most computationally expensive parts of the PIC cycle are the solver for Poisson's equation, the particle charge projection onto the mesh and the particle-pushing algorithm [14]. The computational time consumed by these three modules depend on the number of mesh points N m (which is in three dimensions equal to the product of the number of mesh points in all three directions: N m = N x × N y × N z ) and the total amount of particles N p . Besides N m and N p , the total calculation time of the model depends on the time step Δt and the total number of iterations until the simulation is stopped.
Calculations in 3D consume a large amount of CPU time. For example, an ONIX run for one extraction aperture of a NBI source needs ≈150 h to reach steady state using 4096 cores. If parameter variations are investigated, several such calculations have to be done, multiplying the calculation time. In order to keep the total CPU time as low as possible, the used values of N m , N p and Δt should result from a compromise between computational time and numerical accuracy, respecting the restrictions derived in [1] and introduced briefly in the following along with their physical interpretation.

Mesh size and number of particles
The typical distance in a plasma, in which electrostatic particle to particle interaction is dominant, is the Debye length , where 0 is the vacuum permittivity, k B the Boltzmann constant, T e the electron temperature, n e the electron density and q e the electron charge. Beyond the Debye length, any point charge is screened off by the plasma and only global effects are relevant. It has been theoretically proven in [26,28] that in PIC codes a mesh size larger than 3 · λ D with a first order projection of the charge onto the mesh produces numerical instabilities resulting in non-physical results.
The Debye length decreases with increasing plasma density. The number of mesh points in one dimension is related to the plasma density by N ~ L/λ D ~ n 1/2 e , with L the domain length in the respective dimension. Therefore in a 3D code the total number of mesh points (and consequently also the calculation time) will increase as N m ~ n 3/2 e . An additional stability criterion for PIC codes regards the (average) number of particles per mesh-cell. This number should be large enough in order to reduce the statistical noise and thus ensure that the physics is reproduced accurately. Consequently, for increasing the plasma density and-as described aboveincreasing N m , also the number of macro particles has to be increased for a constant number of particles per cell.

Time step
One of the highest frequencies of global behavior in a plasma is the one of electron density oscillations, ω pe = neq 2 e me 0 , where m e is the electron mass. The time step used in PIC codes should be smaller than the inverse of this frequency. Additionally, the time step should be chosen such that the convergence of a second order particle pusher is ensured as shown in [1].

CFL criteria
The Courant-Friedrichs-Lewy condition (CFL) can be expressed as v max ·Δt < Δx [29], where v max is the velocity of the fastest macro particle in the model. This means that a particle cannot cross over one mesh division in one time step. The physical background of this condition is that the macro particles should see all changes in the electric field while moving through the mesh.

Implementation of the PIC model for NBI sources
The ONIX code, implementing a 3D electrostatic PIC-MCC model, was developed to study the negative ion extraction from negative hydrogen ion sources for NBI. A detailed description of ONIX can be found in [30][31][32].
The plasma parameters (T e ~ 1-10 eV, n e ~ 10 16 -10 17 m −3 ) result in the necessity of using a small mesh size (Δx ~ 10 −5 m), resulting in a very high computational cost. Consequently, it is not possible to model the full ion source in 3D unless some scaling rules for the particle densities or the mesh cell size are applied [33,34]. Reducing the domain to the close vicinity of the PG and restricting to a group of apertures or to even only two apertures is not feasible. Instead, a reduced computational domain around one extraction aperture of the PG is used as schematized in figure 2.
This domain gives an insight into the physics of the plasma sheath, electron and ion dynamics close to the PG as well as the influence of surface produced negative ions and magn etic fields on the meniscus. It is used to test different configurations of parameters in order to suppress the co-extracted electron current or enhance the extraction of negative ions. Additionally shown in figure 2 are the orientation of the predominant components of the FF and the DF. Even using this reduced domain, a calculation using 4096 cores needs the aforementioned ≈150 h computational time for reaching equilibrium.

The used computational domain
For the investigations presented here, the aperture geometry of the large area grid (LAG) system [17], used up to 2016 in BATMAN is applied. The aperture diameter is 8 mm and up to 126 apertures were used in the experiment.
The ONIX domain extends axially between 19 mm upstream of the PG and the upstream surface of the EG, as shown in figures 2 and 3. Along the beam direction the upstream and downstream boundaries of the domain (the prior being equivalent to the transition of the plasma bulk, the latter to the EG surface) are treated as plates while in the directions parallel to the PG periodic boundary conditions are applied. The latter approximation in principle represents the sequential arrangement of an infinite number of apertures and is not capable of correctly describing effects taking place close to the sidewalls of the ion source. The presence of an aperture in the EG is neglected by the model as this simplification does not affect the plasma behavior.
Used at the upstream boundary of the domain, the surface of the PG and the EG, are Dirichlet boundary conditions. The upstream boundary and the PG are set to 0 V, while the EG is set to the extraction potential, typically V ex = 10 kV, which is close to the maximum extraction potential used at the test facilities. While the ITER-like extraction system of ELISE [35] is designed for V ex = 10 kV, this extraction potential is too high for obtaining a good beam optics (i.e. a low divergence) with the LAG system [36]. However, at the test facilities for both the ITER-like and the LAG extraction system extraction potentials up to 10 kV are used in order to extract as many negative ions as possible.
In between the upstream boundary and the PG, a source region is defined. This region is used to approximate the particle fluxes from the thermal bulk plasma into the calculation domain. In order to ensure the desired flux, two conditions have to be satisfied: (i) the outflow of particles from the source region is constant over time, (ii) particles inside this region have a Maxwellian energy distribution. If these conditions are met, it is possible to compare the results from the PIC simulations with analytical models that use a constant Maxwellian particle flux as boundary condition.
A constant outflow from the source region implies that the number of particles in this region should be constantly repopulated over time. Consequently, an injection or re-injection of lost particles has to be implemented in the model. Used by PIC codes investigating the extraction of negative ions [13, 21,   Test region for physical tests (dotted blue). Space used to represent the plasma source (dashed red). Negative ions emitted from the PG are separated according to their origin from the flat or the conical part of the aperture. The EG aperture (dotted green) is not included in the model. 23, 24] are different particle injection/re-injection schemes. Up until now, the effectiveness of the different schemes in reproducing the properties of the plasma sheath was never investigated. Results of such an investigation are presented in section 4.2.
In order to ensure a Maxwellian energy distribution, a thermalization routine is applied to the particles in the source region i.e. the particle velocities in the source region are reset to keep the energy distribution constant over time. Results of validating this thermalization routine are described in section 4.3.
For performing basic tests or verifying and validating the code, the test region (shown in blue in figure 3) is used, typically without magnetic fields, in order to reduce the CPU time needed. This simpler domain consists of a plasma in between two absorbing plates and the source region in between. Dirichlet boundary conditions for the electrostatic potential are used, and the plates are set to 0 V. The dimensions parallel to the plates are treated with periodic boundary conditions. This model represents a uniform plasma in the directions parallel to the plates and it is capable of predicting parameter dependencies of the properties of the plasma sheath (such as potential and width).

Re-injection of particles
The different re-injection schemes used by PIC codes for negative ion sources are described in the following and for the case of a positive ion-electron plasma they are illustrated in figure 4: Single particle re-injection [24] • If a charged particle hits a wall, a particle of the same species is re-injected into the source region. Pair re-injection [21,23] • If an electron hits a wall, it is removed.
• If a positive ion hits a wall, it is removed and an ionelectron pair is injected into the source region. • If a surface produced negative particle hits a wall, it is removed.
Constant flux [13] • If any particle hits a wall, it is removed.
• Pairs of electrons and positive ions are injected, with a constant rate, into the source region at each time step in order to keep a constant density in this region.
In order to benchmark the physical correctness of the different injection/re-injection schemes, calculations have been done using the test region. This test is done for a positive ionelectron plasma, an electron temperature and positive hydrogen ion temperature of 1 eV, and a plasma density of 10 16 m −3 .
For the initial state, a quasi-neutral plasma in the source region is used, i.e. total electric charge equal to zero. The electric charge difference in the plasma sheath should then be selfconsistently generated by the model along with the applied injection/re-injection schemes. Figure 5 shows the 1D distribution of the electrostatic potential along the x-axis calculated using the three injection/ re-injection schemes. The pair re-injection and flux injection methods, at steady state, give similar results. For these schemes, the spatial structure of the electrostatic potential drop agrees with the well-known structure of the plasma sheath. In contrast, the single particle re-injection method shows no clear structure in the potential.
These differences between the results for the different injection/re-injection schemes can be understood as follows: the pair re-injection scheme keeps the number of ion particles present in the calculation domain constant over time, while the number of electrons varies until an equilibrium is found by the system. This steady state is found when the number of electrons remains constant over time, meaning that the number of lost electrons and ions at the surface facing the plasma is the same, i.e. both fluxes are equal. The flux injection scheme reaches a steady state when the number of particles of each species is constant over time, which means that the number of pair of particles leaving the source region and the pairs reaching the surface are equal. This implicitly implies that the flux of electrons and ions to the surface are equal as it is the case for the flux injected into the source region. The single particle re-injection scheme preserves the total number of particles, i.e. the total charge. Since the initial condition consists of an overall null electric charge, the creation of a plasma region with electric charge accumulation, i.e. the sheath, is not possible.
For the different injection/re-injection schemes the sheath potential (in the steady-state) was compared with the predicted one of 2.56 V by the analytical formula ϕ = − Te 2 ln(2π me mi + (1 + Ti Te )) from [37], where m i and T i are the ion mass and the ion temperature, respectively. The sheath potential from ONIX using the pair re-injection scheme is 2.53 V, which is in good agreement (discrepancy of only 1%) with the analytical formula. Besides the good agreement of the sheath potential, the pair re-injection method has the advantage that the plasma density in the system can be defined easily, as the number of ions is kept constant over time. Using the flux injection scheme, a good agreement with the analytical form ula was also found, having a difference of 4% with respect to the analytical description. The flux injection scheme has the advantage that it gives more freedom in the kind of sheath that the system can recreate, as it does not impose any condition on the total electric charge of the system. The calcul ation time needed for reaching the equilibrium state, however, can be significantly higher than for the pair re-injection scheme.
For further validation and application of ONIX, the pair re-injection is chosen since it allows easy definition of the plasma density in the system.

Thermalization process
Injecting or re-injecting particles into the source region should result in constant particle densities and particle temperatures in this region. Non-constant densities or temperatures can affect the sheath potential and thus particle fluxes onto the surface, but also the fluxes of surface emitted charged particles into the plasma.
Particle injection/re-injection is performed using a Maxwell distribution for the particle velocities (defined by the particle temperature) and a uniform distribution of the particle positions. The particle flow out of the source region is proportional to the particle velocity and thus depends on the particle energy. Faster particles escape with a higher flow and they are replaced by the same number of particles with a Maxwell distribution. This procedure results in a reduction of the average particle energy in the source region. To maintain the desired velocity distribution requires either some kind of artificial heating or the use of a more appropriate distribution function for injecting/re-injecting the particles, as suggested in [3].
Introducing an artificial heating process is the simplest option, and it can be implemented by resetting the velocity distribution of the particles with a sufficiently high frequency [23]. This is done in ONIX for the electrons, using a rate of 2 × 10 9 s −1 [23] for resetting the particle velocities based on a Maxwell velocity distribution function. Calculations with and without this thermalization procedure have been done for the test region with T e = 1 eV and n e = 10 16 m −3 . The parameters have been chosen to be as close as possible to the ones of the NBI sources. The electron density, however, is reduced by a factor of ten in order to reduce the computational cost.
The electron temperature in the source region with and without thermalization is shown in figure 6. Without thermalization a reduction of the temperature takes place, reducing the temperature from 1 eV to 0.7 eV, i.e. by about 30%, after 0.5 µs of simulated time. According to the analytical formula given by [37] this temperature reduction causes a reduction of the sheath potential from 2.5 V to 1.7 V. Therefore, thermalizing the source region plasma is mandatory for PIC models treating the production and extraction of negative ions.

Plasma sheath benchmark
For simplified scenarios, analytical formulae can predict the sheath potential in a 1D collisionless plasma. The sheath potential is found to be dependent of the electron temperature and the masses of the plasma particles [37,38]. However, more generally the properties of the sheath, such as the potential drop and the spatial extension, depend on the fluxes of charged particles coming from the plasma towards the wall.
In the following, simplified scenarios for which analytical solutions for the plasma sheath exist are used as test cases for ONIX. In a second step, variations of N m , N p and Δt are performed in order to define the best compromise between computation time and numerical accuracy. For all calculations, parameters as close as possible to the ones of the NBI sources have again been used and the electron density has been reduced in order to reduce the computational time.
In order to validate the dependence of the sheath potential on the electron temperature, the electron temperature is varied between 0.5 eV and 5 eV, at a fixed plasma density of 10 16 m −3 in the source region. The ONIX results are compared against the analytical plasma sheath models by Stangeby [37] and Schwager [39] and results from a 1D PIC code [13]. For low electron temperatures (T e < 3 eV), the results are in very good agreement with the analytical solution and the 1D PIC model, as can be seen in figure 7. For higher electron temperatures, the disagreement slightly increases, up to about 10% for T e = 5 eV.
Additionally, the dependence of the sheath length on the plasma density was analyzed. The sheath length is in the scale of a few Debye lengths, which depends on the plasma density. Since no analytical formulae for the sheath length are known, the results are compared with those of the 1D PIC code [13] and the numerical solution of the Schwager analytical model [39]. Figure 8 shows a good agreement between ONIX results and the two other models.

Secondary particle emission
Correct treatment of secondary particle emission from walls into the plasma is crucial for models describing the low pressure, low temperature plasma of negative ion sources for NBI with a high relevance of negative ion production on the caesium covered low-work function surface of the PG.
Two cases are studied in order to validate ONIX for surface particle emission: the emission of secondary electrons due to impinging electrons, and the emission of negative ions from the wall after a neutral atom hits the wall. For both casesemission of secondary electrons and emission of negative ions-the presence of additional surface produced negative particles results in a reduced sheath potential and eventually to the creation of a potential well.
The emission of secondary electrons caused by primary electrons impinging a surface has been deeply studied and analytical formulae describing the dependence of the sheath potential on the secondary emission probability γ have been deduced by Stangeby [37] and Schwager [39]. Figure 9 presents the sheath potential calculated by ONIX for a variation of γ for a positive ion-electron plasma. The ONIX results are compared with results of the analytical models and the 1D PIC code from [13]. The comparison is done for an electron and proton temperature of 2 eV since the 1D PIC code calculations have been done using this temperature. The secondary electrons are launched at the surface with a temperature of 1 eV.
The secondary electrons increase the negative space charge in the plasma close to the surface. With increasing γ this space charge increases, resulting in a decreasing sheath potential. When increasing γ above a critical value (critical SEC), the negative space charge reaches a value that creates a virtual cathode. Above the critical SEC, the analytical models are no longer valid: in order to describe the virtual cathode, separate models are needed to describe each region as described in [40]. The ONIX results are in good agreement with the analytical results and the 1D PIC code as can be seen in figure 9, presenting a maximum deviation of 8%. This indicates that the pair particle re-injection can correctly describe the emission of charged secondary particles at a surface and their impact on the plasma sheath.
Next, the effect on the sheath potential of negative hydrogen ion emission from a surface is validated. The test region is used and physical parameters are the same as former  [37], Schwager [39] and the 1D PIC code from [13]. Numerical parameters for the simulation are Δx = 0.5 · λ D and Δt = 0.05 · ω −1 p . Figure 8. Sheath length versus the plasma density for a positive ion-electron plasma. ONIX results are compared with the numerical solution of the Schwager analytical model [39] and the 1D PIC code from [13]. Numerical parameters for the simulation are Δx = 0.5 · λ D and Δt = 0.05 · ω −1 p .
ONIX simulations [30], with the exception again of a reduced electron density due to the computational cost. A fixed surface production rate of 55 Am −2 for the negative ions is implemented, estimated from [13] for typical parameters of the NBI test facilities. A fixed production rate is needed because the dominant surface production of negative ions by conversion of hydrogen atoms cannot be implemented self-consistently into PIC codes (treating charged particles only). Shown in figure 10 is the spatial profile of the positive ion density, the electron density and the negative ion density calculated by ONIX, the model by Schwager [39] and the 1D PIC code from [13]. Similar as in the case of secondary electron emission, the surface produced negative ions deposit an additional negative space charge close to the surface, modifying the structure of the plasma sheath and thus affecting the transport of charged particles from and to the surface. If this space charge is large enough, a virtual cathode can evolve.
The profiles are on top of each other, showing again the capability of the pair particle re-injection in general, and of ONIX in particular, to reproduce the physics related to the emission of negative particles from the wall.

Numerical parameters efficiency
A plasma density of about 10 17 m −3 is measured in the negative ion test facilities close to the PG, ten times higher than the density used for the test calculations presented above. Additionally, calculations for the NBI sources need to use the full domain instead of the test region. Consequently, the computational cost is significantly higher (the calculation time increases by a factor of ≈700) compared to the test calculations. Simulations for n plasma ~ 10 17 m −3 are only possible using high performance computing (HPC) systems with a well-parallelized code as ONIX [14].
Due to the limited amount of available computation time at HPC systems, the amount of code runs can be strongly limited. In order to gain as much physical insight as possible, it is crucial to identify the optimal set of numerical parameters for ONIX calculations, i.e. reducing the computation time of the simulations without altering the simulated physics. In order to do so, a variation of three numerical parameters is preformed: the mesh size, the time step and the number of particles per cell.
Increasing the mesh size reduces the number of mesh points (N m ) and consequently decreases the number of computational operations during the charge projection procedure and solving Poisson's equation.
With a first order projection scheme [1] numerical convergence and physical accuracy are only ensured for a mesh size lower than λ D . Nevertheless, it has been numerically proven that higher order projection schemes allow to use a larger mesh size [28], but the model results should be tested in each case. ONIX uses a linear second order scheme for the particle charge projection [41]. The physical accuracy was tested by comparing the structure of the plasma sheath calculated for a positive ion-electron plasma, using different mesh sizes while keeping a high number of particles per cell (>70) and a small time step: Δt = 0.05 · ω −1 p . For this calculation, plasma parameters close to the ones in the NBI sources have been used, and again the electron density was reduced by a factor of ten in order to reduce the computational effort.
The results of these simulations are shown in figure 11. The Debye length is 5 × 10 −5 m. Significant spatial oscillations occur in the plasma potential along the whole domain for a mesh size equal 8 · λ D . Additionally, in this case the spatial width of the plasma sheath is significantly larger when compared to smaller mesh cells, and it is described by two mesh points only. For a mesh size of 4 · λ D , the plasma sheath still has a different spatial structure compared to the calculations Figure 9. Sheath potential versus the secondary electron emission probability. ONIX results are compared against the analytical models of Stangeby [37], Schwager [39] and the 1D PIC code from [13]. Numerical parameters for the simulation are Δx = 0.5 · λ D and Δt = 0.05 · ω −1 p . Figure 10. Density of species in the plasma with emission of negative ions from the wall. ONIX results are compared against the analytical model of Schwager [39] and the 1D PIC code from [13]. Numerical parameters for the simulation are Δx = 0.5 · λ D and Δt = 0.05 · ω −1 p .
done with smaller mesh sizes, but the number of mesh point defining the plasma sheath has increased to 3. For mesh cell sizes below 2 · λ D , no considerable differences in shape of the sheath can be noticed, and only the number of mesh points defining the sheath increases. Therefore, a mesh size of 2 · λ D can be used in ONIX without numerically perturbing the solution of the plasma sheath. It has to be kept in mind, however, that the spatial extent of structures created by additional physical effects, such as a virtual cathode resulting from surface emission of charged particles, can be smaller than the typical extent of the plasma sheath. If such phenomena are included in simulations, it has to be checked specifically if their typical length scale is resolved by the used mesh. A similar study is done for the number of particles per cell. Using 20 particles per cell gives no difference in the average value of the electrostatic potential in comparison to a significantly higher number of particles per cell. Only the statistical noise is slightingly increased, as expected.
Finally, to satisfy the CFL criteria, a maximal velocity for the particles has to be defined. Due to the much smaller mass of electrons, in the plasmas under investigation, the thermal velocity of electrons is higher than the one for all other particle species. In order to consider the majority of electrons, a velocity five times the electron thermal velocity is chosen to calculate, based on the mesh size, the upper limit for Δt as defined by the CFL criteria.
Based on the present results, a set of numerical parameters is chosen that should result in physically correct results while reducing as much as possible the consumed CPU time. Compared to the determined limits for physically correct results described above, these parameters are relaxed by 25% to 50% in order to ensure an additional safety margin.
In order to prove that using this parameter produces correct results, a second set of parameters (set 2) is used with a factor of three smaller mesh size, a factor of three smaller time step and (in average) 70 particles per cell instead of 30. This second parameter sets over-fulfills the numerical restrictions derived for PIC models in [1]. The two sets of parameters are compared in table 1.
It can be assumed that simulation results based on parameter set 2 are correct. This assumption is supported by a negligible numerical noise of the model results. Comparing results based on the relaxed parameters from parameter set 1 with the correct ones from parameter set 2 indicates if the results obtained using set 1 are numerically perturbed or not.
In order to perform such comparisons, the benchmark calculations on the sheath potential and the sheath length performed for the test region, presented in section 5.1, are repeated using the two different parameter sets. The results for the sheath potential and the sheath length are shown in figures 12 and 13, respectively. Additionally, the analytical results from the model by Schwager [39] are shown. No significant discrepancy between the two numerical parameter conditions and the analytical results can be seen.
Consequently, parameter set 1 is applied for future calculations, resulting in a reduction of the computation time up to 20% compared to parameter set 2.

Application to a negative hydrogen ion source for fusion
ONIX is applied for basic investigations on the production and extraction of negative hydrogen ions from an ion source test facility for ITER NBI. The pair injection procedure was used Figure 11. Simulation of the plasma sheath varying the mesh size. Using a mesh size above 3 · λ D results in a numerically perturbed electrostatic potential.
Particles per cell (#PPC) 30 70 Figure 12. Sheath potential versus the electron temperature in simulations using parameter set 1 and parameter set 2 and from the analytical model.
together with particle thermalization in the plasma source region. The mesh size, time step and the number of particles per cell were taken from parameter set 1 derived in section 5.3 (table 1). The simulation domain starts ~2 cm upstream of the PG, i.e. the axial distance from the PG at that plasma parameters are typically measured in the test facilities [42,43]. Thus, measured values of particle densities and temperatures can be used as initial conditions for the model. The domain ends at the upstream surface of the EG, as shown in figures 2 and 3. An extraction potential of 10 kV is used. As mentioned above, for the LAG system [17] this potential results in a high extracted negative ion current, but at non-optimal beam optics [36]. A list of the used parameters taken from measurements or from previous ONIX simulations [24] is presented in table 2.
Surface production of negative hydrogen ions is implemented assuming a constant and uniform negative ion flux from the surface of 550 Am −2 [13]. This flux is the result of multiplying the influx of hydrogen atoms, most of the negative ions are produced by surface conversion of atoms [44,45], with the energy conversion probability [43,46]. A cosine distribution is used for the velocity direction [47].
In the experiment high negative ion densities up to n NI ≈ 10 17 m −3 are observed in a few centimeters distance to the PG [9], however, none of the existing 3D PIC models for negative ion production and extraction can reach self-consistently these negative ion densities measured in the plasma bulk [33]. The cause of this inconsistency is not yet known. One possible explanation is a too high depth of the virtual cathode, created by surface produced negative ions and limiting the negative ion emission from the surface into the plasma volume. It is not possible, however, to prove or disprove this explanation: measurements (e.g. using Langmuir probes) of the virtual cathode are impossible due to its extreme proximity (a few λ D ) to the PG surface. Additionally, other effects taking place on a global scale, and thus being neglected when using a reduced calculation domain around one aperture of the PG, may play a role [40].
The too low transport of negative ions into the plasma and possible counter-measures will be further investigated in the future. For the simulations done within the scope of this publication, an initial amount of negative ions is added into the source volume according to experimental results, as shown in table 2. These negative ions are treated differently from those emitted from the PG: single particle reinjection condition is used in order to keep the negative ion density constant over time in the plasma source region. In contrast, surface produced negative ions are treated with the pair re-injection scheme. The complex 3D magnetic field, composed by the FF and DF, is calculated by the PerMag code [48].
As initial condition, the particles are uniformly distributed in the source region. The plasma expands into the whole domain, and the meniscus is formed close to the extraction aperture by the equilibrium of the applied extraction potential with the plasma pressure. As a result, the domain is separated into two regions: the quasi-neutral plasma, and the region in that the electrostatic potential is no longer screened and only extracted negative particles are present. This separation can be seen in the top part of figure 14, showing the spatial distribution of positively charged particles in the x-y plane (the main component of the FF is along the y axis, the DF along the z axis, see figure 2). A strong drop of the positive ion density towards the meniscus is observed. The meniscus, being defined as the first isoline of null positive ion density, is indicated by a white line in figure 14.
Shown in the bottom part of figure 14 are profiles along the beamlet axis (indicated in by the horizontal green dashed line in the upper part of figure 14) of the electrostatic potential and the densities of positive ions, negative ions and electrons. Additionally, these profiles indicate the position of the meniscus: it is located at the point where the positive ion density approaches zero.
A strong decrease of the electron density, from ~10 17 m −3 to ~10 14 m −3 , takes place between the source region and the meniscus. The reason is a small coefficient for diffusion of the magnetized electrons perpendicular to the magnetic FF and the DF. The electron diffusion is, however, not symmetric over the 3D space. As can be seen in figure 15, showing the electron density for two perpendicular planes, a 3D distribution of the electrons over the simulation domain is developed. This is due to the different topologies of the two components of the  Due to their much higher mass, the negative ions are not magnetized and their density is not affected by the magnetic fields. The positive ion density follows the decreasing profile of the electron density in order to keep the quasi-neutrality of the plasma. The combination of the described effects results in the formation of an ion-ion plasma near the meniscus. Additionally, the co-extracted electron current is significantly reduced compared to the case without magnetic fields. This means that varying the intensity or topology of the magnetic fields can be a key factor for further reducing the co-extracted electron current in the experiment.
The potential in the source plasma has an almost constant value of ~2.7 V, being in close agreement with the sheath potential expected (for an electron-ion plasma) for the electron temperature of 1 eV. Close to the meniscus the potential increases, representing some kind of pre-sheath of the applied 10 kV positive extraction potential. Figure 16 shows   The general shape of the electrostatic potential looks like a normal plasma sheath. Close to the PG surface the virtual cathode can be seen, limiting the number of negative ions generated at the PG that can be transported into the plasma bulk. Responsible for the creation of the virtual cathode is a peak in the negative ion density that strongly exceeds the positive ion density. The depth of the virtual cathode varies between 1.05 V and 1.4 V at different locations of the PG. These values are much lower compared to results of former ONIX calcul ations [24]; these showed a depth of the virtual cathode of up to around 7 V. The reason for the now much shallower virtual cathode is that values for the time step, the mesh size, and the number of particles per cell are used for that the correctness has been verified (section 5.3). The virtual cathode is defined by four mesh points, which is sufficient for a proper description of the physics of this region.
As a result of the less deep virtual cathode, the negative ion density, coming from the surface, in the plasma volume is increased from ~10 13 m −3 in previous calculations [49] up to ~10 16 m −3 . Although this increase is significant, the value is still one order of magnitude below the value of ~10 17 m −3 measured in the experiment. This disagreement is the subject of future investigations. Figure 17 shows the temporal behavior of the total extracted negative ion current density and the density of the co-extracted electron current, measured at the EG. Additionally, the fractions of the negative ion current density coming from the source region (coming from the artificially increased negative ion density in the source region), from the conical part, and the flat part of the PG are shown.
In the first microsecond of the simulation, strong fluctuations of the co-extracted electron current density appear, correlated with a peak of ~400 A m −2 , electron-ion ratio ≈2. During this time, the plasma does not extend up to the PG, the extraction potential is not yet fully screened by the plasma, and the electrons are accelerated directly towards the EG.
A steady state of the simulation can be distinguished from the total extracted negative current density, saturating after ≈1.4 µs. After the plasma reaches the PG, the coextracted electron reduces until equilibrium is reached at an electron-ion ratio of ≈0.46. This value corresponds to typical values for a well-caesiated source [9].
The simulation indicates that the majority of extracted negative ions are created at the conical part of the PG. The effect of this fact on the beamlet optics is illustrated in figure 18, showing in the left part three density maps, again in the x-y plane, for the negative ions according to their origin (plasma volume, flat part and conical part of the PG); additionally indicated by white lines is the meniscus. Shown in the right part of figure 18 are the respective density profiles of negative ions at the surface of the EG.
Negative ions from the plasma volume (top row of figure 18) reach the meniscus predominately in its central (and more flat) part, as well as, with a preferred direction, resulting in a well-focused part of the beamlet. The same holds in principle also for negative ions generated at the flat part of the PG surface (middle row). Due to the virtual cathode, only a small number of these negative ions can reach the plasma from where they can be extracted. In order to illustrate their density profile at the EG surface, the inlet in the right column of figure 18 shows a zoom (factor of 100) of the density.
The meniscus penetrates deeply into the plasma, and a fraction of the negative ions from the conical part of the aperture (bottom row) is directly extracted. This means that these ions do not spend enough time inside the plasma to be affected by its collective behavior before extraction. The directly extracted negative ions cross the meniscus very close to its edge with a  high curvature. Consequently, the directly extracted negative ions form a halo in the extracted beamlet.
Qualitatively, this effect can be measured in the BATMAN test facility: for the LAG system and an extraction potential of 10 kV the beamlet is far away from the perveance optimum and has a high divergence. When reducing the extraction potential to 5 kV the optics of the extracted beamlet improves as expected, as shown in figure 19. Table 3 summarizes the results of the extracted negative ion current and the co-extracted electron current for the extraction voltages 5 kV and 10 kV.
Up until now, however, these results could not be compared with experimental results: at the test facilities, the beam quality can be determined only on a global scale, not on the level of the individual beamlets. At the recently started BATMAN upgrade test facility [50], a tungsten wire calorimeter located 19 cm downstream the extraction system offers the possibility of roughly investigating the shape and structure of single beamlets.
Since ONIX can distinguish the origin of the negative ions and their effect on the extracted beamlet, the code can give hints, which parameters lead to a higher extracted current and better beam optics. In addition, measures for reducing the coextracted electrons can be investigated. Free parameters to be varied during such investigations are the geometry of the grids, the topology of the magnetic fields or the physical parameters such as particle densities in the plasma or the extraction potential. Such variations can help improving the performance of negative ion sources towards the requirements for the ITER NBI system and are planned for the future.

Summary
Self-consistent modelling in 3D is mandatory for improving the understanding of the plasma dynamics close to the extraction aperture in negative hydrogen ion sources for fusion. The 3D PIC-MC code ONIX was validated for application to    0.14 such ion sources. An appropriate particle re-injection scheme was identified. Validations were successfully performed for the behavior of the plasma sheath when varying the electron temperature, the plasma density and the emission rate of secondary emission of particles from a surface. The mesh size, time step, and number of particles per mesh cell used by the model were optimized in order to reproduce the plasma sheath physics and simultaneously reduce the computation time as much as possible. The chosen set of parameters reduces the computation time by more than 20% with respect to the one used during the validation tests. This reduction is quite relevant, taking into account that even using several thousand CPU cores the calculation time is higher than 100 h.
Calculations for a negative ion source for fusion show extracted negative ion currents and electron-ion current ratios that are in good agreement with those measured in the experiments for a well caesiated source. The depth of the virtual cathode at the PG surface is now much lower compared to results of former ONIX calculations. While this represents a step towards describing the transport of surface produced negative ions into the plasma bulk self-consistently, the present results still underestimate the amount of these ions reaching the plasma, and the negative ion density in the plasma has to be increased artificially. This will be the topic of future investigations.
The majority of extracted negative ions are produced on the conical part of the extraction apertures. Since the meniscus penetrates deep into the plasma, some of the negative ions from the conical parts are extracted directly or very close to the meniscus edge. These negative ions create a halo around the extracted beamlet. This demonstrates that the position and the shape of the meniscus play a critical role in the halo formation. Further calculations on the correlation between PG geometry and the meniscus shape are in preparation.
The present work gives valuable insight into the physics of the plasma sheath in the case of wall emission of negative hydrogen ions. Additionally, a 3D PIC code has been validated for physical effects of the plasma volume close to one extraction aperture of negative ion sources for fusion. This code can describe physical effects not accessible by diagnostics, but of high relevance for further improving the performance of such ion sources. Thus, future application of this code can provide essential insight for reaching the target values of the ITER NBI system.