Numerical simulations of global Alfvén eigenmodes excitation and stabilization in NSTX-U

Two types of subcyclotron frequency Alfvén eigenmodes (AEs) of compressional and shear polarization, namely, Compressional Alfvén eigenmodes (CAEs) and Global Alfvén eigenmodes (GAEs) have been frequently observed during neutral beam injection (NBI) in the National Spherical Torus Experiment (NSTX/NSTX-U),^1–3 as well as other beam-heated devices such as MAST^4,5 and DIII-D.^6,7 These modes have frequencies of a fraction of the ion cyclotron frequency (0.1–1ω_ci), and are driven unstable through the Doppler shifted cyclotron resonance with the super-Alfvénic NBI ions.^1,2,8 Subcyclotron frequency AEs can be excited by fast ions in other laboratory devices if the fast ion velocity is large enough to Doppler-shift the mode frequency into the ion cyclotron range. Thus, these modes can be excited in ITER due to super-Alfvénic velocities and strong anisotropy of the beam ions. They can also be excited by alpha particles near the outer edge of ITER plasma due to anisotropies in the alpha particle distribution. In NSTX, observations link these modes to flattening of electron temperature profiles and anomalously low central temperature at high beam power;^9,10 therefore, the ability to control them will have significant implications for NSTX-U, ITER, and other fusion devices where super-Alfvénic fast ions might be present.

The numerical study presented in this paper focuses on properties of counterpropagating Global Alfvén eigenmodes (GAEs) for the NSTX-U experimental parameters. It has been motivated by the exciting experimental discovery of a strong stabilizing effect that large pitch (⁠ $v | | / v ∼ 1$⁠) beam ions from the new beam sources have on these modes.^11,12 Thus, the NSTX-U results show that the use of any of the three additional neutral beams with large tangency radii provides an excellent technique to control the counterpropagating GAEs excited for $| n | =$ 8–13 and frequencies up to 0.5ω_ci. Moreover, the complete stabilization of all unstable GAEs in NSTX-U has been robustly achieved by a relatively small population of resonant high-pitch ions.¹¹ 

The injection velocities of the NBI ions in NSTX were large compared to the Alfvén velocity, $V 0 > 2 − 4 V A$⁠, the beam beta was comparable with thermal plasma beta, while a strong anisotropy in the fast-ion pitch-angle distribution could serve as the energy source; therefore, a large number of beam-driven modes could be excited. In the early stages of the NSTX experimental campaign, all modes observed in the frequency range between 0.2 and 1.2 ω_ci were identified as compressional AEs due to 1) large frequencies, and 2) agreement between the CAE frequency spectra and some observations.¹ However, initial 3D linearized simulations^8,13 using the HYM code predicted instability of counterpropagating AEs with the shear Alfvén polarization and characteristic features consistent with those of GAEs.¹⁴ In simulations, the unstable GAEs were shown to have relatively large parallel wave numbers, $k | |$ (sufficient to satisfy the Doppler shifted cyclotron resonance condition), and frequency below the minimum of the Alfvén continuum. Later, the experimental observations on NSTX also found subcyclotron frequency modes with dispersion more complicated than the CAE dispersion and characteristic features consistent with those of the GAEs.² Since then, GAEs have been routinely observed in both the NSTX and NSTX-U.^11,12,15

General properties of GAEs have been studied previously in cylindrical and toroidal geometries.^14,16–19 A GAE is a global MHD mode with frequency just below the minimum of the Alfvén continuum $ω < min ( ω A )$⁠, and owes its existence to the coupling between the shear and compressional waves. For a typical NSTX q profile, the minimum of ω_A occurs close to the magnetic axis, resulting in the mode localization in the core. Due to their localization, the main features of GAE modes observed in NSTX simulations are generally similar to cylindrical GAE modes. Early studies of energetic particle excitation of GAEs were performed for fast ion velocities $V 0 ≲ V A$⁠, and considered a resonant excitation for regular Cherenkov resonance $ω = k | | v | |$⁠.¹⁶ In this case, the GAE growth rates were found to be very small γ/ω < 0.1%, and modes strongly damped by the resonant electrons (in cylinder) and due to sideband (mostly m + 1) coupling to continuum (in toroidal geometry).¹⁶ In contrast, low-frequency toroidal Alfven eigenmodes (TAE) were found to be more dangerous, with larger growth rates, due to small $k | |$ and $γ ∼ 1 / k | | 2$ scaling.

For NSTX, the beam ion injection velocities are large, V₀ > 3 V_A, and the cyclotron resonance condition $ω − k | | v | | − ω d = l ω ci$ can be satisfied for l = 1 (normal cyclotron resonance) for counterpropagating GAEs, where ω_d is drift frequency. The cyclotron resonance instability has different dependence on $k | |$⁠, and larger $k | |$ modes can be strongly unstable, including the GAE modes. Beam-driven GAEs for different NSTX discharges have been studied using the HYM code,^8,13,20 where growth rates were found to be of the order γ/ω ∼ 0.2%–10%, with frequencies ω = 0.1 − 0.5ω_ci. The calculated range of the unstable toroidal mode numbers, frequencies, and mode polarizations have been compared with experimental observations.²¹ A recent extensive numerical study of GAEs for a wide range of beam ion parameters also revealed that strong energetic particle modifications to GAEs can occur in NSTX-like conditions.²² 

Previous analytical studies of GAE^2,23 (and CAE^2,24) excitation by energetic ions in NSTX have assumed either delta-function or very narrow distribution in pitch parameters, and failed to take into account finite injection velocity of the NBI ions, effectively assuming that $v | | , res / v 0 ≪ 1$⁠, where $v | | , res$ is the resonant velocity. Calculations of the GAE linear growth rate based on the local dispersion relation, given in the  Appendix of this paper, demonstrate why this is a poor assumption, and correct previously obtained instability condition. This also suggests different interpretation for the GAE stabilization mechanism by new beams in the NSTX-U, compared to the one suggested initially.¹² 

This paper presents results of 3D nonlinear numerical simulations that have been performed to study the excitation and stabilization of counterpropagating GAEs in the NSTX-U shot #204707 right before (t = 0.44 s) and shortly after (t = 0.47 s) the additional off-axis beam injection. Simulations using the HYM code closely match the experimental observations and confirm a robust physical stabilizing mechanism via modest off-axis beam injection.¹¹ The numerical model and results of simulations of beam driven GAEs for different toroidal mode numbers are described in Secs. II and III. The calculated range of the unstable toroidal mode numbers, frequencies, and amplitudes are compared with experimental observations^11,12 in Sec. III. Section IV focuses on GAE stabilization by the injection of high-pitch beam ions from new beam sources. Numerical results are compared with a new analytical study described in the  Appendix. The summary and conclusion are given in Sec. V.

