In negative ion sources for the ITER Neutral Beam Injection system, the co-extraction of electrons is one of the main limiting factors. The current of co-extracted electrons can be decreased by applying a positive bias voltage to the Plasma Grid (PG) with respect to its source walls. Simulations using three-dimensional Particle-in-Cell Monte Carlo Collision (3D-PIC MCC) model are a powerful tool for studying the extraction region of such ion sources. However, the inclusion of both PG and source walls in the simulation domain is difficult due to numerical constraints. This study uses the 3D-PIC MCC code ONIX to explore the effects of particle injection models on plasma characteristics, using a flux injection model to regulate particle influx for a flat transition in potential from the bulk plasma to the simulation domain. Biasing of the PG above floating potential is possible using the flux injection scheme and results in a notable reduction in co-extracted electrons, corroborating with established experimental observations.

The ITER Neutral Beam Injection (NBI) requires the acceleration of H/D ions up to 1 MeV.1 The negative ions (NIs) for the ITER NBI will be produced in large ceasiated radio frequency (RF) driven ion sources and extracted through a multi-grid, multi-aperture system of 1280 apertures arranged in 16 groups of 5 × 16 apertures. The NIs will be extracted from the plasma through the apertures of a plasma grid (PG) by the application of an extraction potential of ≈10 kV between the PG and the extraction grid (EG). Along with NIs, electrons are co-extracted and dumped onto the EG by applying a magnetic deflection field (BDF) perpendicular to the beam direction. One of the main limiting factors for the achievable source parameters is the amount and temporal stability of the co-extracted electrons; ITER will therefore require the current of co-extracted electrons to be below that of NIs. To limit the destruction of NIs by electron detachment (ED) collisions with electrons, a horizontal magnetic filter field (BFF) with a typical field strength of several mT is applied.2 In addition, the BFF reduces the co-extraction of electrons by decreasing the electron density and temperature near the PG.

Two prototype NI sources are in operation at Max-Planck Institute for Plasma Physics: BATMAN Upgrade (BUG) and Extraction from a Large Ion Source Experiment (ELISE), which have a source body area ( width × height) of 0.6 × 0.3 and 1.0 × 0.9 m2, respectively. Experimental investigations show that applying a positive voltage to the PG relative to the other surfaces, such as the source walls and bias plate (BP), can produce a positive potential difference between the plasma and the PG, i.e., an electron-attracting sheath and, consequently, a substantial decrease in the co-extracted electron current je. Depending on the potential difference, the extracted negative ion current jex can be decreased.3–5 

In the absence of bias and negative ions, the potential difference between plasma and wall is determined by the electron temperature Te, the mass ratio of electrons to positive ion (PI) m e / m PI, and the ratio of PI and electron temperature T PI / T e. The plasma potential ϕpl is given by the analytical models by Stangeby6 and Schwager and Birdsall7 
ϕ pl = ϕ w k B T e 2 · ln [ 2 π m e m PI · ( 1 + T PI T e ) ] ,
(1)
where ϕw is the wall potential and kB is Boltzmann's constant.

The extraction system in ITER-like ion sources has a complex geometry, and the two magnetic fields BDF and BFF are crossing each other in the extraction region.8 During extraction, an equipotential surface, called meniscus, forms between the quasi-neutral plasma on the upstream side of the PG and the extracted NI beam. The meniscus acts as an emission boundary; every negatively charged particle passing through the meniscus is accelerated out of the extraction region while PIs are reflected.

The plasma close to the PG (the extraction region) has a density of approximately 1017 m−3 and an electron temperature Te of 2 eV, which both are one order of magnitude lower than in the driver region.9 The negative charges are equally composed of NIs and electrons, leading to the formation of an ion–ion plasma near the meniscus. Additionally, protons have a mean free path of approximately 10 cm in the expansion region, which decelerates the proton flux from drivers toward the PG.

To model the plasma close to the PG (the extraction region), the meniscus, the extraction of NIs, and co-extraction of electrons, the topology of the magnetic field, the interaction between particles, source walls, and the grid system need to be considered. These physical phenomena are characterized by intricate and inherently three-dimensional (3D) dynamics (such as the particle interaction with the magnetic field), making analytical methodologies unsuitable for comprehensively examining the extraction mechanisms.

Fluid models can simulate the entire ion source but assume a full Maxwellian distribution, missing details such as single particle trajectories, the extraction of negative ions, and co-extraction of electrons. As a result, deploying predictive numerical modeling using 3D Particle-in-Cell (3D-PIC) Monte Carlo Collisions (MCC)10 is indispensable for studying NI extraction. Simulations with 1D Particle-In-Cell (PIC) using a Neumann condition on the plasma-facing side of the domain have been performed;11 however, PG biasing has not been incorporated in 3D-PIC simulations. Previous 3D-PIC MCC studies12–17 have neglected the effect of biasing the PG with respect to the source walls. A 2D Comsol-based fluid simulation found that BP biasing affects the electron flux in the ion source, but the effect of PG bias on electron co-extraction was not explored.18 

This article examines injection methods and boundary conditions in 3D-PIC MCC modeling using the code Orsay Negative Ion eXtraction (ONIX),16 focusing on its ability to simulate PG bias in large ion sources for ITER NBI and in a simplified slab-domain. In particular, the impact of two injection methods on plasma density, potential, and the extraction of ions and electrons is investigated near the extraction region in the ELISE ion source. Additionally, the paper reports the first application of flux injection for PG biasing in 3D-PIC MCC simulations, which is an essential technique for reducing the co-extraction of electrons in experiments. The insights contribute to advancing knowledge in plasma physics and optimizing NBI systems.

3D-PIC modeling is a numerical method used extensively in plasma physics to self-consistently solve problems involving charged particles and electromagnetic fields. In the PIC model, the computational domain is discretized into a grid or mesh on which Poisson's equation for the electrostatic potential field is solved numerically. The particles in the PIC method represent not single electrons or ions but rather a collection (macroparticles), which move in the electromagnetic fields calculated on the grid points. The Lorentz equation governs the motion of the particles.

