Tutorial: Modeling of the Extraction and Acceleration of Negative Ions from Plasma Sources using Particle-based Methods L. Garrigues 1 and G. Fubiani 1 LAPLACE, Université de Toulouse, CNRS, 31062 Toulouse Cedex 9, France (*Electronic mail: laurent.garrigues@laplace.univ-tlse.fr, gwenael.fubiani@cnrs.fr) (Dated: 24 January 2023)

In this Tutorial we consider plasma sources with applications to fusion devices and high energy accelerators. These ion sources typically produce negative ions from hydrogen-isotope gases which are extracted through one or multiple apertures and accelerated to high kinetic energies. Next, they are either double stripped of two electrons to form positive ions used as precursors in accelerator devices or neutralized to produce a neutral beam injected in Tokamak reactors. Contrary to the working conditions of most ion sources where volume production prevails, the mechanism of negative ion production by dissociative electron attachment on vibrationally excited molecules inside the plasma volume of fusion-type hydrogen-fueled high power discharges is mostly balanced by their destruction by detachment before being extracted rendering this mean of producing negative ions rather inefficient. Surface production through the transfer of electrons from low work function metallic materials to the impacting atoms is the alternative solution to fulfill the requirements for the applications concerned. Negative ions are produced close to the aperture from which they are extracted. As a result, the analysis and understanding of the extraction mechanisms through experimental diagnostics is rather difficult due to the lack of accessibility and can only give a partial view. In addition, most of the experimental work is focused on the validation of requirements for the applications and not to the investigation of the fundamental processes that take place inside these types of sources. This Tutorial is focused on the description and understanding of the physical mechanisms behind the extraction and acceleration of negative ions from hydrogen plasma sources through modeling methods. We describe the numerical techniques of particle-based methods with a specific emphasis on particle-in-cell Monte Carlo collision algorithms. An analysis of the physical processes involved in driving the negative ions from the plasma source, across the apertures and inside the accelerator as reported in the literature is presented in detail. This Tutorial concludes with additional and future works to be addressed in the coming years.

I. INTRODUCTION

After collisions, atoms and molecules can form positive ions (cations) by losing one or more electrons. Some elements, based on nonmetal and halogen categories of the periodic table, have in parallel the capability to capture one or more electrons (so called electron affinity [START_REF] Greenwood | Chemistry of the Elements[END_REF] ) to form negative ions (anions). In the literature, discharges occurring in this kind of gases are called electronegative plasma discharges or negative ion sources when the ions are extracted and accelerated. In several low-temperature plasma applications, it can be preferable to generate plasmas with extremely few electrons if any. Highly electronegative plasma discharges based on or including halogen atoms/molecules can be composed quasi-exclusively of negative and positive ions (see Ref. 2 for a review on ion-ion plasmas). In microelectronics, plasma etching and deposition of thin films may generate surface charging due to hot electrons, leading to damages in the printed circuit boards (Refs. 3 and 4). In gridded ion thrusters for space applications, a flux of positive ions induces a thrust as a result of their extraction and acceleration from a plasma source using a system of grids negatively polarized with respect to the plasma potential. An external thermionic cathode is used as a neutralizer. Its role consists in generating an electron flux to counteract the positive electric charging by the ion flux preventing the charge accumulation on the spacecraft body (Refs. 5 and 6). A novel concept of gridded ion thruster that accelerates alternately positively and negatively charged ions to provide thrust and to eliminate the need for an additional neutralizer has been proposed [START_REF] Aanesland | The PEGASES gridded ion-ion thruster performance and predictions[END_REF][START_REF] Aanesland | A review on ion-ion plasmas created in weakly magnetized electronegative plasmas[END_REF] . In highly electronegative gases, to render the generation of negative ions by attachment processes efficient, the discharge must be rich in low energetic electrons in the plasma volume. They can be produced either in the temporal afterglow of pulsed discharges [START_REF] Economou | Fundamentals and applications of ion-ion plasmas[END_REF] or by magnetic filtering (trapping) of electrons [START_REF] Levko | Chemical composition of SF6 low-pressure plasma in magnetic field[END_REF] .

Hydrogen negative ion type sources are used in a variety of applications with common goals of generating a large negative ion current (as high as tens of amps in fusion applications) and accelerating them to very high energies (keV to MeV range). Depending upon applications, negative ion beams are either neutralized before entering in fusion machines, or serve as precursors to generate positively charged ions in accelerators that later form neutrons or elementary particles. Tokamaktype devices are based on the fusion of deuterium (D) and tritium (T) that produces a helium nucleus, a neutron and energy. Fusion plasmas need high particle temperatures and must be stable for long periods of times. Power from the alpha particles will be re-injected to the plasma but external heating means are indeed necessary to sustain the discharge (neutral beams and high frequency electromagnetic waves which transfers their power resonantly through the electron and ion cyclotron motion). Neutral-beam injection (NBI) is one of the main techniques to inject external power to the plasma (highenergy neutral particles can enter in the magnetic confinement field). These neutrals are ionized by collisions, and the ions formed become trapped by the magnetic field and transfer most of their remaining energy by further colliding with the plasma. By tangential injection in the torus, neutral beams also provide momentum to the plasma and current drive, the one essential feature for long pulses of burning plasmas. High current and a high power neutral deuterium beam is foreseen for ITER. Two injectors are used, each with 17A of equivalent neutral current at 1 MeV (Refs. 10-12 and references therein). 1280 beamlets are extracted from the ion source over a large surface area of 2 × 0.9 m 2 . The necessary choice of negative ions was made regarding the neutralization capability of positive and negative deuterium ions. The neutralization efficiency drops to almost zero for positive ions for energies larger than 100 keV, while it plateaus to 60% for negative ions [START_REF] Hemsworth | Positive and negative ion sources for magnetic fusion[END_REF][START_REF] Tsumari | A review of diagnostic techniques for high intensity negative ion sources[END_REF] . The detailed characteristics of ITER and DEMO neutral beam injectors can be found in Refs 13, 15-17 for instance.

40 mA beams of negative hydrogen ions at 160 MeV are produced in the Linac4 accelerator. Negative ions after being stripped of two electrons are injected in the Proton Synchrotron Booster (PSB) to be accelerated to 2 GeV and later introduced in the Large Hadron Collider (LHC) [START_REF] Briefi | Experimental benchmark of the NINJA code for application to the linac4 Hion source plasma[END_REF][START_REF]Linac4 design report[END_REF] . The same basic concept is used for the Spallation Neutron Source (SNS) at the Oak Ridge National Laboratory (ORNL in Oak Ridge, Tennessee), where the negative ion sources produce 35 mA of H -each [START_REF] Moehs | Negative hydrogen ion sources for accelerators[END_REF][START_REF] Plum | SNS stripper foil development program[END_REF][START_REF] Stockli | Record productions establish RF-driven sources as the standard for generating high-duty-factor, highcurrent H-beams for accelerators[END_REF] .

Compared to other low temperature plasma applications, fundamental differences exist: 1) the need for large energy hydrogen/deuterium negative ions together with high beam brightness (the ratio between current and emittance), which requires the use of cesium coated surfaces inside the ion source and 2) the use of hydrogen-isotope gases, which have low electronegativity compared to that of a halogen gas or molecules, making the production in volume through attachment processes insufficient to meet the requirements of fusion and accelerator machines. The production of negative ions has to be enhanced with surface production (hence the addition of cesium).

This Tutorial is devoted to the analysis using numerical simulations of the physical mechanisms driving the extraction of negative ions from high power hydrogen/deuterium plasma sources and acceleration of the latter inside an electrostatic accelerator. We will focus primarily on particle-based methods, i.e., particle-in-cell, ray-tracing, test particle algorithms including collision processes implemented numerically via the Monte-Carlo method. The physical-chemistry inducing the generation of negative ions inside the discharge volume as well as from particle-surface interactions will not be described in detail in this Tutorial. We refer instead to depth reviews [START_REF] Bacal | Volume production negative hydrogen ion sources[END_REF][START_REF] Bacal | Negative hydrogen ion production mechanisms[END_REF] , recent Tutorials [START_REF] Bacal | Negative ion sources[END_REF] and articles [START_REF] Cartry | Alternative solutions to caesium in negative-ion sources: a study of negative-ion surface production on diamond in H2/D2 plasmas[END_REF][START_REF] Laporta | Vibrational excitation and dissociation of deuterium molecule by electron impact[END_REF][START_REF] Wada | Measurements of work function and surface conditions in cesiated negative ion sources[END_REF] . The reader interested to have a complete overview of recent results in the context of negative ion sources is invited to read a collection of articles published in a special issue of the New Journal of Physics in 2018 with a dedicated focus on Sources of Negatively Charged Ions 29 .

II. TYPICAL NEGATIVE ION SOURCE CONCEPT

Negative ion sources used for fusion and accelerator machines are based on similar principles. A schematic view of the tandem-type BATMAN (BAvarian Test MAchine for Negative ions, in operation at IPP Garching) test-bed (e.g., Refs 30 and 31) in its upgraded version is shown in Fig. 1. The ITER device combines eight sources similar to BAT-MAN into one [START_REF] Fantz | Development of negative hydrogen ion sources for fusion: Experiments and modelling[END_REF] . Gas used is either H 2 or D 2 . An inductively coupled discharge (ICP) is generated in a cylindrical chamber called the "driver" of 8 l volume at the rear of the source. The coupling is made using an RF power supply with a power frequency of 1 MHz generated through an antenna (six turns) wrapped around the dielectric walls (made in quartz or aluminum oxide). The RF generator feeds high alternating current through the coil antenna, which creates an alternating magnetic field in the plasma chamber. According to Maxwell's equations, the oscillating magnetic field will engender an oscillating electric field within the plasma chamber. Finally, this electric field accelerates and heats the electrons generating the plasma through ionization processes [START_REF] Lieberman | Principles of Plasma Discharges and Materials Processing[END_REF] . A Faraday screen (visible in Fig. 1) protects the dielectric from ion bombardment and suppresses the capacitive coupling. In the first versions of the drivers hot filaments were used but they have been abandoned for lifetime duration and maintenance reasons. The evaporation of tungsten from the filaments and its re-deposition on the grid surfaces has a deleterious impact on the production of negative ions [START_REF] Krylov | Caesium and tungsten behaviour in the filamented arc driven kamaboko-III negative ion source[END_REF] . The power injected in the BATMAN source is on the order of 50-100 kW and the pressure is around 0.3 Pa. One of the critical issues is to increase the power coupling efficiency that is actually limited to 60% [START_REF] Fantz | Negative hydrogen ion sources for fusion: From plasma generation to beam properties[END_REF] . The RF driver for the SNS machine is very similar (but with a smaller volume), the pressure is around 3.2 Pa (24 mTorr), the frequency is 2 MHz, and the injected power is in the range of 55-65 kW [START_REF] Stockli | Record productions establish RF-driven sources as the standard for generating high-duty-factor, highcurrent H-beams for accelerators[END_REF] . The volume of the discharge is 0.4 l in the Linac4 source, the injected RF power is ∼ 40 kW (1 MHz frequency), and the gas pressure is about 3 Pa [START_REF] Briefi | Experimental benchmark of the NINJA code for application to the linac4 Hion source plasma[END_REF] . Note that the plasma in the driver region is mainly an electropositive plasma, and the electron temperature is too high for generating negative ions through attachment processes [START_REF] Fantz | Negative hydrogen ion sources for fusion: From plasma generation to beam properties[END_REF] . The H 2 gas is used in these ion sources.

