Nonlinear simulations of GAEs in NSTX-U

Alfvén eigenmodes (AEs) in the sub-cyclotron frequency range have been frequently observed during neutral beam injection (NBI) in the National Spherical Torus Experiment (NSTX/NSTX-U)^1–4 as well as other fusion devices such as MAST, DIII-D, JT-60U, and ASDEX Upgrade tokamak.^5–10 These modes have either predominantly shear Alfvén or compressional polarization and are identified as global Alfvén eigenmodes^11,12 (GAEs) or compressional Alfvén eigenmodes^13,14 (CAEs), respectively. Strong activity of AEs in this frequency range is most often observed for the beam injection velocities significantly exceeding the Alfvén velocity when both co- and counterpropagating AEs have been reported.^1,4,15 Counterpropagating sub-cyclotron frequency modes are driven unstable by the Doppler-shifted cyclotron resonance, whereas co-propagating AEs can be excited by both the Cherenkov resonance and by the anomalous cyclotron resonance with the energetic ions.^16–18 In NSTX, a correlation between strong activity of sub-cyclotron frequency AEs at high beam power and flattening of the electron temperature profile has been reported.^19,20 Since the reduced heating of the plasma core by NBI can significantly limit plasma performance, these instabilities can have important implications for future fusion devices, especially with the super-Alfvénic NBI velocities.

Previous theoretical attempts to explain the flattening of the electron temperature profile by the high-frequency AEs focused on an enhanced electron transport due to the interaction of Alfvén eigenmodes with bulk electrons via a parallel electric field, as well as stochasticity of the electron orbits in the presence of multiple unstable GAE modes.²¹ Energy channeling mechanisms were also proposed, whereby core-localized GAEs or CAEs can channel the energy of the beam ions from the injection region to the location of continuum damping closer to the edge.^22–25 All proposed mechanisms require several unstable modes of sufficiently large amplitudes to explain the observed reduction of the efficiency of central plasma heating by the NBI. It was shown, in particular, that the calculated energy transport is very sensitive to the number²¹ and the amplitudes²⁴ of the unstable modes. Therefore, the ability to predict the number of unstable sub-cyclotron frequency AEs as well as their saturation amplitudes is important for the transport studies in the spherical tokamaks (STs) and other devices.

Earlier analytical and simulation studies^16,26 suggest that the counterpropagating GAEs are more unstable than CAEs in the NSTX and NSTX-U conditions, which is in agreement with the experimental results.² This paper focuses on the numerical study of nonlinear evolution of counterpropagating GAEs for the NSTX-U experimental parameters. Nonlinear 3D simulation studies of unstable GAEs have been performed using the HYM^24,26 code for the NSTX-U shot #204707 (t = 0.44 s), where several unstable AEs were experimentally observed.^2,26 Two sets of nonlinear simulations have been performed. In the first set, a separate simulation run was performed for each toroidal mode number, following the nonlinear evolution of a single toroidal Fourier harmonic (single-n simulation). In the second case, a fully nonlinear simulation run was performed by following the evolution of all resolved toroidal mode numbers in the simulation simultaneously, including $n=0$⁠. A comparison of the results demonstrates the importance of the fully nonlinear simulations for correct predictions of the nonlinear amplitudes of AEs when many modes are unstable.

Numerical and analytical studies of the energetic particle-driven Alfvénic instabilities in fusion devices and their nonlinear evolution is an active area of research, but with most work focusing on the low-frequency modes including the toroidal Alfven eigenmodes (TAEs), reversed shear Alfven eigenmodes (RSAEs) and energetic particle modes (EPMs).^27–29 Saturation due to resonant particle trapping and the scaling of the saturation amplitude with the linear growth rate was studied.^30,31 Excitation of multiple unstable AEs and the resonance overlap is shown to lead to stochasticity, significant fast ion transport, and the profile flattening.^29,32 Previous studies on the nonlinear simulations comparing the single-n vs multiple-n results for the low-frequency AEs have shown a significant difference. In particular, an enhanced radial transport of energetic ions in multiple-n simulations has been reported, as well as a stronger drive of the subdominant modes compared to the single-n results.³³ 

Unlike low-frequency Alfvénic instabilities, the counterpropagating GAEs have frequencies of a fraction of the ion-cyclotron frequency, $ω∼0.1−1ωci$⁠, and are driven unstable by the cyclotron resonance with the NBI ions. The source of the free energy for the GAEs and CAEs is the velocity space anisotropy of the beam ion distribution,^16,26 whereas the TAEs, for example, can be driven by the configuration space gradients. This leads to the differences in the nonlinear evolution of the unstable GAEs in single-n and multiple-n simulations compared to the low-frequency AEs.

This paper is organized as follows. Section II describes equilibrium plasma configuration and simulation parameters. Section III presents the results of several singe-n simulation runs. Fully nonlinear simulation results including all toroidal harmonics are presented in Sec. IV. Section V discusses the effects of GAEs in energetic particle transport, and Sec. VI includes the summary of the results.

