Hybrid simulations of sub-cyclotron compressional and global Alfvén Eigenmode stability in spherical tokamaks

For a chosen set of beam parameters, our simulations find co-propagating CAEs for $n\geq 3$ with $0.28\omega_{ci}<\omega<0.70\omega_{ci}$, marked by blue circles on Fig. 1. For all toroidal mode numbers, the CAE frequencies are larger than those of any unstable GAEs. This is reasonable since CAEs have frequencies proportional to $k$ whereas the GAE frequencies scale with $\left|k_{\parallel}\right|<k$. Moreover, Fig. 2 shows that the CAEs have $\left|k_{\parallel}/k_{\perp}\right|=0.3-1.5$, with larger $\left|k_{\parallel}/k_{\perp}\right|$ corresponding to the modes with higher $n$ numbers. As will also be the case with the GAEs, it is important to note that $\left|k_{\parallel}/k_{\perp}\right|$ can range from small to order unity, violating the $k_{\parallel}\ll k_{\perp}$ assumption that is often made in previous theoretical work that had large aspect ratio tokamaks in mind. This observation motivated in part the reexamination of the instability conditions for CAEs and GAEs in Ref. Lestz et al., 2020a and Ref. Lestz et al., 2020b, which will be used to interpret these simulation results in Sec. IV.1.

This group of modes was identified as CAEs since they are high frequency sub-cyclotron modes with magnetic fluctuation dominated by the compressional component near the core, where $\delta B_{\parallel}/\delta B_{\perp}\gg 1$. Qualitatively, $\delta B_{\parallel}$ for the low to moderate $n$ modes ($n<10$) usually peaks on axis with low poloidal mode number $(m=0-2)$. A typical example is shown in Fig. 3, which shows an $n=4$ co-propagating CAE with beam distribution parameterized by $v_{0}/v_{A}=5.0$ and $\lambda_{0}=0.7$ – the parameters most closely matching the conditions of the NSTX discharge which this set of simulations are modeling. The core-localization of these modes agrees with previous nonlinear HYM simulations Belova et al. (2017), contrasting with the analytic studies of CAEs under large aspect ratio assumptions, which predict localization near the edge Smith et al. (2003).

These modes also exhibit a substantial $\delta B_{\parallel}$ fluctuation on the low field side beyond the last closed flux surface, which is also a generic feature of counter-propagating GAE (see Fig. 5). This similarity has complicated previous efforts to delineate between high frequency AEs in experiment Crocker et al. (2013). While $\delta B_{\perp}$ is very small in the core for linear simulations dominated by a single CAE, there can still be large $\delta B_{\perp}$ closer to the edge. This structure is most prominent on the high field side, as shown in Fig. 3, but is also visible to a lesser degree on the low field side. The feature has been previously identified as a kinetic Alfvén wave (KAW) Belova et al. (2015, 2017), located at the Alfvén resonance location $\omega=k_{\parallel}v_{A}(R,Z)$ where CAEs undergo mode conversion. The KAW appears whenever a CAE is unstable in these simulations.

The modes with higher $n$ are localized closer to the edge on the low field side and have somewhat higher poloidal mode number ($m=2-4$). The full range of poloidal mode structures of CAEs observed in linear HYM simulations is shown in Fig. 4. The first two modes, with $n=4$ and $n=6$, are more concentrated in the core, whereas the last three modes ($n=12$, 10, and 12, respectively) are more localized near the edge.

The labeling of these modes with poloidal mode numbers is somewhat arbitrary, as the modes do not lend themselves well to description with a single poloidal harmonic $m$ along the $\theta$ direction, as the structure can peak on axis or have structures that are poorly aligned with $\psi$ contours. This dilemma is also discussed in Ref. Smith and Verwichte, 2009, where the spectral code CAE3B is used to solve for CAEs at low aspect ratio within the Hall MHD model. Moreover, the CAE has a standing wave structure, indicating the presence of multiple signs of $m$, which in turn broadens the spectrum of $k_{\parallel}$ and $k_{\perp}$ values for each mode.

Regardless of the poloidal mode numbers ascribed to them, the CAE mode structures from the HYM simulations presented here qualitatively match the CAE eigenmodes found by the CAE3B eigensolver for a separate NSTX discharge 130335 Smith and Fredrickson (2017). The similarity between the HYM and CAE3B mode structures provides further evidence that the instabilities seen in HYM and heuristically identified as CAE due to large $\delta B_{\parallel}$ in the core are indeed CAE solutions. Comparison against the CAE dispersion relation further supports the classification of these modes. In a uniform plasma, the magnetosonic dispersion including finite $\beta$ may be written

Here, $u\equiv v_{S}/v_{A}=\sqrt{\gamma_{a}\beta/2}$ where $\gamma_{a}=5/3$ is the adiabatic index and $\beta=2\mu_{0}P/B^{2}$. The finite $\beta$ corrections are important because the large fast ion pressure can make $v_{S}\sim v_{A}$ in the vicinity of the mode. In lieu of well-defined poloidal mode numbers, the mode structures from simulations are Fourier transformed to determined $k_{R}$, $k_{Z}$, and $k_{\phi}$. A quantitative agreement within $10\%$ is found between Eq. 7 and the mode frequencies in the simulations.