The hybrid code HYM^21,25,26 has been used to investigate properties of beam ion driven global Alfvén eigenmodes in NSTX-U. The HYM code is a 3D nonlinear, global stability code in toroidal geometry, which treats the fast ions using δf particle simulations. The HYM code is unique in that it includes a full-orbit kinetic description of thermal and/or fast ions, allowing for the ion-cyclotron resonances. Several different physical models of thermal plasma are implemented, including kinetic description of thermal ions.²⁵ The code version used in this work utilizes the one-fluid resistive MHD model to represent the background plasma. The physical model is briefly described in this section.

The two plasma components, i.e., fast ions and thermal plasma, are coupled using a current coupling scheme.^27,28 In this scheme, the momentum equation for the thermal plasma is

where ρ, V, and p are the thermal plasma density, velocity, and pressure; n_b and J_b are the beam ion density and the beam ion induced current, $B = B 0 + δ B$ is the total magnetic field, E is the perturbed electric field, $J = c / 4 π ∇ × B$ and $δ J = c / 4 π ∇ × δ B$ are the total and perturbed plasma currents, and ν is a viscosity coefficient.

Equation (1) can be obtained by adding momentum equations for the thermal ions with density n_i and thermal electrons, neglecting the electron inertia, and using a quasineutrality condition: n_e = n_b + n_i. The rest of the fluid equations are

Here, the pressure equation includes Ohmic and viscous heating, and γ = 5/3. It is assumed that the fast ion pressure can be comparable to that of the thermal plasma, but the beam ions have a low density $n b ≪ n e$⁠. In this case, the MHD Ohm's law (2) applies.

The above model does not include the electron Landau damping for GAE modes, which is found to be negligible;²⁹ the thermal ion Landau damping is small due to the high frequency of the modes. The continuum damping is included via the dissipative terms in Eq. (1). The current-coupling scheme implemented in the HYM code possesses the physically relevant conservation properties, which enabled, in particular, the proof of the existence of global-in-time weak solutions for the resistive equations of the HYM code.²⁸ In the dissipationless limit, the current-coupling scheme is shown to be Hamiltonian,³⁰ therefore avoiding nonphysical instabilities²⁷ that could be associated with typical pressure coupling schemes.

The model and numerical scheme implemented in the HYM code conserve the total energy,²¹ which is the sum of the fluid energy (magnetic field and thermal plasma) and the beam ion energy

and F is the beam ion distribution function. A second-order, time centered, explicit scheme is used for time stepping, with smaller time steps for field equations (subcycling). A fourth-order finite difference and cylindrical geometry are used to advance fields and apply boundary conditions, while a 3D Cartesian grid is used for the particle pushing and gathering of fast ion density and current density. For a typical nonlinear run, a 144 × 60 × 256 cylindrical (Z, R, ϕ) grid is used, with a 144 × 101 × 101 (Z, X, Y) grid for particle pushing, and the particle time step dt = 0.05/ω_ci with 32 substeps for MHD equations, and a total of 2–10 M simulation particles. The total energy is conserved in fully nonlinear runs within 10% of the wave energy, provided that the numerical resolution is sufficient for the mode of interest.

The beam ions are described as kinetic particles using the particle-in-cell (PIC) simulation method and full-orbit equations of motion. The delta-f method³¹ is used to reduce numerical noise in the simulations. In this method, the equilibrium distribution function of NBI ions needs to be known analytically, and the equation for the perturbed distribution function δF = F − F₀ is integrated along the particle trajectories

where w = δF/P is the perturbed particle weight, and F₀ is the equilibrium distribution function, taken to be a function of the particle integrals of motion³²  $F 0 = F 0 ( ε , λ , p ϕ )$⁠, where ε is the particle energy, $λ = μ B 0 / ε$ is the pitch variable, $p ϕ = − ψ + R v ϕ$ is the normalized toroidal angular momentum, and μ is an adiabatic invariant μ = μ₀ + μ₁. Expression for μ includes first-order and some of the second-order corrections³² in ρ_i/L, where L is the equilibrium magnetic field scale length. In Eq. (3), P is the distribution function of simulation particles (markers) and p = F/P is the “equilibrium” particle weight. At the beginning of simulation, the simulation particles are loaded with distribution $P 0 ( ε , p ϕ )$⁠, which is also a function of the integrals of motion, so that both the distribution function of simulation particles and the physical distribution function satisfy: $dP / dt = dF / dt = 0$⁠, and $p = F 0 ( t = 0 ) / P 0 ( t = 0 ) = const$ along the trajectory.

For the simulations presented in this paper, the equilibrium distribution function is taken to be of the form:³²  $F 0 = F 1 ( v ) F 2 ( λ , v ) F 3 ( p ϕ , v )$⁠, where $v = 2 ε / m i$ is the particle velocity, and functions $F 1 , 2 , 3$ are defined by

where $F 0 → 0$ for v > v₀ or $p ϕ < p min$⁠; v₀ is the injection velocity, and we assumed $v * = v 0 / 2$⁠. The coefficient C in slowing-down distribution F₁(v) is set to 1 for v < v₀, but $C = exp [ − ( v − v 0 ) 2 / Δ v 2 ]$ for v > v₀, with $Δ v = 0.05 v 0$⁠, i.e., there is an exponential “tail” for velocities over the injection velocity, modeling the high energy tail.³³ This is done to keep both F₁ and $∂ F 1 / ∂ v$ continuous, and to avoid dealing with derivatives of a step function in the δf scheme. The pitch parameter dependence on energy (v) is chosen for good match to TRANSP results: $λ 0 ( v ) = λ 00 [ 1 − χ ( 1 − v / v 0 ) ]$ and $Δ λ ( v ) = Δ λ 0 [ 1 + d ( 1 − v / v 0 ) ]$⁠, with χ = 0.3 and d = 2. The typical parameters for the pitch-angle distribution used in this paper are Δλ₀ = 0.3 and λ₀₀ = 0.65. The function $F 3 ( p ϕ , v )$ is used to match the TRANSP profiles of the beam ion density, where α = 3 is a numerical parameter, the condition $p ϕ > p min$ describes a prompt-loss boundary, and $p max = R 0 v − ψ 0$⁠, where R₀ and ψ₀ are the major radius and poloidal flux at the axis (it is assumed that ψ₀ < 0, and ψ = 0 at the plasma boundary). Because of the difficulty in inverting the pitch parameter distribution, the simulation particles are loaded at the beginning of the simulation with $P 0 = F 1 ( v ) F 3 ( p ϕ , v )$⁠, so that the particle equilibrium weights are proportional to F₂(λ, v).

A generalized form of the Grad-Shafranov equation solver has been used, which includes, nonperturbatively, the effects of the beam ions with anisotropic distribution.³² The beam ion beta in the NSTX(-U) can be relatively large, and the beam ion current density can be comparable to that of the thermal plasma. As a result, significant modifications of the equilibrium can occur due to self-consistent inclusion of the beam ions: more peaked current profile, anisotropic total pressure shifted relative to the flux surfaces, and increase in Shafranov shift—which all can have an indirect effect on stability properties. The HYM equilibrium solver is a free boundary solver, which calculates the equilibrium in a cylindrical (Z, R) domain, but it has an option to make an iterative fit to a given plasma shape using the FREE_FIX code.³⁴ Some of the plasma profiles from the simulations described in Sec. III are shown in Fig. 1. Due to a relatively small beam injection velocity in the NSTX-U shots considered here, $v 0 ≲ 2 V A$⁠, the profile fit to the TRANSP results is very good. However in the NSTX cases, with larger normalized beam velocity, i.e., $v 0 ∼ 5 V A$⁠, the kinetic modifications to the calculated equilibrium profiles can be significant,²¹ complicating comparison with TRANSP, which does not account for these effects.