Numerical simulations using the HYM code^24,26 have been performed to study the excitation and nonlinear evolution of counterpropagating GAEs. Initial conditions were taken to correspond to plasma conditions in the NSTX-U shot #204707 (⁠ $t=0.44$ s), where several counterpropagating GAEs have been observed² with toroidal mode numbers $n=−8$ to $n=−11$ (here the sign convention $n<0$ for the modes propagating counter to the plasma current and beam injection directions is used).

In this shot, the plasma was heated by 3.3 MW of two deuterium beams with tangency radii $Rtan≈$ 0.7 and 0.5 m, and the injection energies 76 and 87 keV, respectively. The beam and thermal plasma density were $nb=7.7×1018m−3$⁠, $ne=3.3×1019m−3$⁠; the toroidal field value $Bt=0.546$ T at the axis $R0=1.08$ m, and total current $Ip=0.61$ MA. The corresponding normalized parameters used in the simulations are $nb/ne=0.236$⁠, $βbeam=0.135$⁠, $βplasma=0.08$⁠, the beam injection velocity $v0=2VA$⁠, and the mean pitch parameter $λ0=0.65$⁠. Thermal and beam ion density profiles and thermal pressure and q profiles used to calculate the initial conditions for the 3D simulations are shown in Fig. 1.

The 3D nonlinear simulations have been performed using the hybrid MHD-kinetic option of the HYM code, which is described in Ref. 26 in detail. In this model, the thermal plasma is described using one-fluid resistive MHD equations, and the energetic beam ions are modeled using particle-in-cell simulations (PIC) and the delta-f scheme. Quasineutrality is assumed with $ne=ni+nb≫nb$⁠. Full orbit kinetic description is used for the beam ions to allow the cyclotron resonances. The beam ion equilibrium distribution function is taken to be a function of the integrals of motion $F0=F0(E,λ,pϕ)$⁠, where $E=v2/2$ is the normalized particle energy (zero equilibrium electric field is assumed, v is normalized to $vA$⁠), $λ=μB0/E$ is the pitch parameter, $μ$ is the magnetic moment [ $μ=miv⊥2/(2B)+μ1$⁠, where $μ1$ includes all first-order corrections in $ϵ=ρb/LB$—the ratio of the beam ion Larmor radius to the equilibrium magnetic field scale length Belova, Gorelenkov, and Cheng (2003)], and $pϕ=Rvϕ−ψ$ is the canonical toroidal angular momentum, normalized to $mivA2/ωci$⁠. In the rest of the paper, normalized units are used, where length, time, and frequency are normalized by $vA/ωci0,$ $1/ωci0$⁠, and $ωci0=qiB0/(mic)$⁠, respectively, and $B0$ is the magnetic field at the axis.

The beam ion distribution function is slowing down in energy, has a power law dependence on $pϕ$⁠, and is a Gaussian distribution for the pitch parameter centered at $λ=λ0$ with width $Δλ$⁠. The equilibrium beam ion distribution function and parameters are described in Sec. II of Ref. 26, where the same NSTX-U shot was modeled, but focusing on the linear stability. In the simulations presented in this paper, more realistic values of the plasma resistivity, $η$⁠, and viscosity, $ν$⁠, were used corresponding to Lundquist number $S=RcvA/η=105$ and Reynolds number $Re=RcvA/ν=105$⁠, where $Rc$ is the radius of the simulation region. Initial conditions were calculated using a generalized form of the Grad–Shafranov equation,³⁴ which includes, nonperturbatively, the effects of the beam ions. The thermal plasma and beam profiles were fitted to the corresponding profiles from the TRANSP code. The model used in this paper does not include collisions, sources, or sinks for the NBI distribution.

The HYM code is an initial value, fourth-order finite-difference code, which uses a second-order time-centered explicit time stepping scheme. A cylindrical grid (⁠ $Z,R,ϕ$⁠) is used for the fields advance, and a 3D Cartesian grid (Z, X, Y) is used for the particle pushing. In all simulations discussed in this paper, a cylindrical grid of $120×60×256$ and a Cartesian grid of $120×101×101$ were used. The number of simulation particles in the single-n runs was 2M, and 10M simulation particles were used in the fully nonlinear, all-n runs. Random initial perturbations were used in all cases.

To follow the nonlinear evolution of the unstable GAEs, first, a separate simulation run was performed for each toroidal mode number from $|n|=6$ to $|n|=12$⁠. Longer simulations were performed for lower n modes, because they have smaller growth rates and reach peak amplitude at a later time. For all cases considered, the unstable modes were propagating in the negative toroidal direction (opposite to the plasma current) and were of the shear Alfvén polarization. In the single-n simulations, discussed in this section, the perturbed fields were filtered out, keeping only the $±n$ Fourier harmonics. Therefore, in this model, the MHD nonlinearities were excluded, but the energetic particle evolution was fully nonlinear, including the nonlinear equation for the delta-f weight evolution equation, and the particles following orbits in the total fields (equilibrium plus perturbed).

