Modeling of plasma transport and negative ion extraction in a magnetized radio-frequency plasma source

Negative ion sources for fusion are high densities plasma sources in large discharge volumes. There are many challenges in the modeling of these sources, due to numerical constraints associated with the high plasma density, to the coupling between plasma and neutral transport and chemistry, the presence of a magnetic filter, and the extraction of negative ions. In this paper we present recent results concerning these different aspects. Emphasis is put on the modeling approach and on the methods and approximations. The models are not fully predictive and not complete as would be engineering codes but they are used to identify the basic principles and to better understand the physics of the negative ion sources.


Introduction
Negative ion sources are used in a variety of research fields and applications [1] such as in tandem type electrostatic accelerators, cyclotrons, storage rings in synchrotrons, nuclear and particle physics (for instance to produce neutrons in the Spallation Neutron Source [2]) and in magnetic fusion devices (generation of high power neutral beams [3]). High brightness negative ion sources (i.e., which produces large negative ion currents) use cesium vapor to significantly enhance the production of negative ions on the source cathode surface. Cesium lowers the work function of the metal and hence facilitates the transfer of an electron from the metal surface to a neutral hydrogen atom by a tunneling process. The main types of devices that use cesium are magnetrons, Penning and multi-cusps ion sources. The former have applications in accelerators and the latter are often large volume ion sources like those developed for fusion applications. The plasma in large volume devices can be generated by hot cathodes (heated filaments) or radio-frequency (RF) antennas (inductively coupled-plasma, ICP discharges) standing either inside or outside the discharge [1]. Ion sources for fusion are tandem type devices with a so-called expansion chamber juxtaposed next to the discharge region. The expansion chamber is often magnetized with magnetic field lines perpendicular to the electron flux exiting the discharge. The magnetic field strength is typically of the order of ∼100 G and is generated either by permanent magnets placed along the lateral walls of the ion source or via a large current flowing through the plasma electrode (which is also called 'plasma grid'). The plasma grid (PG) separates the ion source plasma from the accelerator region, where the extracted negative ions are accelerated to high energies. The axial electron mobility is strongly reduced by the magnetic field inside the expansion chamber and the electron temperature is hence significantly reduced as electrons loose a large amount of energy through collisions. In ion sources for fusion, the background gas pressure (either hydrogen or deuterium type) is ∼0.3Pa and the electron temperature is of the order of 10eV in the discharge region. The magnetic filter reduces the electron temperature down to the eV level in the extraction region, close to the PG. The role of the magnetic filter field in the expansion chamber is threefold: (i) a large versus low electron temperature between the discharge and the extraction region allows the production of negative ions through the dissociative impact between an electron and an hydrogen (or deuterium) molecule n ( )  H 4 2 , where ν is the vibrational level. The vibrational excitation of the hydrogen molecule is maximized at high electron temperatures (typicallyT 10 e eV) while the cross-section for the dissociative attachment of H 2 and hence the production of a negative ion is the largest forT 1 eV vicinity of the PG significantly increases the survival rate of the negative ions and (iii) the magnetic filter more or less lowers the electron flux onto the PG but this is not sufficient and a suppression magnetic field is used to limit the co-extracted electron current. Co-extracted electrons have a damaging effect inside the electrostatic accelerator [4]. The electron beam is unfocused and induces a large parasitic power deposition on the accelerator parts. Note that in fusion-type, high power, large volume and low pressure ion sources, negatives ions produced via dissociative attachment of the background gas molecules (so called 'volume processes') range between 10%-20% of the total amount of extracted negative ion current [5,6], the remaining part corresponds to ions generated on the cesiated PG surface through neutral atom and positive ion impacts. In magnetic fusion applications, negative ion sources are a subset of a neutral beam injector (NBI) producing high power neutral beams which are injected into the Tokamak plasma. Neutrals are insensitive to magnetic fields and can hence penetrate into the hot plasma core. The neutral beam provides power to the plasma, current (which is necessary to sustain the poloidal field) and are helpful to minimize the buildup of some type of instabilities. In the future International Thermonuclear Experimental Reactor (ITER), NBIs are designed to inject 33MW of power (split over two beam lines) with an energy of 1MeV into the Tokamak plasma [7]. The ITER project is the first fusion device which will mainly be heated by alpha particles ( + H e 2 ). The plasma will consist of Deuterium and Tritium ions providing 500MW of fusion power. 50MW of additional external power will be necessary in order to heat and control the plasma during the operating phase while the alpha particles will re-inject 100MW of power to the fusion plasma (the total heating power is 150 MW). The remaining 400 MW is carried by the neutrons toward the wall of the Tokamak [8]. The external heating system for ITER also includes 20MW of electron cyclotron heating at 170GHz and 20MW of ion cyclotron heating in the - 35 65 MHz frequency range [9]. Total power is consequently 73MW (including neutral beams), slightly above the required 50MW for ITER.
In this paper we illustrate and analyze, on the basis of new results, the work performed in our group in the last ten years on the modeling and simulation of the negative ion source and negative ion extraction. The modeling of the negative ion source in all its complexity (power absorption, plasma chemistry, coupling with neutral transport and chemistry, transport across magnetic field, negative ion production and extraction) is a formidable task and we address the different questions separately with dedicated models using simplifying assumptions. In most cases the models are applied to the ITER prototype negative ion source BATMAN [5,[10][11][12] (BAvarian Test MAchine for Negative ions) developed at the Max-Planck-Institut für Plasmaphysik, Garching, Germany. The source is a tandem type device similar to the ITER configuration but with one ICP discharge (driver) and a smaller expansion chamber volume, accordingly. The driver dimensions are a cylinder of diameter 24.5cm and length 16cm [11,13]. An external cylindrical antenna confers to the background gas (either molecular hydrogen or deuterium) about 100kW of RF power, which generates a high density plasma of the order of´-4 10 m 17 3 (averaged over the whole ion source volume). The expansion chamber, which is connected to the driver, has a larger volume and is magnetized; its size is approximately 57.9cm in height, width of 30.9 and 24.4cm in depth.
Most results presented and discussed in this paper have been obtained with particle-in-cell Monte Carlo collisions (PIC-MCC) methods. In section 2 we describe in details the method and discuss different points such as (1) parallelization, (2) use of scaling to make the simulations more tractable, (3) 3D versus 2D and 2.5D calculations, (4) a method for injecting power into the plasma (electron heating), (5) implementation of collisions (including physical-chemistry) and lastly (6) the modeling of the production of negative ions on the PG surface.
The negative ion source is a high density plasma in a large volume, i.e. the Debye length is much smaller than the dimensions of the source. The strong constraints on the grid spacing and time steps of a PIC-MCC simulation make it difficult to deal with the real values of the plasma density or dimensions. In several recent publications we have been using scaling of the plasma density (or, equivalently, of the vacuum permittivity) to keep the computation time within reasonable limits. Section 3 is devoted to the study of the effect of density scaling on the simulation results in a simplified problem of plasma transport across a magnetic filter where the Hall effect contributes to the cross-field transport and leads to plasma asymmetry (as in the ITER negative ion source).
We have previously published the results of a fluid model of the negative ion source where the plasma properties were analyzed as a function of power and pressure [14,15]. The plasma fluid model was coupled to a fluid description of the neutral transport and chemistry. Interesting outcomes of the simulation were the strong neutral depletion due to gas heating and ionization, and the high temperature of hydrogen atoms with respect to molecules. The Knudsen number can be close to one in the negative ion source and this may have consequences on the velocity distribution functions of neutral particles. In section 4 we look at the possible effects of the low gas density on neutral transport and on the velocity distribution of hydrogen atoms and molecules. The particle transport in this section is described by a direct-simulation-Monte-Carlo (DSMC) model.
In section 5, we provide a detailed description of the plasma transport across the magnetic filter field inside the expansion chamber of the ion source (potential profiles, electron density and temperature). We discuss the incidence of the Hall effect on the plasma dynamics (induction of a transverse asymmetry), the kinetics of negative ions and the role of positive ions in the production of negative ions on the cesiated surface of the PG. Section 6 compares the model predictions with experimental measurements. We show that the plasma asymmetry (and hence the incidence of the Hall effect) is observed in both cases. Section 7 analyses the numerical issues associated with the modeling of negative ion extraction from the PG surface. The simulation geometry is restricted to a zoom around a single aperture. PIC simulations are difficult to use with the real values of the plasma density. We present calculations obtained for plasma densities lower than in a real negative ion source and discuss the possible scaling laws that can be use to extrapolate the results to more realistic conditions.

PIC simulations
The principles of the PIC-MCCs method are described in textbooks [16,17] and in various publications. In this section we discuss some specific issues associated with the simulation of ITER-type negative ion sources, i.e. (1) need to use parallel computing, (2) difficulty to simulate the high plasma density conditions of these sources (and consequences of performing density scaling), (3) 2D versus 2.5D and 3D calculations (4) implementation of a simplified model for the external power absorption inside the ICP discharge, (5) hydrogen plasma chemistry and (6) on the method to model the generation of negative ions onto the cesiated PG surface.

General features and parallelization of PIC-MCC simulation
We have developed and used 2D and 3D parallel Cartesian electrostatic explicit PIC-MCC models [16,17] to study the plasma of negative ion sources for fusion. Due to he strong constraints on the grid spacing and time step in the high density plasma of negative ion sources, the PIC method must be optimized and parallelized. In this section we recall the principles of PIC-MCC simulations and briefly describe the technique of parallelization that is used in our models and its performance in term of computer time as a function of number of computer cores (or threads).
In an explicit algorithm, the particle trajectories are calculated based on the fields evaluated at the previous time step. The (self) electric field is derived self-consistently from the densities estimated on the grid nodes of the simulation domain. The magnetic fields, filter and suppression fields (the latter is generated by permanent magnets embedded in the first grid of the accelerator), are prescribed in this work. The time step must be a fraction of the electron plasma period and the grid size close to the electron Debye length, accordingly (both are set by the lightest of the simulated particles). The parallelization is performed in an hybrid manner using OpenMP [18] and MPI libraries. We use a particle-decomposition scheme for the particle pusher where each core (thread) have access to the whole simulation domain (as opposed to a domain-decomposition approach). The number of particles per core is nearly identical. We further implemented a sorting algorithm [19] in order to limit the access to the computer memory (RAM) and boost the execution time, Dt push , of the pusher subroutine. The latter includes electron heating (inside the ICP discharge), field interpolations, update of the velocities and positions together with the charge deposition on the grid nodes. Particles are sorted per grid cell. The field and density arrays are hence accessed sequentially. Dt push is shown in figure 1 normalized to the number of particles in the simulation. The best performance is obtained by attaching a MPI thread per socket and a number of OpenMP thread identical to the number of cores per socket. For the simulations of figure 1, we set the number of OpenMP threads to 10. We sort particles every 10 time steps without any loss of performance. The calculation is performed with a 3D PIC-MCC model and the numerical resolution is either´96 64 128 grid nodes or eight times larger with 80 particles-per-cell (ppc). The time gained in the pusher with the particle sorting is a factor ∼4. The sorting algorithm remains efficient as long as there is on average at least one particle per cell per thread. Beyond this limit Dt push converges toward the value without sorting as shown in figure 1. We define the efficiency of the pusher without sorting as, where N core is the number of cores (threads) and D ( ) t push 1 the execution time of the pusher for = N 1 core . β should be equal to 1 for a perfect parallelization of the pusher. We find b 78  % for 20 cores, 70% for 320 cores and lastly, dropping to ∼60% for 640 cores (i.e., about 23% loss in efficiency with respect to 20 cores).
Poisson's equation is solved iteratively on the grid nodes with a 3D multi-grid solver [20]. The latter is parallelized via a domain-decomposition approach. In multi-grid algorithms, a hierarchy of discretizations (i.e., grids) is implemented. A relaxation method (so-called successive-over-relaxation, SOR, in our case) is applied successively on the different grid levels (from fine to coarse grid levels and vice versa). Multigrid algorithms hence accelerate the convergence of a basic iterative method because of the fast reduction of short-wavelength errors by cycling through the different sub-grid levels. Each sub-domain (i.e., a slice of the simulated geometry) is attached to a MPI thread while the do-loops are parallelized with OpenMP (SOR, restriction and prolongation subroutines [20]). Once there is less that one node per MPI thread in the direction where the physical domain is decomposed then the numerical grid is merged between all the MPI thread. The parallelization for the coarsest grids in consequently only achieved by the OpenMP threads. This is clearly a limiting factor and more work is needed to further improve the algorithm. As an example, using a mesh of 512 3 nodes, the speedup is about ∼30 for 80 cores (b 40  %). The execution time of the Poisson solver (normalized to the number of grid nodes) versus the number of cores in the simulation is shown in figure 2.
Lastly, for the numerical resolution which we typically implement to characterize the plasma properties of the ITER-prototype ion source at BATMAN, that is,´192 128 256 grid nodes with 20 ppc, the fraction of the execution time per subroutine averaged over one time step is, ∼55% for the particle pusher, ∼8% for the Poisson solver, ∼16% for Monte-Carlo (MC) collisions, ∼4% for the sorting. The remaining time concerns both the evaluation of the electric field and the calculation of the total charge density on the grid nodes (which involve some communication between MPI threads).

Scaling of PIC-MCC simulations
The high plasma density and large volume of negative ion sources make it practically impossible to perform multi-dimensional PIC-MCC simulations for real conditions. The ratio of discharge dimension to Debye length is on the order of 10 4 (tens of centimeters versus tens of micrometers) so a simple 2D PIC-MCC simulation in real conditions would involve 10 8 grid points and more than 10 9 super-particles. This is clearly prohibitive for parametric studies with 2D simulations and impossible for 3D. The 3D simulation of negative ion extraction is also difficult although the plasma density there is smaller than in the driver and one generally consider a small simulation domain around a grid aperture.
To overcome this problem, one solution is to perform some 'scaling' i.e. to run the simulations for more tractable conditions (e.g. smaller plasma densities or smaller dimensions) and extrapolate the results to the real conditions by using some scaling laws [21][22][23][24]. The simplest scaling is the scaling on plasma density. The Figure 1. Execution time of the particle pusher (per time step) normalized to the number of macroparticles in the simulation versus the number of cores. The time is shown either with (red and gray lines) or without implementing a sorting algorithm (black-line). We use 80 particles-per-cell (ppc), a numerical resolution of´96 64 128 grid nodes (black and red lines) and´192 128 256 (gray line). The calculation is performed with a 3D PIC-MCC model on a 10 cores Intel Xeon processor E5-2680 v2 (25M cache, 2.80 GHz). There is 2 sockets per CPU, 20 cores in total. Boltzmann equation for electrons or ions is linear in the density of charged particles if the collision term is linear (i.e. when only charged particle collisions with neutrals are considered, with a given, constant velocity distribution function of neutral species), i.e. the equation is invariant if the distribution function (hence the charged particle density) f is divided by a constant α: . , 2 where a is the acceleration of the charged particles due to the Lorentz force: and I is the collision operator. In a quasineutral plasma, the electric field is deduced from current continuity and is invariant when the plasma density n is divided by a constant α. This can be easily seen on the electron and ion momentum equations (from which the total current density can be deduced). A simplified form of the electron current density is: where J e is the electron current density, P e the electron pressure and n en the electron-neutral collision frequency. In this equation, the electric field is invariant when the electron density is multiplied by a given factor. A similar argument can be made for ions. Therefore we can conclude that the properties of a quasineutral plasma are not changed when the plasma density is scaled by a constant factor (for linear collision terms) and a PIC-MCC simulation of a quasineutral plasma can be performed with scaled densities. In a real problem the plasma is bounded which implies the presence of non-neutral Debye sheaths next to the walls. The sheath properties are clearly not invariant with a scaling of the plasma density (the sheath length is generally proportional to the Debye length), but the sheath voltage and hence the plasma potential do not depend on plasma density (the plasma potential depends only of the electron temperature and electron to ion mass ratio). Therefore, when a PIC-MCC simulation is performed with scaled densities (i.e. plasma density smaller that the real density) only the wall sheath thickness is modified. If the sheath thickness is still much smaller than the discharge dimensions the scaled simulation gives an accurate description of the real problem. However, care must be taken in the following situations: • Since the sheath is larger for lower plasma densities it may become more collisional. The charged particle fluxes to the walls may be modified if the scaling is too important (ionization can also take place and be enhanced by secondary emission if present).
• In a magnetized plasma, the ratio of the charged particle gyroradius to the sheath length is also modified by the density scaling. This may also impact the charged particle transport and charged particle fluxes to the walls.
• The properties of the ion beam extracted from the plasma source are modified by the plasma density for a given extraction voltage (Child Langmuir law in the case of a collisionless sheath) and therefore the applied voltage should be scaled accordingly. This is clearly an important issue when scaling is used in a model of negative ion extraction (see section 7.5).
• In the presence of instabilities or turbulence associated with space charge separation (e.g. in a magnetized plasma) the density scaling no longer works. The instabilities and associated anomalous transport across the magnetic field do not scale linearly with the plasma density and density scaling cannot be used.
If Coulomb collisions (or other nonlinear collisions such as electron-ion recombination) are important in the real conditions, the collision module of the scaled PIC-MCC simulation can easily be modified to take into account the real collision frequencies (by using the real value of the target particle densities instead of the scaled value).
Note that in some published papers the scaling used is presented as a scaling of the vacuum permittivity instead of a scaling of the density [21][22][23]. It is easy to see (in Poisson's equation) that dividing the plasma density by a factor α is exactly equivalent to multiplying the vacuum permittivity by the same factor. In that case no scaling needs to be done for Coulomb collisions since the plasma density in the simulation is the real one.
Finally some authors have been using scaling on the discharge dimensions instead of a density scaling (see, eg, [25,26]). In a quasineutral plasma, the system formed by the Boltzmann equations coupled with the generalized Ohm's law is invariant when the dimensions are reduced by a factor β provided that the gas density (i.e. collision operator) and magnetic field are multiplied by β (the time scale is also divided by β). The Boltzmann equation in that case can be written as: . . 5 We see that this equation is invariant if the collision term is multiplied by β, i.e. if the gas density is multiplied by β (the collision term is proportional to the gas density in the absence of Coulomb collisions), and if the acceleration, is multiplied by β. The electric field E being a potential gradient is automatically multiplied by β when the dimensions are reduced by the same factor, while the external magnetic field B must be multiplied by β.
The scaling, defined above, on density, permittivity and dimensions, are all equivalent since in these three cases, the ratios l L c and r L L are the same as in the real problem while the ratio l L D is larger than in the real problem by a factor a 1 2 in the density scaling and by a factor β in the dimensions scaling. The density scaling and dimensions scaling are therefore strictly identical, for a quasineutral plasma, if b a = 1 2 . l c , r L , l D , and L are respectively the charged particle mean free paths, Larmor radii, Debye lengths, and discharge dimensions. The density is about 10 5 times lower than the real density (i.e. a = 10 5 , see previous sub-section). In 2D rectangular PIC simulations, on the other hand, one may simulate a significantly larger plasma density. The plasma in the direction perpendicular to the simulation domain is implicitly supposed to be uniform and infinite and no charged particle losses along this direction are taken into account. A solution that allows to implement approximately the charged particle losses in the direction perpendicular to a 2D simulation domain without performing a 3D calculation is the so-called 2.5D model [23,24]. For magnetized plasmas, the particle transport is simulated in the plane perpendicular to B (i.e. where the magnetized drift motion takes place). We assume that the plasma is uniform along the unsimulated direction, perpendicular to the 2D simulation plane (i.e., parallel to the magnetic field lines), and we use the following considerations to estimates the charged particle losses: • The ion dynamics in the direction perpendicular to the 2D simulation plane is not calculated but we estimate the ion losses from the Bohm fluxes to the walls. The loss frequency at a given location in the simulation plane is obtained from [27] is the local Bohm velocity, L y is the length of the ion source in the third dimension, = á ñ h n n 0.5 s  , n s is the local plasma density at the sheath edge, á ñ n the averaged density, T e (m i ), the local electron temperature (ion mass), respectively.
• The electron and negative ion trajectories are followed in the third dimension assuming that the plasma potential is flat (i.e., no electric field). When a negatively charged particle reaches a wall, it is removed if its kinetic energy along the un-simulated dimension is greater than the difference between the plasma potential and the wall, i.e., for a grounded wall. m i is the particle mass.
Macroparticles are created anywhere between   y L 0 y in the third dimension (via ionization processes). The 2.5D model estimates plasma characteristics which are averaged over L y . Note that this model is restricted to simplified magnetic field maps, where the field lines are straight in the un-simulated direction. The comparison between 2D, 2.5D and 3D models has been extensively discussed in [23].

External RF power absorption and Maxwellian heating in the driver
The ITER-type tandem reactors have an ICP discharge which couples a high RF power (typically 100 kW at 1 MHz frequency) to a hydrogen or deuterium plasma. We do not simulate directly the interaction of the RF field with the plasma but assume instead, a given absorbed power. Every time step, macroparticles which are found inside the region of RF power deposition are heated according to some artificial heating collision frequency. Electrons, being the lightest particles, are assumed to absorb all of the external power. Redistribution of energy to the heavier ions and neutrals is done through collisions (both elastic and inelastic) and the ambipolar potential. Electrons undergoing a heating collision have their velocities replaced by a new set sampled from a Maxwellian distribution with a temperature calculated from the average specie energy (inside the power deposition region) added to the absorbed energy per colliding particles, i.e., where T h (eV) is the heating temperature in electron-Volts (eV), á ñ E k h is the average electron energy, P abs (W) is the absorbed power, n h the heating frequency and N eh the number of electrons, respectively. For a given time step, colliding macro-electrons are chosen randomly where N em is the total number of macroparticles inside the heating region.
The simulated electron temperature profile is constant inside the discharge region. This is consistent with plasma conditions which are found in devices running at low pressure, low RF frequency (∼1 MHz) where electrons have typically a large thermal velocity in the driver (T 10 e eV), and a low collision frequency (a mean-free-path of the order of the discharge radius). These non-local heating conditions allow electrons, which can move freely, to deposit energy over the whole driver volume. Kinetic effects such as collisionless heating (leading to a so-called anomalous skin depth) and a non-negligible ponderomotive force (due to the high power and the low frequency of the RF antenna) are also expected inside the ICP discharge. The electron distribution function is typically non-Maxwellian in these conditions. A complete self-consistent model of the energy coupling in the ICP would be required to obtain a better estimation of the electron distribution in the driver. A non-self-consistent but much simpler and useful approach would be to study the influence of a non-Maxwellian distribution on the plasma parameters by imposing in the driver electron distributions that are deduced from experiments.

Implementation of collisions in a particle model-MC and DSMC methods
In a PIC-MCC algorithm, the Boltzmann equation, is solved numerically in two steps [28,29].
T r the total differential cross-section (summed over all the collision processes between the incident and the target particles) and, lastly, Ω the solid angle. Primes denote the distribution function after the collision. For small time steps, equation (7) may be rewritten as, i is known explicitly from the previous time step. This finite-difference analog of equation (7) is second order correct in Dt . The operators D and I are, . , 1 0 Applying the operator ( -DtD 1 ) on the distribution function f i is equivalent to solving the Vlasov equation, The PIC procedure [16,17] is a characteristic solution of equation (11). Once the particle trajectories have been updated, then the second operator ( -DtI 1 ) may be applied on the (updated) distribution function. A macroparticle is equivalent to a Dirac delta function in position-velocity space (Eulerian representation of a point particle) and hence a probability may be derived from equation (9) for each collision processes [28,29]. The probability for an incident particle to undergo an elastic or inelastic collision with a target particle during a time step Dt is with N c corresponding to the total number of reactions for the incident specie, n c the density of the target specie associated with a given collision index and = - is artificially set to its maximum value and hence ( ) P i max is greater than the real probability and is constant over the entire simulation domain. There is consequently a probability, that a particle undergoes a fake collision (dubbed 'null' collision), which will be discarded.
The total number of incident particles which will hence collide during a time step Dt (including a 'null' collision) is, where N i is the number of incident macroparticles in the simulation. N i must be replaced by i for collisions with another particle of the same specie [28]. ( ) P i max is equiprobable for any pairs of incident-target particles and consequently the latter may be chosen randomly inside the simulation domain. In the model, one checks first if the incident macroparticle experienced a real collision, where r is a random number between 0 and 1. The probabilities P c for each reactions (whose total number is N c for a given incident specie) are ordered from the smallest to the largest and a reaction k occurred if, Once a collision type is selected then the macroparticles (both incident and target) are scattered away in the center-of-mass (CM) frame (see next section). In the model, neutrals are either considered as a non-moving background specie with a given density profile or are actually implemented as macroparticles and their trajectories integrated. In the case of the former, collisions between charged particles and neutrals are performed by the so-called MC method while for the latter, actual particle-particle collisions are evaluated using a DSMC algorithm [30]. Both are similar except that in the MC method, one artificially extract a neutral particle velocity from a Maxwellian distribution function. Collisions between charged particles are always performed by a DSMC algorithm in the model. Collisions (both elastic and inelastic), are implemented assuming that particles (incident, target or newly created) are scattered isotropically in the CM. Energy and momentum is conserved and we assume for simplicity that each byproduct partner after the collision have identical momentum in the CM frame. This implies that the lightest particles will equally share most of the available energy. For further details, please refer to [21].

Physical chemistry of charged particles
We describe below the most complete version of the plasma chemistry module embedded in our PIC-MCC model (we often use a simplified sub-set of this module, depending on the purpose of the model). In this model, the plasma consists of electrons, molecular hydrogen (background) gas H 2 , hydrogen atoms H, molecular ions + H 2 and + H 3 , protons and lastly negative ions -H . Collisions between electrons, ions and neutrals are considered; the set of reactions is presented in tables 1 and 2 (66 collision processes in total) and is very similar to the one used by previous authors [15,31]. Table 1 corresponds to the collision processes associated with electrons. Reactions #2, 6, 7, 8 and 14 combine multiple inelastic processes included in the model in order to correctly account for the electron energy loss. Reaction #2 regroups the excitation of the hydrogen atom from the ground state to the electronic level =n 2 5 [32]. Reaction #7 combines the ground state excitation of the hydrogen molecule n S = [32], Rydberg states [34] and lastly rotational levels J=2 [35,36] and 3 [37,38]. Reaction #17 models in a simple manner the generation of negative ions in the ion source volume, which are a byproduct of the dissociative impact between an electron and molecular hydrogen n . The concentration of excited species is not calculated self-consistently in the model. To estimate the volume production of negative ions, we assume that 1% of H 2 molecules are excited in vibrational levels with n  4. This is in accordance with the H 2 vibrational distribution function calculated either with a 0D model [39] or a 3D particle tracking code [40]. Table 2 summarizes the collision processes of heavy ions with neutrals. Reaction #9 corresponds to the excitation of the hydrogen molecule from the ground state to vibrationally excited levels n¢ = -1 2 [41,42] and to the rotational levels = -J 2 3 [43]. To our knowledge there is no reliable data available for the elastic collision between + H 3 and neutral atoms (reaction #2), we consequently use the same crosssection as in reaction #1.

Physical chemistry of neutrals
Cross-sections for collisions between neutrals inside the ion source volume, which are summarized in table 2 (reactions #18-20), as well as backscattering, dissociation or recombination probabilities against the ion source walls are required for the modeling of the neutral particle dynamics (and the associated neutral depletion). Table 3 shows the surface processes and corresponding coefficients. In a low-pressure plasma device such as the one used for fusion applications (ITER or DEMO for instance. DEMO is a concept for the next generation of Tokamaks), the plasma-wall processes have a strong impact on the source characteristics. Low-temperature backscattered molecular hydrogen is assumed to be in thermal equilibrium with the wall. An average  backscattered energy is considered for fast atoms and ions, i.e. a thermal accommodation coefficient g < 1 (g = 1 corresponds to the wall temperature). These estimates are based on Monte Carlo calculations from the code TRIM [44]. Average reflection probability is also taken from the same database. Furthermore, we assume that atoms which are not backscattered will recombine. The interaction of + H 3 and + H 2 ions with the walls and the corresponding coefficients are not well known. The coefficients used in the simulations are reported in table 3. For + H 2 we use coefficients that are consistent with the measurements of [45]. For + H 3 we assume guessed values (the + H 3 flux to the walls is relatively small with respect to the + H 2 and + H , and the results are not very sensitive to these coefficients).

Negative ions
Negative ions are produced on the cesiated PG surface as a byproduct of the impact of hydrogen atoms and positive ions. Our PIC-MCC model has not been coupled to a neutral transport module so the negative ion flux emitted by the surface is not obtained self-consistently from a flux of neutral atoms deduced from the simulation. The magnitude of the negative ion current density due to the impact of H atoms on the cesiated PG is either derived from plasma parameters measured experimentally or from DSMC calculations. Assuming a Maxwellian velocity distribution of the hydrogen atoms, the flux of these atoms on the PG is: where n H is the atomic hydrogen density, m H the mass and e the electronic charge. The negative ion current is deduced from, [46], which was not obtained in a plasma (the experiment produced hydrogen from thermal dissociation in a tungsten oven) and consequently remains approximate for the ITER-type ion sources. For typical BATMAN working conditions, we find [47][48][49]. Negative ions are generated on the PG assuming a Maxwellian flux distribution function with a temperature = T 1eV n in the model. Furthermore, the surface production of negative ions resulting from positive ion impacts is calculated self consistently. For each ion impinging the PG, the yield is evaluated assuming a molecular ion may be considered as an ensemble of protons sharing the incident ion kinetic energy (a + H 3 ion for instance would be equivalent to three protons each with an energy = ). Each of these 'protons' may produce a negative ion. The condition  r Y must be fulfilled for the negative ion to be generated with r a random number between 0 and 1. The yield is taken from Seidl et al [46] for Mo/Cs surface with dynamic cesiation. The negative ions are scattered isotropically toward the ion source volume with a kinetic energy assumed to be = There is experimental evidence that negative ions may capture a large amount of the incident positive ion energy [50]. In addition, for clean metallic surfaces (tungsten) the reflected atomic hydrogen particle energy is numerically evaluated to be around 65% of the impact energy at normal incidence and for = E 1 eV k [51]. Lastly, it has been reported in the experiments that the extracted negative ion current increases only slightly with cesium when the PG is water-cooled [52] while a PG heated to a temperature of~ 100 C-250 • C induces a significant increase of the negative ion current, by a factor~-4 5 in the experimental conditions of [5,52] (the other walls of the ion-source were water-cooled). In the model, we consequently assume that negative ions may only be produced on the cesiated PG surface.

Simulation domain
The simulation domain for the 3D PIC-MCC modeling of the BATMAN device is shown in figures 3(a) and (b). The magnetic filter field is generated in the model by permanent magnet bars which are located on the lateral side of the ion source walls close to the PG. The field is calculated by a third-party code [53]. In 2.5D, solely the XZ plane is considered ( figure 3(b)) but with a higher numerical resolution (or similarly plasma density) than in 3D. We assumed = L 32 cm y for both the driver and the expansion chamber and we implemented a Gaussian profile for the magnetic filter (i.e., mirror effects are neglected [23]), with an amplitude B 0 , width L m and a maximum located at x 0 . The magnetic field generated by the permanent magnets has a shape very similar to a Gaussian profile on the ion source axis [23], as shown in figure 4. Figure 3(c) shows the simulation domain for higher numerical resolution 2D and 3D PIC-MCC modeling of negative ion extraction from the PG surface. B D corresponds to the deflection magnetic field from permanent magnets embedded into the extraction grid (EG). A domain restricted to the vicinity of the PG allows the implementation of plasma densities closer to the real one. This will be discussed in section 7.

Density scaling and Hall effect
As said above, the high plasma density and large volume of negative ion sources for fusion make it practically impossible to run full scale PIC simulations in 3D (and very difficult in 2D) and scaling the plasma density to lower values or the permittivity to higher values allow to perform computationally tractable simulations. As discussed in section 2.1.2 the scaling provides an approximate solution of the problem, the validity of the approximation depending for example on the size and role of the wall sheath in the considered problem.
In this section we provide a quantitative description of the effects of the density scaling in a simplified plasma source with magnetic filter. This example allows us to discuss both the influence of the density scaling on the results and the physics of the Hall effect induced by the presence of the magnetic filter and which contributes to non-collisional charged particle transport across the filter, and to the development of an asymmetry in the plasma properties.   , where G = -e J e e is the electron current density). The force is directed toward the bottom surface of the ion source for the filter configuration schematically shown in figure 3(b). The presence of walls induces a charge separation (polarization) and the creation of an average electric field that opposes the effect of the Lorentz force, as in the Hall effect. This Hall electric field (which is consequently downward-directed), E H , generates in turn anÉ B H drift along the X-axis which significantly increases electron transport across the magnetic filter with respect to an ideal 1D filter without transverse walls [54,55]. We therefore expect that the Hall effect will create a plasma asymmetry with an electric potential and a plasma density higher in the top of the chamber (large Z) than in the bottom (small Z). The electron flux in equation (20) may be expressed as follows [21], The electron motion is consequently dominated by the magnetic drift which is composed of a diamagnetic term Ṕ B e (collective effects), and anÉ B term [54,56]. The electric field is a combination of the Hall and the ambipolar fields.
The Hall effect in low temperature plasmas and its impact on plasma asymmetry have been studied analytically in the simple conditions of a positive column [57,58]. The situation is more complicated in the magnetic filter of the negative ion source because of the non-uniform magnetic field and of the presence of axial plasma density gradients. Furthermore, the general features of the Hall effect (i.e., production of a voltage difference across an electrical conductor, perpendicular to both the direction of the electric current in the conductor and the applied magnetic field) have been clearly observed in other magnetized plasma sources with particle transport properties comparable to those of the ITER prototype ion sources. Experimental measurements have been recently performed in a low power inductively coupled plasma with a magnetic filter and have shown the presence of a strong asymmetry in the collected current density [59]. Note finally that the Hall effect is not present in devices such as Hall thrusters where the electron drift perpendicular to the discharge current is closed and is not impeded by the presence of walls (closed-drift devices).

Transport across a filter in a simplified geometry and influence of density scaling
In order to illustrate the Hall effect in magnetized plasmas in a simplified manner, and, at the same time study the influence of density scaling on the results, we implemented a 2D simulation domain which is a square box of dimensions20 20 cm 2 . The model is a 2D PIC-MCC and there is no particle losses in the plane perpendicular to the simulation domain. We model the XZ plane and the magnetic filter field is along (OY) as in figure 3(b). The magnetic field profile is given by equation . We consider only electrons and + H 2 ions as particle species composing the plasma and therefore we use a subset of the physical-chemistry described in tables 1 and 2. Instead of assuming that an external power is absorbed by the plasma, as described in section 2.1.4, we keep the plasma density constant by re-injecting an electron-ion pair each time a positive ion is lost on the external boundaries of the simulation domain. The latter are absorbing surfaces. The particle re-injection is set inside a magnetic field free region between x = 1.5 and 4.5cm. Furthermore, the electron temperature is maintained constant in that area with = T 10 eV e . The scope is to draw an electron current (flux from left to right) through the magnetic filter and evaluate the Hall effect. For that purpose we assume that there is no ionization processes and hence reaction #5 in table 1 is artificially replaced by an inelastic collision (excitation). The density profile of molecular hydrogen is constant with =´n 5 10 m The profiles are very similar except that the electron flux channels closer to the walls in the higher density case. This is due to the transverse shape of the plasma potential. The size of the Debye sheath is smaller and hence the pre-sheath extends closer to the boundaries. This also indicates that the electron motion across the magnetic filter field occurs mainly in the pre-sheath. This is confirmed by quasineutral fluid calculations. The maximum value of the Hall parameter is h 40 max  in the model and the electron flux is hence well described by equation (23). The electron flux is a combination of a diamagnetic drift, which is a consequence of the particle random motion (i.e., the velocity spread) expressed mathematically in the pressure term and anÉ B drift. The two terms are often of opposite sign, i.e., cancelling each others. The electric field is itself a combination of the Hall (which is downward directed) and ambipolar fields as demonstrated in section 3. In the regions (1) and (2) highlighted in figure 5(a) , i.e., the drift is ofÉ B type, respectively. The electron current density profile depends on the shape and magnitude of the magnetic filter field but the general features described in this section are reproduced in any type of magnetized plasma sources where a current is drawn across the magnetic field (biasing the rhs electrode hence enhance the Hall effect). Figure 5 shows that the plasma density has a small influence on the plasma properties in the considered range (the ratio of the plasma density between the two simulations is equal to 64). Charged particle transport in the plasma occurs mainly inside the quasi-neutral region driven by the pressure gradient, the ambipolar and Hall electric fields, and the effect of the sheath on the electron current density distribution is negligible on the figure. The sheath should be about 8 times larger in the lower plasma density case of figure 5(a). Note that the presence of the magnetic filter can also induce plasma instabilities seeded by charge separation. This happens, in the simulations of the simplified problem considered here, when the length of the transverse direction is significantly increased. In that case the Hall effect is no longer sufficient to ensure cross-field transport and instabilities develop, enhancing electron transport across the filter. The asymmetry of the plasma density resulting from the Hall effect can be seen in figure 5(c). Figure 6 which shows the transverse plasma potential profile versus the average plasma density. The latter is increased from á ñ =´- 3 . The ratio of the densities between the two extreme cases is a = 256. The variations between the potential profiles in figure 6 lie essentially on the size of the Debye sheath. The amplitude of the potential in the quasi-neutral region is similar within ∼10%. The Hall electric field E H is about 15V m −1 (measured between the top and bottom plasma sheath edges). Figure 7 shows the electron current collected on the biased electrode (rhs of the simulation domain) versus the plasma density. The current increases linearly with the plasma density as expected. One

Neutral transport and plasma properties versus power and pressure
In a previous work we studied the properties of the plasma of the negative ion source versus power and pressure based, on a quasineutral fluid description of the plasma, coupled with a Navier-Stokes model of neutral transport [14,15]. It was shown that the neutral density was strongly depleted due to gas heating and ionization and that the temperature of atomic hydrogen was much larger than that of molecular hydrogen. In this section we use the same plasma model but we couple it with a kinetic description of neutral transport based on a DSMC method. The objective is to estimate the consequences of the fact that the gas flow is rarefied (Knudsen number not small with respect to 1) on the model results and on the velocity distribution of hydrogen atoms and molecules.

Modeling of the neutral transport and chemistry in the ITER prototype ion source
In the ITER prototype source BATMAN, the typically working conditions correspond to a low background gas pressure (molecular hydrogen or deuterium) of ∼0.3Pa, together with a high RF power (coupled to the plasma by an external antenna), ∼100kW. Such conditions depletes the neutrals in the experiments [47]. This effect was first described through modeling [15] and then confirmed by experiments [47]. We model neutral depletion by coupling a DSMC algorithm for the neutrals (both molecular and atomic hydrogen in our case) with the 2D implicit fluid model described in details in [14]. In this self-consistent model, the plasma parameters adjust so that charged particle and power balance are satisfied at steady state (e.g. volume ionization is compensated by surface losses). The geometric parameter involved in the charged particle and power balance is the chamber volume over surface ratio therefore the dimensions of the 2D geometry are rescaled accordingly in the model of BATMAN [23]. The dimensions of the ion source in the model is a driver of length 9cm, height 8cm and an expansion chamber of16 16cm 2 . The flow rate for the molecular hydrogen gas injected into the ion source volume is adjusted,  ms in BATMAN). Figure 8 shows the atomic hydrogen energy distribution function in the center of the negative ion source for a 60kW absorbed power and 0.3Pa background gas pressure. The distribution function is highly non-Maxwellian and its properties are mostly controlled by the production and collisions of H atoms against the walls of the ion source (the total mean-free-path is of the order of 1 m, i.e., significantly larger than the dimensions of the device). Reactions which either generate (so-called source term) or remove (loss term) hydrogen atoms from the ion source volume are summarized in table 4. Particle and power (gain or loss) are shown as a percentage of the total. H atoms are mostly created and heated by the wall recombination of protons and molecular ions on the ion source walls (reaction #1 with 45  % of the particle production and 77  % of the energy gain) and by the volume dissociation of H 2 (reaction #2). + H x ions (where = x 1-3) are mainly generated inside the discharge and are accelerated by the plasma potential toward the walls of the ion source. The amplitude of the potential in the driver is about 50 V for 60kW of RF power at 0.3Pa in the experiments [11]   and hence + H x ions impact the walls with a high energy. This explains the origin of the large energy tail in the distribution function of atomic hydrogen as shown in figure 8. Lastly, H atoms loose most of their energy through collisions with the source walls (∼95%, reactions #6 and 7). The H and H 2 temperatures are strongly dependent on the assumptions that are made on the wall reactions (accommodation coefficients provided in table 3).
The energy distribution function for molecular hydrogen is shown in figure 9. H 2 molecules are created uniquely through the recombination of + H x ions and H atoms on the walls. Molecular hydrogen is emitted from the surfaces as a Maxwellian flux at T w =300 K (where T w is the temperature of the surface). The energy distribution function is well fitted by a Maxwellian (up to about T 6 H 2 ). The mean-free-path is ∼10cm, i.e., smaller than the dimension of the ion source. The energy tail is induced by the collisions with the warm H atoms (T 1 eV H  for 60kW and 0.3Pa while T 0.08 H 2  eV). Molecular hydrogen is mainly heated through elastic collisions with atoms (∼65% of power gain, reaction #1) and by electrons (reaction #2) as shown in table 5.
The calculations have been performed either with or without a magnetic filter field in the expansion chamber. The magnetized case corresponds to a maximum field amplitude of = B 15G max close to the PG, which is lower than the field in the actual experiment (∼75 G on axis) but nevertheless m = h B 1 B e  and the electrons are fully magnetized (a smaller magnetic field is used in the simulation because of numerical issues at the time with the fluid model for large magnetic fields). The indirect effect of the magnetic field on the neutral dynamics is that the depletion of H 2 occurs in the area where the electron density is highest (i.e. in the driver when the expansion chamber is magnetized) because molecular hydrogen is dissociated or ionized mainly by electrons (table 5). The density profile of hydrogen atoms is on the contrary quite insensitive to the magnetic field due to the fact that the volume losses (ionization) are significantly smaller than for H 2 (i.e., ∼14% of the total losses, see table 4) and that the mean free path, ∼1m, greatly exceeds the ion source dimensions. The 2D density and temperature profiles for the H atoms are shown in figure 10. Figure 11 displays the electron density and temperature averaged over the negative ion source volume versus the absorbed RF power in the discharge and the background gas pressure. For a pressure of 0.75Pa, the electron temperature is almost independent of the external power while the electron density increases quasi-linearly with power. This behavior may be  explained by a global model. Assuming steady-state conditions, the total amount of particles created through ionization in the ion source volume at a given time is equal to the number of particles lost on the device walls, where n s is the plasma density at the sheath edge, á ñ n is the average plasma density (n e = n i is assumed), is the Bohm velocity, m i is the mass of + H 2 , ( ) k T i e is the ionization rate (which is a function of the  electron temperature), n g is the background gas density and S (V ) is the ion source wall surface area (volume), respectively. Equation (24) is derived from the volume integration of the electron continuity equation. The term on the left-hand side (lhs) corresponds to the electron flux impacting the ion source walls and the rhs term is the volume integrated ionization rate. If the neutral gas density n g is fixed then equation (24) provides an average estimate for the electron temperature that is independent of power. The electron temperature is linked to the ion source geometry through the surface-to-volume ratio S/V. An estimate for the average plasma density may be deduced from the power balance equation, that is, the power absorbed (P abs ) in the ion source volume is equal to the power lost on the walls, where e ( ) T T e is the average energy lost per electron-ion pair lost on the walls [15,21,27]. For a given value of the background gas density n g , the electron density hence increases linearly with the absorbed power. The gas density does vary with power as shown in figure 12 because of increased gas temperature and ionization rate as the absorbed power is increased. For instance at 0.75Pa and 90kW the gas density is depleted by ∼65% compared to the density without discharge.
The charged particle balance in equation (24) provides a relation between the gas density and electron temperature, . This relation can be calculated from the hydrogen ionization cross-section assuming a Maxwellian electron velocity distribution function, and the corresponding curve is shown on figure 13 (dark solid line). The results from the 2D model are also plotted (symbols) on the same figure. The symbols correspond to different gas pressures before the plasma is on, and, for the same symbol, the set of points correspond to several values of the discharge power (leading to distinct values of the averaged gas density due to the depletion associated with gas heating and ionization). We see on this figure that the 2D model results are in excellent agreement with the simple 0D model. High electron temperatures correspond to low gas densities and the asymptotic behavior of the curve at high electron temperatures shows that the discharge cannot be sustained below a given gas density (low gas pressure and high power). The good agreement between the simple 0D model (taking into account only + H 2 ions) and the more complex 2D model with plasma chemistry may seem surprising but is due to the fact that ionization of molecular hydrogen is on the average significantly larger than ionization of atomic hydrogen. Also, in this plot, the neutral depletion is a parameter in the 0D model while it is self-consistently calculated in the 2D model.
The gas density variations with power are displayed in figures 12(b) and 14(b) for different pressures (the plotted densities are averaged over the whole ion source volume). The atomic hydrogen density increases with power because the dissociation rate increases with increasing plasma density and electron temperature while the H 2 density decreases with power because of gas heating and dissociation. However, the average H density reaches a limit when the power increases, because of the increase in the hydrogen atom temperature. The variations with power and pressure of the volume averaged molecular and atomic hydrogen temperatures are shown in figures 12(a) and 14(a). The H temperature (in the 0.1-1.2 eV range) is much larger than the H 2 temperature (on the order of 0.02-0.08 eV). Both temperatures increase continuously with power. It is interesting to note that, for a given input power, the H temperature increases with decreasing pressure while the opposite is true for H 2 . This is due to the larger energy exchange rate (between H and H 2 ) at higher pressure. The large H temperature is due (1) to the fact that H atoms are generated with a large energy during electron impact dissociation of H 2 and (2) to the generation of fast atoms at the walls resulting from the recombination of positive ions.
The wall accommodation coefficients are unknown experimentally. The model may be used to evaluate its influence on the plasma properties. Assuming that the neutral hydrogen particles are backscattered off the walls with the same temperature as the surface (accommodation of 1) for the reactions # 2, 4, 6 and 8 of table 3 instead of g = 0.5, we find that (1) the amplitudes of the molecular hydrogen temperature and density are only slightly modified while (2) the H atom temperature is on average divided by 5 and the density as a consequence is larger by a factor of 2. The relationship between the accommodation coefficient and the neutral atom temperature is almost linear. The calculation was performed for a background gas pressure of 0.3Pa and an absorbed power of 60kW. There is experimental evidence that the temperature of H atoms is significantly larger than that of H 2 but the strong dependence of the H temperature on the accommodation coefficients demonstrated by the simulations shows that systematic comparisons of experiments and model results on the H temperature would be useful to get a better estimation of the accommodation coefficient.   [46] have calculated the averaged conversion yield á ñ( ) Y T H for a Mo/Cs surface with dynamic cesiation assuming that the distribution function of the atomic specie was Maxwellian. The experiments were not performed in a plasma. In our case, the distribution function is highly non-Maxwellian and hence a range for the yield may be estimated either by (1) fitting the energy distribution function of the atoms with a Maxwellian distribution, in the vicinity of the PG in the model (figure 10) or (2) considering the atoms as a 'beam', is assumed [46] (deduced from the experiments). E in is the incident energy of the atom. Using equation (26), we find á ñ Y 11.5  % and hence . j G is the negative ion current density generated on the PG surface.

Plasma transport across the magnetic filter
In this section, we model the plasma properties of the (one driver) ITER-prototype negative ion source at BATMAN with a 2.5D PIC-MCC algorithm (section 2.1.3) using a density scaling of a = 400 (plasma density 400 smaller than in the real source). We describe not only the plasma transport across the filter, but also the extraction of negative ions, considering a limited number of grid apertures (slits). Although care must be taken in the interpretation of the results on negative ion extraction with such a large density scaling, we think that these results are interesting because plasma source and negative ion extraction are described in the same model (extraction models are generally performed on a simulation domain around a single aperture, and the correctness of the boundary condition inside the simulated plasma region is currently an issue [60]).
The magnetic filter field profile is Gaussian following equation ( . We model the XZ plane, as shown in figure 3(b) and the magnetic field is directed along (OY). The length of the third, un-simulated dimension, is = L 32 cm y (also for the discharge). The numerical resolution iś 1024 1536 grid nodes with ∼100ppc and we model a lower plasma density, that is, á ñ =´n 7.5 10 m p 14 3 . We simulate 7 slit apertures in the PG, each with a diameter of 1.5cm and length = L 32 cm y . The deflection magnetic field B D is calculated with a third-party code [53]. Lastly, we consider an absorbed power of 60kW and a background gas pressure of 0.3Pa. The external RF power is coupled to the plasma in the model by artificially heating macroparticles in the driver region following the method described in section 2.1.4 (n =  and T 1 eV H  ( figure 14), respectively. The latter are consistent with experimental observations [47]. Negative ions are produced inside the ion source volume (reaction #17 of table 1) and on the cesiated PG surface either as a byproduct of positive ion impacts or atomic hydrogen. We assume a negative ion current density of a = j 600 A m G 2 generated by H atoms and an ion temperature of = T 1eV n . Figure 15 shows the electron density in (a) and (b), electron temperature (c) and plasma potential profiles (d) in the plane perpendicular to the magnetic field lines (XZ plane) of the BATMAN negative ion source. The PG bias voltage is 25 V in (a), 20 V in (b)-(d). The Hall effect (section 3) generates an electric field directed downward and hence the transverse potential profile is asymmetric. The plasma density is highest in the vicinity of the top wall of the ion source, i.e., where the plasma potential is maximum (the maxima is shifted up compared to a situation without a magnetic field). The extent of the plasma asymmetry is strongly related to the PG voltage [12,24] (biased positively with respect to the ion source walls). A larger electron current is drawn from the driver through the filter when the bias is increased and hence the Hall electric field is strengthened. The temperature profile is also asymmetric (oblique isothermals) due to the magnetized electron drift dynamics which evolve into an oblique electron flux across the magnetic filter. T 10eV e  in the discharge region and ∼1.5eV in the vicinity of the PG. Lastly, the densities of the plasma particle species averaged over the ion source volume are a á ñ´-  Figure 16 shows the negative ion flux (a) and density (b) profiles in the XZ plane of BATMAN. The PG bias voltage is 20 V. The flux corresponds solely to negative ions produced on the PG by neutral atom impacts. The numerical resolution is1024 1536 nodes with ∼35ppc. The scaling factor is a = 400. Streamlines are displayed in figure 16(a) to indicate the direction of the negative ion flux. Except close to the PG, the flux is directed toward the ion source volume. Negative ions which are extracted originate consequently from the PG surface surrounding the apertures. Negative ions are somewhat magnetized, which significantly enhance the skewness of the density profile displayed in figure 16(b) for the elastic collisions with neutrals (reactions #12 and #13 of table 2). The negative ion density averaged over a line-of-sight (LOS) parallel to the PG (from the top to the bottom wall) and 2cm from the latter is a á ñ´n 6.5 10 m n 16 3  which is of the order of experimental measurements [12]. At the same distance from the PG, the negative ion density  for a Maxwellian flux distribution function [49]. Note that the probability for a negative ion produced on the PG to be extracted from the ion source is on the order of 50% for a PG bias voltage of 20 V.

Negative ion dynamics
The asymmetry in the plasma parameters has important consequences on the (extracted) negative ion beamlet profiles. The beamlet currents may deviate by as much as 50% from the average current in the model. The ITER accelerator has a ±10% acceptance and this may likely translate into some beam interception on the accelerator grids. Figure 17 shows the electron and negative ion current extracted from the 7 slit apertures of diameter 1.5cm versus the PG voltage. The EG potential is set to ¢ = V 210 EG V. This value is obtained by assuming that the extracted currents scale with the Child-Langmuir law a m = j V n EG 3 2 in order to estimate the EG voltage for a plasma density a = 400 times smaller than the value encountered in the actual ITER prototype ion source at BATMAN. This approximation will be further discussed in section 7. μ is the perveance and V EG is the potential for a = 1. We find a ¢ = V V EG EG 2 3 assuming = V 11.4 kV EG on the EG, which is located 9mm from the back of the PG in the model (this corresponds to a potential V 7.5 kV EG  for the 6 mm gap of the ITER accelerator). The other parameters of the simulation are identical to the ones of the preceding section.

Electron and negative ion extraction versus the PG bias voltage
The positive ion flux on the PG decreases with a larger bias voltage [22] which has two consequences: (1) the negative ion current produced on the PG by positive ion impacts decreases as well and (2) a smaller positive ion V, the amplitude of the plasma potential is gradually getting lower than the applied bias voltage everywhere parallel to the PG. Negative ions are hence increasingly trapped near the PG surface and the extracted current is dropping as shown in figure 17. The potential in the presheath 1cm from the PG is 1.2 V below the electrode voltage for figure 18. This behavior is confirmed by experimental measurements [12]. Lastly, the gap between the plasma potential at the edge of the Debye sheath in front of the PG and the bias voltage (V PG ) decreases with a larger value of the bias ( figure 18). The Electron temperature in the extraction region is T 1.5 eV e  in the model and consequently more electrons may cross the sheath barrier (and be collected on the PG surface). This explains the continuous drop of the coextracted electron current versus the PG bias voltage shown in figure 17.

Summary
Although they have been obtained with a large scaling factor, the results presented in this section raise some questions about the plasma asymmetry in the source and its possible consequences on the non-uniform distribution of the extracted negative ion current along the PG surface. They show that drawing an electron current across the magnetic filter in a fusion-type negative ion source induces a plasma asymmetry due to the Hall effect (the plasma asymmetry is confirmed by experiments as described in the next section). As a consequence, the extracted negative ion beamlet current density is asymmetric as well in the model, with values exceeding the ±10% current spread foreseen for the ITER NBI electrostatic accelerator. In addition, we found that the PG bias voltage (1) enhances the plasma asymmetry (figures 15 and 16), (2) increases the extracted negative ion current up to (approximately) a floating PG and lastly, (3), induces a monotonous decrease of the co-extracted electron current (figure 17). The extracted particle currents and the plasma potential profiles were derived from a fixed value of the plasma density (scaling factor a = 400). Only the PG bias potential was varied  in the simulations. The profiles strongly resemble experimental measurements [12,61] which seems to indicate that the conclusions drawn from this section may be extrapolated to higher plasma densities.

Comparison with experiments
In this section, we compare the model to experimental measurements. We simulate the conditions reported in Schieko et al [62]. The model is a 3D PIC-MCC algorithm with a numerical resolution of´128 96 192 grid nodes, a scaling factor a =5 10 4 and 20ppc. The magnetic filter field is generated by permanent magnets, positioned against the lateral wall of the BATMAN prototype source and 9cm for the PG. The simulation domain is displayed in figures 3(a) and (b). The discharge is approximated as a rectangular box in the model instead of a cylinder. The filter field map is calculated by a third party software [53]. The background hydrogen gas pressure in the experiment was 0.6Pa and the RF power = P 40 RF kW. The properties of the neutrals are unknown experimentally and we implemented the values derived from the DSMC model ( figures 12 and 14 (this is about 12 times smaller than for a pressure of 0.3 Pa). The latter corresponds to an optimal cesiation of the PG. Lastly, we assumed that T T n H  . Figure 19 shows the axial profile for the electron density, temperature and plasma potential for two Langmuir probe positions along (Oz), that is, =z 10 cm (bottom) and = z 10 cm (top). Both probes are positioned at =y 5 cm. The experimental data [62,63] are plotted in figure 19 for comparison. Both the experiment and the 3D PIC-MCC model exhibit similar features. The magnetic filter field generates an asymmetry in the plasma parameters (the gap in the plasma potential is the hallmark of the Hall effect). The main discrepancy between the experiments and the model comes from the external power (we assumed 15 kW of absorbed power versus 40 kW of RF power in the experiments). The peak for the electron density near the exit of the driver is also more pronounced in the Figure 19. Axial electron density, temperature and plasma potential profiles for two Langmuir probe positions along (Oz), that is, =z 10 cm (bottom) and z=10cm (top). Both probes are positioned at =y 5 cm. The experimental data of Schiesko et al [62,63] are indicated by the symbols (square and triangles). The experimental conditions correspond to a background hydrogen gas pressure of 0.6Pa, a RF power of 40kW and a magnetic filter field generated by permanent magnets on the side walls of the ion source, 9cm from the PG.
model. This may be due to the oversimplifying assumption of implementing a rectangular geometry for the discharge. The scaled electron Debye length is l 6.5 mm De  on average in the calculation. The fraction of the volume occupied by the plasma sheath isL L

Extraction of negative ions
In this section, we simulate the extraction of negative ions from a fusion-type ion source using a model which is restricted to a single aperture [60]. The simulation domain is shown in figure 3(c).

Numerical issues
Recently published results from PIC-MCC models [64][65][66] have led to a counter-intuitive and unexpected description of negative ion extraction. Using chamfered apertures in the simulations, the models [64,66] show that only those negative ions emitted from the tip of the chamfered aperture can be extracted, which is rather surprising and does not seem to correspond to a proper operation of the extraction system. A very small negative ion current is emitted from the rest of the grid surface due to space charge saturation associated with very large values of the potential drop in front of the emitting surface (i.e., a virtual cathode). The numerical grid spacing used in the simulations was much larger than the Debye length (typically by a factor between 5 and 10). Experiments, on the other hand, have shown that extraction of a negative ion beam from a plasma electrode with a flat surface around the aperture is actually possible [5,67,68]. We show in section 7.3 that a grid spacing smaller than the Debye length is required for a proper description of the plasma in the vicinity of the PG (including the shape of the virtual cathode) [69]. We describe in section 7.4 the properties of the plasma meniscus predicted by the model. Lastly, we discuss the use of scaling laws in a 2D model of negative ion extraction (slit apertures) in sections 7.4 and 7.5, i.e., we compare the real plasma density (a = 1) to lower densities (a > 1) and analyze the correlations. An extrapolation to circular apertures (3D PIC-MCC modeling) is also discussed.

Model
The simulation domain is described in figure 3(c). The model is a zoom around a single aperture of dimensionś 3.2 1.6 cm 2 in 2D. The aperture (slit) is not chamfered and with a diameter of 8mm as in the prototype ion source at BATMAN described in [5]. The top and bottom boundaries are periodic while all the others are of Dirichlet type. The PG is set at a given reference potential (0 V in our case) and the left-hand side (lhs) boundary voltage V LB is adjusted in order to simulate the effect of a bias potential. The plasma is numerically sustained by re-injecting an electron-positive ion pair for each positive ions lost on the walls of the simulation domain.
Particles are re-injected on the lhs of the domain between = x 2 and 4cm and the electron temperature is maintained by replacing the macroparticle velocity by a new one sampled from a Maxwellian distribution at a preset temperature = T 2 eV e . Negative ions are uniquely produced on the PG surface (we specifically want to assess the dynamics of negative ions which are produced on the electrode). The physical chemistry is simplified. We only consider negative hydrogen ions, + H 2 and electrons (reactions #14 and #15 of table 2 have a negligible contribution). We implemented a background gas density =´n 4 10 m is the position of the EG, = d 1.6 cm is the distance between magnet bars and lastly = B 600G 0 for the BATMAN configuration. Figure 20 shows the transverse virtual cathode potential profile (at its minimum), parallel to the PG, versus the numerical grid spacing (D = D , numerical heating increases the electron temperature which in turn modifies the plasma parameters. The virtual cathode depth drops in the model and the negative ion saturation current escaping the PG is hence also significantly reduced.   when the EG potential is increased and below a curvature radius  r a c (where = a 2 8mmis the diameter of the aperture), aberrations appears in the extracted negative ion beam.

Scaling laws
Both the extraction (EG) voltage and the plasma density affect the shape of the virtual cathode (the capability of the plasma to screen off the external potential is related to the electron Debye length). The EG potential modifies the amplitude of the plasma potential in the extraction region (it acts similarly to a bias voltage applied on the PG). The trajectories of the charged particles in the vicinity of the aperture are hence altered and the virtual cathode depth is in turn impacted. The latter is increasing (in absolute value) with the extraction potential. This induces a drop of the extracted negative ion current as well. In order to derive scaling laws, i.e., to correlate the plasma properties between high versus low plasma densities, we must preserve the curvature radius of the meniscus. One approximate solution is to calculate the extraction potential with the Child-Langmuir law, where the extracted negative ion currents have hence been retained (primes denote the voltage and current for a > 1). Figure 22 shows the negative ion current density profile along (OY) on the EG versus the EG potential and the plasma density. ¢ = V 200 V EG in (a) and 900 V in (b) correspond to plasma meniscus with nearly identical curvature radius (r a 1 c  ). In both cases, doubling the extraction potential induce the appearance of aberrations in the extracted ion beam profile (r a 1 c ). The profiles displayed in figures 22(a) and (b) are very similar beside a factor 16 in plasma density. This shows the capability of the model to reproduce correctly the beam dynamics with lower densities. This is also visible in figure 22(c) which shows the beam profiles for the same meniscus curvature radius (r a 1 c  ) and for a density ranging from a = 1 to 64. The numerical resolution is2048 1024 (256×128) grid nodes for a = 1 (a = 64), respectively, and 40ppc. = V 0 LB V. The extracted negative ion current varies by a factor 2 (scaling factor) between a = 1 and 64 while the virtual cathode depth by a ratio of ∼1.6. The latter may be used to estimate the negative ion current density extracted from circular apertures instead of slits. The only modification in the model is the geometry of the extraction aperture and hence we posit that the scaling factors are preserved. The scaled negative ion current derived from a V. The numerical resolution waś 256 96 2 nodes with a = 64. We did not model the mechanisms to either produce negative ions in volume via the dissociative attachment of the background gas molecules or from the impact of positive ions on the cesiated PG. Speth et al [5] reports that ∼20% of the extracted ions (gas pressure of ∼0.4 Pa for a RF power of ∼100 kW) in BATMAN originated from volume processes. Adding this small correction to the numerical estimates gives . Note that a slight disparity in the virtual cathode profile between the model and the experiments translates into a significantly larger difference for the amplitude of the extracted negative ion current. The properties of the virtual cathode are unknown experimentally. In addition, the current model (a zoom around a single aperture with lateral periodic boundary conditions) neither account for the complex dynamics of the electrons resulting from the magnetic drifts inside the expansion chamber nor for the properties of the positive ions which are produced in the discharge and impact the PG with a high average energy [22]. Both will modify the characteristics of the virtual cathode. This question should be addressed in future work.

Conclusion
This paper proposes a synthesis, illustrated with new results, of the work that has been carried at LAPLACE on the modeling and simulation of the negative ion sources for the neutral beam injection systems used in fusion applications (ITER, DEMO). We can summarize the important results as follows.
• Fluid models coupled with the plasma chemistry and neutral flow in the negative ion source, show that the gas density is non-uniform in the chamber and is strongly depleted due to gas heating and ionization. Due to the low gas density, the properties of the plasma and gas phase are controlled not only by collisions in the volume but by interaction with the walls. For example, the recombination of high energy ions at the walls contributes to the generation of fast neutrals in the discharge chamber (the presence of high energy, i.e. up to 50 eV ions is due to the high plasma potential resulting from the low gas density). Moreover kinetic (DSMC) simulations of the neutral gas transport show that the velocity distribution function of neutral species is strongly non-Maxwellian, with a highly populated tail. This can certainly affect the production of negative ions on the PG surface. The fluid models of the plasma developed at LAPLACE are very useful and efficient to study the plasma properties with complex chemistry and neutral transport. They are however more difficult to use at high magnetic fields (problems with numerical accuracy and convergence) and are not adapted to the description of negative ion extraction.
• Particle (PIC-MCC) simulations are very powerful but due to constraints on the grid spacing and time step, cannot be used for the high plasma densities and large volume of the negative ion source. We have discussed in this paper how simulations performed under conditions of smaller plasma densities ('scaling'), although not 'exact', can provide a very useful insight in the physics of the source, and that the results can often be linearly extrapolated to the real densities with a good approximation. Care must however be taken when the sheaths play an active role in the discharge or in the simulation of negative ion extraction (and when instabilities and turbulence are present). One important conclusion of the 'scaled' PIC-MCC simulations is the demonstration that the presence of the magnetic filter can induce a strong asymmetry in the plasma properties in the direction parallel to the PG. This is because the diamagnetic andÉ B drift in the magnetic filter region are not closed, as in closed-drift sources (Penning source, magnetron discharge, Hall thruster etc...) but are directed toward a wall. This induces a Hall effect which generates the plasma asymmetry and enhances electron transport though the filter. Calculations performed with a large scaling factor tend to show that the plasma asymmetry leads to a significant non-uniformity of the negative ion current density along the extraction apertures. The validity of linear extrapolation of this result to the real conditions of plasma density is difficult to prove but nevertheless, these simulations raise the question of beam uniformity.
• The description of negative ion extraction by PIC-MCC simulations is also difficult because of the high plasma densities and even if the simulation domain is limited to a small region around a grid aperture, with periodic boundary conditions. We have shown that erroneous and unphysical results can be obtained if the constraints imposed by the PIC method are not strictly taken into account. We have shown that the conclusion, drawn in some recent papers, that only a small negative ion current can be extracted from the PG surface facing the plasma due to space charge current limitation is highly questionable and is a consequence of a misuse of the PIC-MCC technique. 3D simulations respecting the constraints inherent to an explicit PIC simulation seem difficult to perform or to use for systematic parametric studies even for parallel computations on a large number of nodes. We believe that interesting and useful insight could be obtained from simulations with