Numerical study of stripping and stray particles for the one-RF-driver NBI negative ion source prototype of CFETR

The prototype for a negative hydrogen ion source for neutral beam injection of China Fusion Engineering Test Reactor is being developed at the Southwest Institute of Physics. To study the physics of negative ion beam transport and to optimize the design of the source, the stripping loss and the stray particles’ impacts on the one-RF-driver prototype are analyzed. Collision simulation, including both the beam-gas collisions and the particle-grid collisions, is carried out basing on the results of gas flow evaluation and particle tracing. The stripping loss, the distribution of stray particles and the heat loads are calculated, comparing two configurations of grounded grid (GG) (multiaperture or multislot). At the source filling pressure of 0.3 Pa and the vessel pressure 0.05 Pa, the extraction voltage being 8 kV and the acceleration voltage 200 kV and the extraction grid (EG) magnet peak of ±45 mT, the stripping loss of the 200 A m−2 H− beam can be reduced from 25% to 20% by changing GG from multiapertures to multislots. The H− proportion in the total current at 40 mm after GG, however, shows smaller change than the reduction of the stripping loss possibly because the multislot GG’s larger transparency increases the chance for the stray particles to pass through GG. The total heat load on EG in the two cases with different GG configurations are both around 66 kW, while the GG heat load is reduced from 45 kW for multiapertures to 17 kW for multislots. The study provides good comprehension of the transport process and useful guidance for practical operations.


Introduction
Neutral beam injection (NBI) is applied for the auxiliary heating of plasma in fusion. The under-plan China Fusion Engineering Test Reactor (CFETR) [1] employs the negative-ion-based NBI because of the retaining high neutralization efficiency of negative ions with the energy up to 1 MeV. The prototype of CFETR NBI negative hydrogen ion source is being developed in Southwest Institute of Physics [2]. For now, the first setup is equipped with only one RF plasma driver. The target is to produce a 200 A m −2 /200 keV H − beam. The prototype will help preliminarily understand the physics of negative ion beam and provide experience for the future development of the NBI negative ion source for CFETR.
For the transport of negative ion beam, one important issue is the 'stripping loss'. When the ion source operates, the negative ions have the chance of colliding or reacting with the gas molecules or atoms. Due to the low affinity, the electrons on the H − ions can be easily stripped and will cause the loss of H − . The stripping loss should be limited in the design phase.
Another issue for negative ion source is the impacts of the stray particles. This involves the co-extracted electrons from plasma and the products (electrons, positive ions and atoms) from the beam-gas collisions as mentioned before, which are generated in the volume. It also concerns the secondary emission electrons from the grid surfaces, triggered by the particlegrid collisions. These stray particles could spoil the beam quality when flowing out together with the negative ions. They could also cause heat loads if they collide on the grids. In particular, the positive ions flowing backwards may damage the source [3,4]. For the efficiency and the safety of the ion source, it is important to evaluate the impacts of the stray particles.
Stripping and stray particles in the negative ion source have been studied using Monte Carlo programs like Krylov's 3D program [5] and EAMCC [6], which were based on the electromagnetic field results that had been calculated previously by different particle tracing programs. They have also been studied using suites of multiphysics models developed in COMSOL [7][8][9], where the relevant factors (gas profile, electric and magnetic field, beam trajectory and collision) were calculated systematically in one program. In this work, the model in [9] is adopted and updated to analyze the stripping loss and the impacts of the stray particles for CFETR NBI H − source's one-RF-driver prototype. Collision simulation, including both the beam-gas collisions and the particle-grid collisions, is carried out basing on the results of gas flow evaluation and particle tracing. The stripping loss, the distribution of stray particles and the heat loads are evaluated, comparing two configurations of the grounded grid (GG) (multiaperture or multislot).
The paper develops as follows: An overview of the one-RF-driver prototype is given in section 2. Then, the basic principles of the calculations are illustrated in section 3. The study of the stripping loss and that of the stray particles' impacts are discussed separately in sections 4 and 5. Summarization is given in section 6.

