Three-dimensional simulation of electron extraction process and optimization of magnetic field based on multi-aperture structure of negative hydrogen ion source

In most of the simulations of the extraction region of negative hydrogen ion sources, the single-aperture simulation is often adopted by researchers to study the plasma phenomenon due to its small simulation domain and short calculation time. However, due to the complex three-dimensional magnetic field structure in the extraction region of the negative hydrogen ion source, the single aperture often does not meet the periodicity. In this paper, the complex three-dimensional magnetic field topology is established. The magnetic field includes the magnetic filter field and the magnetic deflection field. The influence of the plasma sheath is taken into account. The electron extraction process in the multi-aperture structure of the extraction region of a negative hydrogen ion source is numerically calculated using the PIC method. Besides, the magnetic field structure is optimized. Ultimately, the electron beam uniformity near the plasma grid is improved effectively, which has certain guiding significance for engineering application.


Introduction
Ion source, as the starting device of neutral beam injection (NBI) system, has become the research target of many scholars in recent years. Among ion sources, especially radiofrequency (RF) negative hydrogen ion source, it was listed as the reference scheme for plasma generation in the International Thermonuclear Experimental Reactor (ITER) Project in 2007 due to its advantages of simple structure, almost maintenance-free and low cost [1]. In region, namely, the plasma grid (PG), the extraction grid (EG) and the grounded grid (GG). Since electrons are also electronegative and will be extracted together with the beam, EG is often embedded with deflection magnets, and the resulting magnetic deflection field will dump the co-extracted electrons onto EG, reducing the co-extracted electron current and making the ratio of the co-extracted electron current density to negative ion current density less than 1, which is also one of the requirements of ITER [9].
The quality of the beam is determined by the plasma generated in the driver region and the design of the extraction region. Experimental attempts are often made to optimize the geometry of apertures in PG to improve the quality of the beam. For example, in the experiments of two different extraction systems (aperture diameters of 8 and 14 mm respectively) of the BATMAN device, preliminary results show that there is no major correlation between the extraction current density and aperture diameter [1]. However, the relevant experiments of the 14 mm aperture diameter extraction system are still hampered by the weak power load processing ability, and no clear conclusions can be drawn [10]. To determine the influence of various factors on beam quality by experiment, the technical difficulties need to be overcome. Therefore, it is necessary to make reasonable usage of simulation to understand the factors of beam quality.
Many relevant simulations of the extraction region mainly focus on the extraction of single aperture [11,12], because the geometric structure of apertures is periodic, and secondly, due to the limitation of Debye length (generally on the order of 20 μm), the grid step size should be selected to be slightly smaller than Debye length to avoid the numerical "grid heating" phenomenon caused by excessive grid step. In this way, the calculation of single aperture can save a lot of simulation domain and calculation time. However, for the extraction region of negative hydrogen ion sources, the complex three-dimensional magnetic field structure makes the single aperture no long meet the periodicity, and the positions of the apertures are not always symmetrical. For example, the large area grid (LAG) systems of BATMAN and MANITU are composed of the identical grid aperture structures, and apertures are staggered into 10/11 and 14/15 holes [10,13]. The topological structure of the magnetic field, especially the magnetic deflection field in EG, plays an important role in the numerical results. Small changes in magnetic field intensity or topological structure can greatly change the extracted electron current, and the correlation of simulation results can only be proposed based on a correctly established magnetic field structure [14].
The three-dimensional Particle-In-Cell (PIC) simulation method is adopted based on the self-developed CHIPIC code [15][16][17][18] in this paper. A complex three-dimensional magnetic field topology is established, which includes a magnetic filter field and a magnetic deflection field. The electron extraction process in multi-aperture structure of negative hydrogen ion source is numerically calculated based on the recent research [15], and the magnetic field structure is optimized to improve the electron beam uniformity near PG plate.

Simulation setup
As shown in Figure 1 16,20] are adopted on CHIPIC code. The particle motion, force, charge, and current density updates adopt the PIC method, while the particle collision part adopts the MCC method. Since the emitted electron in the following is determined based on the emitted electron The process cycle diagram of PIC/MCC and SMPM methods on CHIPIC code.