The use of PIC modeling reduces the computational complexity compared to calculating the particle–particle forces directly. Given the need to adequately resolve the Debye length with a substantial number of macroparticles within each cell, 3D-PIC simulations demand both a significant particle count and a large grid size. Solving the Poisson equation on such a large grid and moving each macroparticle require the use of highly parallelized programming and substantial RAM usage. As a result, these simulations are primarily feasible on supercomputers and necessitate strong limitations on domain size. The 3D simulation domain typically only covers the plasma in the vicinity of one of the PG apertures and the extracted beamlet up until the EG due to these numerical constraints. The simulation boundaries are thus far from the source walls. Since a simulation of the entire source volume is not feasible due to constraints in computational cost, the source and the extraction aperture cannot both be directly included in a 3D-PIC simulation without the scaling of physical parameters.

ONIX was developed to simulate negative ion extraction in sources that are dominated by surface production of negative ions, i.e., for caesiated H/D sources.16,19–21 The ONIX code self-consistently simulates particle transport in a hydrogen or deuterium plasma, including dissociation, the effect of collisions, and the influence of the external magnetic field. Negative ions in caesiated ion sources for fusion are predominantly (≈90%) produced via surface production by the conversion of neutral atoms and PIs.3 The remaining ≈10% are produced via volume production. In ONIX, both Volume Produced H ions ( H VP ) and Surface Produced H ions ( H SP ) are included in the model. However, to reach the negative ion density obtained by Cavity Ring-Down Spectroscopy (CRDS) measurements in ELISE during operation with Surface Production (SP),22 additional negative ions are injected in the simulation. When using pair injection, these ions are injected at the start of the simulation, and when using flux injection, they are introduced as an influx from the injection plane (the pair and flux injection are described in Sec. II A). These negative ions are labeled as H VP ions, although in the real ion source, the majority of them are surface produced and have accumulated in the bulk plasma volume. A detailed description of ONIX can be found in Refs. 16 and 23–25 and validation in Refs. 16 and 26. ONIX is a massively parallelized code using Message Passing Interface (MPI) and a domain and particle composition technique.23 

Studies that neglect the bias of the PG12–17 treat the PG boundary and the upstream plasma-side boundary as walls with Dirichlet conditions, where a Debye sheath is formed toward each boundary. In those setups, a change in the PG bias also shifts the plasma potential, such that the PG is still on floating potential.

Different injection methods for PIC models allow for the modeling of various physical phenomena. Typically, a pair injection model [see Fig. 1(a)] is used. In pair injection, the simulation starts with a uniform distribution of electrons, negative ions, and PIs (H+, H2+, H3+, Cs+) injected in a limited volume, a so-called injection region. H VP are re-injected when exiting the simulation domain, and positive ions are re-injected together with an electron when exiting the simulation domain. Electrons and H SP are lost when exiting the simulation domain.

FIG. 1.

The ONIX simulation domain for (a) pair injection, (b) two-sided flux injection, and (c) one-sided flux injection. Γ i is the flux of particle species i (PIs, NIs, or electrons).

FIG. 1.

The ONIX simulation domain for (a) pair injection, (b) two-sided flux injection, and (c) one-sided flux injection. Γ i is the flux of particle species i (PIs, NIs, or electrons).

Close modal

The pair injection scheme allows the user of the code to define the plasma density in the injection region at the start of the simulation. Since the number of positive ions is kept constant while the plasma spreads out in the simulation domain, the overall density in the simulation domain will be lower than the value used as the initial density in the injection region. However, as a direct consequence of the pair injection model, the fluxes across its boundaries cannot be imposed, which restricts the ability of the user to simulate plasmas with biased surfaces. In addition, to avoid simulated electron cooling in the pair injection model, i.e., the nonphysical loss of high-energy electrons to one of the boundaries of the simulation domain, electrons are artificially thermalized in the injection region to keep Te stable over the course of a simulation, as described in Ref. 26.

The second method of particle injection involves injecting electrons, negative ions, and positive ions from an injection plane, adapting the 1D implementation in Ref. 27. Two versions of the flux injection model are used: the two-sided flux injection model [see Fig. 1(b)], where the particles are injected in both perpendicular directions to the injection plane, and one-sided flux injection, where the particles are only injected in one direction, and the injection plane boundary coincides with the boundary of the simulation domain [see Fig. 1(c)]. To keep the charge density continuous across the boundary of the one-sided flux injection scheme, mirrored particles on the left side of the injection plane are used when calculating the charge density. Positive ions passing through the boundary are deleted, while electrons are reflected and thermalized. A Neumann condition can be set on the injection plane, but this does not allow for simulations where the potential difference between the plasma and wall is defined. To simulate a condition near floating potential with the one-sided injection scheme, one boundary of the simulation (ϕw) is set to 0 V and the other to ϕpl, as given in Eq. (1). This approach is valid even in the presence of negative ions as the ϕpl matches well with Eq. (1) in 3D-PIC simulations.16,26,28

The flux injection method allows for a larger self-consistent region of the simulation domain since there is no injection region where electrons are thermalized. Additionally, the fluxes of charged particles across the boundary plane can be altered to ensure the plasma conditions required by the user so that the transition from the bulk plasma to the simulation region is more realistic.

The influx of positive ions Γ PI ( t ) in the flux injection model is regulated by a Proportional-Integral-Derivative (PID)29 
Γ PI ( t ) = P · e ( t ) Γ P + I · 0 t e ( t ) d t Γ I + D · d e ( t ) d t Γ D ,
(2)
where P, I, and D are constants and e(t) is the error at time t, given by
e ( t ) = n PI ( t ) n PI , S S ,
(3)
where n PI ( t ) is the density of positive ions in the regulation zone at time t and n PI , S S is the input value for the steady-state density of positive ions. The densities are taken in a regulation zone located 5 λ De from the injection point (see Fig. 1). The influx of negatively charged particles is equal to that of PIs.
Although the particles are injected from an injection plane, as ONIX is a 3D-3V PIC code, each particle has velocity components in three directions. The velocity components parallel to the injection plane are sampled from an ordinary Maxwellian distribution, while the axial velocity of each species i is selected by random sampling from a Maxwellian flux velocity distribution
f i ( v ) = | v | · ( m i k B T i ) · exp [ m i v 2 2 k B T i ] ,
(4)
where v is the velocity, and m i and T i are the mass and temperature of the species, respectively. In the real ion source, the flux direction of charged particles, created mainly in the volume of the drivers, toward the upstream simulation boundary can be affected by effects such as plasma drifts caused by the BFF. However, as such drifts in the ion source vary strongly by position and are not yet quantitatively determined, no additional drifts are taken into account in the model.

