Relationship between Surface-Produced H⁻ Density and the Potential Structure of an Ion-Ion Plasma near the Extraction Aperture in Negative Hydrogen Ion Sources

Abstract

The so-called “ion-ion” plasma, which consists primarily of negative and positive ions, has been observed near the extraction aperture in negative hydrogen ion (H⁻) sources with a substantial amount of surface-produced negative ions, according to several experiments. In this study, to elucidate the characteristics of ion-ion plasmas with substantial surface H⁻ production, the relationship between surface-produced negative ion density and electric potential is investigated using three-dimensional Particle-In-Cell simulations with the KEIO-BFX code. The results reveal that the plasma in a negative ion source can be classified into four distinct regions based on the relationship between H⁻ density and electric potential: Region (I), where H⁻ follows the Boltzmann relation; Region (II), where H⁺ follows the Boltzmann relation; Region (III), a transition region between the plasma region and the beam acceleration region; and Region (IV), the beam acceleration region. In negative ion sources with surface-produced H⁻ ions, plotting the H⁻ density against the electric potential offers a useful method for identifying the plasma meniscus, which is crucial for optimizing extracted beam quality and ensuring proper operation of beam extraction systems.

1.  Introduction

Negative hydrogen ion (H⁻) sources are used in a range of high-energy technologies, including particle acceleratos [1, 2], plasma heating for nuclear fusion [3–6], and medical applications [7, 8]. In the development of H⁻ sources, high-power beams and good beam optics are required.

To increase H⁻ density and beam intensity, negative ion sources use surface production [9]. Although beam intensity has continued to increase, further improvement in beam convergence is required [10] to avoid potential damage from the halo component and to prevent reductions in operation time.

Recent studies have highlighted several effects of surface-produced H⁻ ions, such as the formation of a virtual cathode in the sheath near the wall where H⁻ ions are generated [11] and the presence of an electrostatic lens near the extraction aperture. These phenomena are closely related to beam optics [12]. Specifically, these studies demonstrate that surface-produced H⁻ ions influence the electrostatic potential structure and the beam-emitting surface, often referred to as the plasma meniscus.

Moreover, when a large number of H⁻ ions are present inside the ion source owing to surface production, the H⁻ ions themselves have a considerable impact on the plasma. For example, surface-produced H⁻ ions not only contribute to the divergent component of the extracted beam through the curvature of the meniscus [13, 14], but also form an ion-ion plasma in which H⁻ ions, rather than electrons, dominate the negative charge [15–19]. Such a unique plasma, composed predominantly of negative and positive ions (H⁻ and H⁺) with no mass difference, can form a potential structure distinct from that of conventional electron-ion plasmas, potentially affecting beam extraction.

Therefore, the aim of this study is to explore the relationship between the surface-produced H⁻ density and the electric potential in an ion-ion plasma, particularly one with a substantial amount of surface H⁻ production. For this purpose, the three-dimensional (3D) Particle-In-Cell (PIC) code KEIO-BFX [18, 20, 21] is used.

2.  Simulation Model

The KEIO-BFX code is a three-dimensional PIC code that self-consistently solves the equations of motion and the Poisson equation to simulate plasma dynamics and the electric field near the extraction hole. The basic equations of the code are given as follows:

  

$$\frac{d\mathit{x}}{dt}=\mathit{v},\frac{d\mathit{v}}{dt}=\frac{q}{m}(\mathit{E}+\mathit{v}\times \mathit{B}),$$(1)

  

$${\nabla }^2\phi =-\frac{\rho }{{\epsilon }_0},$$(2)

where $\mathit{x}$, $\mathit{v}$, $q$, and $m$ denote the position, velocity, charge, and mass of a particle, respectively, while $\mathit{E}$ and $\mathit{B}$ represent the electrostatic and magnetic fields. In this simulation, the magnetic field is kept constant and is determined by the filter magnet used to trap high-energy electrons. In addition, $\phi$ and $\rho$ denote the electric potential and charge density, respectively. In the PIC method, $\rho$ at the simulation grid points is evaluated using a linear weighting scheme based on particle positions obtained from the equations of motion.