Source description
The one-RF-driver prototype is composed of the plasma driver, the expansion chamber, the grid system and the vessel. The sketch is shown in figure 1. Plasma is generated in the driver by RF power. Ions are extracted and accelerated into multiple beamlets through the grid system. The grids are all divided into two segments and each segment has eight rows of 17 apertures. It is expected to extract H − beam in 200 A m −2 (∼8 A in total).
The grid's configuration is shown in figure 1 too. The system consists of four grids: the plasma grid (PG), the extraction grid (EG), the steering grid (SG) and the GG. PG is made of molybdenum and the others are made of copper. PG apertures are chamfered at both sides. The plasma facing side of PG is chamfered to enlarge the surface area for plasma reactions so that the production of H − ions is augmented. Also, it enhances the extraction probability of the ions with a better starting angle [10]. The downstream side of PG is chamfered to better focus the beam [11]. EG is followed by SG with a reduced aperture radius and these two grids are at the same potential. This design is learnt from JT-60U [12] to decrease the leakage of electrons. In addition, permanent magnets are inserted into the E.g. to deflect the co-extracted electrons, with their pole orientations altering row by row. Two configurations of GG are considered: one is the multiaperture grid and the other is the multislot grid. Referring to [13], the multislots on GG are thought favorable for improving the vacuum and lowering the stripping loss. This is to be verified later.

Principles of the study
The study follows the flowchart in figure 2 updated from [9]. Preconditions are firstly evaluated by the molecular flow evaluation and particle tracing. The gas profile, the electromagnetic field and beam trajectory are acquired and the coextracted electrons' trajectories are also calculated. Then, the beam-gas collisions are simulated basing on the preconditions results. The stripping loss can be estimated and the products from beam-gas collisions are traced. To take the secondary emission electrons also into consideration for the stray particles' impacts, the particle-grid collision simulation is added into the study.

