Simulation of neutral beam current drive on EAST tokamak

Neutral beam injection (NBI) is widely used in tokamaks and stellarators for heating plasma.^1–6 In addition to heating, the fast ions resulting from NBI can also drive electric current in plasmas.^7–10 To model this current, one needs to calculate the steady-state distribution of fast ions, and then integrate it to get the current. One method of doing this calculation is to use the Monte Carlo method to sample NBI fast ion source and following the guiding-center/full trajectories of fast ions, taking into account of their collisions with the background plasmas. Fusion community has developed many computer codes doing this kind of simulations, among which are NUBEAM,¹¹ OFMC,¹² ASCOT,¹³ ORBIT,¹⁴ SPIRAL,¹⁵ and many others.^16–19 An advantage of this method is that it can readily take into account the real space effects, such as finite orbit width and edge loss, which are usually approximately treated in analytical models^20–22 and some velocity grid based Fokker–Planck codes.

The fast ion steady-state distribution is of interest not only to current drive problem but also to research topics such as the interactions between fast ions and MHD modes and turbulence. Some mechanisms in the interaction may depend on the delicate phase-space structure of fast ions. Therefore, an accurate fast ion distribution function is of practical importance and the phase space structure needs to be more thoroughly studied and visualized than that has been done previously. In this paper, we carefully examine the phase-space structure of the steady-state fast ion distribution. Conventional wisdom has it that trapped particles do not carry toroidal current. This is roughly correct. However, in some circumstances, the asymmetry of the trapped particle distribution about the parallel velocity $v ∥$ can be so significant that they can carry significant toroidal current. This was demonstrated by Hager and Chang for edge plasmas with steep density gradient.²³ In the present work, we found that the asymmetry of trapped fast ion distribution about $v ∥$ is weak, so the trapped fast ion carry very small current, but still larger enough to be reliably observed in simulations. Specifically, we found that trapped ions carry co-current (relative to the main plasma current direction) current near the edge and countercurrent current near the core. The magnitude of the trapped ion current is one order smaller than that of the passing ions. Therefore, their contribution to the fast ion current is negligible (1% of the fast ion current).

We consider co-current NBI (all the four beams on EAST tokamak²⁴ are now in co-current direction). In order to operate EAST for longer pulse, there are interest in getting better neutral beam current drive (NBCD) efficiency by operating in optimized parameter regimes. In this work, based on realistic EAST configuration and plasma profiles, we examine the dependence of NBCD on two basic plasma parameters, namely, the plasma current I_p and plasma density n_e, using a Monte Carlo test particle code TGCO.^25,26

To the authors knowledge, there is no previous work discussing the dependence of NBCD efficiency on the plasma current. The results of our work indicate that, with the plasma current increasing, the fast ion current first increases and then decreases. This interesting non-monotonic dependence is not well known and need some explanations. We note that the drift orbit width of a fast ion is inversely proportional to I_p.²⁷ This means that the fast ion confinement improves with the increase in I_p. This may imply that NBCD efficiency also improve with the increase in I_p. However, the simulations in this paper indicate this is not the complete picture: the NBCD efficiency turns out to decrease with the increase in the plasma current in the larger I_p regime. This trend is found to be due to the fast ion trapped fraction increasing with the increase in I_p. As mentioned above, the trapped particle carries nearly zero current. Larger fraction of trapped particles implies smaller fast ion current.

With n_e increasing, the fast ion current decreases. The n_e dependence is related to the variation of the slowing-down time with n_e, which is already well known and is confirmed in our specific situation.

To get the net current, one needs to take into account the electron shielding effect, i.e., the current carried by electrons due to their response to the presence of fast ions. This generally requires to solve the steady-state Fokker–Planck equation for electrons with additional collision term corresponding to the electron collision with the fast ions. For current drive problem, we are interested in the first (parallel) moment of the electron distribution function. Then, making use of the self-adjoint property of the linearized collision operator, the electron response to arbitrary fast ion sources can be obtained by using the Green function method.^28–31 The final results of these studies are usually some fitting formulas for the ratio of net current to the pure fast ion current.^28,30 The present work uses these fitting formulas to include the electron shielding effect. The electron shielding model used in this work is a general model applicable to arbitrary collisionality regime and general tokamak flux surface shapes.^28–31 Our results show that the I_p and n_e dependence of the net current is similar to that of the fast ion current.