The KEIO-BFX code is used to analyze the Linac4 radio frequency (RF) ion sources at CERN, as illustrated in Fig. 1. The simulation domain represents the region near the extraction hole, surrounded by the plasma grid (PG). The voltage at the right-hand boundary (z = 3 mm) is set to the value precalculated from the extraction voltage of 9.7 kV, applied to the extraction grid (EG) located outside the simulation region (z = 7 mm). This precalculated potential serves as the boundary condition for solving the Poisson equation.

Fig. 1.  Schematic of the Linac4 RF H⁻ source at CERN (left) and a cross-sectional view of the modeled extraction region near the extraction hole (right).

The calculation conditions are the same as those in the previous study [11] and are summarized in Table 1. Parameters such as density, temperature, and density ratio are selected based on results from simulations using the NINJA code [22]. Key considerations include ensuring that the mesh size is small enough to resolve spatial scales on the order of the Debye length and to avoid numerical plasma heating. Additionally, the timestep width is chosen to be sufficiently small to satisfy the Courant–Friedrichs–Lewy condition [23]. However, performing PIC simulations on a realistic system with such a fine mesh comes with considerable computational costs. To address this challenge, a reduced-size-scaling approach is adopted. In this approach, the simulated system length is scaled as $L_{\text{sim}}=s{L}_{\text{real}}$, where $L_{\text{sim}}$ and $L_{\text{real}}$ denote the lengths in the simulation and the real system, respectively. The scaling factor is set to $s$ = 3.5 × 10⁻². The mean free paths for various collision processes are also scaled as B_sim = $s$^−1/3 B and λ_sim = $s$λ. The exponent of this scaling factor is determined by equating the perveance between the real and simulated systems [24]. Similar scaling approaches have been applied in PIC simulations to study the Scrape-off Layer and divertor plasmas in fusion devices [25], as well as negative ion sources [24, 26–29].

Table 1. Plasma and PIC parameters.

Parameter	Value
Electron density	1.0 × 1018 m−3
Electron temperature	3.6 eV
Temperature for positive hydrogen ion and volume-produced H− temperature	1.6 eV
Surface-produced negative ion temperature	1.0 eV (Half Maxwellian)
Density ratio of initial plasma ($n_e:n_{{\mathrm{H}}^{+}}:n_{{\mathrm{H}}_2^{+}}:n_{{\mathrm{H}}_3^{+}}:n_{{\mathrm{H}}^{-}}$)	0.985 : 0.740 : 0.074 : 0.186 : 0.015
Debye length (λDe)	1.41 × 10−5 m
Electron thermal velocity	7.96 × 105 m/s
Plasma frequency (ωp)	5.64 × 1010 rad/s
Surface production rate SH−	1,000 A/m2
Extraction voltage	9.7 kV
Real size	31.75 × 48 × 48 mm3
Scaling factor	3.5 × 10−2
Number of superparticle	19,200,000
Mesh	127 × 192 × 192
Mesh size	0.625λDe
Timestep	0.4/ωp = 7.09 × 10−12 s
Simulation time	100,000 timesteps = 0.7 μs
Magnetic field strength	10–18 mT

The simulation is performed using 16 parallel processes. The total number of particles across all processes is 19,200,000, which is constrained by available computational resources. Approximately 20% of the mesh cells correspond to regions where no particles are present. As a result, the average number of particles per mesh cell is approximately 6. This particle number may be insufficient for statistical accuracy. To address this, statistical accuracy is improved by performing time-averaging over 2,500 timesteps. This approach leverages the concept of the equivalence between ensemble and time averages [30].

In this study, both volume H⁻ production and surface production are considered [11]. The surface production rate, denoted as S_H⁻, which represents the average current density of negative ions emitted from the PG, is set based on simulation results from the CERN RF ion source [31]. The destruction reaction of H⁻ ions is simulated using the null-collision method [32].

3.  Results

3.1  Spatial profile of ion density and electrostatic potential