Numerical simulations using the HYM code²¹ have been performed to study the excitation and stabilization of counterpropagating GAEs in the NSTX-U shot #204707 right before (t = 0.44 s) and shortly after (t = 0.47 s) the additional off-axis beam injection [Figs. 2(a) and 2(b)]. In this section, the simulation results for plasma and beam parameters corresponding to t = 0.44 s are described and compared with the experimentally observed unstable modes.

In this shot, at t = 0.44 s, the plasma was heated by 3.3 MW of two deuterium beams with tangency radii inboard of the magnetic axis (R₀ = 1.08 m) at R_tan ≈ 0.7 and 0.5 m, and the injection energies 76 keV and 87 keV, respectively (beams labeled 1a and 1c in Ref. 11). The main plasma parameters at t = 0.44 s were $n b = 7.7 × 10 18 m − 3 , n e = 3.3 × 10 19 m − 3$ (peak values), B_t = 0.546 T at the axis R₀ = 1.08 m, and I_p = 0.61 MA. The corresponding normalized parameters as used in the simulations are $n b / n e = 0.236 , β beam = 0.135 , β plasma = 0.08$⁠, the beam injection velocity v₀ = 2 V_A, and λ₀ = 0.65.

A set of linearized simulations has been performed to investigate the stability of different toroidal Fourier harmonics with toroidal mode numbers $| n | = 6 − 13$⁠. The summary of numerical results is presented in Fig. 2(c), which shows the calculated growth rates of unstable modes, and their real frequencies normalized to the ion cyclotron frequency at the axis, f_ci = 4.15 MHz. HYM simulations reproduce the experimental finding [Fig. 2(a)], namely, before additional beam injection the simulations show unstable counter-rotating GAEs with toroidal mode numbers n = −11 to −8, as well as a weakly unstable n = −7 mode, and frequencies that match the experimentally observed unstable GAEs [Fig. 2(c)]. Numerical simulations show that all unstable modes are counter-rotating GAEs, which have shear Alfvén wave polarization in the core with small $δ B | |$ (Fig. 3). The unstable modes in simulations also have small main poloidal mode numbers with $m ≤ 3$⁠, and are localized near the minimum of Alfvén continuum, at $R ∼ 120$ cm for the n = −10 mode. The most unstable toroidal mode number, n = −10, is consistent with the experiment. Frequencies for unstable GAEs, ω/ω_ci = 0.35–0.4, are calculated in the plasma frame (no bulk plasma rotation included in the numerical model). When the bulk plasma toroidal rotation is taken into account via Doppler-shift correction to frequencies, the calculated values [blue line in Fig. 2(c)] match experimental ones (t = 0.44 s) reasonably closely. The calculated growth rates [Fig. 2(c)] for the GAEs with n = −10 and n = −11 are γ/ω_ci = 2.2% and 1.6%, respectively, which is higher than those estimated from the experimental data¹² (γ/ω_ci = 0.84% and 0.6%), probably, due to inaccuracy of the beam ion distribution function fit, and an underestimated damping by the bulk plasma MHD model.

The linear growth rates are sensitive to distribution function parameters, especially since the resonant particles are located in the “tails” of the $F ( ε , λ )$ distribution, as can be seen in Figs. 4(b) and 4(c), where the HYM fast-ion distribution function from $n =$−11 GAE simulations (t = 0.44 s) is shown on the energy vs pitch contour plots. Colored contour plots show the beam ion distribution function as calculated by TRANSP [Fig. 4(a)] and as used in HYM [Figs. 4(b) and 4(c)]. The contour plots of $F ( ε , v | | / v )$ from simulations are overlaid with scatter plots, indicating the location of resonant particles. The location of resonant particles roughly follows the resonant lines, shown as solid black lines for the Doppler-shifted cyclotron resonance corresponding to $v | | = v res = 1.4 V A$⁠, and sideband (m ± 1) resonances. The scatter plots (color dots) show resonant particles (defined as large-weight particles in δf simulations, where w/p = δF/F), with particle color corresponding to beam ion energies: from 35 keV (green) to 100 keV (red) (this color scale is used for particles in Fig. 5). Note that for a comparison with the TRANSP plot, Fig. 4(b) shows a projection of the beam ion distribution into the (⁠ $ε , v | | / v$⁠) plane, where $v | |$ is the instantaneous parallel velocity. In contrast, Fig. 4(c) shows the beam ion distribution in terms of orbit-averaged parallel velocity, showing, in particular, a population of trapped particles with $⟨ v | | ⟩ / v ≈ 0$⁠, and that all resonant ions are passing (which can also be inferred from Fig. 7). An apparent “gap” between lower and higher energy resonant particles at 80 KeV in Figs. 4(b) and 4(c) corresponds to $v | | / v ∼$ 0.6 or λ = λ₀ = 0.65, i.e., where the largest term in the equation for δf weight is zero, because $∂ F 0 / ∂ λ = 0$ (these resonant particles have nearly zero weights, and therefore are not shown).

Since scatter plots in Fig. 4 (and Fig. 5) show simulation (i.e., marker) particles, which have different equilibrium weights (p = F/P), their density on the plots might misrepresent the corresponding density of physical resonant particles. However, these plots do show the location of the resonances in the phase space. In addition, in Figs. 4 and 5 the simulation particles are labeled as “resonant” when their value of w/p = δF/F is large compared to average. Therefore, all resonant regions are shown, including those from regions of phase-space where F (and, possibly, δF) is very small.

For counterpropagating GAEs, the local resonant condition for passing particles can be written as $ω − k | | v | | − ω d − l ω ci = 0$⁠, where ω_d is the drift frequency and l = 1. Due to poloidal angle dependencies of $ω ci ∼ B 0 ( 1 + ϵ cos θ )$⁠, and drift frequency $ω d ∼ sin ( ψ + θ )$⁠, multiple resonances are possible, and the general resonant condition can be written approximately as²⁴ 

where n(m) is the toroidal (poloidal) mode number, s and r are integers related to ω_ci and ω_d dependence on θ, respectively, $q ≈ q 0$⁠, and the orbit-averaged particle parallel velocity and cyclotron frequency are denoted with bars. Multiple resonances can be seen in Fig. 5(a), which shows orbit-averaged cyclotron frequencies and parallel velocities of the resonant particles. Lower-energy ions (green in Fig. 5) satisfy the resonant condition (4) for the same value of (s + r) and their averaged parallel velocity is $v | | , res ≈ 1.4 V A$⁠, whereas resonant ions from the tail of distribution (orange-red) satisfy at least 5 sideband resonances.