The simulations are performed for the EAST tokamak, which is a superconducting tokamak with a major radius $R 0 = 1.85 m$⁠, minor radius $a ≈ 0.45 m$⁠, typical on-axis magnetic field strength $B axis ≈ 2.2 T$ and plasma current $I p ≈ 0.5 MA$⁠.^24,32 Figure 1 plots the magnetic configuration and plasma profiles used in this work, which were reconstructed by the EFIT code³³ from EAST tokamak discharge #101473 at 4.5 s with constrains from experiment diagnostics.³⁴ The toroidal plasma current is in $+ ϕ ̂$ of cylindrical coordinate $( R , ϕ , Z )$⁠. The radial coordinate ρ_p used in this paper is the square root of the normalized poloidal magnetic flux: $ρ p = ( Ψ − Ψ 0 ) / ( Ψ b − Ψ 0 )$⁠, where $Ψ = R A ϕ$ is the poloidal flux function, $A ϕ$ is the toroidal component of the magnetic vector potential, $Ψ 0$ and $Ψ b$ are the values of $Ψ$ at the magnetic axis and last-closed-flux-surface, respectively.

In this work, we assume a Deuterium plasma with Carbon impurities, with the effective charge number of background ions being $Z eff = 2.23$ across the entire plasma. Simulations in this work are performed by using TGCO code,^25,26 which models neutral beam ionization, slowing-down, collision transport, and edge loss of the resulting fast ions. We consider Deuterium NBI with full energy $E full = 65 keV$⁠, and the particle number ratio between full, half, and 1/3 energy being $80 % : 14 % : 6 %$ after the beam leaving from the neutralization vessel. The beam power after leaving the neutralization vessel is fixed at 1 MW. We consider a reference case where NBI tangential radius is $R tan = 1.26 m$ and injects in the co-current direction. As a comparison, we also consider a modified scenario where the tangential radius is changed to $R tan = 0.731 m$⁠. The results presented below are for the $R tan = 1.26 m$ beam unless specified otherwise.

The neutral beam ionization is modeled by the Monte Carlo method.^11,25 Typical number of Monte Carlo markers initially loaded in the simulations is $1 × 10 5$⁠. The beam stopping cross Sec. data used in the simulation are from Ref. 38, which includes the charge exchange with thermal and impurity ions, impact ionization by electrons, thermal and impurity ions, and the multi-step ionization involving excitation states of neutrals. Beam ionization outside of the last-closed-flux-surface (LCFS) is ignored in the simulations.

Figures 2(a) and 2(b) plot the fast ion density in the poloidal plane (averaged over the toroidal direction) and in the toroidal plane (averaged over the vertical direction), respectively. Figures 2(c) and 2(d) plot the 1D histogram of the fast ions along R and $ϕ$⁠, respectively. The results indicate most fast ions are born near the low-field-side. The particle and power shine-through fraction in this case is 3.8% and 4.1%, respectively. Figure 2(e) plots the deposition density profile along the minor radius, which shows that the density reach its maximum at the magnetic axis, i.e., the beam deposition is on-axis.

It is often useful to use $( P ϕ , Λ )$ coordinates to describe the phase space, where $Λ = μ B axis / E$⁠, μ is the magnetic moment, E is the kinetic energy, and $P ϕ$ is the canonical toroidal angular momentum. In these coordinates, it is easier to classify orbit types. For example, whether a particles is passing or trapped can be readily identified. Figure 3 plots the fast ion birth distribution in $( P ϕ , Λ )$ plane, which indicates that most of the fast ions are passing particles (the region within the dashed line triangle is the trapped region). The fraction of trapped fast ions is 0.6% for this case. Figure 4 plots the birth distribution over the pitch $v ∥ / v$⁠, which shows that the dominant pitch is about –0.60 for the $R tan = 1.26 m$ beam, and as a comparison, –0.35 for the $R tan = 0.731 m$ beam.