Figure 2 shows the spatial profile of the electronegativity $\alpha$, which is the ratio of negative ion density to electron density ($\alpha =n_{{\mathrm{H}}^{-}}/n_{\mathrm{e}}$), near the extraction hole. In the extraction region, where the two PGs are closely spaced, α exceeds 10. This high electronegativity indicates the formation of an ion-ion plasma, primarily composed of positive ions and negative ions, owing to the large surface production of negative ions and the reflection of electrons by the sheath. Additionally, the ion-ion plasma extends approximately to z = −10 mm, mainly owing to the tapered structure of the PG that closes toward the extraction hole. In the beam region (z > 0 mm), a high α region is also observed outside the beam aperture because surface-produced H⁻ ions are directly extracted.

Fig. 2.  Spatial profile of electronegativity, $\alpha =n_{{\mathrm{H}}^{-}}/n_{\mathrm{e}}$ near the extraction hole.

The two-dimensional (2D) density profile of surface-produced H⁻ ions is shown in Fig. 3. For z < −10 mm, the H⁻ density increases towards the bulk plasma. This indicates that the surface-produced H⁻ ions, which are not directly extracted, are accelerated by the sheath electric field and transported to the upper source region. On the other hand, in the range −10 mm < z < −2 mm, the contour lines corresponding to 4, 5, and 6 × 10¹⁷ m⁻³ are widely spaced. This indicates that the spatial variation of the H⁻ density is small in this region. For z > −2 mm, the density decreases to 1 × 10¹⁷ m⁻³ owing to the acceleration from the extraction voltage. This spatial variation in density is also evident in the one-dimensional ion density distribution along the extraction axis (z-direction), as shown in Fig. 4. Moreover, the surface-produced H⁻ density is much larger than the volume-produced H⁻ density. Therefore, in this study, the focus is on the surface-produced H⁻ ions.

Fig. 3.  2D density profile of surface-produced H⁻ near the extraction hole. The red line labeled 10¹⁷ indicates the outline of the beam.

Fig. 4.  Ion density profiles along the extraction axis. The subscripts H⁺, H⁻_SP, and H⁻_VP denote the positive hydrogen ion, the surface produced H⁻ ion, and the volume produced H⁻ ion, respectively.

A 2D density profile of H⁺ ions is shown in Fig. 5. For z < −2 mm, the H⁺ ion density profile closely matches the H⁻ ion density profile shown in Fig. 3. In this region, the positive and negative charges are balanced between H⁺ and H⁻ ions, confirming the formation of the ion-ion plasma. However, around z = −1 mm, the H⁺ ion density decreases sharply, while the variation in H⁻ density remains relatively small. This suggests that the quasi-neutrality of the ion-ion plasma no longer holds for z > −1 mm. Therefore, this region is closely associated with the plasma meniscus.

Fig. 5.  2D profile of H⁺ density.

These ion density profiles are closely related to the electric potential profile near the extraction hole, including the sheath in front of the PG and the extraction voltage. The spatial profile of the electric potential φ is shown in Figs. 6 and 7. In this simulation, despite assuming a high surface production rate, a clear virtual cathode is not observed, as shown in Fig. 6. This could be due to the surface-produced H⁻ ions being drawn out along the PG by the extraction potential, which reduces their role in forming the virtual cathode. As shown in Figs. 6 and 7, the positive extraction voltage penetrates up to approximately z = −2 mm, where φ changes steeply, resulting in a concave-shaped φ profile. The surface-produced H⁻ ions are accelerated by this concave-shaped φ profile, leading to the formation of a converging beam, as shown in Figs. 3 and 4. In addition to the penetration of the extraction voltage, the sheath effect causes the potential to decrease from the plasma center toward the PG, as shown in Figs. 3 and 4. As a result, the potential profile exhibits a local minimum at z ~ −8 mm along the extraction axis (y = 0), forming a potential well. This potential well is closely related to the density profiles shown in Figs. 3 and 5. Accordingly, in the next section, we examine the dependence of the negative ion density on the potential along the extraction axis (y = 0).

Fig. 6.  2D profile of the electric potential φ.

Fig. 7.  1D profile of the potential φ along the extraction axis (y = 0).

3.2  Dependence of surface-produced negative ion density on electrostatic potential