In the case of RF powered plasma sources for fusion, charged particles are mostly produced inside the driver and flow out toward a second larger volume rectangular chamber referred to as the "expansion" region (for this reason, these two segments are often called a "tandem source") while in accelerator-type devices, the volume of the discharge is a cylinder with the RF coil attached externally on one end near the unique extraction aperture. In the former, a grid with a network of apertures on the chamber wall facing the driver connects the plasma to the accelerator vessel. In both cases, a static magnetic field barrier is generated with external magnets to decrease the flux and average kinetic energy of the electrons. A 2D plot of the magnetic filter field map on the plasma grid (PG) of the BATMAN fusion-prototype source is displayed in Fig. 2(a). The field lines are essentially directed parallel to the grid surface. The corresponding profile along the x direction in front of one aperture is shown in Fig. 2(b). The magnetic field strength varies from ∼ 1.5 mT at the entrance of the expansion region to about 7 mT near the grid meaning that only the electrons are magnetized while the ions may be considered either non-magnetized or somewhat mag- netized depending on their mass. A set of permanent magnets located outside the chamber body generates the magnetic filter field; the maximum amplitude is found 20 mm above the plasma grid apertures inside the plasma. As stated above, the role of the filter field is to induce a strong reduction of the electron temperature, i.e., from 10 eV in the driver down to 1 eV in the extraction region to (1) generate negative ions through electron dissociative attachment processes on vibrationally excited molecules (the cross section is very sensitive to the electron temperature) and (2) to significantly decrease the flux of electrons onto the grid and consequently their coextraction alongside the ions. High energetic electrons exiting the driver region are trapped in the region of a large magnetic field and loose energy colliding with the background neutral gas. Nevertheless, as underlined in Sec. I, the generation of negative ions in volume is not sufficient to fulfill the fusion machine current requirements. This has been largely documented in the literature (see e.g. Ref. 34). The main reasons are the confinement of the negative ions by the plasma potential and their relatively short mean-free-path before being destroyed by collisions. Note that in the Linac4 and SNS plasma sources, additional multicusp magnets are mounted on the cylinder outer surface to reduce electron losses and, hence, increase the plasma density for given input power [START_REF] Briefi | Experimental benchmark of the NINJA code for application to the linac4 Hion source plasma[END_REF] .

An additional source of negative ions is necessary which is performed through the deposition of cesium on the wall chamber (including the bias plate and the plasma grid in blue and violet colors in Fig. 1, respectively). It is evaporated from a Cs oven and injected at the rear of the expansion chamber as shown in The last part of the source is the extraction system. The electrodes must be polarized positively with respect to the plasma potential for extracting negative charges. In BATMAN upgrade, the plasma grid (PG) is at -45 kV, the extraction grid (EG) is between 5 and 10 kV above the PG voltage, and finally the grounded grid (GG) at ground potential [START_REF] Fantz | Negative hydrogen ion sources for fusion: From plasma generation to beam properties[END_REF] . Extraction and grounded grids are made of copper while the PG is in molybdenum [START_REF] Hemsworth | Neutral beams for ITER[END_REF] . The apertures have a diameter of 10 mm (14 mm for the specific requirements of ITER [START_REF] Fantz | Physical performance analysis and progress of the development of the negative ion RF source for the ITER NBI system[END_REF] ) together with a chamfered shape to increase the surface area on the plasma side where negative ions are produced. Intense research activities have been conducted regarding the design of the hole geometry in the plasma grid (seen for example, Refs. 34, 37-41). Electrons, being negatively charged, are inevitably coextracted from the plasma alongside negative ions [START_REF] Wünderlich | NNBI for ITER: status of long pulses in deuterium at the test facilities BATMAN upgrade and ELISE[END_REF] . They induce non-negligible power deposition on the extraction grid when accelerated to the keV energy range and deflected by the local magnetic field [START_REF] Fubiani | Modeling of secondary emission processes in the negative ion based electrostatic accelerator of the international thermonuclear experimental reactor[END_REF][START_REF] Fubiani | Analysis of the two accelerator concepts foreseen for the neutral beam injector of the international thermonuclear experimental reactor[END_REF] . To reduce the co-extracted electron current, a suppression magnetic field with a cusp profile surrounding each apertures is generated using permanent magnet bars embedded in the EG with alternating north and south poles along (Oy) as shown in Fig. 1. The profile of the suppression field in BATMAN along the x direction shown in Fig. 2(b) presents a sharp peak with more than 125 mT at the EG on the aperture axis and a magnetic field of 25 mT at the PG. The suppression field, hence, penetrates inside the plasma and does deflect the electrons away from the grid apertures. In addition, experimental campaigns have also highlighted a reduction of the co-extracted electron current with an increase of the magnetic filter field strength [START_REF] Franzen | Commissioning and first results of the ITER-relevant negative ion beam test facility ELISE[END_REF][START_REF] Franzen | On the electron extraction in a large RF-driven negative hydrogen ion source for the ITER NBI system[END_REF] and when a positive voltage with respect to the bias plate is applied to the PG (typically +20 V) [START_REF] Heinemann | Towards large and powerful radio frequency driven negative ion sources for fusion[END_REF] . The latter case is caused by a reduction of the sheath potential height in front of the PG surface (see Fig. 18 of Fubiani et al. [START_REF] Fubiani | Modeling of plasma transport and negative ion extraction in a magnetized radiofrequency plasma source[END_REF] ).

The large driver and expansion chamber volumes allows the implementation of several plasma diagnostics providing a measure of a comprehensive set of plasma parameters. In the case of the BATMAN device, intrusive diagnostics im-FIG. 3. Electron temperature and density obtained for varying filter field strengths from OES at different lines of sight : axially through the driver (labeled "driver") or horizontally at a distance of 1 cm after the driver exit ("expansion") and 2.6 cm from the grid surface ("PG"), respectively. For comparison, the values determined from Langmuir probe measurements also at a distance of 2.6 cm from the PG are included (open symbols). Note that the filter field strengths are given for an axial position of 2 cm above the PG surface and both vertically and horizontally centered [START_REF] Fantz | Negative hydrogen ion sources for fusion: From plasma generation to beam properties[END_REF] . Reproduced with permission from Fantz et al., Front. Phys. 9, 09651 (2021). copyright 2021 Frontiers Media SA. mersed in the plasma (Langmuir probes [START_REF] Fantz | Plasma expansion across a transverse magnetic field in a negative hydrogen ion source for fusion[END_REF][START_REF] Neely | A langmuir probe system for high power RF-driven negative ion sources on high potential[END_REF][START_REF] Schiesko | Magnetic field dependence of the plasma properties in a negative hydrogen ion source for fusion[END_REF] ); non-intrusive ones, which are laser-based, [START_REF] Berger | Cavity ringdown spectroscopy on a high power RF driven source for negative hydrogen ions[END_REF][START_REF] Christ-Koch | Laser photodetachment on a high power, low pressure RF-driven negative hydrogen ion source[END_REF] and, finally, optical emission spectroscopy (OES) [START_REF] Fantz | Spectroscopy-a powerful diagnostic tool in source development[END_REF] have been extensively used to measure among others densities and temperatures of charged particles either locally or integrated along line of sights. The influence of the magnetic filter field on electron properties is illustrated in Fig. 3 after the source was cleaned of cesium (with a production of negative ions solely in volume).The electron temperature drops from the driver down to the extraction region as expected. The electron density also decreases by an order of magnitude from 7 × 10 17 to 7 × 10 16 m -3 for a maximum magnetic field strength of ≃ 6 mT, for instance. Negative ion densities measured with OES and cavity ring-down spectroscopy a few centimeters from the PG are about a factor of 10 larger for a cesiated (a density of ∼ 10 17 m -3 ) vs. a noncesiated source [START_REF] Christ-Koch | Laser photodetachment on a high power, low pressure RF-driven negative hydrogen ion source[END_REF] (∼ 10 16 m -3 ), respectively. Electron and negative ion densities are typically similar a couple centimeters above the PG. The maximum strength of plasma potential is on the order of 45 V (corresponding to a few times the local electron temperature as expected in an electro-positive discharge [START_REF] Lieberman | Principles of Plasma Discharges and Materials Processing[END_REF] ). The plasma potential profile along the axis (x direction) is maximum in the center of the driver and decreases down to the potential value of the PG electrode (20 V). The electron pressure gradient is slightly greater than the electric field in the driver and, as a consequence, the electrons flow from the driver toward the expansion region where the magnetic filter perpendicular to the source axis (along the z direction) impedes their motion. A diamagnetic ∂ x P e × B z electron drift (Hall current) takes place along (Oy) directing the electron flux toward the top side of the source (or the bottom side of the source if the filter field is reversed), while the positive ions remain un-magnetized. P e = n e T e is the electron pressure. Excess charges accumulate near the corresponding surfaces and, as a result, the plasma polarization induces an electric field in the y direction (this is called the Hall effect) in order to cancel out the Hall current. This plasma response generates a top-bottom asymmetry in the plasma parameters, which has been largely studied and documented in the literature through numerical simulations (fluid [START_REF] Hagelaar | Plasma transport across magnetic field lines in low-temperature plasma sources[END_REF][START_REF] Boeuf | Physics of a magnetic filter for negative ion sources. ii. exb drift through the filter in a real geometry[END_REF][START_REF] Lishev | Fluidmodel analysis on discharge structuring in the RF-driven prototype ion source for ITER NBI[END_REF] and kineticbased [START_REF] Fubiani | Modeling of plasma transport and negative ion extraction in a magnetized radiofrequency plasma source[END_REF][START_REF] St | Physics of a magnetic barrier in low-temperature bounded plasmas:insight from Particle-In-Cell simulations[END_REF][START_REF] Fubiani | Plasma asymmetry due to the magnetic filter in fusion-type negative ion sources: Comparisons between two and three-dimensional Particle-In-Cell simulations[END_REF] ) and experiments [START_REF] Schiesko | Magnetic field dependence of the plasma properties in a negative hydrogen ion source for fusion[END_REF] in the context of neutral beam injectors for Tokamaks. The consequence of the observed asymmetry on the plasma properties is of first importance regarding the operation of the source. The nonuniform ion flux that bombards the caesiated PG results in a production yield of negative ions that varies along the top-bottom direction and which can in turn impact the uniformity of the extracted ion current density profile beyond the acceptance value of the accelerator. Numerical studies both with fluid and PIC-MCC algorithms on the effect of the magnetic filter field and/or the PG bias voltage on the plasma behavior are reported in Refs. 56 and 58.

A conceptual design of the accelerator for the future DEMO fusion device is shown in Fig. 4. Beyond the PG and EG, four acceleration grids (AG1-AG3 plus a grid at ground voltage) with a potential difference of 200 kV between each grids providing a total of 810 keV energy to the ion beam. This concept is a heritage of the MAMuG (MultiAperture MultiGrid) accelerator (Refs. 12, 59-62) used for ITER. Another depreciated concept called SINGAP (Single Gap, Single Aperture) [START_REF] Kashiwagi | Analyses of high power negative ion accelerators for ITER neutral beam injector[END_REF][START_REF] Hemsworth | Status of the ITER heating neutral beam system[END_REF][START_REF] De Esch | Results of the SINGAP neutral beam accelerator experiment at JAEA[END_REF][START_REF] Hemsworth | Status of the development of the SINGAP accelerator for ITER[END_REF] proposed to first pre-accelerate the negative ion beamlets to 55 keV and within one single step to a kinetic energy of 1 MeV. The accelerator configurations used at LHC and SNS are based on similar electrostatic multi-stage designs [START_REF]Linac4 design report[END_REF][START_REF] Stockli | Record productions establish RF-driven sources as the standard for generating high-duty-factor, highcurrent H-beams for accelerators[END_REF] .

Global performance of the accelerator has been assessed mainly through measurements of electron currents collected by the grids, power deposited on surfaces using calorimeters and beam optics with infrared cameras (see Refs. 12, 47, 59, 61, 62, 64, and 65). The driver and expansion chamber of the ion source have multiple diagnostic ports allowing the collection of a wide range of plasma parameters such as the plasma potential, electron density, temperature, neutral, negative ion properties, etc. using, among others, movable Langmuir probes and optical devices (optical emission spectroscopy -OES and laser photo-detachment). Measurements coupled with numerical simulations have greatly improved the understanding of the source operation. Production on surfaces and extraction of negative ions are significantly more difficult to study experimentally due to the fact that the particle transport occurs within or near the plasma sheath highlighting the preponderant role of numerical simulations in providing a precise description of the physics at play.

III. MODELING PRINCIPLES