Fast ions born from the beam ionization are transformed to guiding-center space and the guiding-center drifts are followed by using the 4th order Runge–Kutta scheme in the cylindrical coordinates $( R , ϕ , Z )$⁠. A fast ion is considered lost/thermalized when it touches the first wall or when it slows down to the energy of $2 T i ( 0 )$⁠, where $T i ( 0 )$ is the thermal ion temperature at the magnetic axis. Charge exchange loss of fast ions³⁹ is not included in the simulations. The loop voltage is well controlled to be near zero during the flat-top phase, indicating fully non-inductive. Therefore, the toroidal electric field is not included in the simulation. The collision model used in this work includes the effect of slowing-down, pitch-angle scattering, and energy diffusion (the details are provided in  Appendix B).

The finite Larmor radius (FLR) effect is included when (1) checking whether a fast ion touches the wall and (2) depositing guiding-center markers to spatial gridpoints to compute the distribution moment. The FLR effect is not included when pushing guiding-center trajectories since it has negligible effects.

In order to have a steady state on the slowing-down timescale, we need to include a continuous beam source. The method of including the continuous beam source is given in  Appendix A.

Figure 5 plots the steady-state fast ion density distribution in the poloidal plane (left-panel) and in the toroidal plane (right-panel). The results indicate that the fast ions density is roughly poloidally uniform and toroidally uniform for the 1 MW beam power. Note that the fast ion source is neither poloidally uniform nor toroidally uniform. The steady state fast ion distribution is not guaranteed to be poloidally uniform or toroidally uniform. It depends on the magnitude of beam power. For the 1 MW beam power considered here, the non-uniformity in the poloidal direction and toroidal direction is negligible. The poloidal uniformity is further confirmed in Fig. 6, which plots the radial profiles of various poloidal harmonics of the fast ion density. The results show that the m = 0 harmonic is dominant (m is the poloidal mode number).

Figures 6 plots the radial profiles of fast ion density, where we distinguish between trapped and passing particles. The radial profile of the trapped particle fraction is also shown is Fig. 6(c). The total trapped particle fraction in the steady-state distribution is 17%, which is significantly different from that in the birth distribution $( 0.6 % )$⁠. The pitch angle scattering can scatter particles between passing and trapped region of the phase space. Hence, it is expected the trapped fraction may be changed by collisions.

Figures 7(a) and 7(b) plot the fast ion toroidal current density distribution in the poloidal plane. Here we distinguish the current carried by passing ions and that carried by trapped ions.

The results indicate that current density of trapped fast ions is nonzero. This is the well known “banana current,”⁴⁰ which arises when the trapped ions' radial density profile has nonzero gradient. The banana current is an analog of the diamagnetic current (the latter is due to the gyro-orbits and is in the perpendicular direction, whereas the former is due to drift-orbits and is in the parallel direction⁴⁰) The banana current can be considered part of the bootstrap current.⁴⁰ 

The banana current mechanism can explain that there exists co-I_p current at a fixed radial location if the trapped ion density decreases radially outward. The results in Fig. 7 show that the trapped particle current density near the edge is in the co-I_p direction, which is consistent with the banana current mechanism. However, near the core region, the trapped current density is in the counter-I_p direction. The current direction reverse happens at $ρ p ≈ 0.6$⁠. The volume integrated trapped particle current is in the co-I_p direction, although being very small (only 1% of the total fast ion current).

Figures 7(c) and 7(d) plots the poloidal harmonics of currents carried by passing ions and trapped ions. For the passing ion current, the m = 0 harmonic is dominant, indicating poloidally uniform of the current. For the trapped ion current, the m = 1 harmonic becomes dominant in the out region $( ρ p > 0.6 )$⁠, indicating poloidal nonuniform.

Figures 8(a)–8(e) plot the steady-state fast ion distribution in $( E , v ∥ / v )$⁠, where E is the kinetic energy and $v ∥$ is the parallel (to the magnetic field) velocity. Both the 2D distribution and 1D distributions are plotted. We also distinguish between the passing particles and trapped particles. The results indicate that the trapped particle distribution is roughly symmetrical in $v ∥ / v$ around $v ∥ / v = 0$⁠, indicating it carries nearly zero parallel current. The results also indicate that trapped fast ions are more localized in the low energy region when compared with passing fast ions which have nearly uniform energy spectrum. There are three jumps in Fig. 8(d), which correspond to the NBI source at full, half, and 1/3 of $65 keV$⁠. Note that collisional energy diffusion makes some particles exceed the full energy, as is shown in Fig. 8(d).