A summary plot of the simulation results is shown in Fig. 2, where the time evolution of the square root of the perturbed magnetic energy is plotted from six simulation runs with $|n|=6−11$⁠. The $n=−12$ mode was also unstable, but with a much smaller growth rate, and it is not shown in Fig. 2. The linear growth rates and frequencies of the most unstable mode for each n are plotted in Fig. 3, showing that the $n=−10$ is the fastest growing mode with $γ/ωci=0.024$⁠, and the $n=−9$ and $n=−11$ have smaller, but comparable growth rates $γ/ωci=0.019$ and 0.016, respectively. Calculated frequencies of the most unstable modes from $n=−8$ to $−11$ are in the range $ω/ωci=0.32−0.39$⁠, which is consistent with the experimentally observed frequencies, if the plasma toroidal rotation is taken into account.^2,26

The most unstable GAE frequency (and, therefore, $k||$⁠) is determined by the normalized beam injection velocity $v0/vA$⁠, as has been demonstrated in previous publications, see, for example, Figs. 7 and 8 from Ref. 35, where the numerical results are compared with theoretical predictions and three different experiments. It was also shown²⁶ that instability favors smaller $k⊥$⁠, and therefore, small poloidal mode numbers. Consequently, the value of the normalized injection velocity determines the most unstable toroidal mode numbers for the GAEs. Linear stability for the same NSTX-U shot was discussed in Ref. 26, where the GAE stability theory has been developed. The simulation results were found to be consistent with the theoretical predictions and the experimental data.

The linear mode structure for the $n=−10$ GAE is shown in Fig. 4, where contour plots of the perturbed electric and magnetic field components are shown. Figures 4(a) and 4(b) show two components of $δE⊥$⁠. Note that the parallel component of $δE$ is negligible in the hybrid MHD-energetic particle (EP) model, where the MHD Ohm's law is used.²⁶  Figures 4(c) and 4(d) show two components of $δB⊥$⁠, while the parallel component of $δB$ is shown in Fig. 4(e). It can be seen that the peak value of $δB||$ is at least a factor of 5 smaller than the amplitude of $δB⊥$⁠. The mode rotates in the $Bpol$ direction and has $m=0,1$ as the main poloidal mode numbers.

The HYM code is an initial value code, and it can follow the evolution of the fastest growing mode with a given value of the toroidal mode number in a single-n simulation. However, in many cases more than one unstable mode can be seen in the nonlinear phase of the simulation. For example, the spectral plot from the $|n|=9$ simulation shown in Fig. 5 indicates the growth of the most unstable mode with $ω/ωci=0.34$ in the earlier phase of the simulation (⁠ $tωci ≲600$⁠). At the later time (⁠ $tωci≳800$⁠), the growth of at least two additional modes with frequencies $ω/ωci=0.365$ and 0.39 can be seen. The weaker modes might have different main poloidal mode numbers and/or different radial mode numbers compared to the dominant $n=−9$ GAE, which is mostly $m=1,2$ mode. A frequency shift due to change in the poloidal mode number, $m→m+1$⁠, corresponds to $Δω=vA/(qR0)≈0.05ωci$ for $q≈1$⁠. Weak dependence of the GAE frequency on the radial mode number is also possible due to non-ideal MHD contribution to the shear Alfvén wave (SAW) dispersion by the adiabatic contribution of the beam ions. The beam ion finite Larmor radius effects, for example, will result in a fine frequency splitting for the SAW. The presence of weaker modes with frequencies different from the most unstable one was a typical feature for most of the single-n simulation runs, except the weakest $|n|=12$ instability.

The time evolution plots in Fig. 2 show similar behavior for the two most unstable modes, $n=−10$ and $n=−9$⁠. Their amplitudes reach comparable peak values at $t$ $ωci≈600$ and 800, respectively, with subsequent slower decay, without reaching a steady state. This is consistent with particle trapping saturation mechanism for a relatively weak damping rate (when $γdrive≫γdamp$⁠) in the absence of collisions or particle sources, which can restore the original slope of the distribution.²⁹ Previous linear studies¹⁸ have shown that the main damping mechanism for the GAEs is the continuum damping, which is included in the HYM model. Weaker modes with $n=−11$ and $n=−8$ show nonlinear oscillations after initial peaks in the amplitude and appear to reach a steady state (or slow decay) after several oscillations. The different nonlinear behavior can be related to differences in continuum damping for these GAEs, or, possibly, to growth of the subdominant modes with different poloidal/radial mode numbers. For example, the apparent saturation of the weakest $n=−6$ and $n=−7$ modes at $tωci>2500$ is, in fact, related to the weak growth of the low-frequency AEs with frequencies $ω/ωci∼0.07$⁠, much lower than the GAEs range. It is interesting that the $n=−7$ and $n=−8$ modes reach peak amplitudes higher than the much more linearly unstable $n=−11$ mode.

Some works related to ion cyclotron resonance heating (ICRH) suggest that $C1$ might be a “better” conserved constant compared to $C2$⁠, the value of which might change during the pass through the localized cyclotron resonance.^36,37 A check of the conservation of both constants can also serve for verification and/or as a measure of accuracy of the nonlinear scheme in the code, therefore the conservation of $C1,2$ was studied for resonant and non-resonant particles in the single-n GAE simulations.