The co-propagating CAEs are driven by the Landau resonance with fast ions: $\omega-n\left\langle\omega_{\phi}\right\rangle-p\left\langle\omega_{\theta}\right\rangle=0$. This condition has been verified in previous HYM simulations Belova et al. (2017) by examining the fast ions with largest weights, which reveal the resonant particles due to Eq. 3c. Such examinations generally show that only a small number of resonances (i.e. integers $p$) contribute nontrivially to each unstable mode – a primary resonance $p$ and often two sub-dominant sidebands with $p\pm 1$.

The global Alfvén eigenmodes in these simulations appeared in two flavors: co- and counter-propagating relative to the direction of the plasma current. The counter-propagating modes were excited for $\left|n\right|=4-10$ in the frequency range $0.05\omega_{ci}<\omega<0.35\omega_{ci}$, and are indicated as red squares on Fig. 1. They exhibit a very wide range of wave vector directions, with unstable modes having $\left|k_{\parallel}/k_{\perp}\right|=0.2-3$ in the simulations.

The cntr-GAEs were distinguished from CAEs due to their dominant shear polarization, e.g. $\delta B_{\perp}\gg\delta B_{\parallel}$ in the core where the fluctuation was largest, and their dispersion which generally scaled with the shear Alfvén dispersion $\omega\propto k_{\parallel}v_{A}$. Counter-GAEs were routinely observed in NSTX Crocker et al. (2013) and NSTX-U Fredrickson et al. (2018) experiments with these basic characteristics, and previous comparisons between HYM simulations and experimental measurements have revealed close agreement between the frequency of the most unstable counter-GAE for each $n$ Belova et al. (2019). All counter-propagating modes which appeared in these simulations were determined to be GAEs, even though cntr-CAEs have also been observed in NSTX experiments Crocker et al. (2013). This may be attributed to the initial value nature of these simulations, which are dominated by the most unstable mode. As discussed in Ref. Lestz et al., 2020a, cntr-CAEs and cntr-GAEs are unstable for the same fast ion parameters, but the growth rate for the cntr-GAEs is typically larger.

Similar to the CAEs, the counter-GAEs did not feature poloidal structure with well-defined mode numbers. Using an effective mode number $m$ loosely corresponding to the number of full wavelengths in the azimuthal direction, the GAEs ranged from $m=0-5$, with lower $\left|n\right|$ modes typically having larger $m=3-5$ and modes with $\left|n\right|>7$ preferring smaller poloidal mode numbers of $m=0-2$. Some representative examples are shown in Fig. 6. The counter-GAE mode structure is typically more complex than that of the co-CAEs, likely due to the close proximity of the GAEs to the Alfvén continuum, which introduces shorter scale fluctuations on a kinetic scale that modulates the slower varying MHD structure. Comparing the mode structures in Fig. 4 and Fig. 6, one can see that while the CAEs are trapped in a potential well on the low field side Smith et al. (2003), the GAEs can access all poloidal angles. The sheared mode structure present in Fig. 6 may be due to nonperturbative energetic particle effects, which have previously been documented in experimental observations Tobias et al. (2011) and simulations Spong (2013); Wang et al. (2015) of shear Alfvén waves in tokamaks.

The most interesting property of the GAEs in these simulations is how significantly the beam distribution influences the frequency of the most unstable mode. This discovery was thoroughly investigated in a recent publication Lestz, Belova, and Gorelenkov (2018), and will be briefly summarized here. It was found that almost uniformly for each toroidal mode number, the frequency of the most unstable GAE scales linearly with the normalized injection velocity of the fast ion distribution. This effect is shown in Fig. 7, where the frequency of the most unstable mode can change by hundreds of kHz as $v_{0}/v_{A}$ is varied in separate simulations. In most cases, the mode structure of the most unstable mode remains relatively stable despite these large changes in frequency – yielding small quantitative changes in the mode’s width, elongation, or radial location, but not changing mode numbers. This behavior is in contrast to what is seen for the CAEs, where the frequency of the most unstable mode is nearly unchanged for large intervals of $v_{0}/v_{A}$, then switching to a new frequency at some critical value, which also coincided with the appearance of a new poloidal harmonic. In this sense, the GAEs behave more like energetic particle modes (EPMs), whose frequency fundamentally depends on fast ion properties, than conventional, perturbatively-driven MHD eigenmodes. Due to their sub-cyclotron frequency and $n<0$, the cntr-GAEs can only interact with fast ions through the ordinary Doppler-shifted cyclotron resonance. It may be written as $\omega-n\left\langle\omega_{\phi}\right\rangle-p\left\langle\omega_{\theta}\right\rangle=\left\langle\omega_{ci}\right\rangle$.