For a standard case of comparison between the pair injection and flux injection model, a 3D-PIC simulation was performed for the ELISE PG geometry.30 The simulations are done both to show the impact of the injection method on the particle distributions in the plasma and that the flux injection model can be used to model the biasing of the PG.

Figure 2 shows the 3D simulation domain. The simulation domain is 20 × 20 × 30 mm3 (i.e., the domain width equals approximately 2% of the ELISE ion source width) and covers one aperture. The aperture is taken as the center aperture in one of the aperture groups, far from the source walls, with periodic boundary conditions in the horizontal and vertical directions. The mesh size is equal to the Debye length of the electrons and the number of macro particles per cell is ≈25. The negative ions are extracted through the meniscus of an aperture of the PG in the axial direction. In this paper, the meniscus is defined by the contour surface where the PI density is three orders of magnitude lower than in the plasma. The simulation domain begins in the bulk plasma and ends at the bordering plane of the EG.

FIG. 2.

(Top) Schematic view of the ELISE ion source and (bottom) ONIX simulation domain for the simulation of one PG and EG aperture in ELISE.

FIG. 2.

(Top) Schematic view of the ELISE ion source and (bottom) ONIX simulation domain for the simulation of one PG and EG aperture in ELISE.

Close modal

A 3D map of the magnetic field configuration is taken into account. This configuration includes the BFF, generated by an electric current of 2.5 kA through the PG and oriented in the horizontal direction, and the BDF is generated by SmCo magnets (VACOMAX 225 HR) with a remanence of 1.10 T embedded in the EG.31 The maximum vertical field in the extraction gap is 61.6 mT. At a distance of 5 mm from the PG into the plasma, it is ≈7.7 mT and thus the electrons are magnetized with a Larmor radius of ≈0.6 mm.

The direction of the BFF and BDF is shown in Fig. 2. The extraction voltage U ext = 10 kV w.r.t. the PG potential. The boundary plane of the downstream exit of the simulation domain is interpolated for a vacuum condition.

Figure 3 shows the axial-vertical slices of the 3D simulation domain for the pair and flux injection simulations, based on the ELISE geometry. The origin of the coordinate system is marked with green arrows.

FIG. 3.

Axial-vertical cross section of the ONIX simulation domain for simulating the extraction of negative ions in the ELISE ion source using (a) pair injection and (b) flux injection. The origin or the coordinate system is located at the center of the downstream boundary and marked with green arrows.

FIG. 3.

Axial-vertical cross section of the ONIX simulation domain for simulating the extraction of negative ions in the ELISE ion source using (a) pair injection and (b) flux injection. The origin or the coordinate system is located at the center of the downstream boundary and marked with green arrows.

Close modal

For the pair injection case, the particles are injected in a 7 mm wide injection region. The simulation boundary at −30 mm and the PG potential ϕPG are set to 0 V; this is equivalent to assuming that there is no bulk plasma for the axial position <−30 mm. In the flux injection case, there is an influx of each particle species denoted by Γ i. The simulation boundary toward the bulk plasma is set to 0 V and the PG potential ϕPG to −5 V, which is close to the potential difference obtained in simulations using pair injection with an electron temperature of Te = 2 eV and expected from theory [see Eq. (1)].

The particle parameters for both cases are based on recent experiments in the ELISE and BUG test facilities and described in detail in Ref. 28. The isotope used for the simulations is Protium, denoted as H. The total density of positive ions (H+, H2+, H3+) is 1.2 × 10 17 m−3 with a ratio of 40:40:20 of each positive ion species, respectively. The electron to negative ion ratio is approximately 1:1 in the bulk plasma, i.e., n e = n H = 6 × 10 16 m−3 and the electron temperature T e = 2 eV. The emission rate of H SP , Γ H SP is 200 A∕m2 and their initial temperature is 0.6 eV.

The point of the slab domain simulation is to try the three particle injection methods in a simple case without the impact of the 3D-magnetic field.

To study the effect of the particle injection methods in a simple case without the impact of the 3D-magnetic field, initial results are presented in a slab-domain of length 0.5 mm (x-direction). The width and height (y- and z-direction) of the domain are 0.052 mm. Since the simulation domain is periodic in the y- and z-directions, the simulation is effectively 1D. Some basic results of a simulation of electrons and H+ ions (H ions are neglected) based on the two-sided flux injection are depicted in Fig. 4. The temperatures T e and T H + are 1 and 0.8 eV, respectively, and the steady-state positive ion density n + , S S is 5 × 10 17 m−3. The figure shows time traces over the first 100 000 iterations, equivalent to 0.5 μs for the PID components: ΓP, ΓI, and ΓD (top), the error in the density e(t) for both electron and H+ species (middle), and total number of macroparticles for both electrons and H+ (bottom).

FIG. 4.

Visualization of the flux injection scheme over 100 000 iterations, equivalent to 0.5 μs. The subfigure (top) shows the PID components' evolution, ΓP, ΓI, and ΓD. The middle subfigure depicts e(t) for both species. The bottom subfigure shows the total number of macroparticles for both species.

FIG. 4.

Visualization of the flux injection scheme over 100 000 iterations, equivalent to 0.5 μs. The subfigure (top) shows the PID components' evolution, ΓP, ΓI, and ΓD. The middle subfigure depicts e(t) for both species. The bottom subfigure shows the total number of macroparticles for both species.

Close modal

The impact of the derivative component in Eq. (2) is set to zero (D = 0) to avoid noise amplification. However, the integral term is needed to avoid oscillations. The proportional term tends to zero as n + ( t ) approaches n + , S S, after 0.05 μs, meaning Γ I Γ P.

At the start of the simulation, the system begins from an unpopulated state, resulting in e ( t ) = 100 % for both species. After approximately 0.05 μs, the error stabilizes, approaching near 0%, implying that the simulated state converges toward the desired system state. This stabilization of the error function is mirrored in the ΓP component, which also approaches zero after 0.05 μs.