The plasma is composed of various positive ion species in the case of hydrogen-isotope gases used in fusion-type ion sources, such as (taking hydrogen as an example) H + , H + 2 , and H + 3 (plus a small amount of Cs + , their density being one or two orders of magnitude below -see e.g. Refs. 67 and 68), H -negative ions, H and H 2 neutral atoms, and electrons. Two different categories of numerical methods exist to model the transport of particles. The first one is derived from the three first moments of the Boltzmann equation through continuity, momentum and energy equations [START_REF] Lishev | Fluidmodel analysis on discharge structuring in the RF-driven prototype ion source for ITER NBI[END_REF][START_REF] Boeuf | Fluid and hybrid models of non equilibrium discharges[END_REF][START_REF] Economou | Hybrid simulation of low temperature plasmas: A brief tutorial[END_REF] . This Eulerian approach considers each of the plasma charged species as a fluid and this technique is adequate when the distributions of charged particles are closed to Maxwellian equilibrium. As already discussed in Sec. II, this approach has been successfully employed to model the ion source operation. The second one is a Lagrangian approach based on the so-called particlein-cell (PIC) method often chosen when modeling beam extraction due to the relative simplicity of developing the algorithm compared to fluid models. Note that the plasma may be assumed as a fluid except for the negative ions when they are produced via surface processes. In this case inertia should be retained for positive ions to properly calculate the characteristics of the space charge limited collision-less sheath induced by the large negative ion current density on the surface. Last, depending on the relevance of electron-electron collisions, the electron distribution function may not be Maxwellian due to the presence of the magnetic field. All these arguments render particle-based methods more practicable to model the transport of charged particles in the extraction region of the ion source. The particle's Newton equations of motion are coupled to Poisson's equation which has the charge density as a source term. The latter is calculated using the particle positions. The current densities produce a negligible magnetic field and, hence, only those generated by external means (permanent magnets, for instance) are included in the simulations. PIC algorithms follow so-called macroparticle trajectories in time (each containing a large number of real particles) interacting both with self and external fields. The self-field is deduced from Poisson's equation and the method is iterative. Fields are calculated on a grid while macroparticles represent moving patches of phase-space (position-velocity space) which can overlap each others. This method cancels out large angle scattering, i.e., particle-particle collisions occurring within a sphere of electron Debye radius (finite-size particles have lost their short-range interactions). The PIC technique, hence, provides a solution of the Vlasov equation but for electric and magnetic fields which are computed as the weighted average over the shape of the macroparticle [START_REF] Lapenta | Particle simulations of space weather[END_REF] . The inclusion of collisions must be done separately through a Monte Carlo collision -MCC -algorithm [START_REF] Nanbu | Probability theory of electron-molecule, ion-molecule, molecule-molecule, and coulomb collisions for particle modeling of materials processing plasmas and cases[END_REF] . Note that in low pressure ion sources (∼ 0.3 Pa), collisions are mainly between charged particles and neutrals, and numerically, the latter are often simply considered as a fixed background profiles.

One of the advantages of PIC-MCC models is their capability to be efficiently parallelized in order to run on high performance computing machines. One computer node is generally composed of two CPU sockets, each of them containing many cores in itself. The objective is to run one simulation in parallel on more than one core to reduce the computational time. Two strategies of parallelization exist. The first option is to split the computational domain in sub-domains (so-called domain decomposition), each sub-domains being associated with a different core. The computation within each core is performed independently in parallel. This strategy is the best choice for 3D simulations with a large number of grid points (typically greater than, say, 512 3 ) when many computer nodes are used (beyond ∼ 4). The 3D codes called ONIX (Orsay Negative Ion eXtraction) [START_REF] Mochalskyy | 3D modelling of negative ion extraction from a negative ion source[END_REF] and ONAC (Orsay Negative ion ACelerator) [START_REF] Revel | Transport of realistic beams in ITER neutral beam injector accelerator[END_REF] developed at LPGP (Laboratoire de Physique des Gas et des Plasmas), Orsay, France, are based on that principle. An attempt to use the ONAC code was made on a GPU machine [START_REF] Revel | Transport of realistic beams in ITER neutral beam injector accelerator[END_REF] . Recent simulations of the BAT-MAN extraction region have been achieved using the ONIX code with 4096 cores (Intel Sandy Bridge) of the HYDRA supercomputer at the Max Planck Computing and Data Facility (MPCDF), Garching, Germany [START_REF] Wünderlich | Review of Particle-In-Cell modeling for the extraction region of large negative hydrogen ion sources for fusion[END_REF][START_REF] Wünderlich | Effects of the magnetic field topology on the co-extracted electron current in a negative ion source for fusion[END_REF] . The second option is called particle decomposition and based on the assignment on each of the cores of a subset of the initial population of particles. This method is very efficient for 2D or 3D simulations when the number of grid points is, hence, not too large because the whole spatial grid is assigned to each core. The numerical resolution is limited by the memory available per CPU. The drawback of this approach is the communication cost between MPI threads for the concatenation of the charge density (source term for Poisson's equation) which may become prohibitive as the number of grid nodes is increased. This option has been implemented in the Electrostatic Parallel PIC (EP-PIC) codes of the LAPLACE laboratory (Laboratoire Plasma et Conversion d'Energie), Toulouse, France. Up to ∼ 200 cores are typically used for 2D/3D simulations on the Olympe supercomputer machine at the Toulouse University and the CINECA supercomputer, Bologna, Italy, [START_REF] Fubiani | Modeling of plasma transport and negative ion extraction in a magnetized radiofrequency plasma source[END_REF][START_REF] Garrigues | Appropriate use of the Particle-In-Cell method in low temperature plasmas: Application to the simulation of negative ion extraction[END_REF][START_REF] Garrigues | Negative ion extraction via particle simulation for fusion: Critical assessment of recent contributions[END_REF][START_REF] Fubiani | Modeling of negative ion extraction from a magnetized plasma source: Derivation of scaling laws and description of the origins of aberrations in the ion beam[END_REF][START_REF] Garrigues | Application of sparse grid combination techniques to low temperature plasmas Particle-In-Cell simulations. I. capacitively coupled radio frequency discharges[END_REF] . Additional information to speed-up the computational time using vectorization and particle sorting methods can be found in Refs. 48, 77, and 80. Different strategies regarding the communication between cores have also been implemented, either exclusively using MPI (message passing interface) libraries to exchange information between cores [START_REF] Mochalskyy | 3D modelling of negative ion extraction from a negative ion source[END_REF][START_REF] Revel | Transport of realistic beams in ITER neutral beam injector accelerator[END_REF][START_REF] Wünderlich | Review of Particle-In-Cell modeling for the extraction region of large negative hydrogen ion sources for fusion[END_REF][START_REF] Wünderlich | Effects of the magnetic field topology on the co-extracted electron current in a negative ion source for fusion[END_REF] , or adopting a hybrid approach by limiting the use of MPI to transfer data between sockets or nodes, while shared memory capabilities are implemented between cores via OpenMP (Open Multi-Processing) [START_REF] Fubiani | Modeling of plasma transport and negative ion extraction in a magnetized radiofrequency plasma source[END_REF][START_REF] Garrigues | Appropriate use of the Particle-In-Cell method in low temperature plasmas: Application to the simulation of negative ion extraction[END_REF][START_REF] Garrigues | Negative ion extraction via particle simulation for fusion: Critical assessment of recent contributions[END_REF][START_REF] Fubiani | Modeling of negative ion extraction from a magnetized plasma source: Derivation of scaling laws and description of the origins of aberrations in the ion beam[END_REF] . The latter technique has the advantage of reducing the cost of communication between threads saving computational time. We will here recall the basics of the PIC algorithm and its implementation in the context of negative ion extraction and acceleration. The PIC-MCC method is detailed in text books (see Refs. 81-84) and review papers [START_REF] Birdsall | Particle-In-Cell charged-particle simulations, plus Monte Carlo Collisions with neutral atoms, PIC-MCC[END_REF][START_REF] Verboncoeur | Particle simulation of plasmas: review and advances[END_REF][START_REF] Lapenta | Particle simulations of space weather[END_REF] . Figure 5 shows a series of numerical steps during one time lag ∆t of the electrostatic PIC-MCC algorithm. In phase (1), the charged density profile ρ is assigned on the nodes of the computational grid from the projection of the charge carried by the macroparticles q i at positions x i using a linear interpolation. In step (2), Poisson's equation is solved to calculate the electric potential profile φ on the nodes of the grid. After discretization using a central finite difference method onto the computational grid nodes, a system of linear equations Aϕ=b has to be solved. The matrix A has constant coefficients coming from the spatial discretization while the vector b contains the space charge and boundary conditions. The resolution of the system can be carried out using an iterative or a direct method (see Ref. 88 and references therein for a panel of existing methods). Iterative methods are the best option over direct ones when the number of grid nodes typically exceeds 10 6 . In addition, the former can be efficiently parallelized to run on multiple nodes at once. Nowadays, many compilers include software libraries for numerical algebra with routines for the resolution of Poisson's equation (e.g. PARDISO [START_REF]PARDISO -parallel direct solver project[END_REF] and PETSc 90 are accessible in the Math Kernel Library of INTEL Fortran and C++ compilers). Other software are also available on a GitHub open-source platform [e.g. FISHPACK [START_REF]FISHPACK, a fortran library for solving of separable elliptic partial differential equations[END_REF] , HYPRE [START_REF]HYPRE, scalable linear solvers and multigrid methods[END_REF] , UMFPACK [START_REF]UMFPACK, unsymmetric multifrontal method package[END_REF] , and PittPack 94 , the last one being dedicated to Graphic Processor Unit (GPU) architecture]. Garrigues et al. [START_REF] Garrigues | Appropriate use of the Particle-In-Cell method in low temperature plasmas: Application to the simulation of negative ion extraction[END_REF][START_REF] Garrigues | Negative ion extraction via particle simulation for fusion: Critical assessment of recent contributions[END_REF] and Taccogna et al. [START_REF] Taccogna | Modeling of a negative ion source. III. Two-dimensional structure of the extraction region[END_REF] are employing a direct method using PARDISO and UMFPACK solvers, respectively, based on an LU decomposition [START_REF] Demmel | A supernodal approach to sparse partial pivoting[END_REF][START_REF] Petra | An augmented incomplete factorization approach for computing the schur complement in stochastic optimization[END_REF] to inverse the matrix. Note that the preconditioning of the matrix A, which is necessary to reduce the computational time for convergence is only performed once at the beginning of the simulation since the coefficients of A remain unchanged, and only the source term (space charge) of Poisson's equation evolves in time. Furthermore, an iterative gradient conjugate algorithm [START_REF] Kaasschieter | Preconditioned conjugate gradients for solving singular systems[END_REF] is implemented by Mochalskyy et al. [START_REF] Mochalskyy | 3D modelling of negative ion extraction from a negative ion source[END_REF] and a bi-conjugate gradient stabilized method [START_REF] Van Der | BI-CGSTAB: a fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems[END_REF] in the 3D code of Nishioka et al. 100 . Fubiani et al. [START_REF] Fubiani | Modeling of plasma transport and negative ion extraction in a magnetized radiofrequency plasma source[END_REF] has implemented a parallel geometric multi-grid approach using a successive over relaxation (SOR) scheme 101 (home-made algorithm).

In step (3), the electric field profile E is computed on the grid nodes from the electric potential with a standard central difference method. The electric field must be calculated with caution on surfaces (Dirichlet boundary conditions) since the electric potential is not computed inside the materials. Garrigues and Fubiani et al. [START_REF] Fubiani | Modeling of plasma transport and negative ion extraction in a magnetized radiofrequency plasma source[END_REF][START_REF] Garrigues | Appropriate use of the Particle-In-Cell method in low temperature plasmas: Application to the simulation of negative ion extraction[END_REF] use a second order linear extrapolation of the fields at forward points inside the simulation domain to deduce the electric field on the surface itself 102 . In the ONIX 3D code [START_REF] Mochalskyy | 3D modelling of negative ion extraction from a negative ion source[END_REF] , at the border, the component of the electric field normal to the walls is calculated from the two first half-mesh points using a correction method 103 . In the ONIX code, particular attention is also made regarding the calculation of the electric potential around the cylindrical aperture inside the grid by smoothing out the electric potential profile calculated on the Cartesian grid to suppress the non pure radial effect [START_REF] Mochalskyy | 3D modelling of negative ion extraction from a negative ion source[END_REF] . In step (4), the field forces including the magnetic field F i are linearly interpolated at the particle locations from the fields calculated on mesh nodes. In step (5), particle position x i and velocity v i are updated. Newton's equations of motion are integrated using the Boris-Buneman leap-frog scheme 104 which has the advantage not to require the storage in memory of positions and velocities at intermediary steps. Electrons can be accelerated to energies up to 1 MeV inside the accelerator and the relativistic factor γ = 1/ 1v 2 e /c 2 , where v e and c are electron and light velocities, respectively, is on the order of 3. Implementation of the relativistic factor in the Boris-Buneman algorithm is described in Ref. 83 and by Revel et al. [START_REF] Revel | Transport of realistic beams in ITER neutral beam injector accelerator[END_REF] in the context of the calculation of negative ion transport inside an electrostatic accelerator using the ONAC code. Note that for large γ factors, new particle pushers, including implicit algorithms to better conserve the kinetic energy, have been proposed and benchmarked; see for, instance, the work of Ripperda et al. 105 ).