Previous simulation studies of GAEs^26,35 have shown that, typically, there are two groups of particles resonant with the GAE: lower energy particles, which are driving the instability due to $∂F0/∂λ>0$ condition, and high energy particles, which have a stabilizing effect on the mode due to large negative derivative of $∂F0/∂E$ near the injection energy. These two groups are separated by the gap seen in Fig. 6(a) at $E∼2$⁠, where the lower energy particles (blue-green-yellow) with $λ<λ0$ contribute to the instability drive, whereas the high-energy particles (red) with $λ∼λ0$ are stabilizing. The gap, separating these groups of resonant particles, corresponds to a region of the phase-space where the change in the distribution function is small (particles have small weights $w∼δF$⁠). It can be shown³⁵ that location of the gap can be approximately described by a condition: $E∂F0/∂E+ωci/ω∂F0/∂λ=0$⁠. The phase-space plot of the resonant particle in (⁠ $λ,pϕ$⁠) plane [Fig. 6(b)] shows that all ions driving the instability are co-passing.

Figure 7 shows time evolution of normalized energy, angular momentum, magnetic moment, and changes in $C1$ and $C2$ along the trajectory of one of the resonant particles from the simulation of the $n=−10$ GAE (with the initial parameters corresponding to one of the green dots in Fig. 6). The initial normalized particle energy is $E≈1.5$⁠, pitch $v||/v≈0.8$⁠, and $λ≈0.3$⁠. It can be seen that both constants of motion $C1,2$ are conserved with good accuracy well into the nonlinear phase of the simulation. The range of changes in the particle energy, $ΔE≈±0.25$⁠, is much larger than the changes in both $C1$ and $C2$ values. Also shown is the time evolution of the particle weight, $w=δF/P$⁠, and magnitude of the magnetic field along the trajectory and projection of the particle equilibrium orbit in the poloidal plane. It can be seen that in the linear phase of the simulation (⁠ $tωci≲500$⁠), $w(t)$ is growing exponentially, indicating that the particle is in resonance with the mode. At the same time, the particle is losing its total kinetic energy, while its perpendicular energy (⁠ $∼μ$⁠) is decreasing and the parallel energy is growing. Large amplitude oscillations in all particle parameters occur in the nonlinear phase of the simulation (⁠ $tωci≳500$⁠), when the particle is trapped in the field of the mode. The largest relative change occurs in the value of $μ$⁠, which varies between $μB0≈0$ and $μB0≈1$⁠. Fast oscillations in the value of |B| are due to the cyclotron motion, and their amplitude correlates with changes in $μ$⁠, whereas longer timescale oscillations correspond to the particle poloidal orbit frequency.

Time evolution of the same parameters, that is $E$⁠, $pϕ$⁠, and $μ$ along a non-resonant, high-energy (⁠ $E≈2.1$⁠, $v||/v≈0.5$⁠, and $λ≈0.7$ at t = 0) passing particle orbit is compared to the changes in $C1$ and $C2$ in Fig. 8. The initial parameters for this orbit correspond to one of the red dots in Fig. 6. In this case, both constants of motion are also well conserved. Somewhat larger oscillations in the value of $ΔC1$ compared to those in Fig. 7 are related to the accuracy of the conservation of magnetic moment $μ$ along the equilibrium orbit, rather than the effects of the perturbed fields. In the HYM code, the expression for $μ$ is used, which includes all first-order corrections in $ε=ρb/LB$— the ratio of the beam ion Larmor radius to the equilibrium magnetic field scale length.³⁴ This accuracy is sufficient for the lower energy particles, but for the particles with the energies close to the injection energy, the first-order corrections to $μ$ are not sufficient, and the higher order corrections might be needed for a good conservation. For the same reason, the small scale oscillations can also be seen in the time plot of $μ$ for this orbit. In any case, the plot of $ΔC1(t)$ does not show any secular changes, in contrast to $E(t)$ and $μ(t)$ behavior. This particle is getting accelerated to an energy above the injection energy due to its nonlinear interaction with the GAE. This can be evidenced in the $E(t)$ plot, but also by the large particle weight in the nonlinear phase of the simulation. As the particle is moving into the initially vacant region of the phase space, where $F0=0$⁠, its weight becomes comparable to the “equilibrium” weight $w=(F−F0)/P≈F/P=O(1)$⁠.

Strictly speaking, the expressions for $C1$ and $C2$ should include the perturbed fields. Specifically, the exact expression for $C2$⁠, for example, includes the perturbed energy and toroidal angular momentum: $C2=E+δφ−(ω/n)(Rvϕ−ψ+RδAϕ)$⁠, where $δφ$ is the perturbed electrostatic potential, $ψ$ is the equilibrium poloidal flux, and $δAϕ$ is the toroidal component of the perturbed vector potential. The exact conservation of the above expression can be verified analytically by calculating $dC2/dt$⁠. A similar expression can be derived for $C1$ using the gyrokinetic formalism for corrections for $μ$⁠.³⁸ Simulation results, however, demonstrate that both constants of motion are sufficiently well conserved even when the equilibrium expressions are used, without including the perturbed fields in the definitions of $C1,2$⁠. This is especially true for the resonant particles that experience secular changes in their parameters due to strong interaction with the mode. However, even for the particles very far from the resonance, which experience only small oscillations in their energy, the amplitudes of the oscillations of the equilibrium expressions for $C1,2$ are several times smaller than changes in their energy (not shown).