Preconditions evaluation
Angular coefficient method is adopted to evaluate the gas profile under the molecular flow assumption. Only the interactions between particle and surface are considered. The basic principle of the method is illustrated in the subplot (a) of figure 2. The gas molecule flux, defined as the molecule flow rate per unit area (unit: 1/(m 2 s)), is assumed to be adsorbed and reemitted diffusely from surfaces. The angles of the emitted particles follow the cosine law. With a flux ψ ′ emitted from the surface S ′ , the flux ϕ received by the surface S is given by where A is known as the angular coefficient and is only determined by the geometry. The profile of the gas flow throughout the system can thus be learnt by calculating the molecular flux arriving at the wall of the chambers. Macroscopic variables like the gas density n and the pressure p are then derived according to the kinetic theory. Since the geometry's angular coefficient does not change with vacuum conditions, the method makes it fast to study the gas profile in a given system under different conditions.
To solve equation (1), boundary conditions should be defined. The emitted flux ψ of the inlet or outlet surface of the vacuum system can be derived from the pressure p s at the surface, according to the relation known as: where m is the particle mass, k B Boltzmann constant and T the gas temperature. The outlet can also be defined as a pump, Figure 2. Flow chart of the study [9]. Subplots are: (a) basic principles of angular coefficient method and (b) algorithm of particle tracing. whose effect is equalized with the transport coefficient f that indicates the probability for a particle to pass the boundary. f is derived from the ratio of the pump speed C and the ideal maximum pump speed C max that can be realized on the same surface area S: f = C/C max = √ 2πm/(k B T)C/S. The relation between the flux of the emitted and the incident particles at the pump boundary is: ( Walls are set to have the emitted flux equal with the incident flux: The outgassing on the wall is not considered. The electromagnetic field in the system and the trajectory of beam are learnt by particle tracing. The algorithm is illustrated in the subplot (b) of figure 2. The magnetic field B and the electric field E are evaluated by solving each's Laplace equation with finite element method: where Λ represents the scalar magnetic potential φ m for magnetic field and the electric potential φ for electric field. The beam's movement (velocity v) is calculated by solving the motion equation basing on the electromagnetic field results: where q is the particle's charge. The magnetic field is determined only by the arrangement of EG magnets because the beam's magnetic effect is negligible, while the electric field should include the contribution of the beam space charge. So the calculations are iteratively executed between solving the electric field's Poisson equation and the beam's motion equation starting from the second turn. The magnetic field does not change during the iteration. The distribution of the space chargeρ is derived from the beam trajectory and used as the source term of Poisson equation: where ε 0 is vacuum permittivity. The iteration reaches convergence when the difference of electric field between two consecutive turns is within limit.

Beam-gas collision simulation
The beam-gas collisions are examined after having reached a stable beam in the previous calculations. It is assumed that the influence of the collisions on the electric and magnetic field is negligible. So, by employing the calculated gas profile and simulating the beam's movement in the stable electric and magnetic field again, the detection of collision with Monte Carlo method is started in the simulation. The occurrence of a collision between beam and gas is examined at each time step. The collision frequency α (t) is given by: where α i (t) is the frequency at which an individual reaction type i happens in one collision and is determined by the gas density n, the collision cross section σ and the relative velocity v r between particle and gas: The probability P that a collision happens at time t is then derived: Employing a uniformly distributed random number U sampled from the interval (0, 1), the collision is detected by comparing U and P. The collision happens when: Equation (11) can be simplified using the equivalence between the probability distribution of U and 1−U, and becomes: When the collision is accepted, the specific reaction type is determined according to the probability P i given by: Only the reactions between H − and H 2 are considered in this work. The products of the reactions are traced. In the program, the collision between H − and H 2 is modeled as the binary collision. The pre-collision velocity of H − is known from the previous particle tracing simulation. The precollision velocity of H 2 molecule follows Maxwell distribution at room temperature. When a collision is detected, the program calculates H − 's post-collision velocity v ′ n and H 2 's post-collision velocity v ′ g . However, because H − or H 2 should be removed after the reactions, their velocity is given to their products. For the stripping reactions, v ′ n is set as the initial velocity of the ion products (H, H + ) and v ′ g is for the stripped electron. For the gas ionization reaction, v ′ n is given to the electron and v ′ g is for the positive ion H 2 + .

Particle-grid collision simulation
Updated from the original process of [9], the particle-grid collisions are simulated additionally. The secondary emission induced by the collisions between the electrons (from plasma or from beam-gas collisions) and the grids is considered. The simulation learns the manner of [6]. Materials of the grids are set as molybdenum for PG and copper for the others. For the incident electrons with energy around tens kilo-electron-volt, the secondary emission mainly involves the true secondary electrons (indicated as j = 1) and the backscattered electrons (j = 2). The emission probability P se,j (j = 1, 2) follows: where θ e is the acute angle between the incident velocity and the surface's normal direction. P 0,j and A 0,j are the coefficients varying with the incident particle's energy and are taken from figure 3 and equation (9) of [6]. The emission is assumed isotropic. Though the initial energy of the emitted electron depends on the incident particle's state, it is set as constant here for both kinds of the emitted electrons due to the limit of the program. As the typical energy spectra of the true secondary electrons is low (0-50 eV) [14], they are emitted with energy 10 eV following [6]. The backscattered electron's energy can be comparable to the incident electron's energy, depending on the collision angle [15]. Their initial energy is set as half of the incident electron's energy in the simulation for simplicity.
With the statistics of all the stray particles, the particle distribution and the heat loads on the grids are evaluated. The original process of [9] has been applied for the beam loss analysis of the HUST negative ion source whose design takes BATMAN [16,17] as reference. In [9], the calculation results of the pressure corresponded with the measurements of the HUST source, and the beam loss estimation was close to the experiment results of BATMAN under similar conditions. The reference manner [6] for the updated particle-grid collision model has also been proved reliable in many cases [18][19][20]. It is reasonable to say that the whole process here can be trusted.

Stripping study
In the following calculation (both sections 4 and 5), only the bottom half of the grid system is modeled considering the symmetry between the two segments of the grids. Both configurations of GG are considered. For simplicity, the grid system with the multiaperture GG is referred as 'the aperture-GG case' and that with the multislot GG as 'the slot-GG case'.

Gas profile and field
To calculate the gas profile, the boundary at 10 mm in front of PG is set as the inlet with a given pressure p driver . The boundary at 40 mm after GG is set as pump boundary with f = 0.3 (equivalent to a pumping speed about 11 000 l s −1 for H 2 ). It is assumed that the gas does not leak through lateral structures, so the grid surfaces as well as the lateral boundaries are set as wall boundary. The boundary conditions for the molecular flow evaluation are illustrated in figure 3.
The pressure in the driver without plasma discharge is recommended to be 0.3 Pa for the NBI negative ion source in International Thermonuclear Experimental Reactor (ITER) in order to keep the stripping loss below 30% [21]. Given p driver = 0.3 Pa (gas density of 7.4 × 10 19 m −3 ) at the inlet, the corresponding gas profile along the central aperture axis l at room temperature of both cases are presented in figure 4. It shows that the pressure drops quickly at the grids. The pressure after GG in both cases are slightly over 0.05 Pa (gas density 1.2 × 10 19 m −3 ). However, the multislot GG is proved to allow a better vacuum at the acceleration gap. This predicts a smaller value for the stripping loss in the slot-GG case.
The beam density is set 200 A m −2 in the calculation. The electrical conditions for the optimal beam optics given in [22] are employed: the extraction voltage is set 8 kV and the acceleration voltage is 200 kV. Ions are emitted from a semiellipsoid surface with the same radius of PG aperture and the depth of 1 mm. The ions have the initial energy of 1 eV and their initial velocity is normal to the emitter surface. The emitter is set in the PG aperture near the chamfer tip and is at a potential 0.5 V higher than PG's. Each EG magnet has the remanence of 0.9 T and their pole orientations are parallel with the axis and opposite between rows. The calculated electric field and EG magnets field along the central aperture's axis l are given in figure 5. Only small difference in the electric field near GG is observed between the two configurations. EG magnets produce a pair of magnet peaks of ±45 mT in the vertical (y) direction.

Stripping loss
The three reactions: single stripping (a), double stripping (b) and ionization (c) indicated in figure 6 are the main reactions for the collision between H − and H 2 gas. The cross section of each reaction varies with H − ion's energy and is interpolated from the data in [23]. As can be seen, single stripping has the largest cross section and plays the most important role all along the beam route. Double stripping usually takes the smallest part among the three reactions and the H 2 ionization becomes much significant mainly in the energy range from 10 keV to 100 keV.
The evolution of the stripping loss s can be estimated by integrating the product of the gas density n and the cross section σ s of stripping along the route: s =ˆnσ s dl (15) where σ s = σ a + σ b with σ a the cross section of reaction a and σ b the reaction b. Basing on the gas profile in figure 4 and the electromagnetic field in figure 5, the integration along the beam route l (figure 3) for the two cases is shown in figure 4. It is seen that the rate of the loss varies along the route. In both cases, the place where the loss increases the most rapidly is at  PG-EG gap because the gas there is dense. At EG-GG gap, the loss rate has slowed down, as the gas density decreases and the collision cross section is smaller under higher particle energy. The estimated final stripping loss in the slot-GG case is around 21% compared to that around 27% in the aperture-GG case. The multislot GG is verified to lower the stripping level by 6% than the multiaperture GG.
In fact, the more precise values of the reactions' occurrence can be given by the beam-gas collision simulation. Reactions are examined with Monte Carlo method. According to the simulation, the proportion of the particles involved in stripping is 25%: reaction (a) at 24% and (b) at 1% for the aperture-GG case and 20%: (a) at 19% and (b) at 0.9% for the slot-GG case. The simulated absolute value of the stripping loss in each case is lower than that estimated by equation (15), while the simulated relative difference between the two cases is close to equation (15)'s result. The simulation also gives the ionization (c) level: 3% for aperture-GG and 2.2% for slot-GG.
The stripping loss in the slot-GG case is further simulated under different pressure conditions to find the proper ones that can keep the loss within limit, using the same beam and the same electromagnetic field conditions as previous. Figure 7 shows the results. Under the same driver pressure, the vessel pressure is higher at smaller pump (lower f ). The stripping loss increases with the pressure in the system as can be expected. The dash lines are fitted from the simulation results. Their nonchange slope is due to the system's fixed flow conductance. It is learnt that if the stripping loss here is to be limited below 30% with f = 0.3 at the outlet, the driver pressure must be kept under 0.5 Pa.

Stray particles study
Stray particles are studied, including the co-extracted electrons from plasma, the products from beam-gas collisions and the electrons from secondary emission. The calculation is still based on the precondition results in figures 4 and 5 and the two GG configurations are considered. Typical results of the stray particles' trajectories, their distribution statistics, the heat loads on the grids are shown in figures 8-11. Unless indicated, the results are all doubled from those calculated with the half model.

Co-extracted electrons.
Assuming that the coextracted electrons are released in the same way and having the same current as the ions, figures 8(a) and 9(a) presents their trajectories. They do not show much difference between the two cases: EG magnets effectively deflect the co-extracted electrons at the extraction gap. This avoids them from accelerating to the full energy. However, a large number of collisions are induced on the EG upstream surface. By calculation, the power deposition of the co-extracted electrons in the two cases are close, about 64 kW. Due to their large amount, the co-extracted electrons are the main source of the EG heat load.

Beam-gas collision products.
The product particles (electrons, positive ions and atoms) from the beam-gas collisions are traced in the simulation and their trajectory in the two  In general, the movements of the products in the two cases are similar. Electrons are generated by all three reactions. Due to the presence of EG magnets, the product electrons at the extraction gap can largely be prevented from passing through EG. However, the influence of the magnetic field is less significant at the acceleration gap. The product electrons generated at the exit of EG are deflected less and some could be transmitted out with the beam. Those generated at the acceleration gap are mostly accelerated out directly.
In the aspect of the positive ions (H + , H 2 + ), they are mostly accelerated in the opposite direction to the negative ions and the electrons and they flow to the upstream components of the source. The flow's max energy is up to 200 keV due to the positive ions generated near GG. It carries the power about 25 kW at PG upstream in the aperture-GG case and 15 kW in the slot-GG case. Apart from the back flow, the rest positive ions fly downstream. They are either generated after GG where the field is negligible, or they are generated from double stripping before GG with enough energy to overcome the decelerating field. As for H atoms, they are not shown in the figures for clearness. Their trajectories are not influenced by the fields. They are either sent out or intercepted by the grids.
The difference between the two cases is mainly noticed in the proportion of the particles distributed in different areas. For each product particle category, the proportions of those flowing out are monitored at 40 mm after GG and are summarized in figure 10. It is seen that GG with multislots allows more electrons and atoms to pass through. About 47% (1 A) of the current of the product electrons flow out in the aperture-GG case and it is 50% (0.85 A) in the slot-GG case. The flowing-out atoms are 74% and 80% of their total numbers for aperture-GG and slot-GG respectively. The slots also benefit for the transmission of positive ions in the backward direction, so the flowing-out positive ions' proportion in the slot-GG case (about 11%, 27 mA) is smaller than the aperture-GG case (13%, 40 mA). In addition, the positive ions flowing to the source take up 80% (245 mA) of their total current and 7% (22 mA) are intercepted by the grids in the aperture-GG case. For the slot-GG case, they are 85% (220 mA) and 4% (10 mA) respectively.
The beam purity is calculated as the H − proportion in the total beam current (sum of the absolute value of H − , electron and positive ion currents) at 40 mm after GG and is 85% for the aperture-GG case and 88% for the slot-GG case. The difference is not as large as that of the stripping loss. This can be explained by the fact previously observed that the larger transparency of the slots increases the chance of passing through the grid and makes up for the reduction of the stray particles' generation. As a result, the multislot GG provides the lower stripping loss but not necessarily improve the beam purity by an equal level.

Secondary emission electrons.
The co-extracted electrons and the electrons from the beam-gas collisions induce secondary emission on the grids. Due to the large amount of the co-extracted electrons, EG surface is where the secondary emission mainly takes place, as figures 8(c) and 9(c) show. Some of the secondary emission electrons are depressed back to the surface, while others can be accelerated downstream. The current of the secondary emission electrons from EG is about 6 A in both cases, while the current of those from GG is ∼100 mA in the aperture-GG case and ∼60 mA in the slot-GG case. The flowing out current is about 32 mA and 16 mA for aperture-GG and slot-GG respectively. Figure 11 shows the heat load on EG and GG. It is observed that the load of EG concentrates on a narrow part beside each aperture and is opposite between rows because of the altering arrangement of magnets' pole orientations. The load on GG is point-like and more distributed. The maximum power density on GG is about 7.5 MW m −2 in the aperture-GG case and 4.5 MW m −2 in the slot-GG case, both higher than their maximum on EG (∼4 MW m −2 ). Although the amount of the incident particles on GG is much smaller than that on EG, the particles' energy after acceleration can make GG's localized heat load considerable. The total heat load on each grid is acquired. The heat load on EG is around the same value ∼66 kW in both cases, caused by the co-extracted electrons (∼64 kW) and the product particles from the beam-gas collisions (∼2 kW). Their secondary emission electrons take up 18 kW on EG. GG's heat load is all caused by the product particles. It is about 45 kW in the aperture-GG case and 17 kW in the slot-GG case. The secondary electrons take up 8 kW and 3 kW respectively. It is learnt that the slots reduce a large quantity of the load on GG.

Summary
Beam stripping and stray particles in the one-RF-driver prototype for CFETR NBI H − source are studied with a suite of multiphysics models. Half of the grid system is modeled, considering two configurations of GG. At the source filling pressure of 0.3 Pa and the vessel pressure 0.05 Pa, the voltages being 8 kV/200 kV and EG magnet peak of ±45 mT, the stripping loss is evaluated to be 20% for the 200 A m −2 H − beam when using the multislot GG, reduced from 25% when using the multiaperture GG. For the same voltage and magnet conditions in the case with the multislot GG, the stripping loss can be kept below 30% with the equivalent pumping speed 11 000 l s −1 at 40 mm after GG if the filling pressure is under 0.5 Pa.
Stray electrons and atoms are either absorbed by the grids, or transmitted out together with the negative ion beam. Positive ions from the beam-gas collisions mostly (80%-85%) flow backwards to the source. On the same preconditions as mentioned above, the H − proportion in the total beam current at 40 mm after GG shows smaller change between the two cases than the reduction of the stripping loss. A possible explanation is the larger transparency of slots increases the chance for the stray particles to pass through at the same time when decreasing the stripping loss. The total heat load on EG is around 66 kW in both cases, while the GG heat load is about 45 kW for multiapertures and is reduced largely to ∼17 kW for multislots.
The study gives a preliminary view on the stripping loss and the stray particles' impacts for the one-RF-driver prototype considering the reactions between H − and H 2 . However, H − can also interact with H atoms produced in the dissociation of plasma and it is not included in the model. This should be considered in the future. Other impacts, such as the thermal mechanics issue of the grids and the influence of the power deposition of positive ions on the driver, should also be learnt. Despite the decrease of stripping loss, the effects of the multislot GG on beam optics remains to be studied.

Data availability statement
The data generated and/or analysed during the current study are not publicly available for legal/ethical reasons but are available from the corresponding author on reasonable request.