Figure 8 shows the plot of the surface-produced H⁻ ion density ($n_{{\mathrm{H}}^{-}}$) as a function of the potential along the extraction axis. The log₁₀($n_{{\mathrm{H}}^{-}}$)-φ plot reveals four distinct regions as follows: Region (I)—a region with a positive linear trend, Region (II)—a region with a negative linear trend, Region (III)—a region where the density remains almost constant despite an increase in φ, and Region (IV)—a region where the density begins to decrease monotonically. The boundaries between Regions (I) and (II), (II) and (III), and (III) and (IV) are located at z = −8, −4.5, and −2.25 mm, respectively.

Fig. 8.  Surface-produced H⁻ ion density as a function of electric potential along the extraction axis (y = 0). The dependence of $n_{{\mathrm{H}}^{-}}$ on φ varies across Region (I) (z < −8 mm), (II) (−8 mm < z < −4.5 mm), (III) (−4.5 mm < z < −2.25 mm), and (IV) ( −2.25 mm < z).

The physical characteristics of each region can be understood as follows.

Regions (I) and (II) are closely related to the potential structure with a local minimum. In Region (I), the surface-produced H⁻ ions flow upstream owing to the positive potential and are trapped between the source plasma region and the potential minimum, as shown in Fig. 9. As a result, H⁻ ions settle into a steady-state distribution in Region (I). In this case, where the H⁻ ions are trapped by the potential for a sufficiently long time, the H⁻ density follows the Boltzmann relation. The Boltzmann relation is generally written using the density constant $N_1$ as follows:

Fig. 9.  Schematic diagram of the potential profile along the extraction axis and the behavior of ions in Regions (I) and (II).

  

$$n_{{\mathrm{H}}^{-}}=N_1\exp \left(\frac{q\phi }{k{T}_{{\mathrm{H}}^{-}}}\right),$$(3)

where $q$, $k$, and $T_{{\mathrm{H}}^{-}}$ are the elementary charge, Boltzmann constant, and negative ion temperature, respectively. Taking the logarithm of this equation and converting it to base 10 for clarity in the numerical results, we obtain

  

$${\log }_{10}{n}_{{\mathrm{H}}^{-}}\propto {\log }_{10}e•\frac{q}{k{T}_{{\mathrm{H}}^{-}}}\phi .$$(4)

From this equation, the H⁻ ion density log₁₀ ($n_{{\mathrm{H}}^{-}}$) exhibits a positive linear dependence in a log₁₀($n_{{\mathrm{H}}^{-}}$)-φ plot in Region (I). Using Eq. (4) and the slope in Region (I), the H⁻ ion temperature is evaluated to be 1.0 eV. This temperature closely corresponds to the initial temperature of the surface-produced H⁻ ions, which supports the idea that the H⁻ density follows the Boltzmann relation.

In contrast to Region (I), the H⁺ ions in Region (II) are trapped around the potential well, as shown in Fig. 9. In this case, the trapped ions are H⁺ ions, and therefore their density follows the Boltzmann relation. Using the density constant $N_2$ and the positive ion temperature $T_{{\mathrm{H}}^{+}}$, the Boltzmann relation for ions can be written as follows,

  

$$n_{{\mathrm{H}}^{+}}=N_2\exp \left(\frac{-q\phi }{k{T}_{{\mathrm{H}}^{+}}}\right).$$(5)

Taking the base 10 logarithm, we obtain

  

$${\log }_{10}{n}_{{\mathrm{H}}^{+}}\propto -{\log }_{10}e•\frac{q}{k{T}_{{\mathrm{H}}^{+}}}\phi .$$(6)

Therefore, the H⁺ density exhibits a negative linear trend with φ. As shown in Fig. 3, an ion-ion plasma is formed in this region, meaning that quasi-neutrality is maintained by H⁺ and H⁻ ions as follows:

  

$$n_{{\mathrm{H}}^{-}}\approx n_{{\mathrm{H}}^{+}}.$$(7)

Therefore, the H⁻ density follows the H⁺ density profile, meaning that the H⁻ density also exhibits a negative linear trend on the log₁₀($n_{{\mathrm{H}}^{-}}$)-φ plot. Based on Eq. (6) and the slope in Fig. 10, the positive ion temperature is determined to be 1.2 eV, which is 0.4 eV lower than the initial temperature. This decrease is due to the high-energy tail in the velocity distribution of H⁺ ions escaping from the potential well.