Despite an equal injection rate of electrons and H+ ions, the total number of H+ ions surpasses that of electrons. This discrepancy arises due to the positive space-charge in the Debye sheath, leaving an excess of H+ ions in the simulation domain. Stability is reached after 0.75 μs, with a slight delay after the flux is stable since more time is needed for the particles to spread out across the full simulation domain. In summary, the flux injection scheme reaches a converged particle density quickly.

The impact of biasing of the right wall on the plasma potential was investigated by simulating a plasma with T e = 1 eV, T H + = 0.8 eV, and n + , S S = 10 17 m−3 for varying bias of the right wall ϕrw between −10 and +10 V. The potential of the left wall ϕlw was 0 V in all cases.

Figure 5 shows the potential profiles for three simulations with ϕrw set to −10, 0, and +10 V, respectively. The injection plane, the regulation zone, and the detection zone are indicated. On each side of the simulation domain, a Debye sheath is formed toward the wall. When both walls are set to ϕ lw = ϕ rw = 0 V, the sheath is equally deep on both sides, approximately 2.5 V, as given in Eq. (1). The ϕrw has a large impact on ϕpl, in particular, when biasing the wall positively with respect to the floating potential, as this increases ϕpl almost by an equal amount.

FIG. 5.

Electrostatic potential for simulations with varying ϕrw. The injection plane is marked with a dashed red line, the regulation zone is marked with an orange box, and the detection zone where the plasma potential is evaluated is marked with a red box.

FIG. 5.

Electrostatic potential for simulations with varying ϕrw. The injection plane is marked with a dashed red line, the regulation zone is marked with an orange box, and the detection zone where the plasma potential is evaluated is marked with a red box.

Close modal
Figure 6 shows the x-, y-, and z-components of the kinetic energy ( E kin x , E kin y , E kin z) for electrons and H+ ions. The values are taken in the center of the “detection zone,” located between the injection plane and the right-hand side wall for the ϕ lw = ϕ rw = 0 V case. The sign of each kinetic energy component is preserved in the distribution, such that E kin x , E kin y , E kin z are defined as follows:
E kin dir = sgn ( v dir ) · m v dir 2 ,
(5)
where sgn ( x ) is the sign function and v dir is the velocity component in one direction (x, y, or z). For the electrons, the y- and z-components of the velocity distribution function (VDF) correspond well to a Maxwellian distribution with 1 eV, indicating that neither numerical heating nor cooling has occurred. In the x-direction, the distribution is cut at −2.5 eV; the electrons with higher energy in the x-direction will overcome the potential barrier in the plasma sheath and be collected on the right side. For the H+ ions, the horizontal and vertical components of the VDF also correspond well to a Maxwellian distribution with 0.8 eV. A drifting Maxwellian of 1 eV is observed in the x-direction, which is caused by the so-called “source sheath” between the plasma and the injection plane.7,27 Since all H+ ions on the right-hand side of the injection plane will be collected at the wall independent of their energy in the x-direction, no ions are traveling in the opposite direction, and thus E kin x > 0 for all H+ ions to the right of the injection plane. In summary, the flux injection scheme correctly reproduced a Maxwellian velocity distribution in the plasma.
FIG. 6.

Distribution of the x, y, and z-components of kinetic energy ( E kin x , E kin y , E kin z) for (top) electrons and (bottom) H+ ions in the center of the plasma simulation volume.

FIG. 6.

Distribution of the x, y, and z-components of kinetic energy ( E kin x , E kin y , E kin z) for (top) electrons and (bottom) H+ ions in the center of the plasma simulation volume.

Close modal

The differences between plasma potential ϕ plasma and the potentials on the left wall and right wall are shown in Fig. 7 as a function of ϕ rw. For a right-hand side, bias below −2.5 V, ϕpl − ϕlw is unaffected. As the right-hand side, bias approaches and passes 0 V, the ϕpl − ϕlw compensates the potential and follows the increase in the ϕrw bias. In simulations where the particles are injected into the plasma, and not from a boundary, the fluxes toward the boundaries adjust to ensure flux balance and quasi-neutrality. Thus, the bias between the PG and plasma observed experimentally cannot be simulated with a two-sided injection scheme. Therefore, to simulate the bias of the PG in 3D, the one-sided flux injection scheme is used.

FIG. 7.

Difference between the plasma potential and the potential on the left and right walls as a function of ϕrw.

FIG. 7.

Difference between the plasma potential and the potential on the left and right walls as a function of ϕrw.

Close modal

1. Particle densities

For the 3D domain around one aperture of ELISE, a steady-state is reached after a simulation time of ≈ 2μs. The extracted currents are stable in time, and a meniscus has formed. Figure 8 shows the negative ion density at the end of the simulation in the central axial-vertical plane, indicated in Fig. 3, for the pair injection and flux injection simulations, both with T e = 2.0 eV. The injection region is marked in the pair injection case, where negative ions are clearly accumulating. The meniscus, here defined by the contour plane where the positive ion density is reduced by three orders of magnitude compared to the injection zone, is indicated for both simulations in yellow. In the pair injection case, negative ions are accumulated in the injection zone, while in the flux injection case, the highest density is found closer to the meniscus.

FIG. 8.

Negative ion density in the axial-vertical central plane (see Fig. 3) for the (top) pair injection and (bottom) flux injection. The injection zone is marked with orange dashed lines, the meniscus with a yellow line, and the plotting line along which the quantities in Fig. 9 are plotted is marked with a purple line.

FIG. 8.

Negative ion density in the axial-vertical central plane (see Fig. 3) for the (top) pair injection and (bottom) flux injection. The injection zone is marked with orange dashed lines, the meniscus with a yellow line, and the plotting line along which the quantities in Fig. 9 are plotted is marked with a purple line.

Close modal

Figure 9 shows the potential and particle density along the axial direction, along the plotting line 2 mm from the bottom of the simulation domain for the pair injection simulation and one-sided flux injection simulation. In both cases, 200 A/m2 of H SP is emitted from the PG, causing the formation of a local minimum in the potential, a so-called virtual cathode (VC) of approximately 1 V and a corresponding peak of negative ion density.

FIG. 9.