In addition to the counter-GAEs, high frequency co-propagating GAEs were also found to be unstable in these simulations with $0.15\omega_{ci}<\omega<0.60$ across $n=8-12$. Almost uniformly these modes have $m=0$ or 1, similar to the low $m$ cntr-GAEs shown in Fig. 5 and Fig. 6. Due to the large $n$ values and small $m$, co-GAEs tend to have large $\left|k_{\parallel}/k_{\perp}\right|>1$ in simulations. The co-GAEs may simultaneously resonate with the high energy fast ions through the “anomalous” Doppler-shifted cyclotron resonance, $\omega-n\left\langle\omega_{\phi}\right\rangle-p\left\langle\omega_{\theta}\right\rangle=-\left\langle\omega_{ci}\right\rangle$, as well as with fast ions with $v_{\parallel}\sim v_{A}$ through the Landau resonance $\omega-n\left\langle\omega_{\phi}\right\rangle-p\left\langle\omega_{\theta}\right\rangle=0$, as shown on Fig. 8. Although there may be comparable numbers of particles participating in these two resonances in the simulations, the energy exchange is dominated by the high energy fast ions (interacting via the $\ell=-1$ resonance). The inclusion of realistic two-fluid effects may increase the efficiency of the Landau resonance for GAEs in experiments relative to what is present in the single fluid hybrid simulations Lestz et al. (2020b). Just as with the cntr-GAEs, the high frequency co-GAEs behave more like EPMs than MHD eigenmodes, exhibiting large changes in frequency in proportion to changes in $v_{0}/v_{A}$, without any significant changes to the mode structure (also shown on Fig. 7).

Whereas cntr-GAEs are frequently observed experimentally and have been the subject of recent theoretical studies, high frequency co-GAEs are less commonly discussed. The early development of GAE theory did involve co-propagating modes, but these were restricted to very low frequencies and consequently only considered the $\ell=0$ Landau resonance Ross, Chen, and Mahajan (1982); Appert et al. (1982); de Chambrier et al. (1982); Fu and Dam (1989); Dam, Fu, and Cheng (1990). This conventional type of low frequency co-GAE driven primarily by the Landau resonance was not found to be the most unstable mode for any set of fast ion parameters in these simulations. As explored in Ref. Lestz et al., 2020b, this may be due to the fact that they have a similar instability condition to the co-CAEs, but with growth rate reduced by a typically small factor $(\omega/\omega_{ci})^{2}$. High frequency co-GAEs were never observed in NSTX, nor have they been documented in other devices. As detailed in Sec. IV.2, this is at least partly because their resonance condition favors very tangentially injected, super-Alfvénic beams which was not possible with the beam lines available on NSTX. However, the additional neutral beam installed on NSTX-U allows for much more tangential injection Gerhardt, Andre, and Menard (2012), which should be able to excite high frequency co-GAEs for discharges with sufficiently large $v_{0}/v_{A}$ (experimentally achievable with a low magnetic field).

This section presents the linear stability trends from simulations, which will be subsequently interpreted with analytic theory. Unstable CAE/GAE modes are found for a variety of beam parameters, and all simulated toroidal mode numbers $|n|=3-12$. In general, co-GAEs have the largest growth rate, followed by cntr-GAEs, and then co-CAEs. For realistic NSTX beam geometry $(\lambda_{0}=0.5-0.7)$, cntr-GAEs usually have the largest growth rate.

The summary of results from a large set of simulations is shown in Fig. 9. Each subplot corresponds to simulations restricted to a single toroidal harmonic $\left|n\right|$. Within each plot, each circle represents a separate simulation with the fast ion distribution in Eq. 4 parameterized by injection geometry $\lambda_{0}\approx v_{\perp}^{2}/v^{2}$ and injection velocity $v_{0}/v_{A}$. The white circles indicate simulations where all modes were stable (a small random initial perturbation decays indefinitely), while the color scale indicates the growth rate of the most unstable mode. The most unstable mode in each simulation is identified as a co-CAE, cntr-GAE, or co-GAE based on the descriptions given in Sec. III.

The normalized growth rates span three orders of magnitude in the simulations, ranging from $\gamma/\omega_{ci}=10^{-4}$ to $10^{-1}$. When instead normalized to the mode frequency, growth rates for CAEs are in the range $\gamma/\omega=0.005-0.05$, while most of the GAEs have $\gamma/\omega=0.01-0.2$, with a few more unstable cases having $\gamma/\omega$ up to 0.4. While there is reason to believe that GAE growth rates are typically larger than those of CAEs, the nearly order of magnitude difference seen in these simulations is primarily due to the different resonant interactions, as will be explained shortly. The colored boxes on Fig. 9 indicate the most unstable type of mode in each simulation: co-CAE, cntr-GAE, or co-GAE, as determined based on the properties discussed in Sec. III.