Frontiers in Physics
frontiersin.org distribution after the collision, this paper does not consider the influence of collision factor, and focuses on tracking the distribution and motion trajectory of electrons after the spatiotemporal iteration of simulation. Figure 2A shows the schematic view of negative hydrogen ion source. The aperture structure in the extraction region of negative hydrogen ion source is established by referring to the LAG extraction system [10] of BATMAN testbed of German IPP. The extraction domain is composed of bias plate (BP), PG and EG, and the simulation domain is cut off from behind the EG to reduce calculation time. Figure 2B is the distribution diagram of the aperture tunnel on the EG, showing the longitudinal section of the EG plate in the Y-Z direction; Figure 2C is the schematic diagram of the simulation domain in the X-Z direction. The Port structure in Figures 2B, C acts as the feed-in voltage and a boundary condition to truncate the simulation domain respectively. There are 126 apertures on the PG and EG plates. The inner radius of the aperture is 3 mm, the outer radius is 4 mm, and the extraction area is about 0.0063 m 2 . The movable magnetic boxes are placed 9 cm in front of the PG. Each box contains the "2X4" permanent magnet. The formed long-range weak magnetic filter field can reduce the electron density near the PG. The EG is embedded with permanent magnetic rods with alternating magnetization, forming a short-range strong magnetic deflection field in a direction perpendicular to the magnetic filter field [21]. In this paper, SMPM method is adopted to calculate the three-dimensional static magnetic field generated by permanent magnets. The bias current is measured on the PG [22], and the co-extracted electron current is measured on the EG [23]. The BP and PG are located 1 cm apart axially. The BP is connected to the chamber and at the same potential [24]. The feed-in voltage parameters are shown in Table 1. Since the time step is on the order of 1 picosecond, the calculation amount of 3D model is huge. In this paper, non-uniform grid is used to reduce the number of grids in the simulation domain. The maximum grid step is 4 mm and the minimum one is 1 mm. More simulation parameter settings are shown in Table 2.
The electron velocity conforms to the three-dimensional Maxwell distribution function, i.e., the electron velocity obeys Eq. 1, where, v is the electron velocity, m is the electron mass, kT is the electron temperature, and e is the natural constant.
References [26,27] show the different plasma potentials measured at different positions when electrons move from the expansion region to the extraction region, which indicates that the electron density is obviously non-uniform. The main reason for this phenomenon is the magnetic drift of electrons caused by the magnetic filter field in the expansion region. This phenomenon has also been seen and explained on the basis of our recent work [15]. Therefore, the distribution of electrons at the BP also presents a non-uniform phenomenon. Based on the previous calculation [15], the electron current density distribution on the BP is re-assigned as the initial setting for electron beam emission in the extraction region in this paper. The current density distribution of electron beam emission on the BP is shown in Figure 3. The emission area is 180 mm × 140 mm. The maximum electron current density

Ion source component Potential
Bias plate (connected with source body) −15kV Frontiers in Physics frontiersin.org can be up to 1,040 mA/cm 2 at non-uniformity and 136 mA/cm 2 at uniformity. Due to the density of the electron beam up to 1 × 10 16~1 × 10 17 m −3 , the interaction between electrons is ignored to avoid the emission surface being affected by the space-charge-limited current. Besides, grid step is obviously much larger than the Debye length (λ D ε0k b Te e 2 ne ), which can lead to artificial numerical heating [28,29]. That is, the energy of particles rises, especially for electrons, which will lead to an increase in the mobility of electrons and may cause an increase in local electron density. In our model, because the massive simulation domain, the calculation amount and calculation time have become unbearable, the interaction between electrons and electromagnetic fields is not considered. After the electron is emitted from the emission surface, it is a unilaterally forced state. Although such a setting avoids the effect caused by grid heating and roughly calculates the motion of the electron, it will inevitably lead to the difficulty in the simulation of the plasma sheath.
Moreover, the bias voltage will greatly affect the plasma sheath structure, thus affecting the co-extracted electron current [25].
Therefore, the method in [27] is used in this paper to consider the sheath potential drop near the wall and PG. If electrons reach the upper/side/bottom walls, those electrons with the energy Ek greater than the plasma potential Vpl are absorbed at the wall, while those with the energy Ek lower than Vpl are reflected. The PG plate is often coated with Cs layers to where, T b is the temperature of H − ions, set as 1 eV, j b is the emission rate of H − ions from the PG surface, set as 16 mA/cm 2 , j bmax is the space charge limited negative ion current density from the PG at the virtual cathode, i.e., the transported negative ion flux to the plasma. j bmax is related to the electron current density across the virtual cathode to the PG plate. V pl and U b are given by the experimental data in Ref. [25]. Similarly, if electrons flow to the PG, the electrons with energy Ek greater than the plasma sheath The potential distribution from the emission surface (Distance to emission surface = 0 mm) to the PG (Distance to emission surface = 10 mm). Frontiers in Physics frontiersin.org depth V sd are absorbed at the PG, and the absorbed electrons will form a bias current after reaching the steady state, while those with energy lower than V sd are reflected. The steady state of the whole simulation is determined by the steady state of the bias current. The potential distribution from the emission surface to the PG is shown in Figure 4.