Even in the case considered in this section, that is, of the unstable GAEs driven by the beam ions from the old NSTX-U beam sources, there are two groups of resonant particles, either driving the mode or taking the energy from the mode, depending mostly on their pitch parameter. As described earlier, there is a gap in the resonant particle scatterplot at energies close to 80 keV in Figs. 4(b) and 4(c), occurring where $∂ F / ∂ λ ≈ 0$⁠. This gap (approximately) separates the two particle groups. Analysis of the simulation data shows that the beam ions with initial energies of about half of the injection energy (50 keV–70 keV) and high pitch of $v | | / v ∼$ 0.8–0.9 are responsible for the instability drive, and higher energy resonant particles (over 80 keV) have a stabilizing effect. The energy exchange rate between the beam ions and the mode can be calculated as $Δ K = ∫ ( δ j b · δ E ) d 3 x = ∑ ( v m · δ E ) w m$⁠, where $δ j b$ is the perturbed beam ion current, w is the δf weight, and m is the particle index. Figure 5(b) shows time-averaged values of $( v · δ E ) w$ for the resonant particles. It can be seen that the lower energy particles are giving energy to the mode, while energetic particles from the tail of the distribution are taking energy (particle energy color scales are the same as in Fig. 4). The sum over all particles in the simulation gives ΔK_tot = −0.037 (a.u.) for a total rate of change of the beam ion kinetic energy, which is comparable to ΔK_res = −0.029, which is the change in the kinetic energy of the resonant particles only (defined as large-weight particles, so that the fraction of resonant particles is N_res/N_b = 0.037). As expected, there is little energy exchange between the mode and nonresonant particles in the simulation. Calculation of partial sums for groups of resonant particles with different initial energies gives: $Δ K res ( ε < 2 ) = − 0.104$ and $Δ K res ( ε > 2 ) = 0.075$⁠, and shows that particles from the tail have a strong stabilizing effect (here, the normalized injection energy is ε₀ ≈ 2). This has been directly confirmed in simulations by “turning off” (i.e., setting their weights to zero, w = 0) the contribution from high-energy ions with $ε ∼ ε 0$⁠, which increased the linear growth rate almost by a factor of two for the same unstable mode frequency.

Discussion in the  Appendix compares simulation results with a simple analytical expression for the linear growth rate derived in a local approximation, and compares them with previous studies. It is shown, in particular, that the approximate instability condition: $2 < k ⊥ ρ b < 4$ previously obtained for counter-GAEs² does not apply in general, and that the most unstable modes can have small values of $k ⊥ ρ b < 1$⁠, as has been found in the simulations. In addition, expression for the growth rate Eq. (A2) becomes

for the beam ion distribution $f = f ( v , λ ) = A exp ( − ( λ − λ 0 ) 2 / Δ λ 2 ) H ( v 0 − v ) / ( v 3 + v ∗ 3 )$⁠, where H(x) is a step function, $ε | | = v | | , res 2 / 2$⁠, and A is the normalization. The integrand in Eq. (5) is expressed in terms of λ using the resonant condition $v | | = v | | , res$ and $λ ≈ ε ⊥ / ε$⁠, so that 2 $ε = v 2 = v | | , res 2 / ( 1 − λ )$⁠, etc. The integration limit is $λ m = 1 − v | | , res 2 / v 0 2$⁠, where $v | | , res = ( ω ci − ω ) / | k | | |$ and v₀ is the injection velocity. For typical beam parameters, the expression in square brackets is positive, when $λ ≲ λ 0 − 2 / 3 Δ λ 2 ω / ω ci$⁠; this condition, therefore, defines regions of the phase-space, which drive the instability. As can be seen from Fig. 4(b), for a given resonant velocity, lower energy particles will have small values of λ (driving), and high energy particles have larger values of λ, which can be larger than λ₀, i.e., stabilizing. The additional stabilizing effect comes from $∂ f / ∂ ε$ in Eq. (5), because the distribution function has large gradients near the injection energy.

As shown in the  Appendix, the main driving term in the expression for the growth rate is proportional to $f ( v 0 , λ m ) ∼ exp [ − ( λ m − λ 0 ) 2 / Δ λ 2 ]$⁠. It is also shown that condition λ_m < λ₀ is a sufficient condition for the instability (all resonant particles are driving), and that the growth rate will be largest in cases when $| λ m − λ 0 | < Δ λ$⁠. Strongly unstable modes will have frequencies in the range $( 1 + v 0 / v A ) − 1 < ω / ω ci < ( 1 + v 0 / v A 1 − λ 0 ) − 1$ and small $k ⊥$⁠. This instability regime has been missed by previous studies,^2,23 which effectively assumed that λ_m = 1 (⁠ $v | | , res ≪ v 0$⁠). Higher frequency modes (for which λ_m > λ₀) can also be unstable with smaller growth rates, provided that condition: $2 < k ⊥ ρ b < 4$ is satisfied. Since the linear growth rate peaks at small $k ⊥$ [Fig. 11(b)], the most unstable GAEs will have small poloidal mode numbers, as has been seen in the 3D simulations. A detailed and more general study of the instability conditions for GAEs and CAEs will be published in a separate paper.³⁵ 

The expression in Eq. (5) has been obtained neglecting p_ϕ dependence of the distribution function. Simulation results show that, for the case considered, the contribution from ∂f/∂p_ϕ is destabilizing, that is, the growth rate was decreased by 30% and a different mode (with somewhat larger frequency ω/ω_ci ∼0.38) becomes unstable, when terms proportional to $∂ f / ∂ p ϕ$ are turned off in the equation for the particle weight Eq. (3).

Nonlinear simulations of the n = −11 to −9 GAEs have been performed separately for each toroidal harmonic, but include a full nonlinear beam ion response. In the simulations, the instability saturates due to nonlinear particle trapping, and the amplitude evolution in the nonlinear phase (Fig. 6) exhibits characteristic oscillations consistent with the particle trapping saturation mechanism. The simulations for n = −11 (Fig. 6) show the peak saturation amplitudes of $δ B ⊥ / B 0 ∼ 5 × 10 − 3$ at R ∼ 1.2 m close to the minimum of the Alfvén continuum, and δB/B₀ ∼ 10⁻³ near the edge at the midplane. For comparison, experimental estimates of the peak mode amplitudes in the same NSTX-U shot at t = 0.44 s are $δ B / B 0 ∼ 2.7 × 10 − 3$ for the n = −10 and $δ B / B 0 ∼ 1.8 × 10 − 3$ for the n = −11 modes, which is fairly close to the simulation values. Note that experimental amplitudes, which are based on reflectometer reconstructions and averaged over the time window, include a large uncertainty due to their dependence on accuracy of density profile gradient data.³⁶ 