Electrostatic potential and particle density along the purple plotting line shown in Fig. 8 for the flux injection case and the pair injection case with PG at floating potential. Parameters: n e = 6 × 10 16 m−3, T e = 2 eV, Γ H SP = 200 A/m−2, and U ext = 10 kV.

FIG. 9.

Electrostatic potential and particle density along the purple plotting line shown in Fig. 8 for the flux injection case and the pair injection case with PG at floating potential. Parameters: n e = 6 × 10 16 m−3, T e = 2 eV, Γ H SP = 200 A/m−2, and U ext = 10 kV.

Close modal

The transition from bulk plasma on the left-hand side to the simulation domain is smooth in the flux injection case. In contrast, with pair injection, a sheath toward the left-hand side of the simulation domain is formed. This is because this boundary is treated as any other wall, and thus a Debye sheath is formed. This sheath does not exist in the real ion source and is only a numerical remnant from the injection scheme that has been used. This sheath impacts the simulation in multiple ways. It accelerates positive ions from the injection zone toward the bulk plasma and reflects electrons and negative ions. The latter is valid also for H SP ions produced on the PG that may be reflected by the sheath and not exit the simulation domain, even though they are traveling away from the PG aperture with high energy, accelerated by the Debye sheath formed in front of the PG. The NI density in the axial direction differs between the two simulations. In the pair injection case, the NI density is significantly higher at −25 mm (located approximately in the middle of the injection zone in the pair injection case) than near the PG at the axial position −12 mm. For the flux injection case, the NI density is approximately equal in the two locations. However, as in the 1D simulation, a source sheath is formed near the injection plane, causing an increase in particle density in its vicinity.

While the PG in the pair injection simulation shares the same potential as the upstream wall, the H SP ions' perpendicular kinetic energy to the wall is less than their initial acceleration energy from the PG due to the PG's 50° angle. A majority of the H SP ions will, therefore, be reflected by the unphysical sheath formed on this side of the simulation domain. This is a problem for any pair injection model with an upstream Dirichlet boundary condition since the sheath toward the boundary reflects the approaching H SP ions. The flux injection scheme instead produces a plasma potential closer to what is taking place in the real ion source, with a flat transition from the simulation domain to the bulk plasma.

Figure 10 shows a comparison of the menisci for the two simulations. The penetration of the meniscus into the plasma is deeper in the pair injection case since the particle density is lower in the vicinity of the meniscus compared to the flux injection case. This shows that although the particle densities for both particle injection types are similar in the injection zone, particle fluxes are different, resulting in different meniscus shapes.

FIG. 10.

Meniscus profile for the pair injection and flux injection simulations with PG at the floating potential. Parameters: n e = 6 × 10 16 m−3, T e = 2 eV, Γ H SP = 200 A/m−2, and U ext = 10 kV.

FIG. 10.

Meniscus profile for the pair injection and flux injection simulations with PG at the floating potential. Parameters: n e = 6 × 10 16 m−3, T e = 2 eV, Γ H SP = 200 A/m−2, and U ext = 10 kV.

Close modal

The extracted currents are summarized in Table I. The total extracted negative ion current is higher in the flux injection case (146 A/m2) compared to that in the pair injection case (108 A/m2). The co-extracted electron currents are 108 and 93 A/m2 for the flux and pair injection simulations, respectively. The differences in the extracted negative ion current stem from the differences in the penetration of the meniscus and the distribution of negative ions coming from the bulk plasma. In the flux injection case, the plasma density (including the density of H ions) in the meniscus region is higher than for pair injection, leading to a higher extracted current of H VP . However, the higher plasma density leads to a meniscus shape that penetrates less on the plasma side of the PG aperture and thus fewer negative ions, produced on the PG are directly extracted.

TABLE I.

Extracted currents for the flux injection and pair injection simulations.

j e j ex
Sim. A / m 2 A / m 2 j e / j ex
Flux inj.  108  146  0.74 
Pair inj.  93  108  0.86 
j e j ex
Sim. A / m 2 A / m 2 j e / j ex
Flux inj.  108  146  0.74 
Pair inj.  93  108  0.86 

2. Electron transport

The transport of electrons in the axial direction is reduced due to the perpendicular magnetic fields (both BFF and BDF), leading to a reduction of electron density from the injection zone toward the aperture. This thermalization of electrons in the pair injection scheme has a great impact on the electron transport.

Figure 11 shows sample trajectories of electrons with the thermalization turned on and off; the particle tracking was done using the a-posteriori particle tracking developed for ONIX in Ref. 32 and using the electrostatic potential obtained from the flux injection simulation. The trajectories of electrons follow the typical banana orbits of the BDF in the axial-vertical plane (see Fig. 2). However, the trajectories are interrupted by the thermalization, which leads to an unphysical gradient in the electron density for the pair injection case, corresponding to the bump in electron density observed in the injection zone in Fig. 9. The thermalization can also have an impact on the co-extraction of electrons by introducing an artificial electron transport between field lines.

FIG. 11.

Electron trajectories with (top) and without (bottom) thermalization. The thermalization takes place in the injection zone, marked with a gray box; the meniscus is shown with a red line, and one field line of the BDF is shown with a black line.

FIG. 11.

Electron trajectories with (top) and without (bottom) thermalization. The thermalization takes place in the injection zone, marked with a gray box; the meniscus is shown with a red line, and one field line of the BDF is shown with a black line.

Close modal

To show the effect that each injection method plays on the formation of an ion–ion plasma in front of the meniscus, Fig. 12 shows the fraction of electrons to ions for the flux and pair injection simulations. Due to thermalization in the pair injection case, the electron trajectories are interrupted and the transition to an ion–ion plasma is not smooth. The thermalization region is not placed far enough to give electrons sufficient space for self-consistent redistribution.

FIG. 12.

H to e ratios for pair injection (top) and flux injection (bottom), PG at the floating potential. Parameters: n e = 6 × 10 16 m−3, T e = 2 eV, Γ H SP = 200 A/m−2, and U ext = 10 kV.

FIG. 12.

H to e ratios for pair injection (top) and flux injection (bottom), PG at the floating potential. Parameters: n e = 6 × 10 16 m−3, T e = 2 eV, Γ H SP = 200 A/m−2, and U ext = 10 kV.