Across all types of modes, the most unstable modes are found for $\left|n\right|=6-9$. The growth rate dependence on toroidal mode number is shown inFig. 10. At low $\left|n\right|=3-4$, co-CAEs were the most common mode. For each value of $\left|n\right|=3-12$, there exists at least one set of beam parameters $(\lambda_{0},v_{0}/v_{A})$ such that a co-CAE was the most unstable mode, with the growth rate maximized for $n=4$. In contrast, the GAE growth rates peaked at moderately large values of $\left|n\right|=6$ and $\left|n\right|=9$ for cntr-GAEs and co-GAEs, respectively. Note that the co-GAEs are excited only at large $n$ since the $\ell=-1$ resonance requires a large Doppler shift $k_{\parallel}v_{\parallel,\text{res}}\approx nv_{\parallel,\text{res}}/R$. In simulations restricted to $\left|n\right|=1$ and $\left|n\right|=2$, the only unstable modes had $\omega/\omega_{ci}<0.05$, with high $m$ numbers and mixed polarization. Hence these modes are qualitatively different from the CAEs/GAEs being studied, and their further investigation is left to future work.

The simulation results can be interpreted with linear analytic theory. In Ref. Lestz et al., 2020a and Ref. Lestz et al., 2020b, a local expression for the fast ion drive due to an anisotropic neutral beam-like distribution is derived. Notably, the effect of finite injection energy of NBI distributions is included consistently, yielding a previously overlooked instability regime Belova et al. (2019). Terms to all orders in $\omega/\omega_{ci}$, $\left|k_{\parallel}/k_{\perp}\right|$, and $k_{\perp}\rho_{\perp b}$ are kept for applicability to the entire possible spectrum of modes. Analysis was restricted to 2D velocity space, neglecting the influence of plasma non-uniformities. In particular, this excludes drive/damping due to gradients in $p_{\phi}$, which is addressed in Sec. VI. Moreover, the calculation does not include bulk damping sources, hence it is most reliable when applied far from marginal stability. To supplement this analysis, quantification of the magnitude of damping present in HYM simulations as well as an estimate of the thermal electron damping rate (absent in the simulation model) will be presented in Sec. VII. In a nutshell, the drive (or damping) due to fast ions is a weighted integral over gradients of the fast ion distribution:

Here $\ell$ is the cyclotron resonance coefficient appearing in Eq. 5 and $h(\lambda)$ is a complicated positive function that weights the integrand. The full expression for the fast ion drive determined by Eq. 4 can be found in Eq. 21 of Ref. Lestz et al., 2020a, which can be integrated numerically for given beam and mode parameters in order to quantitatively predict the fast ion drive. However, we can first gain a qualitative understanding of the simulation results by analyzing Eq. 16. For the model distribution, the $\partial f_{0}/\partial v$ term is always negative (damping) since the slowing down function decreases monotonically. Velocity space anisotropy $(\partial f_{0}/\partial\lambda)$ can provide either drive or damping depending on its sign. For $\ell\neq 0$ resonances and the experimental value of $\Delta\lambda\approx 0.3$, the $\partial f_{0}/\partial v$ contribution is much smaller than that from $\partial f_{0}/\partial\lambda$, though they can be more comparable when $\ell=0$.

Now it is clear why the co-CAE growth rates are typically smaller than those of the GAEs. As discussed previously, the co-CAEs are driven by the $\ell=0$ resonance, while the co-GAEs and cntr-GAEs are driven by $\ell=\pm 1$, which leads to the large factor $\omega_{ci}/\omega$ multiplying the GAE growth rates. Cntr-CAEs have also been observed in experiments Crocker et al. (2013); Sharapov et al. (2014), which should also have relatively large growth rates due to this factor for the $\ell=1$ resonance, but they are not seen as the most unstable modes in HYM simulations. As discussed in Sec. V of Ref. Lestz et al., 2020a, local theory predicts the linear growth rate of cntr-CAEs to be somewhat smaller than that of cntr-GAEs driven by the same fast ion distribution, so these subdominant modes would not appear in linear initial value simulations.

One might also ask why the co-GAEs are driven by the $\ell=-1$ resonance, but the co-CAEs are not, since it would enhance their growth rate based on the argument above. Due to the difference in dispersion, fast ions must have a larger parallel velocity to resonate with CAEs than GAEs. For GAEs, the $\ell=-1$ resonance requires $v_{\parallel,\text{res}}/v_{A}\approx(1+\omega_{ci}/\omega)$, while for CAEs, $v_{\parallel,\text{res}}/v_{A}\approx\left|k/k_{\parallel}\right|(1+\omega_{ci}/\omega)$. For the ranges of $\left|k/k_{\parallel}\right|$ and $\omega_{ci}/\omega$ inferred from modes in the simulations, this resonance would require fast ion velocities far above the simulated beam injection velocity, hence the condition can not be satisfied and co-CAEs are restricted to the $\ell=0$ resonance with typically smaller growth rates.