Conservation of two constants $C1$ and $C2$ implies linear relations between the changes in the resonant particle energy, magnetic moment, and toroidal momentum: $ΔE=(ω/ωci)ΔμB0$ (or $ΔE=ωΔμ$ in normalized units) and $ΔE=(ω/n)Δpϕ$⁠. Figure 9 shows scatter plots of resonant particles demonstrating the relations between $ΔE$⁠, $Δμ$⁠, and $Δpϕ$ taken from the $|n|=10$ simulation at $tωci=825$⁠. Here, the particle color corresponds to particle energy as in Fig. 6. For comparison with changes in $pϕ$⁠, in this and subsequent plots the change in the normalized poloidal flux from the magnetic axis to the edge is $|ψ0|=60.$

It can be seen that the linear relations between the particle parameters hold with good accuracy, and the spread in the scatter plots is very small, especially compared to the large variations in the values of $ΔE$ and $Δμ$⁠. At the saturation (⁠ $t=500/ωci$⁠), the amplitudes of the relative changes in resonant particles energy, toroidal angular momentum, and magnetic moment can be estimated as $ΔE/E∼20%$⁠, $Δpϕ/pϕ∼15%$⁠, and $Δμ/μ∼25%−100%$⁠. Change in the pitch parameter can be calculated as $Δλ=(ωci/ω−λt=0)ΔE/E$⁠, showing that for a relatively small frequency, $ωci/ω≫1$⁠, even a small amplitude mode can significantly modify the $λ$ distribution of the resonant ions, especially in the lower energy region.

Conservation of two constants implies that motion of the resonant particle in the phase space is one-dimensional. Due to interaction with the mode, the resonant particles are moving along a straight line defined by $E−ωμ=const$ and $E−(ω/n)pϕ=const$ conditions in a 3D $(E,μ,pϕ)$ space. Figure 10 shows a contour plot of the projection of the perturbed beam ion distribution function, $δF$⁠, into the $E$ and $μ$ plane. Blue/red colors correspond to negative/positive values of $δF$ in the range from −0.005 to 0.006 (when F is normalized to 1). Straight lines show condition $C1=E−ωμ=const$ for three values of $C1$⁠, and arrows indicate the directions of the particle motion along these lines. It can be seen that the beam ions are moving out of the central region, where $δF<0$ into higher and lower energy regions, where $δF>0$⁠. The higher energy particles are moving toward the upper-right corner (increasing both their total energy and their perpendicular energy), while the lower energy particles are moving into the lower $E$ and $μ$ corner. This is consistent with the previous discussion in Sec. III A that there are two groups of the resonant particles: lower energy particles that are driving the instability (and losing energy) and the particles with energies close to the injection energy, which have a stabilizing effect on the GAE (being accelerated by the mode).

These relations suggest that for super-Alfvénic injection velocities, the beam-driven GAEs can serve as a very effective mean of redistribution of the resonant ions in the velocity space. The excited GAEs transfer EP energy from the perpendicular to the parallel component (for the particles driving the instability with $ΔE<0$⁠) or vice versa (for the particles near the injection energy with $ΔE>0$⁠). A smaller fraction of the energy goes into the excitation of the mode itself, when $v0≫vA$⁠.

Overall, the main effect of the GAE on the evolution of the beam ion distribution function is the reduction of the anisotropy driving the instability. The growth of the mode ceases when the slope of the distribution, $∂F0∂λ$⁠, reduces, making the instability drive comparable with the damping. Since the instability saturates at relatively small amplitudes,²⁶ the change of the slope of the total distribution (including all, resonant and non-resonant, particles) is rather small, as can be seen in Fig. 11(b) showing results for the $n=−10$ mode at $t=0$ and $tωci=700$⁠. Conservation of two constants of motion of the wave–particle interaction can be used to limit the analysis in the phase-space focusing on the evolution of the resonant particle distribution, as was demonstrated for the low-frequency modes, for example, in Ref. 31. Figure 11(a) shows the contour plot of the projection of the energy exchange rate, defined as $ΔK=∫δjb·δE d3x=∑m(vm·δE)wm$⁠, in the $(C1,C2)$ plane, where m is the simulation particle index. It can be seen that strong wave–particle interaction happens in the localized region of the $(C1,C2)$ plane. This region is then used to limit the particle phase-space, calculating the evolution of the distribution function in this region only. Figure 11(c) shows the beam ion distribution $F(λ)$ at $tωci=700$ (solid) and $t=0$ (dashed) calculated only for the particles in the limited range of constants of motion: $C1=[1.2,2.0]$ and $C2=[4.0,5.2]$⁠. A significant reduction of the slope of the distribution $F(λ)$ can be seen in a wide range of pitch parameter $0<λ≲0.8$ in that case.