Particle interactions with walls are handled in step (6). The production of negative ions under ion and neutral particle impacts on the caesiated grid surface facing the plasma has been self-consistently calculated in the work of Taccogna and Fubiani et al. [START_REF] Taccogna | Three-dimensional structure of the extraction region of a hybrid negative ion source[END_REF][START_REF] Taccogna | Particle modelling of the hybrid negative ion source[END_REF]106,107 . For an incident neutral flux Γ n , the incoming flux of negative ions is simply expressed as Γ H - = Γ n γ(H -), where γ(H -) is the yield of conversion from a neutral to a negative ion. The yield has been derived from experimental measurements by Seidl et al. 108 . Positive ion impacts can also be precursors to negative ions, but calculations have shown that the positive ion fluxes onto the walls are typically much smaller, their contributions represent less than ∼ 1 % of the total [START_REF] Taccogna | Three-dimensional structure of the extraction region of a hybrid negative ion source[END_REF][START_REF] Taccogna | Particle modelling of the hybrid negative ion source[END_REF] . The induced negative ion current density j H -deduced from these calculations for the BAT-MAN source [START_REF] Fubiani | Modeling of plasma transport and negative ion extraction in a magnetized radiofrequency plasma source[END_REF][START_REF] Wünderlich | PIC code for the plasma sheath in large caesiated RF sources for negative hydrogen ions[END_REF]109 is on the order of 600 A.m -2 (similar to what is expected for the KAMABOKO negative ion source in the JT-60U Japanese Tokamak 110 which is powered by filaments instead of a RF coil). The neutral flux is estimated from the experimental assessment of the neutral density and temperature. In the numerical models, a half-Maxwellian flux [START_REF] Garrigues | Appropriate use of the Particle-In-Cell method in low temperature plasmas: Application to the simulation of negative ion extraction[END_REF] corresponding to a fixed current density of 600 A.m -2 of negative ions is injected at a temperature typically between 0.5 and 1 eV on the grid surface facing the plasma (e.g. Refs. 48, 77, 95, and 111). A large fraction of these negative ions will be extracted across the apertures toward the accelerator and further accelerated to high energies (1 MeV in the case of ITER). Electrons despite the suppression magnetic field will be extracted as well. Inside the accelerator vessel secondary particles will be generated due to mostly the interaction between negative ions and the residual background gas which leaks from the ion source or the neutralizer. These particles, namely electrons and heavy ions, are byproducts of negative ion single and double stripping reactions or ionization of the gas. These particles are themselves accelerated to high energies as well and impacts the accelerator grids, producing in turn other secondary particles. This is a critical issue for neutral beam injection devices in fusion plasmas as secondary particles absorb a substantial amount of the powersupply power and can damage the accelerator parts. A detailed Monte-Carlo model of the complex physics occurring inside the accelerator can be found in the EAMCC 3D code of Ref. 43 (see also references therein for input data).