In the simulations, the unstable n = −11 mode has been identified as a counterpropagating GAE based on a dominant perpendicular component of the perturbed magnetic field in the core (Fig. 3); however, large parallel component $| δ B | | | ∼ | δ B ⊥ |$ has been found at the plasma edge on the low-field side (LFS) for this mode. Figure 6 shows the time evolution of $δ B | |$ and two components of $δ B ⊥$ at the core, and close to the plasma edge in the equatorial plane from the same nonlinear simulation. It is seen that in the core, the magnetic perturbations have shear Alfvén wave polarization with dominant $δ B ⊥$⁠, and the mixed compressional/shear polarization near the plasma edge. Strong coupling between shear Alfvén and compressional perturbations in the NSTX simulations is related to the small aspect-ratio, relatively high beta, and kinetic effects due to energetic particles. Experimental magnetic measurements at the edge also show the mixed compressional/shear Alfvén polarization for both CAEs and GAEs in NSTX.¹⁵ 

At saturation, the GAE structure changes to a much broader one compared to the linear growth phase, possibly due to growth of other unstable modes with different values of poloidal and/or radial mode number. At the same time, the width of wave-particle resonances increases compared to that seen in Figs. 4(b) and 5. Change in the beam ion distribution function projected into the (λ, p_ϕ) phase space is shown in contour plots in Figs. 7(a) and 7(b) at two different times, tω_ci = 630 (growth phase) and tω_ci = 776 (near saturation). The changes of δf occur over a large range of the pitch parameter, λ.

Redistribution of resonant particles at saturation corresponds to reduction of their numbers near the peak of the linear instability drive, i.e., at $λ ≲ λ 0 = 0.65$⁠, and increase in f at λ > λ₀ (stabilizing), and at λ ≈ 0 (which, for fixed $v | |$⁠, corresponds to lower energies). Projection of the ion distribution in parallel velocities, $F ( v | | )$ (dashed line), and its change close to the saturation at tω_ci = 776 (solid line) are also shown in Fig. 7(c), where both functions, F and δF, are scaled to 1. It can be seen that the largest change in $F ( v | | )$ occurs at the resonant velocity $v | | , res ≈ 1.4 V A$⁠, and that the sign of δF corresponds to a reduction of the slope of $F ( v | | )$⁠. The highest order terms $∼ ω ci / ω$ in the expression for the growth rate, Eq. (5), can be rewritten as the integral of expression proportional to $− ∂ F / ∂ v | | | ε$⁠; therefore, the observed change in the beam ion distribution function corresponds to a reduction of the resonant particle drive.

In NSTX-U, an addition of new beam sources injecting fast ions nearly parallel to the magnetic field allowed effective modification of the fast ion distribution.¹¹ It has been experimentally demonstrated, in particular, that all unstable GAEs can be completely stabilized with the injection of a relatively small amount of fast ions with $v | | / v ∼ 1$⁠. Reliable suppression of the counterpropagating GAEs has been observed in most shots, and the measured GAE suppression time on the order of few milliseconds was much smaller than slowing-down time (∼50 ms), suggesting that it takes relatively few high-pitch fast ions to stabilize the GAEs.

In this section, the simulation results are described for plasma and beam parameters corresponding to t = 0.47 s in NSTX-U shot 204707, after one of the new neutral beam sources was added [Figs. 2(a) and 2(b)]. Additional 1.3 MW deuterium beam was injected with tangency radii outboard of the magnetic axis at R_tan ≈ 1.1 m, and the energy 76 keV (beam labeled 2c in Ref. 11). The corresponding normalized parameters as used in the simulations are n_b/n_e = 0.261, β_beam = 0.18, and λ₀ ≈ 0.

HYM simulations (Fig. 8 for n = 11) show that off-axis neutral beam injection strongly suppresses all unstable GAEs even for a relatively weak added beam. The off-axis beam injection has been modeled by adding beam ions with λ distribution as $F add ∼ exp ( − λ 2 / Δ λ a 2 )$⁠, i.e., with small pitch parameter values, λ₀ = 0, and smaller width $Δ λ a =$ 0.2, which roughly reproduces the TRANSP data. Figure 8 shows time evolution of n = −11 GAE magnetic energy from the linear phase to saturation for t = 0.44 s NSTX-U parameters (red line), and decay of initial perturbation for the stable case corresponding to t = 0.47 with additional off-axis beam injection (blue line).

In the experiments, new beam sources injected particles with the approximately same voltage as the old beams, increasing total beam power by ∼25%–30%, with the corresponding increase in the total population of fast ions by the same percentage. In simulations, the fraction of the off-axis beam population of the total beam ion inventory has been varied from 2% to 17%. Figure 9(b) shows the fast ion distribution function for the same beam ion parameters as in Fig. 4(b), but with added 5% of fast ions with F_add distribution, corresponding to a new beam source. Significant modification of the beam distribution function occurs only in a relatively small region of phase space, where $0.8 < v | | / v ≤ 1$⁠, but this is the region responsible for the GAE instability drive. For a case shown in Figs. 9(b) and 9(c), the modes with n = −7 to −10 are completely stabilized, and the n = −11 GAE remains unstable with a factor of 3 reduced growth rate, γ/ω_ci = 0.5%. Complete stabilization of the n = −11 GAE occurs when a fraction of fast ions from the new source increases to 7% (Fig. 10). Therefore, the simulations show that a complete stabilization of all unstable GAEs (⁠ $| n | = 7 − 11$⁠) requires an additional beam power significantly lower than what was used in this shot in the NSTX-U (i.e., 7% instead of experimental 25%), even though the HYM calculated GAE growth rates were relatively large compared to the experimental estimates.

For GAEs with n = −11 to −9, a set of simulations with varying N_add/N_tot parameter has been performed in order to find the stabilization thresholds. Figure 10 shows the growth rate of the n = −11 to n = −9 GAEs vs N_add/N_tot parameter. The unstable n = −11 GAE is stabilized when the fraction of the additional beam ions is larger than 7%. The stabilization threshold is lower for lower $| n |$ modes (Fig. 10), namely, 5% and 4% for n = −10 and n = −9, respectively. Figure 9(b) shows that the cyclotron resonance curves move to higher energies for lower $| n |$ due to the change of the corresponding resonant velocity, $v | | , res$⁠. It is also clear that the modification of total fast ion distribution $F 0 ( v | | / v )$ by a new beam source is stronger for larger energies for a fixed value of N_add/N_tot [Fig. 9(c)], which explains the lower threshold values for lower toroidal mode numbers.

An approximate expression for the GAE growth rate, Eq. (A3), given in the  Appendix, can be used to explain the experimental and simulation results on GAE stabilization. Thus, the main instability driving term includes the derivative $∂ f / ∂ λ | ε$⁠, which for the new beam with $f = F add ∼ exp ( − λ 2 / Δ λ a 2 )$ is always negative and, therefore, stabilizing. The addition of a small fraction of beam ions with high pitch changes the slope of the total beam distribution function at small λ, resulting in the net stability. In order to compare with TRANSP plots, the expression for the growth rate can also be written in terms of the pitch dependence of the distribution function, in which case it is proportional to the integral of $− ∂ f / ∂ ( v | | / v ) | ε$⁠. Negative sign of the derivative is needed for instability, and for the original beam sources this corresponds to pitch values $v | | / v ≳ 0.7$ for energies over 50 keV [Fig. 4(a)]. Addition of the beam ions from new beam sources changes the slope of the pitch distribution at $v | | / v ≲ 1$ [Fig. 9(c)], resulting in stabilization.