Figure 11(d) also shows the projection of the beam ion distribution in orbit-averaged parallel velocities $F(u||)$ (dashed line), where $u||=⟨v||⟩$⁠, calculated for the particles in the same range of constants of motion: $C1=[1.2,2.0]$ and $C2=[4.0,5.2]$⁠. The perturbed distribution, $δF(u||)$⁠, scaled to have the same range as $F(u||)$⁠, is shown by the solid line, and the vertical line indicates the resonant velocity. It can be seen that $F(u||)$ becomes a flat-top distribution in a wide range of parallel velocities $1≲u||/vA≲1.6$ after the instability reaches its peak amplitude at $tωci=700$⁠.

In addition to the flattening of $∂F0∂λ$ in the resonant region of the phase space, the nonlinear interaction of the NBI ions with the GAE leads to smoothing of the sharp gradient of $∂F/∂E$ near the injection energy, $E0$⁠, and accelerates the energetic particles to energies above $E0$⁠. Simulation results demonstrate that the instability can lead to generation of a high-energy tail of the beam ion distribution, analogous to the ICRF heating.^36,39,40 Figure 12(a) shows projection of the beam ion distribution $F(E)$ at $tωci=700$ (solid) and $t=0$ (dashed) calculated for all simulation particles. Small changes in $F(E)$ can be seen only above the injection energy $E0≈2$ (energy is normalized by $mivA2=45 keV$⁠). The same plot [Fig. 12(b)], but including only particles from the resonant region of phase space, i.e., particles with $C1$ and $C2$ values in the fixed range of constants of motion: $C1=[1.2,2.0]$ and $C2=[4.0,5.2]$⁠, shows a significant reduction of the slope of the distribution $F(E)$ in the tail in that case. As a result of the GAE instability, the high energy particles are accelerated to the energies above $E0$⁠, and the number of particles with $E>E0$ increases compared to the initial loading.

The self-consistent effect of the counterpropagating GAEs on the beam ions is, therefore, twofold. First, the GAE serves as a means of the isotropization of the lower energy resonant particles, with their perpendicular energy converting into the parallel energy, as described in Subsection III B. This effect is stronger for more perpendicular injections (with larger $λ0$⁠), because more perpendicular injection leads to stronger GAE instability with larger growth rates and amplitudes.³⁵ Second, the unstable GAEs channel a fraction of the free energy from lower-energy higher-pitch particles to the lower-pitch energetic particles in the tail.

It is worth noting that, according to the recently developed GAE linear stability theory,^16,26 the most unstable GAE frequency corresponds, approximately, to the condition: $λ0=1−(v||,res/v0)2$⁠, where $v||,res=(ωci−ω)/k||$ is the resonant velocity and $v0=2E0/mi$ is the injection velocity. This condition defines the GAE frequency for which most of the resonant ions have $λ<λ0$ and, hence, are located in the region of the phase space where $∂F0∂λ>0$⁠, contributing to the instability drive. The case like that can be seen in Fig. 6, for example, showing that the resonant curve for $v||/vA=1.5$ passes through the point corresponding to $λ=λ0$ at the energies close to $E=E0$⁠. Therefore, for the most unstable mode, the resonant curve passes through the region of the peak density of the high-energy particles, suggesting strong wave–particle interaction with the tail ions. There could be exceptions to this rule, particularly for a very weak instability, i.e., for weak anisotropy (small values of $λ0$ and/or large $Δλ$⁠) or sub-Alfvénic injection velocities.¹⁶ In the general case, however, the above argument suggests that the development of the energetic tail in the beam ion distribution could be a generic feature of the counterpropagating GAE instabilities.

The energetic particles do not have to satisfy the orbit-averaged resonance condition Eq. (3) to experience significant net changes in their energy. For example, for a particle with energy close to the injection energy, Fig. 13 shows jumps in all parameters as the particle passes through the localized ion cyclotron resonance. The jumps happen when the particle crosses the midplane on the low-field side (LFS) where the mode amplitude peaks and the magnitude of |B| is the lowest. Both $C1$ and $C2$ are conserved during the resonance crossing, whereas the particle energy $pϕ$ and $μ$ experience jumps at the resonance location, but are conserved between the jumps. A net energy gain in this case occurs due to the finite interaction time, as the mode growth and damping time scales are comparable to the particle poloidal orbit period.

Experimentally, a high-energy feature (HEF) or “spike-on-tail” on the NBI energetic ion spectrum has been observed in some of the NSTX discharges with strong GAE activity.⁴¹ However, a significant increase above the injection energy was not reported, instead an increased beam ion density near the injection energy was observed.⁴¹ The authors of Ref. 41 concluded that the observed spike-on-tail phenomena were related to the quasilinear interaction of the NBI ions with the GAEs. How these observations relate to the HYM simulation results for an NSTX-U shot discussed above is not clear, and a further study and nonlinear simulations for the specific NSTX shots where this effect was observed are needed.