Fig. 10.  H⁺ density as a function of electric potential along the extraction axis (y = 0).

In Region (III), a positive electric field is established owing to the penetration of the extraction voltage. At the same time, an ion-ion plasma is present, as shown in Fig. 2. In other words, H⁻ ions begin to be accelerated by the extraction electric field, while quasi-neutrality is still maintained by the continuous influx of surface-produced H⁻ ions from the PG. The maintenance of quasi-neutrality is confirmed by the charge density distribution shown in Fig. 11 for −5 < z < −2 mm. In positive ion sources, the meniscus is defined as the location where extracted particles (positive ions) begin to accelerate. However, since quasi-neutrality is still maintained in negative ion sources, this definition may not be valid. Therefore, Region (III) is considered the transition region between the plasma region and the beam acceleration region.

Fig. 11.  1D profile of the charge density along the extraction axis.

In Region (IV), H⁻ ions are extracted by the extraction voltage, resulting in a decrease in H⁻ density. Conversely, H⁺ ions are reflected back toward the plasma by the extraction voltage. Therefore, the quasi-neutral ion-ion plasma cannot be maintained for z > −2 mm, as shown in Fig. 11.

As shown in Figs. 7 and 8, the meniscus position is conventionally defined by the condition dφ/dx = 0, typically estimated from the continuity between the ion saturation current and the space-charge-limited current (Child-Langmuir law). This position corresponds to the boundary between Regions (I) and (II). However, as discussed above, quasi-neutrality is still maintained in Region (III) owing to the continuous supply of surface-produced negative ions perpendicular to the extraction axis. Therefore, in negative ion sources that include surface-produced H⁻ ions, it becomes challenging to estimate the meniscus position using the conventional one-dimensional approach, which is commonly applied in positive ion sources. Based on the discussion in this paper, the boundary between Regions (III) and (IV), where quasi-neutrality breaks down, can be interpreted as the meniscus when the influence of surface-produced negative ions is considered.

4.  Conclusion

In this study, the relationship between surface-produced H⁻ ion density and the electric potential in an ion-ion plasma with substantial surface H⁻ production is analyzed using three-dimensional PIC simulations with KEIO-BFX.

In the simulation, an ion-ion plasma was formed owing to the large influx of surface-produced negative ions. Additionally, a potential well (a local minimum) was created near the extraction hole because of the positive extraction voltage (z = 3 mm) and the potential drop caused by sheath formation. The relationship between the surface-produced H⁻ ion density ($n_{{\mathrm{H}}^{-}}$) and the electric potential (φ) in the ion-ion plasma was then investigated. The log₁₀($n_{{\mathrm{H}}^{-}}$)-φ plot along the extraction axis revealed four characteristic regions: Region (I), exhibiting a positive linear trend; Region (II), exhibiting a negative linear trend; Region (III), where log₁₀($n_{{\mathrm{H}}^{-}}$) remains nearly constant despite increasing φ; and Region (IV), where log₁₀($n_{{\mathrm{H}}^{-}}$) decreases as φ increases. The characteristics of each region can be understood as follows. In Region (I), owing to the potential well structure, H⁻ ions are trapped, while H⁺ ions are trapped around the potential well. As a result, in Region (I), the H⁻ density follows the Boltzmann relation and shows a positive linear trend on the log₁₀($n_{{\mathrm{H}}^{-}}$)-φ plot. In contrast, in Region (II), the H⁺ density follows the Boltzmann relation, and the H⁻ density distribution mirrors that of the H⁺ density. As a result, a negative linear trend is observed on the log₁₀($n_{{\mathrm{H}}^{-}}$)-φ plot. In Region (III), the H⁻ ions begin to be accelerated by the extraction electric field, while quasi-neutrality is still maintained. This region can therefore be considered a transition region from the quasi-neutral plasma region to the beam region. In Region (IV), the positive extraction voltage prevents the plasma from maintaining quasi-neutrality, and the negative ions are fully accelerated as a beam. Based on these characteristics, the plasma meniscus can be interpreted as the boundary between Region (III) and Region (IV).