In order to compare the stability properties of each mode against one another as a function of beam parameters, the results contained in Fig. 9 have been condensed into Fig. 11, which includes all simulated toroidal harmonics and examines the existence of modes instead of the magnitude of their growth rate. Clearly, the cntr-GAEs prefer large $\lambda_{0}$ whereas the co-GAEs prefer small $\lambda_{0}$. This is reasonable based on the form of Eq. 16. Cntr-GAEs interacting with beams ions through the $\ell=1$ resonance are driven by $\partial f_{0}/\partial\lambda>0$, while co-GAEs are driven by regions of phase space with $\partial f_{0}/\partial\lambda<0$ due to the $\ell=-1$ resonance. Hence when the distribution peaks at large $\lambda_{0}\rightarrow 1$, a larger region of phase space contributes to drive the cntr-GAEs (and damp the co-GAEs). Conversely when $\lambda_{0}\rightarrow 0$, more resonant fast ions contribute to driving the co-GAEs (and damp cntr-GAEs). A more quantitative comparison will be shown in Sec. IV.2 to elaborate on this intuition.

The co-CAE dependence on $\lambda_{0}$ is somewhat more subtle. From Eq. 16, one would expect that small $\lambda_{0}$ favors their excitation, just as in the previous argument given for co-GAEs. Yet, Fig. 11 shows unstable CAEs across a wide range $0.3\leq\lambda_{0}\leq 0.9$. Although the fast ions do provide drive for co-CAEs for $\lambda_{0}=0.1$ according to theory, the predicted growth rate is too small to overcome the background damping in simulations. As will be shown in Sec. IV.2, numerical evaluation of Eq. 16 predicts that the drive from fast ions peaks at some intermediate injection geometry $\lambda_{0}\approx 0.5$, which is similar to what occurs in simulations which also include damping on the background plasma.

With respect to the injection energy, the cntr-GAEs are unstable at significantly lower beam voltages than the CAEs. Whereas no CAE is found to be unstable for $v_{0}/v_{A}<4$, unstable cntr-GAEs were found with $v_{0}/v_{A}>2.5$. This is consistent with NSTX experiments, where cntr-GAEs were more routinely observed and in a wider array of operating parameters. This result is also qualitatively similar to dedicated MAST experiments performed at a range of magnetic field values $(v_{0}/v_{A}=1.5-2.25)$ which found a transition from cntr- to co-propagating modes as $v_{0}/v_{A}$ was increased Sharapov et al. (2014). For co-CAEs driven by the $\ell=0$ resonance, the fast ion damping from $\partial f_{0}/\partial v$ competes more closely with the drive from anisotropy $\partial f_{0}/\partial\lambda$ than when $\ell\neq 0$, leading to a larger $v_{0}/v_{A}$ for instability. As discussed in Sec. III.B.1 of Ref. Lestz et al., 2020b, net drive for co-CAEs driven by beams with $\lambda_{0}=0.7$ and $\Delta\lambda=0.3$ requires $v_{0}/v_{A}>4.1$, similar to what is found in HYM simulations. The co-GAEs also require relatively large $v_{0}/v_{A}$ for excitation, due to the requirement of a sufficiently large Doppler shift $k_{\parallel}v_{\parallel,\text{res}}$ in order to satisfy the strong $\ell=-1$ resonance which drives them. Note that $v_{0}/v_{A}\gtrsim 2.5$ and $v_{0}/v_{A}\gtrsim 4$ should not be regarded as universal conditions necessary for cntr-GAE and co-GAE excitation, respectively. CAE/GAE excitation also depends on equilibrium profiles, which determine the background damping rate as well as the spectrum of eigenmodes. For instance, cntr-GAEs were routinely excited in early operation of NSTX-U Fredrickson et al. (2018) with $v_{0}/v_{A}\approx 1-2$, while co-CAEs were rarely observed. In addition, cntr-propagating, sub-cyclotron modes have also been observed in dedicated low field DIII-D discharges Heidbrink et al. (2006); Tang et al. (2021) with $v_{0}/v_{A}\approx 1$, in agreement with corresponding HYM simulations Belova et al. (2020).

In this section, we will quantitatively examine the dependence of the growth rate on the beam injection geometry $\lambda_{0}$ and velocity $v_{0}/v_{A}$. Since the growth rate is sensitive to the specific mode properties (frequency $\omega/\omega_{ci}$ and wave vector direction $\left|k_{\parallel}/k_{\perp}\right|$) whereas only the mode with the largest growth rate survives in the simulations, it makes sense to theoretically calculate the growth rate of the most unstable mode across all values of $0<\omega/\omega_{ci}<1$ and $\left|k_{\parallel}/k_{\perp}\right|$, using the growth rate expression in Eq. 21 of Ref. Lestz et al., 2020a. While that derivation contains general two fluid effects, they will be neglected here since our aim is interpretation of the simulations which use a single fluid background plasma. The growth rate dependence is shown in Fig. 12, with separate subplots for co-CAEs, cntr-GAEs, and co-GAEs. In these figures, the background color gives the growth rate calculated analytically, while the individual points show the growth rates of unstable modes in separate simulations.