Results of full nonlinear simulation including all numerically resolved toroidal harmonics from the $n=0$ to $|n|=32$ are shown in Fig. 14, where the time evolution of $δB$ amplitude for different unstable toroidal mode numbers is displayed with different colors. The simulation run was performed for the same beam ion and plasma parameters as that shown in Fig. 2. Comparison of the amplitude evolution with that from the “single-n” runs shows that the single-n results are close to the full nonlinear simulations only for the most unstable mode (⁠ $n=−10$⁠, red), in which case the saturation amplitudes and changes in the fast ion distribution are comparable. In contrast, the peak amplitudes of all subdominant modes in the fully nonlinear simulation are strongly reduced by a factor of 3–10 compared to single-n runs. In particular, the peak amplitude of the $n=−9$ GAE, which is the second most unstable mode, is reduced by an order of magnitude compared to the single-n simulation result.

A numerical convergence check of the fully nonlinear results was performed by increasing the number of simulation particles by a factor of 5. Simulation with larger number of particles shows very similar nonlinear evolution of the $n=−10$ GAE and a strong reduction of the amplitudes of all subdominant modes compared to the single-n runs. The peak amplitudes of the subdominant modes were somewhat lower in the better resolved run due to the reduction of numerical coupling between different Fourier harmonics in the simulation.

The reduction of the amplitudes of all subdominant modes in the simulations where all toroidal harmonics are evolved simultaneously, including all MHD nonlinearities, can be explained by a flattening of the beam ion distribution by the fastest growing $n=−10$ mode. As can be seen in Figs. 2 and 14, the linear growth phase of $n=−10$ GAE ceases at $tωci∼550−600$⁠, which happens before the weaker modes can grow to a large amplitude. The $n=−10$ instability saturation is caused by a significant change in the distribution of the resonant particles, which also reduces the drive for the weaker instabilities.

The frequencies of the three most unstable modes: $ωn=−9,−10,−11=0.34,0.33,0.32ωci$ are very closely spaced, so that the frequency separation $Δω/ωci≈0.01$ is smaller than the corresponding linear growth rates or variation of $ωci$ over the cyclotron orbit. The frequency separation is also much smaller than the change in the Alfvén frequency due to a change in n or m by $±1$⁠, suggesting that for the most unstable mode, the poloidal mode number changes with n to keep $ω$ approximately the same (⁠ $q∼1)$⁠. This is consistent with previously published theory,^16,26 which shows that for $k⊥∼k||$ the GAE growth rate peaks in a relatively narrow range of frequencies and depends weakly on $k⊥/k||$⁠.

Small frequency separation results in an overlap of the resonant regions in the phase space for several GAEs with different toroidal mode numbers. Figure 11(d) also demonstrates that the flattening of $F(v||)$ distribution by the most unstable mode occurs over a rather wide range of parallel velocities $⟨v||⟩/vA≈1.1−1.5$⁠. Therefore, the nonlinear behavior of all modes with resonant velocities in this range can be affected by the changes in the distribution function caused by the most unstable mode.

Since the amplitude of the most unstable mode with $n=−10$ remains much larger than other unstable GAEs, the constants of motion, applicable for a monochromatic case of wave–particle interaction, are still reasonably well conserved in the fully nonlinear simulations. Figures 15(a) and 15(b) show scatter plots of the resonant particles demonstrating the linear relation between $ΔE$ vs $Δμ$ and $ΔE$ vs $Δpϕ$⁠, from the all-n simulation. Here, the particle colors indicate the particle energy as in Fig. 6. Solid black lines correspond to $ΔE=ωΔμ$ and $ΔE=(ω/n)Δpϕ$ relations, respectively, for $n=−10$ and $ω=0.33$⁠. The scatter of the particles from a single line is more noticeable in this case compared to that in Fig. 9 due to the presence of the weaker modes with different toroidal mode numbers and frequencies.

The change in the fast ion distribution function in the resonant region is shown in Fig. 15(c), where the initial (⁠ $t=0$⁠, dashed) and final (⁠ $t=1160$⁠, solid) distributions are plotted for the same range of $C1$ and $C2$ as in Fig. 11(c). Comparison with Fig. 11(c) demonstrates that the changes in the particle distribution function in the resonant region are comparable to that in the single-n run for the $n=−10$ GAE.

In addition to fast particle redistribution in the velocity space by the GAEs, significant changes in $pϕ$ seen, for example, in Figs. 7, 8, and 13 for the resonant and high-energy-tail particles could imply a GAE-induced radial transport. Analysis of the simulation results for the most unstable $|n|=10$ case, however, shows that on the simulation time scales (⁠ $Tsimωci∼2000$⁠, which corresponds to $Tsim≲0.1$ ms) most of the energetic particle redistribution occurs in the velocity space, while the average changes in the radial location of the resonant particles are relatively small. This was confirmed by comparing the initial and final radial profiles of the particles within the resonant range of $C1,2$ values. However, a significant fast ion diffusion by the GAEs is still possible on timescales much longer than the simulation time. In NSTX-U, for example, the sub-cyclotron frequency modes (GAEs and CAEs) can persist for 100 ms with repetitive bursts.