As described above, plotting the surface-produced H⁻ ion density (log₁₀($n_{{\mathrm{H}}^{-}}$)) as a function of φ provides valuable insight into the characteristics of ion-ion plasmas. Moreover, analyzing this relationship may offer a method to estimate the meniscus position using an approach different from the one typically used for positive ion sources. These findings are expected to serve as a foundation for future research on meniscus position estimation in negative ion sources, which is crucial for controlling the quality of extracted beams and ensuring the proper operation of beam extraction systems.

 Acknowledgements

This work was supported by JST HAKASE 229761.

References

• [1]   J.  Lettry et al., Rev. Sci. Instrum. 85, 02B122 (2014).
• [2]   J.  Lettry et al., Rev. Sci. Instrum. 87, 02B139 (2016).
• [3]   R.  Hemsworth et al., Nucl. Fusion 49, 045006 (2009).
• [4]   M.  Kashiwagi et al., Rev. Sci. Instrum. 83, 02B119 (2012).
• [5]   K.  Tsumori et al., Rev. Sci. Instrum. 87, 02B936 (2016).
• [6]   U.  Fantz et al., Rev. Sci. Instrum. 87, 02B307 (2016).
• [7]   H.  Etoh et al., AIP Conf. Proc. 1869, 030050 (2017).
• [8]   M.  Onai et al., AIP Conf. Proc. 1869, 030043 (2017).
• [9]   Y.I.  Belchenko et al., Nucl. Fusion 14, 113 (1974).
• [10]   N. den  Harder et al., Nucl. Fusion 64, 076046 (2024).
• [11]   K.  Hayashi et al., Plasma Fusion Res. 18, 1401008 (2023).
• [12]   R.  McAdams et al., Plasma Sources Sci. Technol. 20, 035023 (2011).
• [13]   M.  Lindqvist et al., Phys. Plasmas 31, 033903 (2024).
• [14]   K.  Miyamoto et al., Plasma Sources Sci. Technol. 31, 105012 (2022).
• [15]   K.  Tsumori et al., Rev. Sci. Instrum.   83, 02B116 (2012).
• [16]   S.  Christ-Kouch et al., Plasma Sources Sci. Technol. 18, 025003 (2009).
• [17]   P.  McNeely et al., Plasma Sources Sci. Technol. 18, 014011 (2009).
• [18]   S.  Nishioka et al., J. Appl. Phys. 119, 023302 (2016).
• [19]   Z.  Zhang et al., Plasma Sources Sci. Technol. 27, 124003 (2018).
• [20]   S.  Nishioka et al., J. Appl. Phys. 123, 063302 (2018).
• [21]   M.  Lindqvist et al., J. Appl. Phys. 126, 123303 (2019).
• [22]   S.  Mattei et al., J. Comput. Phys. 350, 891 (2017).
• [23]   C.K.  Birdsall and  A.B.  Langdon, Plasma Physics via Computer Simulation, (CRC Press, Boca Raton, 2004).
• [24]   K.  Miyamoto et al., Appl. Phys. Lett. 102, 023512 (2013).
• [25]   T.  Takizuka et al., Nucl. Fusion 49, 075038 (2009).
• [26]   K.  Miyamoto et al., Appl. Phys. Lett. 100, 233507 (2012).
• [27]   A.  Hatayama, Rev. Sci. Instrum. 79, 02B901 (2008).
• [28]   S.  Kuppel et al., J. Appl. Phys. 109, 013305 (2011).
• [29]   A.  Hatayama et al., New J. Phys. 20, 065001 (2018).
• [30]   F.  Reif, Fundamentals of Statistical and Thermal Physics, (McGraw Hill Inc., New York, 1965).
• [31]   S.  Mochalskyy et al., New J. Phys. 18, 085011 (2016).
• [32]   R.K.  Janev et al., Elementary Processes in Hydrogen-Helium Plasmas: Cross Sections and Reaction Rate Coefficients, (Springer, Berlin, 1987).
• [33]   P.C.  Stangeby, The Plasma Boundary of Magnetic Fusion Devices, (Institute of Physics, Bristol and Philadelphia, 2000).