Step (7) focuses on the Monte Carlo module for collisions inside the volume. The probability P for one particle to collide during one timestep ∆t is P = 1e -ν∆t , where ν is the collision frequency. The collision frequency can be expressed as ν = n t (r)v r σ T (ε r , θ ) and depends on the target particle density n t at a given location r, the relative velocity v r between the incident and target species, the total collision cross section σ T that is a function of relative energy ε r and angle θ , respectively. In most of the PIC models developed in the context of negative ion extraction and acceleration, the gas density is uniform in space. It is determined from ideal gas law for fixed pressure and temperature. It can also be estimated from experiments or actual implementation of the neutral dynamics in a numerical model of the entire negative ion source [START_REF] Fubiani | Modeling of plasma transport and negative ion extraction in a magnetized radiofrequency plasma source[END_REF]112 . The knowledge of the collision cross section angular dependence for a large range of energy is rare and only available for a few gas types. All collision cross sections are consequently assumed angularly equiprobable (isotropic). A comprehensive physical chemistry database for hydrogen has been compiled by Janev et al. [START_REF] Janev | Collision processes in lowtemperature hydrogen plasmas JUEL-4105[END_REF] , while for the deuterium gas, data are coming from the work of Barnett 114 . The total collision frequency is expressed as σ T (ε r ) = ∑ j σ j (ε r ) where the integration is made on all the processes j. The calculation of the probability for a particle to undergo a collision is not trivial since it changes along its path in the simulation domain. One can calculate at each time step for each of the particles the collision probability, but this method is very time consuming and inefficient as far as P ≪ 1. One alternative and more efficient approach is to consider a maximum collision frequency instead independent of space and energy ν max = n t × max[(v r σ j (ε r )]. This renders the collision artificially equiprobable for all incident-target macroparticle pairs, which implies that the macroparticles may be chosen randomly inside a given grid cell. Implementing a maximum collision frequency instead of the real one implies that a socalled null collision 115 event might occur. A rejection method is, hence, necessary to check whether or not a real collision actually happened and particles are scattered out. If this is the case then new velocities of either the incident only or both the incident and target particles are computed depending on the model. The former corresponds to the situation where the target species is simulated as a fixed background density and temperature profiles assuming the particle distribution function (PDF) is Maxwellian. The target particle velocity is artificially sampled from the PDF using random numbers. In the latter case both the incident and target particle trajectories are calculated including collisions and irrespective of the charge (meaning even neutrals may be considered). This type of model is called DSMC [START_REF] Nanbu | Probability theory of electron-molecule, ion-molecule, molecule-molecule, and coulomb collisions for particle modeling of materials processing plasmas and cases[END_REF] (direct simulation Monte Carlo). The number of particles undergoing collisions (including a null one) is expressed as N c = NP T , with P T = 1e -ν max ∆t and N being the total number of particles of a given species in the simulation domain. Since one must choose ν max ∆t ≪ 1 to ensure that a particle collides with a maximum of one partner during a time step, the number of particles colliding is small, N c ≪ N, and chosen randomly in the list of the N particles. For this subset of particles, collision probabilities (including the null-collision component) are computed and new component of velocities calculated if necessary (see Refs. 86 and 116 for more general details and Refs. 48 and 77 in the context of this Tutorial). This algorithm models the Boltzmann collision operator which is also suited to the calculation of Coulomb collisions as the Landau-Fokker-Planck equation can be derived from the Boltzmann equation under the assumption of dominant small-angle interactions (this is the case in plasmas). Taccogna and Fubiani et al. [START_REF] Fubiani | Modeling of plasma transport and negative ion extraction in a magnetized radiofrequency plasma source[END_REF][START_REF] Taccogna | Particle modelling of the hybrid negative ion source[END_REF][START_REF] Taccogna | Modeling of a negative ion source. III. Two-dimensional structure of the extraction region[END_REF] have implemented a DSMC algorithm for charged-charged species collisions. In the 3D code developed at KEIO University, electron-neutral and electron-electron collisions are treated differently with a random-walk diffusion model 100,110,[117][118][119] .

The last step (8) concerns the introduction of particles in the computational domain. A 2D slab geometry in the (x, y) plane corresponding to a zoom in front of one grid aperture of the plasma grid (PG) is shown in Fig. 6. Periodic boundary conditions are introduced along the y direction which is equivalent to model the effect on the plasma of an infinite number of apertures. Numerically, any particle at a given location y = y old and crossing the y = 0 (y = y max ) boundary, is re-injected at a new position y new = y old + y max (y new = y oldy max ), respectively, while its velocity components remain unchanged. For 3D simulations [START_REF] Taccogna | Three-dimensional structure of the extraction region of a hybrid negative ion source[END_REF][START_REF] Wünderlich | Effects of the magnetic field topology on the co-extracted electron current in a negative ion source for fusion[END_REF]120 , periodic boundary conditions are also considered along the z direction. Negative ions H -are injected from the PG surface which is assumed at ground potential, whereas the EG is at a fixed positive electric potential (to extract negative charges). Electron and positive (and eventually negative) ion currents are either loaded on the left boundary (LB) itself assuming a semi-infinite plasma 121 and a half-Maxwellian flux distribution function [START_REF] Garrigues | Appropriate use of the Particle-In-Cell method in low temperature plasmas: Application to the simulation of negative ion extraction[END_REF] or considering a finite zone of injection inside the simulation domain where particles are generated 122 , as illustrated in Fig. 6 (gray area). In both cases, an artificial region close to the left boundary is formed. The injection of a charged particle flux through a given plane is rather complex since their magnetization, i.e., their Larmor rotation, induces a fraction of these particles to be recollected by the boundary. In this case, a refluxing method is used to re-inject the particles 123 . This technique assumes that each particle escaping the computational domain across the LB is replaced by a new particle coming from the ion source unsimulated volume; its velocity is computed from the plasma properties at that location which are generally deduced from experimental measurements. The re-injection of negative ions following this procedure is questionable since the distribution function of the negative ion specie generated inside the ion source chamber differs greatly from the one produced on the PG surface after being accelerated toward the left-hand-side by the axial electric field (including the plasma sheath). To reduce the appearance of a transition region, a Neumann condition (normal component of electric field equals to zero) is applied on the LB.

The injection of particles in a finite region inside the domain (instead of a flux on the boundary) has been more generally considered in the context of PIC-MCC modeling of the extraction region. A Dirichlet boundary condition fixing the potential to the ground is applied on the LB plane. Particles striking that plane are eliminated from the simulation and a standard sheath is formed in front of the surface. The region of injection appears in gray in Fig. 6. The position and length of that area are chosen such that the artificial generation of particles does not perturb the appearance of the plasma sheath in front of the PG aperture. One can start injecting a given source term into a simulation domain free of charged particles as it has initially been done in the context of scrape-off layer studies 124,125 . This method has been successfully used in the context of 1D PIC simulations of negative ion formation on caesiated surfaces by Wünderlich et al. [START_REF] Wünderlich | PIC code for the plasma sheath in large caesiated RF sources for negative hydrogen ions[END_REF] . The main drawback of that method is that the plasma density is known in retrospect forcing an iterative adjustment of the source term to reach the desired plasma density. One can in-sure also that the plasma density remains fixed at all times by re-injecting an electron-ion pair (or, eventually, a negative ion instead of an electron with a given flux ratio) each time a positive ion impacts a Dirichlet (metal, i.e, absorbing) boundary condition [START_REF] Fubiani | Modeling of plasma transport and negative ion extraction in a magnetized radiofrequency plasma source[END_REF][START_REF] Garrigues | Appropriate use of the Particle-In-Cell method in low temperature plasmas: Application to the simulation of negative ion extraction[END_REF][START_REF] Garrigues | Negative ion extraction via particle simulation for fusion: Critical assessment of recent contributions[END_REF]126,127 . The velocity components are sampled from a Maxwellian distribution function at a given temperature. This method mimics the physical-chemistry inside the injection region preserving current continuity. The choice of the plasma density, the stoichiometric ratio between particles of a similar charge sign, and species temperature are external parameters to the model and taken either from experimental or simulations of the whole ion source (which is feasible in a 2.5D-3V configuration [START_REF] Fubiani | Modeling of plasma transport and negative ion extraction in a magnetized radiofrequency plasma source[END_REF] ). Note that some authors in the literature have used an alternative which, turned out to be incorrect [START_REF] Taccogna | Three-dimensional structure of the extraction region of a hybrid negative ion source[END_REF][START_REF] Mochalskyy | 3D modelling of negative ion extraction from a negative ion source[END_REF]111,120,128 ; i.e., when one charged particle of a given species impacts a surface, it is directly re-injected with a random position in the injection region [START_REF] Garrigues | Appropriate use of the Particle-In-Cell method in low temperature plasmas: Application to the simulation of negative ion extraction[END_REF] . This means that all charged macroparticles are kept in the calculation (their number remains constant throughout time). This numerical technique is incorrect since the self-consistent formation of a plasma sheath results in an overall different number of positive and negative charges in the simulation domain. The two approaches for re-injection have been illustrated in a simple one-dimensional collisionless and nonmagnetized discharge model composed of electrons, H + 2 ions, and grounded walls. The electron and ion temperatures are set to 2 and 0.1 eV, respectively. The electric potential profile is shown in Fig. 7. The re-injection of the same type of particles leads to a nonphysical electric potential profile for an electro-positive discharge with a minimum of potential in the center of the discharge. The explanation is simply related to the fact that to maintain the current conservation (equal amount of charged particle losses on surfaces), electrons have to be accelerated toward the walls. When one electron-ion pair is re-injected each time an ion strikes the walls, an ion sheath is formed in front of the grounded walls and the electric potential is now maximum in the discharge center, as predicted by the theory [START_REF] Lieberman | Principles of Plasma Discharges and Materials Processing[END_REF] .

Another important issue is how to maintain a Maxwellian distribution at a given temperature for the electrons. Injected electrons diffuse toward the walls and those whose kinetic energy is larger than the potential drop escape from the field, and only low energetic electrons are trapped by the electric potential. As a result, the electron temperature is below the one imposed. A thermalization process can be used [START_REF] Fubiani | Modeling of plasma transport and negative ion extraction in a magnetized radiofrequency plasma source[END_REF][START_REF] Wünderlich | PIC code for the plasma sheath in large caesiated RF sources for negative hydrogen ions[END_REF][START_REF] Garrigues | Appropriate use of the Particle-In-Cell method in low temperature plasmas: Application to the simulation of negative ion extraction[END_REF][START_REF] Garrigues | Negative ion extraction via particle simulation for fusion: Critical assessment of recent contributions[END_REF] . This simply can be done considering at one given time step the number of electrons in the injection region N inj e and a frequency of thermalization ν th e , such that ν th e ∆t≪ 1. A good choice for ν th e is such that electrons are thermalized at least once when crossing the injection zone. At each time step, the number of electrons which are randomly chosen from the list of N inj e electrons and that will be thermalized is deduced from N th e =N inj e ν th e ∆t. For these thermalized electrons, old velocity components are discarded and new ones are chosen from a Maxwellian distribution function at the chosen temperature. This method, which acts as a pseudo-collision, hence does not properly describe the transport of the electrons in the presence of a magnetic field. By changing the velocity components, an artificial cross-magnetic field transport is induced. However, the location of the injection zone is far from the region of interest, as we see in Fig. 6 and one can also choose ν th e below the actual physical total collision frequency of a given electron. This method has been used in the example illustrated in Fig. 7. The maximum of plasma potential found in the simulations is 7.4V which is in agreement with the theoretical value for a Maxwellian population of electrons at a temperature of 2 eV [START_REF] Lieberman | Principles of Plasma Discharges and Materials Processing[END_REF] .

A full PIC-MCC description is generally not employed in the literature to model the transport of charged particles in the accelerator, especially when many apertures and beamlets are included. Tentatives to couple both the plasma (modeling only the region in front of a single aperture) together with the accelerator were so far scarce and for some unsuccessful [START_REF] Revel | Transport of realistic beams in ITER neutral beam injector accelerator[END_REF] . Miyamoto et al. 129 studied the origin of halos in the beam optics that way but had to resort to downscaling of the size of the simulation domain. The study of beam-beam interactions and deflections with a 3D full PIC-MCC description in the real geometry would require tremendous computational resources, which is not feasible. In most of the cases, the electric field map is calculated including the beam space charge only (i.e., not the plasma). Macroparticles are used to carry a current element instead of a charge in these socalled "ray-tracing" algorithms. The plasma meniscus which forms in front of the aperture (a plasma sheath) can be approximated by assuming that it corresponds to an iso-potential of a given value (generally 0V). Typical 3D ion optics codes based on the ray tracing technique are OPERA-3D 130

IV. EXTRACTION OF NEGATIVE IONS

To maintain stability and accuracy in explicit momentumpreserving PIC-MCC algorithms, and to avoid numerical heating, some criteria must be respected (see Refs. 81 and 83 for more details). The first criterion is associated with the grid spacing ∆x that must be smaller than the electron Debye length λ D,e = ε 0 k B T e /e 2 n e where ε 0 , T e , k B , e, n e , are the vacuum permittivity, the electron temperature, the Boltzmann constant, the elementary charge, and the electron density, respectively. The second criterion requires the timestep ∆t to be on the order of or smaller than 0.2/ω p,e , where ω p,e = e 2 n e /ε 0 m e is the electron plasma frequency, with m e being the electron mass. The last criterion refers to the number of particles-per-cell (ppc) that must be used to statistically properly sample the charged particle distribution functions. For typical plasma conditions in the extraction region, with densities of about 10 17 m -3 and an electron temperature of, say, 2 eV (see Fig. 3), 3D simulations of a zoom in front of one extraction aperture require tremendous computational resources even on high performance computing machines with ∼ 10 9 cells for a domain of 8 cm 3 and 4 × 10 10 particles (for 40 ppc). These criteria were not always strictly respected by some authors studying the mechanisms underlying negative ion extraction in the literature, most of the time to reduce the computational time. For example, Mochalskyy et al. 111,120,128 with the ONIX code used a very low number of ppc and a high grid spacing-to-electron Debye length ratio (about 5-6 in the injection region, 2 cm upstream the meniscus region). In the 3D simulations of Taccogna et al. [START_REF] Taccogna | Particle modelling of the hybrid negative ion source[END_REF] , the number of ppc is also very low and ∆x/λ D,e ∼ 6. No numerical convergence analysis had been performed in these studies. As a consequence, this triggered a debate in the negative ion extraction modeling community regarding whether or not one must strictly respect the numerical constraints of explicit PIC schemes 140,141 .

Garrigues et al. [START_REF] Garrigues | Appropriate use of the Particle-In-Cell method in low temperature plasmas: Application to the simulation of negative ion extraction[END_REF][START_REF] Garrigues | Negative ion extraction via particle simulation for fusion: Critical assessment of recent contributions[END_REF] have performed 2D PIC simulations around a grid aperture to study the effect of grid spacing-toelectron Debye length ratios (other parameter variations are shown in the two references above). Electric potential and electron temperature profiles are given in Fig. 8. The grid spacing was varied from 21.3 µm (∆x/λ D,e ∼ 0.7 in the injection region) to 250 µm (∆x/λ D,e ∼ 6). The ratio is not strictly proportional because the electron density (and the electron Debye length accordingly) varies between cases. In the injection region (hatched zone in Fig. 8) where charged particles are loaded, the electron temperature is maintained to 2 eV through the thermalization process described in Sec. III. When the electron Debye length is not resolved (for a grid length larger than 31.2 µm), the electron temperature increases together with the grid size. This is a numerical artifact called selfheating and well documented in reference textbooks [START_REF] Hockney | Computer Simulation using Particles[END_REF][START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF] . The artificial heating induces the electron Debye length to be on the same scale than the grid width while when the electron Debye length is resolved (for a grid length of 31.2 µm or smaller in Fig. 8), the size of the grid no longer influences the tem- perature profile which slightly decreases toward the wall. The consequence of using an inappropriate grid size is visible on the electric potential profile as well which increases to values of about 4 V and sharply decreases to -3.5 V near the wall at x = 26.5 mm when ∆x = 250 µm. The potential well in front of the grid surface is called a virtual cathode and is caused by the negative ions, which are space charge limited in the conditions of the simulation (note that their charged density is partially neutralized by the incoming positive ion flux). The well width is, hence, strongly correlated to the negative ion density profile, which are the dominant species there (meaning that the electron Debye length is likely not a relevant parameter, although this is still an open question). The precise numerical modeling of the negative ion trajectories across the well is very sensitive to the choice of the grid size which must be tailored to resolving the potential profile. The depth of the virtual cathode drops by almost a factor 2 in amplitude in Fig. 8 when the grid size is too large. Unresolved virtual cathodes with potential depths as low as -6 V have been reported in the literature [START_REF] Taccogna | Three-dimensional structure of the extraction region of a hybrid negative ion source[END_REF]128 inducing as a result a negligible extracted negative ion current density that was contradictory with the measurements.

A typical 2D simulation domain of 25 × 16 mm 2 , for a plasma density of 2 × 10 17 m -3 and an electron temperature of 2 eV would require 1100 × 700 grid cells to resolve the electron Debye length. This leads to using important computational resources for calculations taking just a couple of days. One alternative is to make simulations for a reduced plasma density using a scaling law derived from the Child-Langmuir (CL) equation. The CL law has been first introduced to calculate the maximum electron current density, which can be extracted due to its own space charge from two parallel plates in vacuum and maintained at a constant potential difference 142 . Larger currents would induce the formation of a virtual cath-FIG. 10. Axial distribution of the charged particle densities (electrons, n e , positive ions, n i , and negative ions, n n ) and plasma potential Φ in the conditions of Fig. 9 at two different transverse positions, (a) at y = 0.8 cm, i.e., along the axis of the grid aperture, and (b) at y = 1.4 cm between the grid surface and the left boundary (along the top dashed line displayed in Fig. 9a). The front PG surface is at x = 2 cm. The maximum of the negative ion density near the grid surface in Fig. 10b ode such that the fraction above the CL limit would be reflected back and, hence, not transmitted across the gap (at the limit the electric field drops to zero in front of the electrode). This description is merely 1-dimensional (1D). The CL law is of paramount importance for ,instance, to calculate the positive ion beam properties extracted from a plasma source. The shape of the plasma meniscus (its curvature radius) is a consequence that the plasma enforce j ext = j CL in a multi-dimensional problem where j CL is the space charge limited current density and j ext is the ion current density at the sheath entrance, respectively. Therefore, j ext = j CL = A U 3/2 /d 2 , where A is a constant accounting for the vacuum permittivity, the ion mass and, the charge state, U is the voltage between PG and EG (equal to EG for PG at ground), and d is the distance between the plasma meniscus and the extraction electrode. The meniscus can take a concave, flat or convex shape depending on the ion current density crossing its boundary. The CL law is still valid in the case of negative ion extraction but the co-extraction of electrons may affect the menis-cus behavior [START_REF] Fubiani | Modeling of negative ion extraction from a magnetized plasma source: Derivation of scaling laws and description of the origins of aberrations in the ion beam[END_REF] . The CL may, hence, be used to scale down the plasma density in numerical simulations of beam extraction. Assuming n r = β n s then j ext,r = β j ext,s where subscript r stands for the real quantities while s for the scaled quantities; β is the scaling factor (greater than 1), respectively. Considering the same distance between electrodes, the EG voltage U EG,s can be expressed as a function of β and the EG voltage applied in the real conditions U EG,r : U EG,s = U EG,r /β 2/3 . In the context of negative ion sources where negative ions are generated from the PG surface, Boeuf et al. 127 have made the approximation that the plasma meniscus profile is not modified if the source operates with n s or n r as long as the extraction potential is consistent with the CL law. The negative ion current density emitted from the surface must also be scaled down by β . A more general 2D formulation of the scaling law in the context of negative ion extraction is presented in Fubiani et al. [START_REF] Fubiani | Modeling of negative ion extraction from a magnetized plasma source: Derivation of scaling laws and description of the origins of aberrations in the ion beam[END_REF] .

Starting from typical plasma conditions n r = 2 × 10 17 m -3 , U EG = 5 kV, and j H-,r = 500 A • m -2 , Boeuf et al. 127 have performed calculations in 2D with β = 20 fixing the initial plasma density to n 0 = 10 16 m -3 . For a simulation domain of 2.7 × 1.6 cm 2 along the x and y directions, respectively, this leads to model the extraction region around one grid aperture resolving the electron Debye length with a reasonable number of grid points (270 × 160 nodes for a Debye length of 100 µm). Following the scaling law described above, the EG voltage U EG is 700 V and j H -= 25 A.m -2 . The electric potential at the LB is fixed to -5V. A negative ion loss mechanism was accounted for in the calculation by introducing a detachment process with a constant frequency of 10 -5 s. This choice of potential value does traps the negative ions in the simulation domain, which is not a realistic situation for a negative ion source as some current might cross the LB. Nevertheless, this model is sufficient to describe in detail the transport of negative ions from the grid surface toward the plasma volume and across the aperture. The negative ions are injected at the PG according to a half-Maxwellian at a temperature T n = 1 eV. 2D plasma properties in front of the aperture are sketched in Figs. 9 and 10. The axial profiles of Fig. 10 correspond to the dashed lines visible in Fig. 9. In the quasineutral injection region close to the LB plane, positive charges compensate negative charges which are composed of two-thirds of electrons and one-third of negative ions. In front of the LB, an electropositive sheath is formed. The electric potential at the sheath entrance is around 1 V (6 V of the potential drop between the LB and the sheath entrance) is consistent with the equation of a collisionless sheath potential drop for an electro-positive discharge with an electron temperature of 2 eV 4 . The space charge is dominated by the negative ion species in the vicinity of the aperture because electrons are magnetized and deflected from the aperture by the cusp field of the suppression magnets. The meniscus that is the boundary above which the plasma is no longer quasineutral, corresponds to roughly 1% of the initial ion density ∼ 10 14 m -3 and is located around ∼ 1 mm upstream the front face of the aperture.

A sharp increase of positive and negative ion densities occurs in front of the grid surface, and an ion-ion sheath is formed. The electric potential presents a nonmonotonic varia-tion resulting first from an excess of positive charges and then negative charges within the sheath in contact with the grid surface. A virtual cathode is formed with a potential well whose depth is around 1 V at y = 1.4 cm. The negative ions emitted at the surface and whose energy is below 1 eV are recollected by the surface; more energetic negative ions cross the potential well and can be transported into the plasma region. The potential well in front of the grid along (Oy) varies from -0.8 to -1.6 V. An estimation of the saturation current 109 j sat is given by j sat = j H -e Φ well /T n , where Φ well can be taken as the averaged potential well depth in front of the PG (in a one dimensional approximation). For Φ well ∼ -1.3 V, the saturation current is on the order of 6.8 A.m -2 , close to the calculated value (6 A.m -2 ). The threshold for the appearance of the saturated regime has been recently derived analytically and extended to the case of multiple charge species in the work of Schiesko et al. 143 . By applying -5 V at the LB plane, negative ions move back and forth trapped by the electric potential before being destroyed by detachment processes or extracted. In the calculations shown here, almost all the negative ions are extracted (5 A.m -2 ). Assuming 0 V on the LB instead, one now includes the possibility for the negative ions, which leaves the PG surface (i.e., cross the virtual cathode) and do not undergo any collisions to be collected by the left-handside surface (roughly 50% of the saturated current density in Ref. 127 ). Calculations show that in the quasi-neutral region the negative ion density significantly drops and the electric potential axial profile between the LB plane and the meniscus is flatter; the extracted negative ion current is halved. These properties are consistent with 2D simulations performed for real plasma conditions (after multiplying the scaled values by the appropriate β ; see the work of Fubiani et al. [START_REF] Fubiani | Modeling of negative ion extraction from a magnetized plasma source: Derivation of scaling laws and description of the origins of aberrations in the ion beam[END_REF] for further details). One may, hence, conclude that, qualitatively, the overall understanding of the physical mechanisms responsible for the extraction of negative ions analyzed through Figs. 9 and 10 is still valid 127 irrespective of the LB potential value.

Results shown above were 2D simulations assuming as a consequence that the extraction hole can be represented by a slit of infinite length in the direction perpendicular to the simulation domain. The use of a simpler 2D approach is sufficient to understand the physics of negative ion extraction but these models do not provide the capability to 1) quantify the actual value of the extracted negative ion current density and 2) describe the transport properties of the co-extracted electrons. 3D simulations are necessary to study these two open problems. Calculations performed by the ONIX code assumed an average plasma density of n 0 = 2 × 10 17 m -3 and a j H -= 550 A.m -2 negative ion current density produced on the PG. The half-Maxwellian negative ion flux has a temperature T n = 0.6 eV. The edge of the hole has been chamfered (to accurately model the aperture geometry presented in Fig. 1), and finally, the EG voltage has been fixed to 10 kV. A comprehensive set of physical-chemistry reactions including electrons, multiple positive ion species (H + , H + 2 , H + 3 , and a fraction of Cs + ), H - negative ions, and H and H 2 neutral atoms has been considered. The grid size over electron Debye length is set to 1.5 with a time step of 0.1/ω pe and 30 particles-per-cell, respectively. The electron-ion pair re-injection method introduced in Sect. III has been used. The computational time to reach convergence is about 150 h using 4096 cores. The negative ion current density extracted from the model 144 j H -, ext is on the order of 250 A.m -2 similar in principle to experimental observations but with (1) a negligible amount of particles coming from the flat side of the PG facing the plasma while it represents ∼ 25% of the total surface area, (2) a higher ratio of ions originating from the plasma volume (∼ 35% instead of ∼ 10% 34 ), and (3) the remaining (∼ 65%) coming from the chamfered surface with a high divergence angle, all being in contradiction with measured data in the experiments [START_REF] Fantz | Negative hydrogen ion sources for fusion: From plasma generation to beam properties[END_REF] . The latter also predicts an increase of the extracted negative ion current somewhat proportional to the aperture surface area in a cesiated ion source [START_REF] Kashiwagi | Development of negative ion extractor in the highpower and long-pulse negative ion source for fusion application[END_REF] . The reasons for the discrepancies between the model and the experiments are currently an open problem but one hypothesis is that the virtual cathode profile in front of the hole must be well described numerically, which impose more stringent constraints on the grid size than just resolving the electron Debye length. In the paper of Montel-lano et al. 144 , for instance (see Fig. 16 therein), the potential well is just described by four points. Figure 12 plots a comparison between a beam extracted from a bevel and chamfered aperture geometry 145 for a configuration in which the beam is in a so-called perveance match condition (a plasma meniscus which is slightly concave providing a laminar beam at the exit plane). The model is a 2D PIC-MCC and the apertures are, hence, infinite slits. The simulation domain is 3.2×1.6 cm 2 with an average plasma density of 3×10 17 m -3 , a 1.5 kV extraction voltage, and a current density of 600 A/m 2 for the negative ions produced on the PG together with a 1 eV temperature (half Maxwellian flux). The number of grid points is 4096 × 2048 associated with a minimum resolution of ∆x = ∆y ≃ 0.5 λ De and providing enough points to model the potential well induced by the space charge limited negative ions. The calculation only includes negative ions produced on the PG surface. One finds an extracted ion current ratio of j ca / j b = 1.5, while the surface ratio is S ca /S b = 1.2, where "ca" stands for the chamfered aperture and "b" for the bevel, respectively, which is more in line with experimental data [START_REF] Kashiwagi | Development of negative ion extractor in the highpower and long-pulse negative ion source for fusion application[END_REF] . The cathode depth has a minimum of about -1.3 V in front of the flat side and -2.1 V for the chamfered one (near the tip). In the case of the chamfered geometry implemented in the model, the tilted over flat surface ratio is S t /S f = √ 2, while the corresponding extracted negative ion current ratio is found to be j t / j f ≃ 1.7; hence, the extraction probability is larger for the tilted segment, 58% vs. 33% once an ion escaped the well, but the saturation current is smaller. Note that it does induce the generation of aberrations in the beam optics as observed also in the simulation of Fig. 11. This is not the case with the bevel design. Last, note that in accordance with experiments (see Fig. 5 of Geng et al. 146 ), we observe in Fig. 12 a shift of the x-point from the axis of the aperture which, is due to slight magnetisation of the negative ions by the magnetic filter field. Ions are no longer extracted beyond the x-point.