Close modal

To minimize the effect of the thermalization region on the electron trajectories, the simulation domain for the pair injection must be enlarged to make room for the self-consistent formation of the plasma, which is highly computationally intensive.

The electron drift, represented by the accumulated horizontal and vertical distances that the electrons travel, can be tracked by saving the global coordinate for each particle throughout its trajectory, including laps across the periodic boundaries of the simulation. To investigate the effect that thermalization, an integral part of the pair injection method, plays in electron transport, the electrons were (a) not thermalized, and thermalized with a frequency (b) f = 10 9 s−1, as is used in Ref. 26, and (c) f = 2 × 10 9 s−1, as is used in Ref. 15.

Figure 13 shows the global coordinates of the electrons when passing through the meniscus in the horizontal–vertical plane. The trajectories of the electrons are again tracked using a-posteriori particle tracking. The electrons are injected with a temperature of 2 eV in the injection zone between −10 and 10 mm, marked in red in Fig. 13. The horizontal and vertical boundaries are periodic, but since the accumulated horizontal and vertical distances traveled by each electron is saved, one can calculate how many aperture rows or columns each electron moves before extraction. The electrons are allowed to move indefinitely in the horizontal and vertical directions, but only the five nearest aperture columns and three nearest aperture rows to the injection zone are shown. For the case without thermalization, where f = 0 s−1, 29.2% of the electrons are extracted through the same aperture row as where they are injected, meaning that a majority of electrons are drifting to neighboring apertures by crossing at least one periodic boundary. For the two cases with thermalization, a higher percentage, 45.7% and 45.8%, respectively, are extracted from the same aperture as the injection zone. In all three cases, the distance traveled is higher in the horizontal direction than in the vertical direction due to the travel along the BFF field lines, which span over multiple aperture rows. In contrast, the BDF spans over only one aperture, causing the electrons to follow the typical banana orbits.

FIG. 13.

Histograms of the global coordinate for the electrons once crossing the meniscus. The electrons were injected in the injection zone of the central aperture, marked with a red square. The particle tracking was done with (a) no thermalization, (b) f = 10 9 s−1, and (c) f = 2 × 10 9 s−1. The apertures are shown with red circles.

FIG. 13.

Histograms of the global coordinate for the electrons once crossing the meniscus. The electrons were injected in the injection zone of the central aperture, marked with a red square. The particle tracking was done with (a) no thermalization, (b) f = 10 9 s−1, and (c) f = 2 × 10 9 s−1. The apertures are shown with red circles.

Close modal

Figure 14 shows the horizontal profiles of the global positions where electrons pass through the meniscus. For the non-thermalized case, there is a clear left–right asymmetry, where more electrons are traveling in the negative horizontal direction, due to the impact of the downward-facing BDF. The profiles for the two cases with electron thermalization are almost equal, indicating that a thermalization rate of 109 s−1 is sufficient to hinder the horizontal transport of the electrons along the filter field lines.

FIG. 14.

Horizontal profile of the global positions where electrons pass through the meniscus with ( f = 10 9 s−1) and without (f = 0 s−1) thermalization.

FIG. 14.

Horizontal profile of the global positions where electrons pass through the meniscus with ( f = 10 9 s−1) and without (f = 0 s−1) thermalization.

Close modal

The flux injection scheme permits a seamless transition of the plasma potential from the bulk plasma on the upstream side of the simulation domain to the simulated plasma. Biasing of the PG is achieved by modifying the potential difference between the bulk plasma and the PG, treating it as an adjustable input parameter for the simulation. The PG bias ranged from −5 to +1 V relative to the plasma potential. The extracted current densities for each simulation are presented in Table II.

TABLE II.

Extracted current for varying PG bias. Parameters: n e = 6 × 10 16 m−3, T e = 2 eV, Γ H SP = 200 A/m−2, and U ext = 10 kV.

ϕ PG j e j ex
V A / m 2 A / m 2 j e / j ex
−5  108  146  0.74 
−3  79  147  0.54 
−1  68  148  0.46 
+1  23  151  0.15 
ϕ PG j e j ex
V A / m 2 A / m 2 j e / j ex
−5  108  146  0.74 
−3  79  147  0.54 
−1  68  148  0.46 
+1  23  151  0.15 

In the specific scenario where the PG bias is −5 V, the sheath close to the PG is strongly electron repelling. Consequently, electrons from the plasma are unlikely to reach the PG and are thus more likely to be extracted. The co-extracted electron current stands approximately at 110 A∕m2 under an extraction voltage of 10 kV. Increased PG potential leads to more e that can overcome the potential barrier, decreasing the electron density in the vicinity of the aperture and, in turn, a decrease in the co-extracted electron current. This is observed in the simulations; when increasing the PG bias toward the plasma's potential, a strong reduction in the co-extracted electron current is observed (see Table II). This reduction exceeds a factor of 4, driving the values significantly below the ITER target value of 1, which follows the trend shown in experimental investigations.3 

Figure 15 shows the potential profiles for varying PG bias along the axial direction, centered in the horizontal direction and near the bottom of the simulation domain in the vertical direction. The plasma potential for simulations with −5 and −3 V are unchanged. However, when increasing the potential to −1 and +1 V, the potential of the plasma also increases. This potential increase is not proportional to the increase in the PG potential, which explains why a further reduction of co-extraction of electrons is observed. This is possible due to the one-sided flux injection scheme, which allows more electrons to be collected on the PG and thus the simulation of PG bias above floating potential. While the Debye sheath toward the PG is affected by the change in PG bias, the meniscus shape is not impacted (not shown). The depth of the VC is equal in all cases (around 1 V).

FIG. 15.

Potential profiles along the axial direction for varying PG bias, taken 2 mm from the bottom of the simulation domain. Parameters: n e = 6 × 10 16 m−3, T e = 2 eV, Γ H SP = 200 A/m−2, and U ext = 10 kV.

FIG. 15.

Potential profiles along the axial direction for varying PG bias, taken 2 mm from the bottom of the simulation domain. Parameters: n e = 6 × 10 16 m−3, T e = 2 eV, Γ H SP = 200 A/m−2, and U ext = 10 kV.