Figure 8(f) plots the steady-state distribution in the $( P ϕ , Λ )$ plane. The fast ions seem to reach their peak density near the passing-trapped boundary. As a comparison, the birth distribution, as is shown in Fig. 3, has no obvious structure near the passing-trapped boundary. The structure of the steady-state distribution near the passing-trapped boundary is not sensitive to the fast ion birth profile. For instance, Fig. 9 considers the $R tan = 0.731 m$ beam and compares the birth distribution and steady-state one, which shows that the latter still reach its peak near the passing-trapped boundary.

Figure 11(a) plots the fast ion current as a function of the plasma current, which indicates that the dependence is not monotonic: in the lower I_p regime, the fast ion current I_f increase with I_p increasing whereas, in the higher I_p regime, I_f decreases with I_p increasing.

The increasing of I_f with I_p increasing in lower I_p regime is probably due to the improvement of fast ion confinement. The drift orbit width is inverse proportional to the poloidal magnetic field.²⁷ Larger plasma current implies stronger poloidal magnetic field, hence smaller drift orbit width and thus better confinement of fast ions. The evidence for this is shown in Fig. 11(b), which indicates that the fast ion loss fraction decreases rapidly with the plasma current increasing in the lower I_p region.

Next, we try to understand why the fast ion current decrease with plasma current increasing in high I_p regime. The only hint that we can identify is that the trapped fraction increase with the plasma current. The increase in the trapped fraction happens for both the birth distribution and the steady-state distribution, as is shown in Figs. 11(e) and 11(f).

Figure 12 plots the radial profiles of fast ion currents, the trapped fraction, and the current carried by the trapped fast ions, for various plasma currents. The results indicate again that trapped particles carry negligible current. Increasing the plasma current makes the trapped fraction increased across a wide radial region except the edge.

We also performed similar simulations as above but with a smaller NBI tangential radius, $R tan = 0.731 m$⁠. This case has larger trapped particle fraction. The results are plotted in Figs. 13 and 14. We observe the same non-monotonic dependence of the fast ion current on the plasma current. Also the results indicate that the trapped particles carry nearly zero current.

In the above scanning of plasma current, the density and temperature profiles are kept fixed. Therefore, the force-balance is not guaranteed. This is a weakness of this work. Similarly, the equilibrium is not re-computed in the density scanning presented in the next section.

Next, we examine the plasma density dependence of fast ion current. The density affects both the ionization process and the fast ion collisional transport. The ionization process affects the shine-through fraction and ion birth location. The latter influences the first-orbit loss fraction and the ratio between trapped and passing fast ions. The collision transport influences the fast ion distribution in both position space (edge loss) and velocity space (slowing-down and passing-trapped transition). These effects may have opposite effects on the fast ion current. Therefore, it is not obvious what is the trend of fast ion current changing with the plasma density. Next, we examine this trend via simulations.

We scan the electron density via scaling the profile in Fig. 1 by a factor of values ranging from 0.2 to 1.8. Figure 15(a) plots the dependence of the volume integrated fast ion current on the density, which indicates that the current decreases with increasing of the density. Figure 15(b) shows that the shine-through loss decrease with the density increasing (as is expected). Also, Fig. 15(c) shows that the ion edge loss fraction decreases with the density increasing. These reduction of loss is beneficial for current drive since more fast ions can stay in the plasma to contribute electric current. On the other hand, Figs. 15(f)–15(h) shows that, with the density increasing, the slowing-down time becomes shorter and the trapped particle fraction becomes larger. These are deleterious effects for current drive. Shorter slowing-down time implies that fast ions can remain energetic for shorter time and thus fewer fast ions can contribute to the fast ion current. Larger trapped fraction means less particles can efficiently contribute to the current. The results in Fig. 15(a) indicates that the deleterious effects turn out to defeat those beneficial effects, making the fast ion current decreases with the density increasing.

Figure 16 plots the radial profiles of fast ion current density. The results indicate that, for all radial positions, the fast ion current density decreases with the increasing plasma density.

Taking into account the electron shielding effect on the fast ion current does not qualitatively change the dependence of the driven current on the plasma current, i.e., net current follows a similar trend as the fast ion current, as is shown in Figs. 17 and 18. For the reference case (⁠ $I p / I p , ref = 1$ and $n e / n e , ref = 1$⁠), the net current is weaker than the fast ion current by a factor of 0.7.