A typical cartography of the electron density calculated with the 3D ONIX code is shown in Fig. 11 (taken from the work of Montellano et al. 144 ). The magnetic filter field is along (Oy) in this calculation, while the suppression field (cusp profile surrounding each apertures) lies in the xz plane, respectively. The electron density, which reaches ∼ 10 17 m -3 in the injection region reduces to less than 10 16 m -3 about 3 mm from the PG surface, while the negative ion density at the same location is 5 times larger. The co-extracted electron current density j e, ext is ∼ 2 times lower than its negative ion counterpart (∼ 250 A/m 2 ). Typical experimental data for similar conditions find an electron current density [START_REF] Fantz | Negative hydrogen ion sources for fusion: From plasma generation to beam properties[END_REF][START_REF] Christ-Koch | Laser photodetachment on a high power, low pressure RF-driven negative hydrogen ion source[END_REF] between 40 and 200 A/m 2 . The electron density profile is somewhat symmetric in the xz plane of Fig. 11(b), while it is asymmetric close to the aperture in the xy plane, which is a telltale sign of a cross B drift induced by the magnetic field of the suppression magnets. This mechanism drives the electron current outward toward the accelerator [START_REF] Wünderlich | Effects of the magnetic field topology on the co-extracted electron current in a negative ion source for fusion[END_REF]147 . The difference in the plasma density on the PG surface leaves an imprint on the depth of the virtual cathode, which as a consequence, is asymmetric as well as shown in Fig. 13(c). 1D axial density profiles, i.e., along (Ox), of H + , H -, and electrons are shown in Fig. 13 , an ion-ion sheath is formed, and the potential well minimum is on the order of -1.2 V, while (2) at the bottom of the aperture (middle figure) and at the location of the mirror point of the cusp field, the electron and negative ion densities in front of the PG are similar and the depth of the potential is about -0.9 V. Furthermore, there is a left-right asymmetry caused by the electron drift. As a consequence the negative ion saturation current density (the actual current density leaving the surface locally) is non-homogeneous. In addition, the slight magnetization of the negative ions by the magnetic filter field renders also the extracted beam current density profile asymmetric (see Fig. 12). Both may have dire consequences for the beam transport inside the accelerator as the latter has been designed for a homogeneous current density profile 132 .

Last, Wünderlich et al. [START_REF] Wünderlich | Effects of the magnetic field topology on the co-extracted electron current in a negative ion source for fusion[END_REF] have carried out parametric studies on the incidence of the magnetic field for reducing the coextraction of electrons and found that increasing the strength of the suppression magnets would be an effective solution. In addition, simulations performed for pure electro-positive discharges (without any negative ions) in a cusp-shape magnetic field topology have revealed (1) that the electrons are diffusing across the aperture due to an E × B drift and (2) that plasma instabilities in the meniscus region induce a time-varying elec-

• Doubled deflection field: The filter field is the same as in the standard case, but the deflection field strength is increased by a factor 2.

The top part of Fig. 2 shows for these three configurations: the ratio of the deflection field strength and the filter field strength along the center axis of the simulation domain. Close to the PG, the deflection field is strongly dominating. Only when using the doubled filter field, the filter field dominates over the deflection field in the first ≈5 mm at the upstream side of the simulation domain.

The plasma parameters present in the experiment close to the PG are the result of the plasma transport from the driver through a large fraction of the magnetic filter into the calculation domain. Consequently, changing the filter field strength may have an impact on the plasma parameters in the calculation domain. Diagnostics results show, however, a weak dependence only of the electron temperature and electron density on the field strength for B > 1.0 mT. [START_REF] Yu | Negative ion source improvement by introduction of a shutter mask[END_REF][START_REF] Fantz | Plasma expansion across a transverse magnetic field in a negative hydrogen ion source for fusion[END_REF] These observed weak changes have been omitted, and identical input values have been used for the plasma parameters in the runs modeling the three magnetic field cases in order to demonstrate purely the effect of the magnetic field on the particle transport in the vicinity of the extraction aperture.

The current densities of extracted negative ions ( j ex ), co-extracted electrons ( j e ), and the electron-ion ratio ( j e /j ex ) are determined in ONIX by the particle fluxes reaching the plane of the EG. The values of j ex , j e , and j e /j ex predicted by the code for the three magnetic field cases are shown in Table III.

The results for the standard case can be compared directly with the experimental hydrogen data: the extracted current density

Cas

Stan Dou Dou