The analytic growth rate is typically larger than that from simulations, so their magnitudes have been re-scaled as indicated in the figures so that a comparison of trends can be made. It is unsurprising that the analytic growth rates are larger than those from simulations for two reasons. First, the simulations include drive/damping from fast ions and also damping of the mode on the single fluid resistive background plasma. As will be discussed in Sec. VII, this damping can be of order $20\%-60\%$ of the fast ion drive, while it is not present at all in the analytic calculation which perturbatively considers the contribution from fast ions alone. Second, the analytic theory being used is a local approximation which does not take into account equilibrium profiles of the background plasma, fast ion density, or spatial dependence of the mode structure. The value of $n_{b}/n_{e}$ that we use in the local calculation is its peak value ($5.3\%$ for the simulations in this section), leading to a calculated growth rate larger than it would be if spatial dependencies were taken into account. Notably, the local calculation does not include the contribution of $\partial f_{0}/\partial p_{\phi}$ to the growth rate, which can be either positive or negative depending on the sign of $n$, as will be discussed. For the high frequency modes studied in this work, this contribution is typically smaller than the terms which are kept in the local approximation, but it nonetheless represents a potential source of error in this simplified calculation of the growth rate. Consequently, the value of comparing this calculation with the simulation results lies in the trends with beam parameters, not the absolute values of the growth rate.

We see that for each type of mode, there is very reasonable agreement between the regions of instability predicted by analytic theory and the beam parameters of unstable modes. For co-CAEs, theory predicts the largest growth rate near moderate values of $\lambda_{0}\approx 0.5$ and large $v_{0}/v_{A}$, which are also the beam parameters preferred by unstable modes in simulations. According to theory, cntr-GAEs generally become more unstable for larger values of $\lambda_{0}$ (more perpendicularly injected beams), while co-GAEs are most unstable for small $\lambda_{0}$ (very tangential injection). This is exactly the trend seen in simulations, where the unstable cntr-GAEs generally have largest growth rate for $\lambda_{0}=0.7-0.9$ and the co-GAEs are most unstable for $\lambda_{0}=0.1$ (the smallest simulated value). Hence, numerical evaluation of the analytic expression for fast ion drive confirms many of the qualitative arguments given when interpreting Fig. 12 through the lens of Eq. 16 in Sec. IV.1. For all types of modes, larger beam velocity leads to a larger maximum growth rate, though that growth rate may correspond to different mode parameters $\omega/\omega_{ci}$ and $\left|k_{\parallel}/k_{\perp}\right|$ than the most unstable mode for some smaller beam velocity.

Lastly, note that the trends from the calculation shown in Fig. 12 rely on the aforementioned search over all mode parameters for the most unstable mode. For instance, when the same calculation is done for cntr-GAEs with fixed values of $\omega/\omega_{ci}$ and $\left|k_{\parallel}/k_{\perp}\right|$, the growth rate is no longer a strictly increasing function of $v_{0}/v_{A}$ and $\lambda_{0}$, but rather it can peak and then decrease, as in Fig. 2 of Ref. Lestz et al., 2020a. This occurs when the beam parameters are varied in such a way that not as many particles resonate with the particular mode of interest. An example of such behavior is given in the $\left|n\right|=8$ subplot of Fig. 9, as the cntr-GAE growth rate for injection geometry $\lambda_{0}=0.7$ first increases, peaks, and then decreases, with the cntr-GAE eventually being replaced by a more unstable co-CAE. Since only the $\left|n\right|=8$ toroidal harmonic is kept in that simulation, the range of unstable modes of a given type is also restricted.

Approximate stability boundaries with respect to beam injection parameters for CAEs and GAEs have recently been derived under conditions relevant to the simulated NSTX conditions Lestz et al. (2020a, b). Such boundaries apply when $0.2\lesssim\Delta\lambda\lesssim 0.8$ (the NSTX(-U) value is $\Delta\lambda\approx 0.3$) and either $\zeta\equiv k_{\perp}v_{\parallel,\text{res}}/\omega_{ci}\lesssim 2$ or $\zeta\gg 2$. Note that this constraint is related to a condition on the finite Larmor radius (FLR) of the fast ions, since $k_{\perp}\rho_{\perp b}\approx\zeta\sqrt{\lambda/(1-\lambda)}$. Moreover, the parameter $\zeta$ can be re-written using the resonance condition as $\zeta=\left|\omega-\ell\left\langle\omega_{ci}\right\rangle\right|\left|k_{\perp}/k_{\parallel}\right|/\omega_{ci}$, allowing the calculation of $\zeta$ from mode properties alone. Using the data shown in Fig. 2, one can calculate that $\zeta=0.5-1$ for the co-CAEs and $\zeta=0.5-1.5$ for the co-GAEs, so all of the simulated co-propagating modes fall within the $\zeta\lesssim 2$ regime. Meanwhile, $\zeta=0.5-3$ for cntr-GAEs, though the most unstable ones do have $\zeta<2$. Hence the approximate instability conditions derived in Sec. IV.B.1 of Ref. Lestz et al., 2020a for $\ell=\pm 1$ resonances and Sec. III.B.1 of Ref. Lestz et al., 2020b for the $\ell=0$ resonance in the small FLR regime are applicable to the unstable modes in NSTX simulations presented here. Using these approximations, and defining $\left\langle\bar{\omega}_{ci}\right\rangle\equiv\left\langle\omega_{ci}\right\rangle/\omega_{ci}$ for convenience, the following condition is found for GAE instability due to beam anisotropy