Thus, theoretical analysis presented here demonstrates why the injection of high pitch, deeply passing particles has a stabilizing effect on the counterpropagating GAEs; however, the details of the stabilization mechanism are different from the original publication.¹¹ The simulation and analytical results here show that both the on-axis resonant ions (driving the mode) and the off-axis resonant ions (stabilizing) have $k ⊥ ρ i < 1$⁠, and the key difference in their contribution to the mode stability is in the way they affect the sign of $∂ f / ∂ λ | ε$ at resonance. In addition, using Eq. (A2), a fraction of new beam ions required for stabilization can be calculated. Comparison of driving and damping contributions corresponding to old and new beam sources for GAE with ω/ω_ci = 0.35, $v | | , res ∼ 1.4 v A$⁠, gives an estimate for the stabilization threshold at about 10%. This is reasonably close to self-consistent simulation results, considering the approximations used in the derivation of Eq. (A2), which does not account for p_ϕ dependence of the distribution function or the continuum damping of GAEs.

The excitation and stabilization of counterpropagating GAEs by energetic beam ions have been studied using 3D nonlinear simulations right before and shortly after the additional off-axis beam injection in NSTX-U shot #204707. The equilibrium plasma shape, profiles, and beam parameters have been carefully matched to those obtained by TRANSP for this shot. Numerical results show very good agreement with the experimentally observed unstable GAEs, reproducing the range of unstable toroidal mode numbers, mode frequencies, and nonlinear amplitudes fairly closely.¹² 

Simulations also confirm that neutral beam injection from the new beam sources, depositing beam ions with large pitch, can very effectively stabilize all unstable GAEs. It has been demonstrated both numerically and analytically how small changes in the fast ion distribution function due to the addition of a relatively small amount of fast ions with $v | | / v ∼ 1$ can suppress the GAEs. Additional off-axis beam injection has been modeled by using a two-beam distribution and varying the fraction of beam ions with large pitch, $v | | / v$⁠. The complete stabilization of all unstable GAEs in simulations has been found at less than 7% of the total beam ion inventory. In NSTX-U, the fraction of beam power from new sources vs total beam power was in the range of 24%–30%, but the GAEs become suppressed at the point where the fast ion population has increased by 6% (based on the neutron rate increase by 6%) in excellent agreement with 3D simulations.^11,12

Numerical results have been compared with analytical studies,^2,23 and it has been found that the previously derived and widely cited instability condition for GAEs: $2 < k ⊥ ρ b < 4$ has limited validity. In particular, the most unstable modes in the simulations have $k ⊥ v ⊥ / ω ci ≪ 1$ for the resonant particles driving the instability. The local dispersion relation derived for two-component plasma (thermal plasma plus the energetic beam ions) also shows that the linear growth rate is largest for small values of $k ⊥ ( k ⊥ ρ b < 1 )$⁠, as long as the slope of the distribution function is positive, i.e., $∂ f / ∂ λ | ε > 0$⁠. For a specific form of beam ion distribution function: $f ( v , λ ) ∼ exp [ − ( λ − λ 0 ) 2 / Δ λ 2 ]$⁠, the sufficient condition for the instability λ_m < λ₀, where $λ m = 1 − v | | , res 2 / v 0 2$ has been obtained. The growth rate is largest in cases with narrow λ distribution, when $| λ m − λ 0 | < Δ λ$ is satisfied. It is shown that, when λ₀ is not too small, the most unstable counter-GAEs have frequencies in the range $( 1 + v 0 / v A ) − 1 < ω / ω ci < ( 1 + v 0 / v A 1 − λ 0 ) − 1$⁠, and their growth rate depends very weakly on $| k ⊥ / k | | |$ and $k ⊥ ρ b$ parameters [up to $| k ⊥ / k | | | ∼$ 3 and $k ⊥ ρ b ≲ 2$⁠, Fig. 11(b)]. For NSTX-U parameters, the predicted unstable range of frequencies, 0.33 < ω/ω_ci < 0.46, is in agreement with experimental and simulation results. The higher frequency modes with $ω / ω ci > ( 1 + v 0 / v A 1 − λ 0 ) − 1$ can also be unstable with smaller growth rates, when the approximate condition $2 < k ⊥ ρ i < 4$ is satisfied (Fig. 11).

New analytical results also show the difference in stability properties of counterpropagating GAEs and CAEs. It is shown, in particular, that for the same beam parameters, the unstable compressional modes will have smaller growth rates for finite values of $| k ⊥ / k | | |$ parameter, and the instability drive will be greatly reduced for $| k ⊥ / k | | | ≳ 2$ [Fig. 11(c)]. These findings might explain why most of the modes seen in a subcyclotron frequency range in the NSTX-U were the GAEs.

Results, presented in the  Appendix, might also explain some of the published DIII-D observations. In particular, for high-frequency modes with $ω / ω ci ∼ 0.6$ (identified as CAEs in Ref. 6), values of the $k ⊥ ρ b$ parameter (⁠ $k ⊥ ρ b ∼$ 0.8) were smaller than those previously suggested for the unstable CAEs or GAEs,² but consistent with results of this paper. In addition, the observed frequency range of most unstable modes (GAEs, or CAEs with $| k ⊥ / k | | | ≲ 1$⁠) for $v 0 / v A ∼ 1$⁠, and its scaling⁶ with B₀, n_e, and λ₀ are consistent with the presented theory.

For NSTX-U relevant beam parameters, theory predicts the peak of the linear growth rate at small $k ⊥$⁠, so that the most unstable GAEs have small poloidal mode numbers, as has been seen in the 3D simulations. Based on analytical results, a new interpretation for the GAE stabilization mechanism by new beams in the NSTX-U has been suggested. Namely, the stabilization occurs as a result of a change of slope of the beam ion distribution function in the resonant region, rather than due to small $k ⊥ ρ b$ values of ions from new beam sources, as has been suggested initially.¹² 

Analysis of the energy exchange between the resonant particles and the GAE in the nonlinear simulations shows that lower energy, high pitch particles (green dots in Fig. 4) contribute to the instability drive, while the higher energy particles with $v | | / v ∼ 0.5$ (orange-red dots) are stabilizing, and are taking energy from the mode. These findings are also consistent with the analytical study. Since lower energy (⁠ $ε ∼ 1 / 2 ε 0$⁠) particles are losing energy, while particles with energies near the injection energy are getting more energetic, the large amplitude GAE can cause significant changes in beam ion distribution, creating an energetic tail and reducing the population of midrange particles. Unstable GAEs will also cause redistribution of resonant beam ions along the resonant curve in the (λ, p_ϕ) plane, as shown in Fig. 7, which might result in the net beam ion transport in the configuration space; however, this is beyond the scope of this paper.