Close modal

Conducting 3D-PIC MCC simulations of the extraction region in fusion negative ion sources poses considerable challenges. Computational limitations necessitate a significantly smaller simulation domain compared to the full ion source. Therefore, an open system is simulated, which strongly depends on the boundary conditions and the fluxes of particles crossing them. Modifying the particle injection model provides several improvements to the study of plasma physics. The improvements include the ability to bias the PG with respect to the plasma potential. This is achieved here for the first time using a 3D-PIC code, making simulations with a realistic potential without an unphysical sheath between the simulation domain and the boundary to the bulk plasma possible. Biasing of the PG was performed in 3D-PIC simulations using flux injection. Both negative biases (from −5 V) and positive biases (up to +1 V) with respect to the plasma were performed. Similar tendencies as observed in experiments were obtained, namely, a substantial reduction in the co-extracted electron current. Additionally, with the pair injection model, which is commonplace in 3D-PIC MCC simulations for negative ion extraction, a thermalization and particle (re-)injection region are needed. Electron transport along field lines is therefore not possible since the electron trajectories are interrupted by the unphysical thermalization. To give the electrons sufficient space to travel along the field lines, the simulation domain would have to be further extended into the bulk plasma when using the pair injection model, which can be unfeasible due to the massive computational cost of 3D-PIC MCC simulations. The flux injection method allows for studying electron losses along field lines. Additionally, it allows for a more realistic simulation of particle fluxes, with positive ions flowing toward the PG and electrons freely following the magnetic field lines without thermalization. For a given input density, the meniscus becomes flatter with the flux injection scheme, showing that the particle injection scheme can have an impact on the beam optics.

Future studies include the simulation of positive bias of the PG. This can be done by increasing the influx of electrons from the bulk plasma.

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 and 2019–2020 under Grant Agreement No. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No. 101052200 – EUROfusion). Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.

We express our gratitude to the Cineca consortium for providing access to the Marconi fusion supercomputer, which was essential for the simulations presented in this paper.

The authors have no conflicts to disclose.