Figure 19 plots the radial profiles of various quantities related to the driven current, namely, the pure fast ion current density J_f, the electron shielding factor F, the net current density $J net = J f F$⁠, the effective trapped electron fraction f_t, the electron collision frequency ν_e, the thermal electron bounce frequency ω_b, and the normalized electron collision frequency $ν e ⋆$⁠. (The details of these quantities are given in  Appendix C.) The formulas for the shielding effect used here are valid for general tokamak equilibria and arbitrary collisionality regime,^30,31 where the equilibrium shaping effects are included via the effective trapped fraction f_t.

As a comparison, we also plot the shielding factor F for the case of $ν e ⋆ = 0$ and the resulting net current. The results show that the $ν e ⋆ = 0$ approximation slightly overestimates the net current.

Simulations of neutral beam current drive on the EAST tokamak were performed, providing detailed information about the fast ion distribution in both real space and velocity space. We distinguish between current carried by passing particles and that carried by trapped particles, and found that trapped ions carry co-current current near the edge and countercurrent current near the core. However, the magnitude of the trapped ion current is one order smaller than that of the passing ions, making their contribution to the fast ion current negligible.

We examine the dependence of the fast ion current on the plasma current I_p. With I_p increasing, the drift orbit width decreases, which helps reduce the first-orbit loss and collisional transport, and thus is beneficial for fast ion confinement. This mechanism makes the fast ion current increase with I_p in the low I_p regime. Yet/However for high I_p regime, the effect of trapped fraction increasing becomes dominant. Since the trapped fast ions carry nearly zero current, the increasing of trapped fraction then implies the decreasing of fast ion current.

The fast ion current decreases with the increase in plasma density n_e. This dependence is mainly determined by the variation of the slowing-down time with n_e, which is already well known and is confirmed in our specific situation.

Youjun Hu thanks Youwen Sun, G. S. Xu, Yingfeng Xu, Lei Ye, Yifeng Zheng, and Baolong Hao for useful discussions, and thanks Xiaotao Xiao for careful reading of the first version of this paper. The manuscript of this paper was written in GNU TeXmacs, a free structured WYSIWYG editor for scientists.⁴¹ Numerical computations were performed on Tianhe at National Super-Computer Center in Tianjin, Sugon computing center in Hefei, and the ShenMa computing cluster in Institute of Plasma Physics, Chinese Academy of Sciences. This work was supported by Users with Excellence Program of Hefei Science Center CAS under Grant No. 2021HSC-UE017, by Comprehensive Research Facility for Fusion Technology Program of China under Contract No. 2018–000052-73–01-001228, and by the National Natural Science Foundation of China under Grant No. 11575251.

The authors have no conflicts to disclose.

Youjun Hu: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (lead); Methodology (lead); Resources (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). Xingyuan Xu: Formal analysis (supporting); Investigation (supporting). Yunchan Hu: Data curation (supporting); Formal analysis (supporting); Resources (supporting). Kaiyang He: Conceptualization (supporting); Data curation (supporting); Formal analysis (supporting). Jinfang Wang: Conceptualization (supporting).

The data that support the findings of this study are openly available in Zenodo at https://doi.org/10.5281/zenodo.8176352, which archives one version of TGCO code in https://github.com/Youjunhu/TGCO.⁴² 

To obtain the steady state of fast ions, we need include the continuous beam source. A straightforward Monte Carlo implementation of this continuous injection would be to introduce new Monte Carlo markers to represent the newly injected physical particles at each time step. This method is computationally expensive. For a time-independent background plasma, there is an efficient method that involves only a single injection and then utilizes the time shift invariant to infer the contribution of all the other injections. This method is illustrated in Fig. 20.

The above method works only for a time-independent background plasma and constant beam power. For time-dependent background plasma, re-injecting new Monte Carlo markers seems to be the only method available. The present work considers a time-independent background plasma with a constant beam power and use the above efficient method.

Both $τ s e$ and v_c have a radial dependence via their dependence on the plasma density and temperature. The radial profile of τ_s for the plasma specified in Fig. 1 is plotted in Fig. 21 for fast ions of $65 keV$ to be slowed down to a cutoff velocity [chosen as $m v 2 2 / 2 = 2 T i ( 0 )$ in this article]. The result shows that typical slowing-down time in the core region is about $60 ms$⁠.