Meanwhile, for co-CAEs, the $\partial f_{0}/\partial v$ terms must also be taken into account, which leads to a more complicated condition (where we also define $x_{0}=\left\langle\bar{\omega}_{ci}\right\rangle\lambda_{0}$ and $\Delta x=\left\langle\bar{\omega}_{ci}\right\rangle\Delta\lambda$ to shorten the expression)

Note that in each of the above equations, $v_{\parallel,\text{res}}$ implicitly depends on $\omega/\omega_{ci}$ and $\left|k_{\parallel}/k_{\perp}\right|$ through the resonance condition and the appropriate dispersion relation for each mode. Hence, these conditions place constraints on the beam parameters ($\lambda_{0}$, $v_{0}/v_{A}$, as well as $\Delta\lambda$ for co-CAEs) and mode properties for a mode to be driven unstable by a neutral beam distribution. Similar to our previous analysis, Eq. 17 and Eq. 18 demonstrate that the cntr-GAEs and co-GAEs have complementary instability conditions resulting from the opposite sign of $\ell$ in the resonance driving each mode. Consequently, co-GAE excitation is favored when $\lambda_{0}$ is small, while cntr-GAEs prefer large $\lambda_{0}$.

A comparison between these conditions and the simulation results can be used to simultaneously verify theory and also understand the simulation results. Such a comparison is shown in Fig. 13, where $v_{\parallel,\text{res}}=(\omega-\ell\left\langle\omega_{ci}\right\rangle)/k_{\parallel}$ is determined from the resonance condition with $k_{\parallel}$ calculated from the mode structure in each simulation as described at the beginning of Sec. III. Each point represents an unstable mode from an individual simulation, with the beam injection velocity used in the simulation plotted against the marginal value of $v_{0}/v_{A}$ necessary for instability, as determined by the preceding equations. The shaded regions – green for co-GAEs, red for cntr-GAEs, and blue for co-CAEs – indicate the regions of instability predicted by the theoretical conditions. For the GAEs, there is quite good agreement between theory and simulations, with only a few of the co-GAEs appearing far from the predicted boundary and all of the cntr-GAEs falling in the predicted region.

The agreement for CAEs is not as strong, though the majority of unstable modes are still in the theoretically predicted region. There are three reasons why this might be the case. First, recall that the local theory which led to the predicted boundaries neglected the $p_{\phi}$ gradient contribution to the fast ion drive/damping. As will be discussed briefly in Sec. VI, this term gives a contribution proportional to $n\partial f_{0}/\partial p_{\phi}$. Hence it provides additional drive for co-CAEs and co-GAEs $(n>0)$, which would further extend the region of instability. Second, the tail of the distribution for $v>v_{0}$ was not accounted for in Eq. 19. However, these particles can provide additional drive for co-propagating modes because they provide more resonant particles with $\lambda>\lambda_{0}$ (corresponding to $\partial f_{0}/\partial\lambda<0$). Lastly, the approximate determination of $v_{\parallel,\text{res}}\approx\omega/k_{\parallel}$ may be inadequate, since it assumes that sideband drift resonances are sub-dominant. Since much of the disagreement between theory and simulation occurs for beams with $\lambda_{0}\gtrsim 0.5$, it may be that trapped particles are playing a more important role and that the sideband resonances need to be properly weighted, in addition to the primary resonance.

Altogether, this comparison suggests that although the local theory is broadly sufficient for understanding the excitation of the different types of modes, some features of nonlocal theory are needed to improve the accuracy for the co-CAEs, which are generally closer to marginal stability where additional corrections could tip the scales.

Another approach for comparison between simulation and theory is to instead fix the beam parameters and consider how the growth rate depends on the properties of each mode, namely its frequency $\omega/\omega_{ci}$ and direction of wave vector $\left|k_{\parallel}/k_{\perp}\right|$. Such analysis can be useful in the interpretation of experimental observations, since we will find that different types of modes are more unstable for different characteristic properties. This comparison is made in Fig. 14. In these figures, the background color is the analytically computed growth rate, with darker red indicating larger predicted growth rate, and blue indicating regions with negative growth rate (damped by fast ions). Gray regions indicate where no resonance is possible, occurring when $v_{0}<v_{\parallel,\text{res}}$. The gold points represent unstable modes in HYM simulations for beam distributions with fixed values of $\lambda_{0}$ and $v_{0}/v_{A}$ as specified on the plots for each type of mode. A small range of $v_{0}/v_{A}$ around a central value is included in each case in order to include enough examples from the simulations to show a trend.