The top part of Fig. 2 shows for these three configurations: the ratio of the deflection field strength and the filter field strength along the center axis of the simulation domain. Close to the PG, the deflection field is strongly dominating. Only when using the doubled filter field, the filter field dominates over the deflection field in the first ≈5 mm at the upstream side of the simulation domain.

The plasma parameters present in the experiment close to the PG are the result of the plasma transport from the driver through a large fraction of the magnetic filter into the calculation domain. Consequently, changing the filter field strength may have an impact on the plasma parameters in the calculation domain. Diagnostics results show, however, a weak dependence only of the electron temperature and electron density on the field strength for B > 1.0 mT. [START_REF] Yu | Negative ion source improvement by introduction of a shutter mask[END_REF][START_REF] Fantz | Plasma expansion across a transverse magnetic field in a negative hydrogen ion source for fusion[END_REF] These observed weak changes have been omitted, and identical input values have been used for the plasma parameters in the runs modeling the three magnetic field cases in order to demonstrate purely the effect of the magnetic field on the particle transport in the vicinity of the extraction aperture.

The current densities of extracted negative ions ( j ex ), co-extracted electrons ( j e ), and the electron-ion ratio ( j e /j ex ) are determined in ONIX by the particle fluxes reaching the plane of the EG. The values of j ex , j e , and j e /j ex predicted by the code for the three magnetic field cases are shown in Table III.

The results for the standard case can be compared directly with the experimental hydrogen data: the extracted current density of 154 A/m 2 and the electron-ion ratio well below one are in good agreement with the results obtained for a good but not perfect cesium conditioning status (see, for example, Ref. 30) taking into account that the temporal dynamics of long plasma pulses (up to 1 h) is not included in the calculations.

The change of the magnetic field topology by only doubling the filter field has almost no effect on the current densities of both species. In contrast, doubling the strength of the deflection magnets reduces the co-extracted electrons by about a factor of four with only a slight reduction of about 10% on the negative ions. Consequently, the electron-ion ratio is significantly smaller than that in the previous two cases.

The slight reduction of the extracted negative ion current observed for the doubled deflection field can be attributed to the absolute field strength, which is higher than 15 mT at the axial position 5 mm before the PG. From this position and up to the PG, the Larmor radius of the negative ions is smaller than 9 mm, which is by a factor 2 smaller than the simulation domain width. Thus, tron transport 147 , which seems to be attenuated when negative ions are present 145 ; further investigations are needed.

V. ACCELERATION OF NEGATIVE IONS

The acceleration systems must fulfill two fundamental requirements in terms of (i) beamlet divergence (for example, ≤ 7 mrad with a dispersion of ± 2 mrad for the ITER reactor) and (ii) minimizing the heat load on the grids caused by beam halos and secondary particles 148,149 . The beam profile, divergence, and power are generally assessed through calorimetric measurements on a carbon target using an infrared camera.

The target is positioned on the beam path [148][149][150] . The ITER accelerator, for instance, will have 1280 beamlets [START_REF] Hemsworth | Overview of the design of the ITER heating neutral beam injectors[END_REF] with an acceptance of ±10% and was designed at first based on ideal conditions using ray-tracing algorithms 130,133,134,136,138 where the beam is assumed extracted from the ion source with a homogeneous profile. Any deviation from these properties will require modification of the grid arrangement, the aperture shape and acceleration voltages. The ITER concept has been designed following the development of prototypes such as the one in Cadarache testing the SINgle GAP (SINGAP) architecture 149,150 or the so-called Multi Aperture Multi Grid (MAMuG) accelerator in Japan's JT-60U Tokamak 148 . To date the latter was chosen for ITER and a scale 1 prototype of the whole Neutral Beam Injector (NBI), MITICA, is currently being tested in RFX, Padova (Italy) 151 . Recent experiments and numerical simulations have shown a number of critical issues related to the beam transport which must be corrected : (i) a beamlet deflection a factor of 2 higher than predicted by the models using ray-tracing algorithms which is likely caused by an asymmetric current density profile within each negative ion beamlets extracted from the ion source 132,152 (see the discussion in Sec. IV), (ii) the beamlet current variation along the lines of aperture perpendicular to the magnetic filter field probably exceeds the accelerator acceptance due to the electron Hall current inside the ion source, which induces an asymmetry in the plasma parameters [START_REF] Fantz | Negative hydrogen ion sources for fusion: From plasma generation to beam properties[END_REF][START_REF] Fubiani | Modeling of plasma transport and negative ion extraction in a magnetized radiofrequency plasma source[END_REF]153 (including the negative ion species), (iii) the latter also generates an uneven repartition of the co-extracted electron current which is mostly collected by the extraction grid (EG) and might result in large power deposition locally, and finally (iv) the estimation of the absorption of the power delivered by the accelerator's power supplies by secondary particles amounts to tens of percents in the calculations with possibly deleterious consequences for the accelerator parts if the adequate cooling is not taken into consideration [START_REF] Fubiani | Modeling of secondary emission processes in the negative ion based electrostatic accelerator of the international thermonuclear experimental reactor[END_REF][START_REF] Fubiani | Analysis of the two accelerator concepts foreseen for the neutral beam injector of the international thermonuclear experimental reactor[END_REF] . The question of secondary particle production is best studied though a comparison between model and experiments 154 . The former can provide a detailed description of the particle kinetics, while the second can measure the power loads on the different grids.

The EAMCC 3D code [START_REF] Fubiani | Modeling of secondary emission processes in the negative ion based electrostatic accelerator of the international thermonuclear experimental reactor[END_REF][START_REF] Fubiani | Analysis of the two accelerator concepts foreseen for the neutral beam injector of the international thermonuclear experimental reactor[END_REF] has been developed to model secondary particle physics in accelerators where the background neutral gas density remains non negligible besides an efficient pumping as is the case in the NBI foreseen for ITER for instance. The neutrals come either from the ion source itself (leaking though the 1280 apertures) or from the gas neutralizer. Their density is estimated to be on the order of ∼ 10 19 m -3 close to the grounded grid of the ITER-MAMuG and as high as 4 × 10 19 m -3 in the gap between the plasma grid (PG) and the EG [START_REF] Fubiani | Modeling of secondary emission processes in the negative ion based electrostatic accelerator of the international thermonuclear experimental reactor[END_REF] . Secondary particles can either be (i) co-extracted electrons as well as (ii) electrons originating from single or double stripping of negative ions producing alongside fast neutral atoms and D + ions or finally (iii) byproducts of the ionization of the background gas. These charged particles are accelerated by the local electric field and deflected by the magnetic field resulting in a highly divergent distribution which have a significant probability to hit the surface of a grid and, hence, producing in turn more secondary electrons. Ions backscatter as neutrals and, in addition, have a larger secondary electron emission (SEE) yield on impact than fast incident electrons. The problem of secondary particle production is modeled via a Monte-Carlo collision algorithm calculating a statistically relevant number of testparticles trajectories. With the space charge of these particles being negligible compared to one of the negative ion beam, the potential distribution in the accelerator may be calculated by a ray-tracing code without including the secondary particles. Models, such as EAMCC 3D take into consideration the complex particle-surface interactions, i.e., the probability for a particle to be backscattered at the given angle, energy loss inside the material, reflection coefficient, SEE, etc. Figure 14 shows an application of the algorithm to the ITER-MAMuG accelerator. From left to right, the first grid is the PG, the second grid is the EG at 9.4 kV, the other five acceleration grids (AG) provide an energy gain of 200 keV each (a total gain of 1009 keV per negative ion). The power is fixed to 40 MW corresponding to 40 A of negative ion current on the grounded grid with an assumption of one co-extracted electron per ion leaving the ion source (≃ 1.5% of these electrons are transmitted past the EG). One finds that the total power necessary to accelerate the beam must amount to ∼ 53 MW as ∼ 13 MV (∼ 25%) corresponds to parasitic power absorbed by secondary particles [START_REF] Fubiani | Analysis of the two accelerator concepts foreseen for the neutral beam injector of the international thermonuclear experimental reactor[END_REF] . The latter amounts to 9.6 MW deposited on the accelerator grids, while power transmitted toward the neutralizer is of the order of 3 MW (mostly carried by neutrals). 840 kW goes back inside the plasma source (exclusively from positive ions) hitting the back plate of the ion source (this has been observed experimentally). Last, note that secondary electron power transmitted toward the neutralizer is significant, 820 kW. Beamlet halos (a highly divergent sub-fraction of the ion distribution) have been measured in the SINGAP accelerator [START_REF] De Esch | Negative ion beam halo mitigation at the 1MV testbed at IRFM[END_REF] with a current ratio carried by the latter of about 10%. A similar fraction was, hence, added in the simulation of Fig. 14 (14%; the ions were loaded at the back of the PG facing the EG) leading to 2.5 MW out of the 9.6 MW caused by direct hitting of negative ions within the halo with the accelerator grids. Furthermore 30% of the beamlet ions are lost due to stripping reactions inside the volume of the accelerator vessel. The EAMCC model has also been used to simulate the SINGAP concept [START_REF] Fubiani | Analysis of the two accelerator concepts foreseen for the neutral beam injector of the international thermonuclear experimental reactor[END_REF] and recently to calculate the power heat load on the grids of the NIO1 ion source's accelerator operating at Consorzio RFX in Padova, Italy 155 .

VI. SUMMARY AND PROSPECTS

Negative ions are used mainly as precursors for producing energetic neutrals in neutral beam injectors (NBIs) for fusion machines as well as for neutron sources and particle accelerators (with application to high energy physics, medicine, etc.). The ions are generally produced by the so-called low temperature magnetized plasma source (with T e of the order of 1 to 10 eV and T neu ≲ T i < T e , where T neu , T i , T e are the neutral, ion, electron temperatures, respectively). In the case of ion sources for fusion, the devices are typically tandem type with a plasma discharge powered by a RF antenna attached to a larger chamber, which is magnetized (it is the design foreseen for the ITER NBI) or, alternatively, simply a single chamber driven by filaments (JT60-U concept). The role of the magnetic filter field is twofold : (i) to reduce the electron temperature in the extraction region, which is beneficial both for the production of negative ions in the plasma volume through the dissociative attachment of vibrationally excited hydrogen molecules and to increase the mean free path of the negative ions before being destroyed by collisions. (ii) To decrease the electron flux onto the plasma grid (PG), which induces in turn a lower amount of these particles co-extracted alongside the negative ions. The latter are mostly generated by neutral impacts on a cesiated PG surface as the volume process is insufficient to produce the required beam current density (these ions represent ∼ 10% of the total current). One of the reasons is the difficulty to extract them due to the electric potential, which confines negative charges. Producing the ions in the direct vicinity of the aperture allows one to greatly increase their extraction probability. Understanding the inner workings of the ion source and accelerator system requires both experiments as well as numerical modeling. The latter are typically central to describe in detail the plasma behavior (especially near the negative ion emitting PG surface), to design the grid ensemble of the electrostatic accelerator, and to model the beam transport there. The numerical techniques to simulate the negative ion kinetics are typically two-and three-dimensional three-velocities particle-in-cell with Monte-Carlo collisions (PIC-MCC) inside the plasma, ray-tracing, and test-particles (including the physical-chemistry and interaction with surfaces) inside the accelerator. Tentatives to use PIC-MCC models coupling the ion source with the accelerator have been discussed in the lit-erature, but to date, using the explicit PIC technique with fixed grid spacing, requires a large number of grid nodes, leading to high computational cost. The limited gain in physical insight over the use of the conventional ray-tracing models for the accelerator does not outweigh the numerical burden.