The simulations reported here were carried out using resources of the National Energy Research Scientific Computing Center (NERSC). This research was supported by the U.S. Department of Energy (NSTX Contract No. DEAC02-09CH11466 and Grant No. DE-SC0011810). The digital data for this paper can be found at http://arks.princeton.edu/ark:/88435/dsp018p58pg29j.

Several analytical studies related to NSTX,^2,23,24 as well as much earlier works on beam ion driven Alfvén modes of shear polarization,³⁷ predicted the instability for sufficiently large fast ion Larmor radius, ρ_b. The approximate condition: $2 < k ⊥ ρ b < 4$ had been obtained for counter-GAEs, based on the growth rate dependence on derivative of the Bessel function in the integral of the form² 

where $g = J 1 2 / ξ 2 , J 1 = J 1 ( ξ )$ is the Bessel function of first order with argument $ξ = k ⊥ v ⊥ / ω ci$⁠, l = 1 cyclotron resonance for the counterpropagating GAE has been assumed, $f b = f b ( v , λ )$ is the beam ion distribution function, and $λ = μ B 0 / ε$ is the pitch parameter. Similar expressions have been obtained for compressional (fast magnetosonic) waves with the instability condition being:²  $1 < k ⊥ ρ b < 2$⁠.

In this appendix, the stability condition for counter-GAEs is reconsidered. Thus, analysis of simulation results, for example, for n = −11 GAE, shows that for this mode $k ⊥ v ⊥ / ω ci ≲ 1$ for the resonant beam ions with energies near injection energy and larger λ. For lower energy resonant particles, with small λ and $ε ∼ 1 / 2 ε 0$⁠, condition $k ⊥ v ⊥ / ω ci ≪ 1$ is satisfied. Moreover, calculation of the energy exchange between the resonant particles and the mode shows that lower energy, high pitch particles (green dots in Fig. 4) contribute to the instability drive, while the higher energy particles with $v | | / v ∼ 0.5$ (orange-red dots) are stabilizing, and are taking energy from the mode, whereas the opposite might be expected from the expression in Eq. (A1), which is positive only for large arguments of Bessel function, i.e., large $v ⊥$⁠.

In order to understand numerical results, a dispersion relation for shear Alfvén type mode can be derived using the same physical model as in the HYM code. In local approximation and neglecting coupling to compressional modes, the dispersion relation becomes³⁸ 

where

and $Φ ⊥ = 1 v ⊥ ∂ f ∂ v ⊥ + k | | ω ( ∂ f ∂ v | | − v | | v ⊥ ∂ f ∂ v ⊥ )$⁠, $ζ l = 1 / ( ω − k | | v | | − l ω ci )$⁠. Keeping only the nonadiabatic contribution of the beam ions, for l = 1 the growth rate can be written as

where $k | | < 0$ for the counter-rotating mode, and f is the beam ion distribution function normalized so that $∫ f d ε ⊥ d v | | = 1$⁠. In the limit $k ⊥ ρ b ≪ 1 , J 1 2 / ξ 2 ≈ 1 / 4$⁠, and for $ω / ω ci ≪ 1$⁠, this expression for γ is the same as derived in Ref. 38 for shear Alfvén wave in $k ⊥ = 0$ case. Using the resonant condition and keeping only the highest order terms in ω_ci/ω in the integrand, for the beam ion distribution $f = f ( v , λ ) H ( v 0 − v )$⁠, the growth rate can be written approximately as

where $v | | , res = ( ω ci − ω ) / | k | | |$⁠, H(x) is a step function, and A is normalization, which, for $f = A exp [ − ( λ − λ 0 ) 2 / Δ λ 2 ] / ( v 3 + v ∗ 3 )$⁠, can be approximated by $A − 1 = π / 3 ln [ ( v 0 / v ∗ ) 3 + 1 ] Δ λ ( 1 − λ 0 ) − 1 / 2$ with good accuracy when $λ 0 < 1$⁠. The integrand in Eq. (A3) is expressed in terms of λ using the resonant condition and $λ ≈ ε ⊥ / ε$⁠, so that 2 $ε = v 2 = v | | , res 2 / ( 1 − λ )$⁠, etc. The integration limit is $λ m = 1 − v | | , res 2 / v 0 2$⁠, where v₀ is the injection velocity. Numerical integration of γ from the complete expression in Eq. (A2) including the energy dependencies of λ₀ and $Δ λ$ described in Sec. II gives: $γ / ω ci ≈$ 0.05 for the beam parameters used in the 3D simulations: $λ 00 = 0.65 , Δ λ 0 = 0.3 , v 0 = 2 v A , v ∗ = v 0 / 2$⁠, and for $ω / ω ci = 0.35$⁠, and $v | | , res ∼ 1.4 v A$⁠. This value is higher than the 3D simulation result: $γ / ω ci ≈$ 0.02, obtained assuming the same distribution function dependence on v and λ, however including in addition the dependence on fast ion canonical angular momentum p_ϕ,³² which was neglected in the derivations of Eq. (A2). The overestimated growth rate can also be attributed to the use of local approximation and continuum damping present in the self-consistent simulations.

Expressions (A2) and (A3) clearly show that beam ions will drive the instability provided $∂ f / ∂ λ > 0$⁠, i.e., when λ < λ₀ for the resonant particles, and will have a stabilizing effect otherwise. Stabilizing contribution from $∂ f / ∂ ε$ is weaker for most cases of interest, i.e., for $ω / ω ci ≪$ 1 and $Δ λ ≪$ 1. It is also clear from Eq. (A3) that this instability condition does not depend on the value of $k ⊥$ in the argument of the Bessel function, and, in fact, the largest growth rate will be obtained for small $k ⊥$ (such that $k ⊥ ρ b ≪$ 1), because function $J 1 2 / ξ 2$ peaks at small values of ξ (effectiveness of wave-particle interaction is reduced for large values of $k ⊥ ρ b$⁠).

For the case of slowing-down distribution $f ∼ 1 / v 3$ (setting v_* = 0), most of explicit pitch parameter dependence in expression Eq. (A3) will cancel, leaving

where $λ m = 1 − ( ω ci − ω ) 2 / ( k | | v 0 ) 2$ is the limit of the integration which depends on mode parameters and the injection velocity. Assuming small λ and integrating by parts, one can obtain the expression for the growth rate including the derivative of $J 1 2$⁠, similar to Eq. (A1). However, in addition, there will be a positive term from the integral of a complete derivative proportional to $f ( v 0 , λ m ) ∼ exp [ − ( λ m − λ 0 ) 2 / Δ λ 2 ] / v 0 3$⁠, which had been neglected in previous studies. This term can be neglected only for very narrow λ distribution and provided that $| λ m − λ 0 | ≫ Δ λ$ is satisfied, in which case it is exponentially small. However, since the most unstable modes have $| λ m − λ 0 | ≲ Δ λ$⁠, this is not a reasonable assumption. Therefore, expression (A1) is valid only in the case of delta-function distribution $δ ( λ − λ 0 )$ (when $λ m > λ 0$⁠), as derived in the original publication.³⁷ 