For co-CAEs, a minimum value of $\left|k_{\parallel}/k_{\perp}\right|$ is needed in order to resonate with the mode at all, and an even larger value is needed in order to have net fast ion drive. This is reasonable since combination of the $\ell=0$ resonance condition and approximate CAE dispersion $\omega\approx kv_{A}$ yields $v_{\parallel,\text{res}}/v_{A}=\sqrt{1+k_{\perp}^{2}/k_{\parallel}^{2}}$. For a resonance to exist, $v_{0}>v_{\parallel,\text{res}}$ must be satisfied, which requires sufficiently large $\left|k_{\parallel}/k_{\perp}\right|$ since this corresponds to sufficiently small $v_{\parallel,\text{res}}$. Meanwhile, the approximate instability condition written in Eq. 19 is of the form $v_{0}>v_{\parallel,\text{res}}g_{c}(\lambda_{0},\Delta\lambda)$ where $0<g_{c}(\lambda_{0},\Delta\lambda)<1$ depends on the beam injection geometry and weakly on the degree of anisotropy. Hence an even larger value of $\left|k_{\parallel}/k_{\perp}\right|$ is necessary for the modes to be unstable than simply for the resonance to be satisfied. A maximum value of $\left|k_{\parallel}/k_{\perp}\right|$ for the existence of co-CAEs can be derived heuristically based on the condition that the CAE is trapped within a magnetic well on the low field side Lestz et al. (2020b), which gives $\left|k_{\parallel}/k_{\perp}\right|\lesssim n/2\pi$, shown as a dashed line on Fig. 14a. The agreement between simulations and theory in the figure is close but imperfect. While the bulk of the unstable modes from simulations cluster around the region of maximum growth rate, some of the modes are very close to or even on the wrong side of the predicted boundary. A likely explanation is the lack of $p_{\phi}$ gradient in the theory calculation.

The same comparison is very successful for both cntr- and co-GAEs. The figures illustrate a key difference between the unstable spectrum of GAEs vs CAEs. For co-CAEs, resonance requires sufficiently large value of $\left|k_{\parallel}/k_{\perp}\right|$, with an even larger value needed for net fast ion drive. In contrast, the GAEs require sufficiently large frequency to satisfy the resonance condition. Again, this can be understood as a consequence of the distinct dispersion relations for CAEs vs GAEs. For GAEs, $\omega\approx\left|k_{\parallel}\right|v_{A}$, so $v_{\parallel,\text{res}}/v_{A}\approx\left|\ell(\omega_{ci}/\omega)-1\right|$. Consequently, a resonance is possible when $v_{0}>v_{\parallel,\text{res}}$, which requires sufficiently large $\omega/\omega_{ci}$. For cntr-GAEs, the approximate condition for instability given in Eq. 17 takes the form $v_{0}<v_{\parallel,\text{res}}g_{g}(\lambda_{0})$, where $0<g_{g}(\lambda_{0})<1$. Consequently, excessively large frequencies can violate this condition (since they correspond to small $v_{\parallel,\text{res}}$), so theory predicts a band of unstable cntr-GAE frequencies. Nearly all of the simulation points fall within this band.

In the case of co-GAEs, there is no maximum frequency for unstable modes since the inequality in their approximate instability condition is flipped relative to that for the cntr-GAEs. Hence, both the resonance condition and the instability condition provide upper bounds on $v_{\parallel,\text{res}}$, or equivalently lower bounds on $\omega/\omega_{ci}$. For the beam parameters shown in Fig. 14c, the lower bound on $\omega/\omega_{ci}$ coming from the instability condition is less restrictive than that coming from satisfaction of the resonance condition, so there are no frequencies where a resonance is possible but the mode is stable. However this is not always the case. For smaller values of $v_{0}/v_{A}$ (for example $v_{0}/v_{A}=4.0$), such a band of stable frequencies exist that are damped by the fast ions. Consistency is found between the properties of unstable co-GAEs in simulations and those predicted by analytic theory, as all of the simulation points appear above the minimum frequency required for resonance, and with values of the wave vector direction $\left|k_{\parallel}/k_{\perp}\right|$ near the region of largest growth rate.

To summarize, analytic theory predicts that instability for co-CAEs occurs for modes with $\left|k_{\parallel}/k_{\perp}\right|$ above a threshold value which depends on beam parameters. A maximum value of $\left|k_{\parallel}/k_{\perp}\right|$ for CAEs can also be determined heuristically. Meanwhile, co-GAEs are unstable for frequencies that are sufficiently large. Lastly, the most unstable GAEs are predicted to occur within a specific range of frequencies. The unstable modes in HYM simulations exhibit these properties in the majority of cases, with outliers likely explained by known limitations of the local analytic theory used to calculate the stability boundaries.