M. Lindqvist: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). D. Wünderlich: Conceptualization (supporting); Data curation (supporting); Formal analysis (supporting); Funding acquisition (supporting); Investigation (supporting); Methodology (supporting); Project administration (supporting); Resources (supporting); Supervision (supporting); Validation (supporting); Writing – review & editing (supporting). S. Mochalskyy: Conceptualization (supporting); Data curation (supporting); Formal analysis (supporting); Investigation (supporting); Methodology (supporting); Supervision (supporting); Writing – review & editing (supporting). N. den Harder: Investigation (supporting); Methodology (supporting); Software (supporting); Visualization (supporting); Writing – review & editing (supporting). A. Revel: Data curation (supporting); Formal analysis (supporting); Methodology (supporting); Software (supporting); Validation (supporting); Writing – review & editing (supporting). T. Minea: Conceptualization (supporting); Formal analysis (supporting); Funding acquisition (lead); Investigation (supporting); Methodology (supporting); Project administration (lead); Resources (supporting); Supervision (lead); Validation (supporting); Writing – review & editing (supporting). U. Fantz: Conceptualization (supporting); Data curation (supporting); Formal analysis (supporting); Funding acquisition (supporting); Investigation (supporting); Project administration (lead); Resources (lead); Supervision (lead); Validation (supporting); Writing – review & editing (supporting).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
International Atomic Energy Agency
,
ITER Technical Basis: ITER EDA Documentation Series No. 24
(
International Atomic Energy Agency
,
Vienna
,
2002
).
2.
U.
Fantz
,
C.
Hopf
,
D.
Wünderlich
,
R.
Friedl
,
M.
Fröschle
,
B.
Heinemann
,
W.
Kraus
,
U.
Kurutz
,
R.
Riedl
,
R.
Nocentini
, and
L.
Schiesko
, “
Towards powerful negative ion beams at the test facility ELISE for the ITER and DEMO NBI systems
,”
Nucl. Fusion
57
,
116007
(
2017
).
3.
C.
Wimmer
,
L.
Schiesko
, and
U.
Fantz
, “
Investigation of the boundary layer during the transition from volume to surface dominated H production at the BATMAN test facility
,”
Rev. Sci. Instrum.
87
,
02B310
(
2016
).
4.
D.
Wünderlich
,
W.
Kraus
,
M.
Fröschle
,
R.
Riedl
,
U.
Fantz
, and
B.
Heinemann
, “
Long pulse, high power operation of the ELISE test facility
,”
AIP Conf. Proc.
1869
,
030003
(
2017
).
5.
D.
Yordanov
,
D.
Wünderlich
,
C.
Wimmer
, and
U.
Fantz
, “
On the effect of biased surfaces in the vicinity of the large extraction area of the ELISE test facility
,”
J. Phys.: Conf. Ser.
2244
,
012050
(
2022
).
6.
P. C.
Stangeby
, “
Plasma sheath transmission factors for Tokamak edge plasmas
,”
Phys. Fluids
27
,
682
(
1984
).
7.
L. A.
Schwager
and
C. K.
Birdsall
, “
Collector and source sheaths of a finite ion temperature plasma
,”
Phys. Fluids B
2
,
1057
1068
(
1990
).
8.
B.
Heinemann
,
U.
Fantz
,
W.
Kraus
,
L.
Schiesko
,
C.
Wimmer
,
D.
Wünderlich
,
F.
Bonomo
,
M.
Fröschle
,
R.
Nocentini
, and
R.
Riedl
, “
Towards large and powerful radio frequency driven negative ion sources for fusion
,”
New J. Phys.
19
,
015001
(
2017
).
9.
U.
Fantz
,
S.
Briefi
,
A.
Heiler
,
C.
Wimmer
, and
D.
Wünderlich
, “
Negative hydrogen ion sources for fusion: From plasma generation to beam properties
,”
Front. Phys.
9
,
473
(
2021
).
10.
C. K.
Birdsall
and
A. B.
Langdon
,
Plasma Physics via Computer Simulation
(
McGraw-Hill
,
1991
), pp.
58
63
.
11.
D.
Matsushita
,
S.
Kuppel
,
A.
Hatayama
,
A.
Fukano
, and
M.
Bacal
, “
Modeling of the plasma electrode bias in the negative ion sources with 1D pic method
,”
AIP Conf. Proc.
1097
,
38
46
(
2009
).
12.
F.
Taccogna
,
P.
Minelli
,
S.
Longo
,
M.
Capitelli
, and
R.
Schneider
, “
Modeling of a negative ion source. III. Two-dimensional structure of the extraction region
,”
Phys. Plasmas
17
,
063502
(
2010
).
13.
F.
Taccogna
,
P.
Minelli
,
M.
Capitelli
,
S.
Longo
, and
R.
Schneider
,
Plasma Grid Shape and Size Effects on the Extraction of Negative Ions
(
AIP
,
2013
).
14.
F.
Taccogna
and
P.
Minelli
, “
PIC modeling of negative ion sources for fusion
,”
New J. Phys.
19
,
015012
(
2017
).
15.
G.
Fubiani
,
L.
Garrigues
,
G.
Hagelaar
,
N.
Kohen
, and
J. P.
Boeuf
, “
Modeling of plasma transport and negative ion extraction in a magnetized radio-frequency plasma source
,”
New J. Phys.
19
,
015002
(
2017
).
16.
A.
Revel
,
S.
Mochalskyy
,
I. M.
Montellano
,
D.
Wünderlich
,
U.
Fantz
, and
T.
Minea
, “
Massive parallel 3D PIC simulation of negative ion extraction
,”
J. Appl. Phys.
122
,
103302
(
2017
).
17.
S.
Nishioka
,
I.
Goto
,
K.
Miyamoto
,
A.
Hatayama
, and
A.
Fukano
, “
Study of ion-ion plasma formation in negative ion sources by a three-dimensional in real space and three-dimensional in velocity space particle in cell model
,”
J. Appl. Phys.
119
,
023302
(
2016
).
18.
D.
Wünderlich
,
D.
Yordanov
,
R.
Riedl
,
A.
Mimo
,
A.
Hurlbatt
,
B.
Heinemann
, and
U.
Fantz
, “
Paving the road towards ITER relevant long deuterium pulses at ELISE by investigating improved operational scenarios
,” in
Talk Presented at 8th International Symposium on Negative Ions, Beams and Sources
(
NIBS
,
2022
).
19.
M.
Bacal
and
M.
Wada
, “
Negative Hydrogen ion production mechanisms
,”
Appl. Phys. Rev.
2
,
021305
(
2015
).
20.
M.
Bacal
, “
Volume production of Hydrogen negative ions
,”
Nucl. Instrum. Methods Phys. Res. B
37–38
,
28
32
(
1989
).
21.
M.
Bacal
and
M.
Wada
,
Negative Ion Production by Plasma-Surface Interaction in Caesiated Negative Ion Sources
(
AIP
,
2013
).
22.
A.
Mimo
,
H.
Nakano
,
C.
Wimmer
,
D.
Wünderlich
,
U.
Fantz
, and
K.
Tsumori
, “
Cavity ring-down spectroscopy system for the evaluation of negative hydrogen ion density at the ELISE test facility
,”
Rev. Sci. Instrum.
91
,
013510
(
2020
).
23.
S.
Mochalskyy
,
D.
Wünderlich
,
B.
Ruf
,
P.
Franzen
,
U.
Fantz
, and
T.
Minea
, “
3D numerical simulations of negative hydrogen ion extraction using realistic plasma parameters, geometry of the extraction aperture and full 3D magnetic field map
,”
Rev. Sci. Instrum.
85
,
02B301
(
2014
).
24.
S.
Mochalskyy
,
D.
Wünderlich
,
U.
Fantz
,
P.
Franzen
, and
T.
Minea
, “
Towards a realistic 3D simulation of the extraction region in ITER NBI relevant ion source
,”
Nucl. Fusion
55
,
033011
(
2015
).
25.
S.
Mochalskyy
,
U.
Fantz
,
D.
Wünderlich
, and
T.
Minea
, “
Comparison of ONIX simulation results with experimental data from the BATMAN testbed for the study of negative ion extraction
,”
Nucl. Fusion
56
,
106025
(
2016
).
26.
I. M.
Montellano
,
D.
Wünderlich
,
S.
Mochalskyy
, and
U.
Fantz
, “
3D-PIC modelling of a low temperature plasma sheath with wall emission of negative particles and its application to NBI sources
,”
J. Phys. D: Appl. Phys.
52
,
235202
(
2019
).
27.
D.
Wünderlich
,
R.
Gutser
, and
U.
Fantz
, “
PIC code for the plasma sheath in large caesiated RF sources for negative hydrogen ions
,”
Plasma Sources Sci. Technol.
18
,
045031
(
2009
).
28.
M.
Lindqvist
,
D.
Wünderlich
,
A.
Mimo
,
S.
Mochalskyy
,
A.
Revel
,
R.
Nocentini
,
T.
Minea
, and
U.
Fantz
, “
Sensitivity of the negative ion beam extraction to initial plasma parameters by 3D particle modelling
,”
Plasma Sources Sci. Technol.
31
,
125001
(
2022
).
29.
A.
Isidori
,
Control Theory for Automation: Fundamentals – Pringer Handbook of Automation
(
Springer Berlin Heidelberg
,
Berlin, Heidelberg
,
2009
).
30.
R.
Hemsworth
,
H.
Decamps
,
J.
Graceffa
,
B.
Schunke
,
M.
Tanaka
,
M.
Dremel
,
A.
Tanga
,
H.
de Esch
,
F.
Geli
,
J.
Milnes
,
T.
Inoue
,
D.
Marcuzzi
,
P.
Sonato
, and
P.
Zaccaria
, “
Status of the ITER heating neutral beam system
,”
Nucl. Fusion
49
,
045006
(
2009
).
31.
R.
Nocentini
, private communication (
2020
).
32.
M.
Lindqvist
,
N.
den Harder
,
A.
Revel
,
S.
Mochalskyy
,
A.
Mimo
,
R.
Nocentini
,
T.
Minea
, and
U.
Fantz
, “
From meniscus formation to accelerated H beam: Coupling of 3D-PIC and ion-optics simulations
,”
Nucl. Fusion
62
,
126068
(
2022
).