Numerical calculation results
Comparison of magnetic filter field and magnetic deflected field with experimental data and FEM The placement schematics of the filter magnets and the deflection magnets are shown in Figures 5A, B, respectively. The filter magnets are magnetized along the +Z direction at a distance of 9 cm from the PG. The deflection magnets are embedded in the EG, and the magnetization direction of magnets is alternately magnetized along the X direction, which is perpendicular to the magnetic filter field. The magnetic field components distribution of By and Bz on PG at X = 14 mm are shown in Figures 5C, D, respectively. It can be found that the magnetic field near the apertures is mainly dominated by the deflection magnetic field. Such a setting is to allow the electrons to be deflected by the short-range strong deflection magnetic field, so as to avoid a large number of impacts on the PG plate, resulting in excessive thermal load on the PG plate, which is also the reason for such setting in the experiment. Figure 6 shows the distribution of the long-range weak magnetic filter field and the short-range strong magnetic deflection field near the PG along the X-axis, and the observation lines go across a beam axis of the PG aperture [21]. SMPM is used to calculate the 3D magnetic filter and deflection fields. The sign of the magnetic deflection field depends on the individual magnet row. As shown in Figure 6, the yellow dotted area represents the BP, which is also the region where electrons start. The magnetic field in the area before the BP (X < 0 mm) is set to zero. The solid green line area represents the PG, and the thickness of both BP and PG is 4 mm. The red and blue solid lines represent the calculated magnetic filter field B z and the magnetic deflection field B y on a beam axis of the PG The distribution of the long-range weak magnetic filter field and the short-range strong magnetic deflection field near the PG along the X-axis. Effect of magnetic field on the spatial distribution of electrons Figures 7A, B shows the spatial distribution of electrons in the X-Z direction and X-Y direction (U b 20.6 V, V pl 20.9V, V sd 0.3V), respectively. As shown in Figure 7A, electrons are relatively uniformly distributed alone the Z-axis through the magnetic field, but nonuniformly distributed along the Y-axis. This is due to the magnetic gradient drift of electrons caused by the magnetic filter field [15] during the diffusion of electrons from BP to PG, which is inevitable. Electrons passing through the PG will be accelerated by the high voltage between PG and EG, which is about 5 kV, and the energy of electrons arriving at EG is also on the order of 5 keV. There are deflection magnets embedded in the EG. The formed strong magnetic deflection field makes the electrons deviate from the original path and hit the EG wall, which makes it difficult to penetrate into the apertures in EG, so almost no electrons are extracted at the end of the EG tunnels. The magnetic field strength at the apertures of the EG (X = 21 mm) is about 86 mT, the electron gyration radius is about 2.8 mm, which is larger than the grid step of 1 mm, and the electron gyration period is 4.15e-10 s, which is much larger than the time step of 1.54 ps. So, in this model, within a time step, the FIGURE 8 Diagram of EEPF as a function of electron energy in front of PG.