The PIC-MCC models of the ion source have shown that the electron Hall current when intercepted by a surface generates an asymmetry in the plasma parameters with repercussions down to the extracted negative ion beam current density profile. This has been observed also experimentally proving the usefulness of the combined experiment/modeling analysis of physical phenomena to understand the mechanisms at play. Experiments cannot assess the complex sheath behaviour near the PG surface where negative ions are emitted. PIC-MCC models have shown that the negative ion density is in a space charge limited regime which induces the generation of a potential well (virtual cathode) in front of the surface. The well depth may be limited by the canceling effect on the charge density of the incoming positive ions from the plasma volume. The mere existence of a well implies that not all the negative ions produced on the surface will escape the sheath, only a fraction that is linked to a saturation current density j sat . The latter is correlated with the residence time of the ions in that area. Numerically it is of paramount importance to have sufficient numerical resolution to describe the potential drop as j sat is an exponential function of the well depth for a given negative ion temperature. The well properties seem to depend mostly on the negative ion density profile and, hence, not on the electron Debye length, although this is currently an open question. Assuming this turns out to be indeed the case (see the numerical resolution in 2D used in Fig. 12, for instance), then it would prohibit any simulations in 3D due to the necessary implementation of a significantly large number of grid nodes. Nevertheless, PIC-MCC calculations restricted to a simulation domain around a single aperture with periodic boundary conditions have demonstrated that the extracted negative ion current density profile within a beamlet itself is asymmetric due (i) to the slight magnetization of the ions near the PG by the filter field and (ii) likely to a nonhomogeneous virtual cathode profile along the grid surface caused by the electron cross B drifts inside the plasma. This asymmetric intra-beamlet profile has been posited to explain the higher deflection inside the accelerator observed experimentally. In addition, the model has shown that the extracted current seems somewhat proportional to the surface area of the aperture facing the plasma (comparison between a chamfered and a bevel grid geometry). This has also been observed in the experiment. Note that the reported simulation results are the sole work of one author to date and should, hence, be confirmed by other researchers before conclusion. The use of a chamfered aperture generates aberrations in the beam optics, which are posited to be at least one of the reasons for the generation of the beam halos. The other one is the migration of cesium to the back of the PG facing the accelerator vessel. Electrons are co-extracted alongside the ions from the plasma source, and 3D PIC-MCC calculations have shown that it is caused by an electron cross B drift near the meniscus driving the current across the aperture beside the implementation of a cusp magnetic field profile to act as a barrier. The drawback of such a 3D numerical configuration is the use of periodic boundary conditions which are equivalent to modeling an infinite amount of apertures. As a consequence, the influence of the extraction potential on the simulated plasma is possibly artificially large and caused by a much higher than real aperture over the wall surface ratio in the calculation. This implies that an electron has likely a larger probability to be extracted in this case (this has implications also on the calculated curvature radius of the plasma meniscus for a given extraction voltage). It is currently an open problem. Last, beamlet transport inside the accelerator using test-particles with Monte-Carlo collisions (assuming at the present time an ideal beamlet without any asymmetry in the current density profile) has demonstrated the paramount impact of the secondary particle production leading to a parasitic absorption of ∼ 25% of the power supply power in the ITER NBI which is then redistributed mostly on the grid surfaces (∼ 10 MW if one assumes a 14% halo ratio which is in line with experimental data). The field map was calculated by a ray-tracing algorithm.

Future directions beside working on the open problems described above would be to combine with a ray-tracing algorithm the negative ion phase space distribution at the meniscus surface calculated by the PIC-MCC models of the extraction region. One of the current weaknesses of the ray-tracing technique is the limited description of the meniscus profile (the plasma is not included in the calculation). This would allow to simulate the effect of an asymmetry in the ion current density profile as well as the transport of the beam halo. The improved field map would be used to model secondary particle production. Another incremental study would be to implement the particle flux distribution functions at a given location in front of the PG deduced from a PIC-MCC model of the whole ion source as a boundary condition for the higher resolution calculation provided by a simulation domain restricted to one aperture.

FIG. 1 .

 1 FIG. 1. Schematic view of the ion source setup and the RF circuit of the prototype source at the BATMAN Upgrade testbed 30 . Reproduced with permission from Fantz et al., Front. Phys. 9, 09651 (2021). Copyright 2021 Frontiers Media SA.

Fig. 1 .

 1 The role of cesium is to lower the work function of the metal (molybdenum in the case of BATMAN) to favor the production of negative ions via atom/molecule or positive ion impacts on surfaces by capturing one or two electrons (see, e.g., Refs. 24-26) by a quantum tunneling process. The main region of interest for the generation of negative ions with a large extraction yield is close to the grid apertures (socalled extraction region).

FIG. 2 .

 2 FIG. 2. (a) Contour plot of the z component of the magnetic filter field at the position of the plasma grid apertures. (b) Profiles of the filter and electron deflection field in the x-direction through a beam axis. The sign of the deflection field depends on the individual magnet row. Dashed line: position of the plasma grid apertures. Solid line: position of the permanent magnet boxes for the filter field 31 . Reprinted with permission from Gutser et al., Plasma Phys. Control. Fusion 51, 045005 (2009).Copyright 2009 IOP Publishing Ltd.

FIG. 4 .

 4 FIG. 4. Sketch of the grids of the conceptual acceleration system of the DEMO NBI 66 . Adapated from Sanoto et al., Nucl. Fusion 57, 056026 (2017), Copyright 2017 IOP Publishing Ltd.

FIG. 5 .

 5 FIG. 5. Flowchart of the electrostatic PIC-MCC algorithm.

FIG. 6 .

 6 FIG. 6. 2D schematic view of the simulation domain. Lengths are indicated in mm. B F and B D are the filter (perpendicular to the simulation plane) and suppression magnetic fields (inside the simulation domain), respectively. The gray region corresponds to the location where particles are artificially injected. The left-hand-side boundary, the plasma grid, and the first extraction grid are shown as LB, PG and EG, respectively 77 . Reproduced with permission from Garrigues et al., J. Appl. Phys. 120, 213303 (2016). Copyright 2016, AIP Publishing LLC.

FIG. 7 .

 7 FIG. 7. Plasma potential profiles obtained with the two injection methods in a simple one-dimensional collisionless and nonmagnetized discharge with grounded walls, with only electrons and H + 2 . Charged particles are injected in the region indicated by the two vertical dashed lines 78 . Reproduced with permission from Garrigues et al., Nucl. Fusion 57, 014003 (2017). Copyright 2017 IOP Publishing Ltd.

FIG. 8 .

 8 FIG. 8. Profiles along the axial direction at y = 0 for grid spacing varying from 21.3 to 250 µm of (a) electron temperature and (b) electric potential. Charged particles are injected in the hatched region. The initial number of electrons and positive ions is 40 particles-percell 77 . Reproduced with permission from Garrigues et al., J. Appl. Phys. 120, 213303 (2016). Copyright 2016 AIP Publishing LLC.

FIG. 9 .

 9 FIG. 9. 2D profiles of (a) the electric potential Φ, (b) the positive ion density, and (c) the negative ion density for an initial plasma density n 0 = 10 16 m -3 uniform in the left half simulation domain, a negative ion current density j H -= 25 A.m -2 along the grid surface, and an EG voltage of 700 V. The voltage on the LB plane is set to -5 V. The plasma potential is in volts and the unit for the positive and negative ion densities, n i and n n is, 10 16 m -3 (the large values of the negative ion density along the grid surface are not visible). The plotted potential contours above 10V are 70, 140, 210 up to 700V 127 . Reprinted with permission from Boeuf et al., Plasma Sources Sci. Technol. 25, 045010 (2016). Copyright 2016 IOP Publishing Ltd.

  FIG. 10. Axial distribution of the charged particle densities (electrons, n e , positive ions, n i , and negative ions, n n ) and plasma potential Φ in the conditions of Fig. 9 at two different transverse positions, (a) at y = 0.8 cm, i.e., along the axis of the grid aperture, and (b) at y = 1.4 cm between the grid surface and the left boundary (along the top dashed line displayed in Fig. 9a). The front PG surface is at x = 2 cm. The maximum of the negative ion density near the grid surface in Fig. 10b (not apparent on the figure) is on the order of 2 × 10 16 m -3 from Ref. 127. Reprinted with permission from Boeuf et al., Plasma Sources Sci. Technol. 25, 045010 (2016). Copyright 2016 IOP Publishing Ltd.

FIG. 11 .

 11 FIG. 11. Electron density (top) horizontal y -axial x plane at z = z max /2, (bottom) vertical z -axial x plane at y = y max /2 obtained with the 3D ONIX code. The white line is the position of the meniscus. The magnetic filter field is directed along (Oy), while the suppression field lies in the xz plane 144 . Reprinted with permission from Montellano et al., J. Phys. D:Appl. Phys. 52, 235202 (2019). Copyright 2019, IOP Publishing Ltd.

FIG. 12 .

 12 FIG. 12. 2D normalized negative ion current density profiles for (a) bevel vs (b) chamfered slit apertures. 2D PIC-MCC calculation. The maximum ion current density j max is 200 A/m 2 . The current streamlines are shown in green and the plasma meniscus in white 145 . From G. Fubiani, private communication.

FIG. 6 .

 6 FIG. 6. (a) Density of the positive ions, negative ions, and electrons as a function of the dis zontal position 1 mm, vertical position 10 mm) and on the bottom of the aperture (horizontal p a function of the position on the PG surface.

FIG. 6 .

 6 FIG. 6. (a) Density of the positive ions, negative ions, and electrons as a function of the distance from the PG along the axial direction on the right of the aperture (horizontal position 1 mm, vertical position 10 mm) and on the bottom of the aperture (horizontal position 10 mm, vertical position 1 mm). (b) Depth of the potential well f well as a function of the position on the PG surface.

FIG. 13 .

 13 FIG. 13. Density of the positive ions, negative ions, and electrons as a function of the distance from the PG along the axial direction (a) on the right of the aperture and (b) at the bottom, respectively. (c) Minimum of the virtual cathode Φ well vs the location on the front PG surface 76 . Reprinted with permission from D. Wünderlich et al., J. Appl. Phys. 130, 053303 (2021). Copyright 2021 IOP Publishing Ltd.

FIG. 14 .

 14 FIG. 14. Geometry of the MAMuG negative ion based electrostatic accelerator. The plot shows a zoom over one beamlet. From left to right: plasma grid (PG), extraction grid (EG) (at an extraction voltage of 9.4 kV), and acceleration grids (AG) 1-5 (total acceleration voltage of 1 MV). Primary and secondary particles are shown; negative deuterium ions (D -) (red color), neutrals (D 0 ) (green), positive deuterium ions (D + ) (blue), positive deuterium molecular ions (D +2 ) (purple), and electrons (black). The negative ion beam aims downward with an average divergence <y'> ≈ 5.5 mrad induced by the PG magnetic filter field (generated by a 4 kA current)[START_REF] Fubiani | Modeling of secondary emission processes in the negative ion based electrostatic accelerator of the international thermonuclear experimental reactor[END_REF] . Reproduced with permission from G. Fubiani et al.,PRAB 11, 014202 (2008). Copyright 2008 American Physical Society.

  used in Refs. 59, 61, 131, and 132, and SLACCAD 133 (see, e.g., Refs. 43, 44, 64, and 132) KOBRA-3D 134 (Ref. 135), BTR 136 used by Dlougach et al. 137 , and IBSIMU 138 . A first utilization of the COMSOL Multiphysics® 139 software has been recently proposed in Ref. 132. Secondary particle production inside the accelerator is treated using MCC algorithms for fixed background gas density profiles 43,44 . The different approaches have been summarized in the article of Denizeau et al. 132 .

TABLE III .

 III Effect of the magnetic field variation on the extracted negative ion current density and the co-extracted electron current density, as well as its ratio. Taken from Ref.47.

	Journal of Applied Physics
	Case	j ex (A/m 2 )	j e (A/m 2 )	j e /j ex
	Standard	154	21.2	0.14
	Double filter field	154	19.5	0.13
	Double deflection field	145	5.2	0.04

J. Appl. Phys. 130, 053303 (2021); doi: 10.1063/5.0054949 © Author(s) 2021• Doubled deflection field: The filter field is the same as in the standard case, but the deflection field strength is increased by a factor 2.

ACKNOWLEDGMENTS

Part of this work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. Support from CEA and from the French Fédération de Recherche sur la Fusion Magnétique is also acknowledged. In addition, some of these projects were granted access to the HPC resources of CALMIP supercomputing center under the allocation 2013-P1125. The authors are very grateful to their colleagues J. P. Boeuf and G. Hagelaar, from LAPLACE, Toulouse and to A. Simonin, CEA, Cadarache for his continued support.

AUTHOR CONTRIBUTIONS

L. Garrigues: Writing -original draft (equal). G. Fubiani: Writing -original draft (equal).