In order to further illustrate this point, Eq. (A2) has been integrated numerically for the beam parameters used in 3D simulations: $λ 00 = 0.65 , v 0 = 2 v A , v ∗ = v 0 / 2$ and for different values of ω and $k ⊥$⁠, assuming a simple relation between mode frequency and resonant velocity: $v | | , res = ( ω ci / ω − 1 ) v A$ (no sideband resonances). Figure 11(a) shows the contour plot of normalized growth rate and plot of $γ / ω ci$ vs $k ⊥ ρ i$ for very narrow pitch parameter distribution with Δλ = 0.03; blue-green contours correspond to negative γ values, and red-orange to positive. The peak value of $γ / ω ci ≈ 2.0$ occurs at $ω / ω ci ≈ 0.446$ for small $| k ⊥ / k | | | ≪ 1$⁠. At higher frequencies, the growth rate is strongly reduced and positive only for $| k ⊥ / k | | | ≳ 2.5$⁠. Figure 11(a) also shows normalized growth rate vs $k ⊥ ρ i$ for two frequencies ω/ω_ci = 0.446 and 0.50, where $k ⊥ ρ i$ is calculated approximately as $k ⊥ / | k | | | ( 1 − ω / ω ci ) λ 0 / ( 1 − λ 0 )$ using the relation $v ⊥ = v | | , res λ / ( 1 − λ )$⁠. It can be seen that for larger frequencies, for example, for ω/ω_ci = 0.5, the growth rate does show the oscillatory behavior related to derivative of the Bessel function as in Eq. (A1), and the approximate instability range $1.7 < k ⊥ ρ i < 3.9$ is similar to that given in Ref. 2. Sharp transition between the strongly unstable low $k ⊥$ regime and the “Bessel” regime occurs at frequency $ω / ω ci = ( 1 + v 0 / v A 1 − λ 0 ) − 1 ≈ 0.46$⁠, corresponding to the λ_m = λ₀ condition. These calculations demonstrate that even in the limit of very narrow pitch distribution, the strongest instability occurs at small $k ⊥$⁠, and it cannot be described by Eq. (A1). For narrow λ distribution, the peak instability growth rate is very large due to the large gradient of f, and because most of the beam ions are resonant. In practice, the band of strongly unstable frequencies around 0.446 will be wider than that shown in Fig. 11(a), because of the possibility of sideband resonances.

For comparison, Fig. 11(b) also shows growth rates calculated for a realistic width of pitch parameter distribution [including λ₀(v) dependence and the widening of $Δ λ ( v )$ at lower energies, as described in Sec. II], i.e., for $Δ λ 0 = 0.3$⁠. In this case, calculated γ values are smaller, with the peak value of $γ = 0.036 ω ci$ at $ω / ω ci ≈$ 0.40 and small $k ⊥ / k | |$⁠. Transition to the “Bessel-like” regime occurs at $ω / ω ci ≳$ 0.45 with smaller growth rates and $k ⊥ ρ i > 2$⁠. Note that there is no instability at $ω / ω ci ≲$ 0.33, because for lower frequencies the cyclotron resonant condition cannot be satisfied for a given $v 0 / v A = 2$⁠, unless sideband resonances are allowed. Calculations based on Eq. (A2) also show the importance of including Δλ(v) dependence for growth rate calculations (and in 3D simulation), because resonant particles can have energies as low as 1/2 of the injection energy, where Δλ is wider than at the injection energy. Larger width of pitch parameter distribution reduces the instability drive, and therefore, results in smaller growth rates.

Figure 11(b) shows that even for a realistically wide pitch parameter distribution, the most unstable modes have small $k ⊥ / | k | | |$ and a relatively narrow range of frequencies given approximately by condition: $( 1 + v 0 / v A ) − 1 < ω / ω ci < ( 1 + v 0 / v A 1 − λ 0 ) − 1$⁠. When the value of λ₀ is reduced, the upper limit of this range and the most unstable mode frequency will reduce as well, and this scaling holds as long as the value of λ₀ is not too small. In the cases of a very weak drive, that is, for smaller λ₀, modes with $k ⊥ / | k | | | < 1$ will remain unstable, but the peak of the growth rate may shift to shorter perpendicular wavelengths with $k ⊥ / | k | | | >$ 4. For v₀/v_A = 2 and $Δ λ 0 = 0.3$⁠, this transition happens at $λ 0 ≤$ 0.45. For larger injection velocities (v₀/v_A > 3) or for narrow λ-distributions (⁠ $Δ λ 0 ≲ 0.2$⁠), the peak instability still occurs at small $k ⊥ / | k | | | ≪ 1$ for a larger range of pitch distributions, i.e., for $λ 0 ≳ 0.25$⁠.

The dispersion relation for compressional modes (CAEs) can be derived in the same way as that for the GAEs. Using the same assumptions as those used for the derivation of Eq. (A2), and neglecting coupling between compressional and shear waves, the expression for the counterpropagating CAE growth rate can be shown to be identical to that of GAEs, except that the Bessel function $J 1 2 / ξ 2$ in Eq. (A2) should be replaced with $1 / 4 ( J 0 − J 2 ) 2$⁠. Therefore, in the limit $k ⊥ → 0$ the CAE has the same growth rate (and same instability condition) as the GAE. However, for finite values of $k ⊥$⁠, the instability condition and the range of unstable frequencies change significantly mostly due to a different relation between the resonant velocity and the mode frequency: $v | | , res = k / | k | | | ( ω ci / ω − 1 ) v A$ (no sideband resonances). Figure 11(c) shows the contour plot of the normalized growth rate and plot of $γ / ω ci$ vs $k ⊥ ρ i$ for the same beam parameters as for GAEs in Fig. 11(b) with $Δ λ 0 =$ 0.3. The main difference from GAEs is that both the growth rate and most unstable mode frequency of CAEs depend strongly on the $k ⊥ / | k | | |$ parameter, and the instability becomes very weak for $k ⊥ / | k | | | ≳ 2$⁠. Since the smallest value of $k ⊥$ is limited by the width of the CAE's potential “well,”²¹  $Δ R$⁠: $k ⊥ ≥ k r ∼ π / Δ R$⁠, and $| k | | | ≤ ω / v A$⁠, it can be estimated that $k ⊥ / | k | | | ≥ π v A / Δ R ω$⁠. For NSTX-U parameters, this implies that $k ⊥ / | k | | | ≳ 1$⁠, and weaker instability drive for CAEs is predicted compared to GAEs. From condition $λ m < λ 0$ follows the range of unstable CAE frequencies: $( 1 + | k | | / k | v 0 / v A ) − 1 < ω / ω ci < ( 1 + | k | | / k | v 0 / v A 1 − λ 0 ) − 1$⁠. Interestingly, since the CAE instability is mostly limited to the $k ⊥ ρ i ≲ 1$ range [Fig. 11(c)], the growth rates and unstable frequencies of counter-CAEs correspond to those of counter-GAEs with reduced effective “injection velocity” $v 0 → v 0 | k | | / k |$⁠.