FIGURE 9
The plasma potential V pl as a function of the bias voltage U b [25] (red solid line) and the relevant setting of the plasma sheath depth V sd in front of the PG in this model (blue solid line) under the BATMAN hydrogen working gas of 0.6 Pa.

FIGURE 10
The diagram of the bias current I bias and the co-extracted electron current density J e as a function of the PG sheath potential drop U b − V pl .

FIGURE 11
The placement diagram of the added magnets.  Figure 8 shows the diagram of electron energy probability distribution function (EEPF) as a function of electron energy in front of PG at 0.6 Pa. It can be seen that the simulated electron temperature is maintained at 1 eV by calculating the slope of the curve. It is in good agreement with 0.9 eV experimental data [34] at 0.5 Pa. It can be verified that the electron temperature is not affected by grid heating and energy distribution of electrons basically agrees with experimental data. Figure 9 shows the plasma potential V pl as a function of the bias voltage U b [25] (red solid line) and the relevant setting of the plasma sheath depth V sd in front of the PG in this model (blue solid line) under the BATMAN hydrogen working gas of 0.6 Pa. Figure 10 shows the diagram of the bias current I bias and the co-extracted electron current density J e as a function of the PG sheath potential drop U b − V pl . It should be noted that the points of the curve in Figure 10 correspond to the related points in Figure 9. If the PG sheath potential drop is set singly, the variation trend of the curve in Figure 10 cannot be obtained. As can be seen from Figure 10, the simulation results are in good agreement with the experiment [25], which proves the validity of the simulation setting.

Optimization of the magnetic field
Due to electron drift effect, non-uniform plasma potential is formed on the PG metal surface. The virtual cathode depth near the PG surface is also non-uniform, so that the non-homogeneous negative ion flux generated on the PG surface enters the plasma or is extracted by the PG apertures  where n e,top and n e,bot represent the electron density upstream and downstream of the apertures respectively. The smaller the value of the electron vertical asymmetry coefficient, the more symmetrical the electron distribution. The optimization of magnetic field is similar to the setting of multicusp magnetic field [27]. Four longilineal magnets are selected, and the size of them is 4 mm × 4 mm × 70 mm. The magnetization direction of the magnets magnetizes alternately along the Y direction. The residual magnetism is 0.1 T. The placement diagram of the magnets is shown in Figure 11. The magnets are placed in the middle of the electron beam from BP to PG. Figures 12A, B show the spatial distribution diagram of electrons in the X-Y direction before and after magnetic field optimization and C, D show the electron density distribution in the Y direction 1 mm before PG before and after optimization. Compared Figures 12A, B, the optimized spatial distribution of electrons is more symmetrical in the Y direction. Compared Figures 12C, D, it can be seen that the symmetry of electron density before optimization is poor, while the optimized electron density distribution is almost symmetric with respect to the center point of the Y-axis. Besides, the electron vertical asymmetry coefficient s ey in the region facing the apertures (−0.075 m < Y < −0.015 m and 0.015 m < Y < 0.075 m, represented by the orange dotted line) is decreased from nonoptimized 0.015 to about 0.002, which is favorable for the formation of uniform plasma potential, so that the negative ion flux generated on the PG surface will be more uniformly into the plasma or be extracted from the PG apertures.
The simulation results show that this magnet placement can effectively reduce s ey values for sets of values in the experimental Ref. [25] and improve the uniformity of electrons in front of the PG plate.

Conclusion
The total magnetic field near the apertures has a complex 3D structure due to the superposition of the magnetic filter and deflection fields. In this paper, based on the multi-aperture structure of three-dimensional RF negative hydrogen ion source BATMAN extraction region, the complex magnetic field topology is calculated by SMPM, and the plasma sheath depth is reasonably considered. The correlation between the plasma sheath potential drop and the bias current and the coextracted electron current density is tellingly reflected. The dynamic characteristics of electrons in the multi-aperture extraction region are successfully simulated under nonperiodic magnetic structure. In this paper, the magnetic field of the extraction region is optimized, and the electron uniformity is effectively improved, which has certain guiding significance for engineering application. However, the optimization of the magnetic field is not yet optimal, and the optimal solution remains to be further studied.

Data availability statement
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.