Simulation of neutral beam current drive on EAST tokamak

The simulations are performed for the EAST tokamak, which is a superconducting tokamak with a major radius $R_{0}=1.85m$, minor radius $a\approx 0.45m$, typical on-axis magnetic field strength $B_{\operatorname{axis}}\approx 2.2T$ and plasma current $I_{p}\approx 0.5\operatorname{MA}$Wan_2017 (30, 22). Figure 1 plots the magnetic configuration and plasma profiles used in this work, which were reconstructed by the EFIT codelao1985 (31) from EAST tokamak discharge #[email protected] with constrains from experiment diagnosticshuyunchan2023 (32). The toroidal plasma current is in $+\hat{\mathbf{\phi}}$ of cylindrical coordinate $(R,\phi,Z)$. The radial coordinate $\rho_{p}$used in this paper is the square root of the normalized poloidal magnetic flux: $\rho_{p}=\sqrt{(\Psi-\Psi_{0})/(\Psi_{b}-\Psi_{0})}$, where $\Psi=RA_{\phi}$ is the poloidal flux function, $A_{\phi}$ is the toroidal component of the magnetic vector potential, $\Psi_{0}$ and $\Psi_{b}$ are the values of $\Psi$ 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_{\operatorname{eff}}=2.23$ across the entire plasma. Simulations in this work are performed by using TGCO codeyoujunhu2021 (23, 24), 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_{\operatorname{full}}=65\operatorname{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 1MW. We consider a reference case where NBI tangential radius is $R_{\tan}=1.26m$ 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.731m$.

The neutral beam ionization is modeled by the Monte-Carlo methodpankin2004 (11, 23). Typical number of Monte-Carlo markers initially loaded in the simulations is $1\times 10^{5}$. The beam stopping cross section data used in the simulation are from Ref. Suzuki_1998 (36), 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.

Figure 2(a) and (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. Figure 2(c) and (d) plot the 1D histogram of the fast ions along $R$ and $\phi$, 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_{\phi},\Lambda)$ coordinates to describe the phase space, where $\Lambda=\mu B_{\operatorname{axis}}/E$, $\mu$ is the magnetic moment, $E$ is the kinetic energy, and $P_{\phi}$ 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_{\phi},\Lambda)$ 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_{\parallel}/v$, which shows that the dominant pitch is around $-0.6$.

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,\phi,Z)$. A fast ion is considered lost/thermalized when it touches the first wall or when it slows down to the energy of $2T_{i}(0)$, where $T_{i}(0)$ is the thermal ion temperature at the magnetic axis. Charge exchange loss of fast ions Kramer_2020 (37) 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 time scale, 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 1MW 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 1MW beam power considered here, the non-uniformity in the poloidal direction and toroidal direction is negligible. The poloidal uniformity is further confirmed in Figure 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).

In Fig. 6, we distinguish between trapped and passing particles. The radial profile of the trapped particle fraction is also shown is Fig. 6c. 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 7a-b plot the fast ion toroidal current density distribution in the poloidal plane. Here we distinguish the current carried by passing particles and that carried by trapped particles. The results indicate that the trapped particle current is negligible. Also note that the trapped particle currents in the core region and near the edge have opposite directions, with the edge current being in the co-current direction. This phenomenon has connection with the so-called “banana current”, which is an analogue 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 directionPeeters_2000 (38)). The banana current can be considered part of the bootstrap currentPeeters_2000 (38). In the case shown in Fig. 7, the current direction reverse happens at $\rho_{p}\approx 0.6$. The volume integrated trapped particle current is in the co-current direction, although being very small (only 1% of the total fast ion current). Figures 7c-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 $(\rho_{p}>0.6)$, indicating poloidal nonuniform.

Figures 8(a)-(e) plot the steady-state fast ion distribution in $(E,v_{\parallel}/v)$, where $E$ is the kinetic energy and $v_{\parallel}$ 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_{\parallel}/v$ around $v_{\parallel}/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\operatorname{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_{\phi},\Lambda)$ 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.731m$ beam and compares the birth distribution and steady-state one, which shows that the latter still reach its peak near the passing-trapped boundary.

Let us examine the dependence of NBCD on the plasma current. We change the plasma current by multiplying the poloidal magnetic flux $\Psi$ (2D data from EFIT) by a factor $\alpha_{p}$. Since the poloidal magnetic field $\mathbf{B}_{p}$ is related to $\Psi$ via

where $\hat{\mathbf{R}}$ and $\hat{\mathbf{Z}}$ are the cylindrical unit vectors, the above multiplication corresponds to multiplying the $\mathbf{B}_{p}$ by $\alpha_{p}$. Similarly, since the plasma toroidal current density $J_{\phi}$ is related to $\Psi$ via

where $\mu_{0}$ is the vacuum magnetic permeability, the above scaling also scales the plasma current $I_{p}$ by $\alpha_{p}$. We keep the values of the toroidal magnetic field $B_{\phi}$ fixed when changing $\Psi$. As a result, the safety factor $q$ is also scaled by a factor of $\alpha_{p}$. Figure 10 plots the safety factor profiles corresponding to different values of $\alpha_{p}=I_{p}/I_{p,\operatorname{ref}}$, where $I_{p,\operatorname{ref}}$ is the plasma current in the original configuration.

Figure 11a 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 fieldwesson2004 (25). 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. 11b, 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 trapped fraction happens for both the birth distribution and the steady-state distribution, as is shown in Fig. 11(e) and (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.731m$. 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. And 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. Fig. 15a plots the dependence of the volume integrated fast ion current on the density, which indicates that the current decreases with increasing of the density. Fig. 15b shows that the shine-through loss decrease with the density increasing (as is expected). Also Fig. 15c 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, Fig. 15(f-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 less 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 plasma density increasing.

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 Fig. 17 and 18.

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_{\operatorname{net}}=J_{f}F$, the effective trapped electron fraction $f_{t}$,  the electron collision frequency $\nu_{e}$, the thermal electron bounce frequency $\omega_{b}$, and the normalized electron collision frequency $\nu_{e\star}$. (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 regimeHonda_2012 (28, 29), 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 $\nu_{e\star}=0$ and the resulting net current. The results show that the $\nu_{e\star}=0$ approximation slightly overestimates the net current.

The authors have no conflicts to disclose.