Despite a weak radial redistribution of the resonant particles obtained on the simulation timescale, a strong GAE instability can lead to increased EP losses. Nonlinear simulation for the most unstable $n=−10$ GAE shows an increase in the beam ion losses compared to the prompt losses calculated separately from the $n=0$ run. These losses are mostly caused by the changes in the perpendicular energy (Larmor radius) and parallel energy of the particles close to the prompt loss boundary. Figure 16(a) shows the initial (⁠ $pϕ,λ$⁠) distribution of the beam ions, which are lost due to GAE by $tωci=1500$⁠. For a reference, the orbit boundaries are plotted for full injection energy. The lost ions have a wide spread in initial energies $E∼20−90$ keV, but a relatively narrow distribution in pitch parameter with $λ∼0.85$⁠. Comparison with Fig. 6 shows that the lost particles are located far from the resonant region of the phase space. The initial (dashed line) and final (i.e., right before the loss) energy distribution of the lost ions is shown in Fig. 16(b). At the moment of the loss, most ions have an increased total energy, and a corresponding larger increase in the perpendicular energy (decrease in the parallel energy) according to Eq. (4). As can be seen in Fig. 16(a), the majority of the lost orbits are of the trapped (potato) orbit type.⁴⁰ An increase in their Larmor radius results in the particle losses to the wall. The total number of the GAE-related losses is about 0.45% of the total beam ion inventory by the time the $n=−10$ mode reaches its peak amplitude and subsequently drops to about one-third of the peak value (Fig. 2), which is a fairly short timescale $Tsim≲0.1$ ms. The simulations, therefore, suggest that multiple GAE bursts of a sufficiently large amplitude can lead to noticeable beam ion losses.

Nonlinear simulations for a single toroidal mode number and fully nonlinear runs including all resolved toroidal harmonics have been performed to study the evolution of the beam-driven GAEs in the NSTX-U.

Two integrals of motion $C1,2$⁠, Eqs. (1) and (2), of a particle in a cyclotron resonance with a single mode have been calculated and very good conservation properties are demonstrated. Conservation of both constants is shown for the resonant particles satisfying the orbit-averaged resonance condition, and for the particles which do not satisfy condition Eq. (3), but experience jumps in their parameters when they cross the local resonance region. Limiting the numerical diagnostics to the range of values of $C1,2$⁠, which shows a strong wave–particle interaction, allows us to follow the evolution of the distribution functions in the resonant part of the phase space. This update allows us to show the changes in the distribution function at the saturation even for very small perturbation amplitudes.

The 3D nonlinear simulations (both single-n and fully nonlinear) demonstrate that GAEs can very efficiently redistribute the beam ions, reducing the anisotropy in the pitch distribution in the resonant region of the phase space, and accelerate the fast ions to energies above the injection energy. It is shown that these effects are especially strong for super-Alfvénic injection velocities. Changes of the beam ion distribution correspond to a wider spread in the pitch parameter, $λ$⁠, distribution, and reduction of the instability drive. For GAEs with $ωci/ω≫1$ (corresponding to $v0>vA$⁠), the changes in the resonant particles' parallel and perpendicular energies can be several times larger than the change in its total energy, increasing as $ωci/ω$⁠. The excited GAEs, therefore, can very efficiently transfer energy from the perpendicular to the parallel component (for the particles driving the instability) or vice versa (for the particles near the injection energy, which are taking energy from the mode). A small fraction of the energy goes into the excitation of the mode itself, when $v0≫vA$⁠. This implies that even a relatively small amplitude mode can significantly modify the beam distribution in the resonant region. Experimentally, the indirect evidence of significant redistribution of fast ions by the GAEs was found, for example, in NSTX, where the GAE bursts were followed by lower frequency TAE avalanches or EPMs.⁴² 

Comparison of the simulation results from the single-n and all-n simulations show that the mode amplitudes in the single-n runs are close to the fully nonlinear case only for the most unstable mode, in which case the saturation amplitudes and changes in the fast ion distribution are comparable. The peak amplitudes of subdominant modes in all-n simulations are smaller by a factor of 3–10 compared to single-n runs due to modification of the distribution of the resonant particles by the strongest mode. The resonance overlap is caused by the close spacing of the unstable mode frequencies with different n and relatively large growth rates of most unstable GAEs. This result has important implications for transport studies, which might rely on linear growth rates and amplitude estimates from the single-n nonlinear evolution. Simulation results presented in this paper demonstrate the importance of including all unstable toroidal mode numbers in the nonlinear simulations for accurate estimates of mode amplitudes in cases when the resonant regions for several modes overlap in the phase space.

The simulations reported here were performed with computing resources at the National Energy Research Scientific Computing Center (NERSC). This research was supported by the U.S. Department of Energy (NSTX Contract No. DE-AC02-09CH11466 and Grant No. DE-SC0021271).

The authors have no conflicts to disclose.

E. V. Belova: Conceptualization (lead); Formal analysis (equal); Investigation (equal); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). E. D. Frederikson: Conceptualization (supporting); Data curation (equal); Formal analysis (equal); Investigation (equal); Writing – review & editing (equal). N. A. Crocker: Conceptualization (supporting); Data curation (equal); Formal analysis (equal); Investigation (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.