Zonal structure effect on the nonlinear saturation of reverse shear Alfven eigenmodes

Alfven eigenmodes in a tokamak plasma can be driven unstable by energetic particles (EP). For near marginal instabilities, it is widely believed that the main nonlinear saturation mechanism is the trapping of resonant particles in the wave field and the local flattening of the energetic particle distribution. This is because the excursion of trapped particles, or the resonant island width in the phase space, scales as $∼A$ (see, for instance, Sec. 7.2 of White¹). Here A is the wave amplitude. If the initial particle distribution function has a locally linear profile in radius or energy, then this island width leads to a perturbed distribution that also scales as $∼A$⁠. In comparison, other nonlinear effects scale as ∼A² and become negligible for near marginal instabilities.^2–4 In a numerical study of such modes, the mode structure and frequency can then be taken from linear eigenmode analyses and assumed fixed, and only the mode amplitude or phase evolves in response to the energetic particle drive.^5–7 Numerical simulation of the nonlinear evolution problem can be greatly simplified. However, the validity of the near marginal assumption should be established with a more self-consistent model. In the kinetic-magnetohydrodynamic (MHD) simulation of n = 0, 1, 2 of Fu and Park,⁴ it is found that the wave-particle trapping is indeed the main nonlinear mechanism for $δBr/B≥10−3$⁠. However, a study by Todo et al.,⁸ also based on the kinetic-MHD model, found that the MHD nonlinear effect is negligible for $δBr/B∼10−3$ but becomes important for $δBr/B≥10−2$⁠. In the Landau-fluid simulation of Spong,⁹ the n = 0 mode plays a crucial role in the coupled evolution of three modes that display repetitive predator-prey phenomena. Zonal structures (flux-surface averaged fluctuations), or the n = 0 perturbations in general, play an important role in the theoretical studies of nonlinear toroidal Alfven waves when mode coupling is considered.^10–12 The n = 0 mode is of special interest because it can either be force generated by the unstable $n≠0$ Alfven eigenmode⁸ or spontaneously generated (see Chen and Zonca¹² and references therein). We use the term “force generated” or “force driven” to refer to a situation where the n = 0 mode is driven by a nonlinear source due to the self-coupling of an $n≠0$ mode and grows at twice the driving mode growth rate. As such, n = 0 fluctuations naturally accompany any n > 0 instabilities.

In this paper, we study the effect of the n = 0 perturbation on the saturation of the n = 4 reverse shear Alfven eigenmode (RSAE), using a gyrokinetic ion/fluid electron hybrid model. The RSAE case studied here is based on DIII-D discharge #159243 at t = 805 ms. RSAEs of mode number n = 2, 3, 4, 5, 6 are excited by injected beam particles. We focus on the n = 4 mode. We first carry out linear simulations to identify the dominant instability as RSAE, then scan over the energetic particle pressure to vary the linear growth rate, and obtain the initial saturation amplitude of the driving n = 4 mode. Our main goal is to answer the question, What is the range of linear growth rate for the validity of a reduced model that treats only the beam particles nonlinearly? The answer to this question depends on the model for the thermal plasma. In a gyrokinetic model for the thermal species, an initial zonal flow is subject to collisionless damping,¹³ which is not present in the MHD model. The hybrid model adopted here is of particular interest, because the gyrokinetic ions provide an accurate description of the collisionless ion response to the zonal fields, while the fluid electron model greatly reduces simulation time and allows us to scan the low growth rates that are impractical for models using kinetic electrons. The main conclusion is a quantitative result on the validity of reduced models that includes only nonlinear beam effect. With the n = 0 perturbations included, the n = 4 saturation amplitude scales with the growth rate as A ∼ γ² at low growth rates, $γ/ω≤3%$⁠. This is the scaling expected from the theory of a single mode nonlinearly saturated by wave-particle trapping. In this range, the n = 0 mode and the thermal plasma nonlinear effects are not important. At higher growth rates, the zonal structure effect becomes significant. Zonal structures are force generated⁸ and saturates with the driving mode. Spontaneous generation of the zonal structures has not been observed. Zonal structures cause nonlinear modification of the mode structure and enhance damping, therefore reduces the saturation amplitude.

A brief description of the hybrid model is given in Sec. II, with emphasis on the description of zonal structures. Linear results are presented in Sec. III to identify the experimentally observed modes as RSAEs. Section IV presents nonlinear results.

In this work, we use the gyrokinetic ion/fluid electron hybrid model of the GEM code.^14,15 The field perturbations are given by $E=−∇ϕ−(∂A||/∂t)b$ and $δB⊥=∇A||×b$⁠, with $b$ being the unit vector along the equilibrium magnetic field. The parallel magnetic field perturbation is neglected. The fluid electron model consists of the electron continuity equation and the isothermal condition for the electron temperature perturbation. The vorticity equation, or the time derivative of the gyrokinetic Poisson equation, is used to obtain ϕ.¹⁵ The potential $A||$ is obtained by time evolving $∂A||/∂t=−∇||ϕ−E||$⁠. The electron continuity equation is

where $u||e$ is the electron parallel flow velocity, $vE=E×b/B, Ωe=eB/me$ is the electron gyro-frequency, and $δp⊥$ and $δp||$ are the perturbed perpendicular and parallel electron pressure, respectively.

The parallel electric field $E||$ is determined by the Ohm's equation

This equation is derived by combining the time derivative of the Ampere's equation with the $v||$-moment of the electron drift-kinetic equations.¹⁶ Here, f_e1 and f_e0 are the perturbed and equilibrium electron distribution, respectively. $vG=v||b+vD+vE+v||δB⊥/B≡v||b+vG⊥$ is the guiding center velocity, and v_D is the $∇B$ and curvature drift. In the fluid model, the mirror force term, the second term on the right-hand-side (RHS), is neglected. Nonlinear terms are contained in the first and the last terms on the RHS. In a mass-less fluid electron model, one takes the limit $me→0$ while holding $mev||2$ and $mev⊥2$ as fixed. The mass-less fluid electron model used here differs from the previous model¹⁴ in that, here nonzero electron mass effect is included in the second term on the left-hand-side (LHS). This Laplace term has the desirable effect of reducing perpendicular grid-scale $E||$ fluctuations. These approximations lead to the following fluid Ohm's equation:

where n_e is the equilibrium electron density, $p||e$ is the equilibrium electron pressure, and $b̃=b̂+δB⊥/B0$ is the unit vector along the total magnetic field.

The main closure relation for this fluid-electron model is the linearized isothermal condition for electrons¹⁷ 

which, for MHD waves with $E||≈0$⁠, is used to derive an evolution equation for $δTe$⁠.¹⁵ Here, T_e is the equilibrium electron temperature. This isothermal condition is derived from the electron drift kinetic equation by assuming $ω≪k||vTe$ (v_Te being the electron thermal speed). It clearly breaks down for the zonal components, which have $k||=0$⁠. This assumption mainly affects the zonal field generation. The zonal flow, or the flux-surface-averaged ϕ, is generated by the gyrokinetic ion response, which is not affected by the electron isothermal condition. For perturbations that contain only the flux-surface-averaged components, Eq. (3) in its linear form gives $〈E||〉=0$ (⁠$〈…〉$ stands for flux-surface-average), which leads to $∂〈A||〉/∂t=−〈E||〉−〈∇||ϕ〉≈0$⁠, i.e., no zonal field generation. By using the above Ohm's law with the isothermal condition, we are further assuming that the n = 0 electron temperature perturbation is zero, and that the zonal field generation is primarily a nonlinear effect, which is captured by the nonlinear terms in Eq. (3). These assumptions are subject to further study with a fully kinetic electron model. We note that some desirable features of the n = 0 mode are included in the hybrid model. The model in linear simulations of the n = 0 mode yields no zonal electron density perturbation, $〈δne〉≈0$⁠, which is reasonable and resembles the often-used adiabatic electron response, $δne=ene(ϕ−〈ϕ〉)/Te$⁠. The geodesic acoustic mode (GAM) and collisionless damping of zonal flows are included because ions are gyrokinetic. The nonlinear mechanism for zonal flow generation is the $E×B$ and the magnetic fluttering nonlinearity in the ion gyrokinetic equation and the electron continuity equation. The zonal field (⁠$〈A||〉$⁠) is generated via the nonlinear terms in the Ohm's equation, which can be shown to be equivalent to the zonal field model of Chen and Zonca¹² [see their Eq. (2.30)–(2.35)].

We will discuss various nonlinear effects in the gyrokinetic equation

Here, $vG0=v||b+vD$ is the guiding center motion in the equilibrium magnetic field, $ϵ=mv2/2$ and

Here, q is the charge. We have suppressed the species index. The two nonlinear terms on the LHS of Eq. (5) are the E × B nonlinearity and magnetic fluttering nonlinearity, respectively. The nonlinear term in Eq. (6), in the case of electrons, gives rise to the last term in Eq. (3). The E × B nonlinearity due to $〈ϕ〉$ gives rise to the zonal flow shearing effect.

The quasi-neutrality condition is

Here, $np=−qi∑k⊥qneTi[1−Γ0(b)]ϕk⊥ exp (ikxx+kyy)$ is the ion polarization density. $δn¯i$ is the perturbed ion density calculated from the gyro-phase independent perturbed ion distribution. Notice that the equilibrium electron density $ne(r)$ is used in n_p. In simulations, we will vary the beam ion density to vary the instability growth rate. In doing so, the electron density is fixed. The polarization density term will also be fixed if the main ion density is replaced by the electron density as in the above equation. In effect, for the purpose of calculating the polarization density only, we assume that the beam ion temperature is the same as the main ion temperature. This is reasonable because the beam density is small. In the long wavelength limit, Eq. (7) also agrees with the MHD model where the polarization term is only proportional to the total ion mass density, independent of temperature. Nonlinear correction to n_p is neglected, so that n_p and $δn¯i$ are formally of the same order in the gyrokinetic ordering parameter.

We use a plasma equilibrium based on DIII-D discharge #159243 at time t = 805 ms. In this experiment, RSAE activity is created with neutral beam injection during the current ramp phase, resulting in multiple frequency sweeping RSAEs, which occur near the minimum in the magnetic safety factor profile (q_min).¹⁸ The density and temperature profiles of the thermal ions, the beam ions, and the electrons are shown in Fig. 1.

Linear simulation of the n = 3, 4, 5, 6 modes will be carried out using both the gyrokinetic ion/fluid electron hybrid model as described in Sec. II and a gyrokinetic ion/drift-kinetic electron model.¹⁹ Both models use $(A||,ϕ)$ to represent the fields and neglect $δB||$⁠, the parallel magnetic field perturbation. While GEM's kinetic electron model has been extensively benchmarked and verified with other gyrokinetic codes for the study of drift-wave turbulence, this is the first time it is applied to low-n toroidal Alfven eigenmodes. A direct comparison of the two models for low-n Alfven eigenmodes is crucial for either model. The simulation domain is chosen to be 0.1 < r/a < 0.9. The number of grids in the field-aligned coordinates is $(Nx,Ny,N||)=(256,32,64)$ for linear simulations and $(Nx,Ny,N||)=(256,64,64)$ for nonlinear simulations. There are two ion species, deuterium thermal ions and deuterium beam ions. The beam ions distribution is assumed to be Maxwellian for simplicity. The mode frequency and structure of the RSAEs are not sensitive to the beam particle distribution. The number of particles is 16 per grid per ion species, and 64 electrons in the kinetic electron model. The time step is $ΩcΔt=5$ in the hybrid model, and $ΩcΔt=3$ with kinetic electrons. Here, Ω_c is the proton gyro-frequency. In linear simulations, the toroidal mode number is first selected, then the domain size in the y-dimension is chosen such that the fundamental mode corresponds to the chosen toroidal mode. In terms of computational efficiency, a simulation with kinetic electrons uses about 10 times the computer time of the fluid electron model.

Figure 2 shows the mode frequency and growth rate from the two models, for toroidal mode n = 3, 4, 5, 6. The n = 1 and n = 2 modes are avoided because a high-n approximation is assumed in the field solvers. For n > 6, the dominant modes are generally driven by thermal temperature or density gradients and are also avoided. All four toroidal modes are unstable, with the n = 4 mode having the largest growth rate. All modes are located near r/a = 0.47 where the safety factor is at minimum with q_min = 2.945. The dominant poloidal mode number is given by m ≈ nq_min, as expected for RSAEs. The radial structure of the m = 11, 12, 13 poloidal components of the n = 4 mode is shown in Fig. 3. The n = 4 poloidal contour plot of ϕ with kinetic electrons is shown in Fig. 4. The mode is identified as RSAEs by the sensitivity of the mode frequency to small changes in q_min, as shown in Fig. 5. Here, the q-profile in Fig. 1 is shifted to vary q_min, while the shape of the profile is unchanged. A Doppler shift of 10 kHz (corresponding to toroidal rotation frequency of 2.5 kHz) is assumed and added to the simulation results. In Figs. 2, 3, and 5, good agreement is seen between the kinetic electron model and the hybrid model. The two models have been compared previously for ion-temperature-gradients (ITGs) and kinetic ballooning modes²⁰ in flux-tube simulations. The agreement seen here for RSAEs is better than in the previous study. The reason is that RSAEs are primarily an MHD mode excited by the energetic particles, and electron kinetic effects are not very important.

These linear results are in broad agreement with the experimental observations as well. The toroidal mode number can be obtained with cross-power spectrograms of ECE channels, and near t = 805 ms sweeping n = 2, 3, 4, 5, 6 modes have been observed. The mode location is found to be near r/a = 0.4 using ECE measurement of the electron temperature fluctuations.

This section focuses on nonlinear simulation of the n = 4 mode including the effects of zonal structures. Zonal structures are always generated whenever there is an n > 0 instability, due to nonlinear coupling of the instability to itself. The important question is whether these n = 0 perturbations play an important role in the nonlinear evolution of the driving mode. In drift-wave turbulence, it is well-known that zonal flows play a crucial role in regulating the turbulence and, in general, reduce anomalous transport. For EP driven Alfven waves, the role of the zonal structures is also of great interest and has been the subject of many theoretical studies. GEM nonlinear simulation of RSAEs has been carried out previously,¹⁵ but the n = 0 mode was not included. The n = 0 mode was included in Spong's study of RSAEs with a Landau-fluid model⁹ and was found to play an important role. Spong's simulation, however, assumes that a beam particle source is present that maintains the beam particle density profile. It does not address the importance of the n = 0 mode relative to the beam nonlinear effects. In gyro-fluid simulation of multiple TAEs,²¹ zonal structures have been found to be the dominant nonlinear saturation mechanism. Our goal here is to quantify the effect of n = 0 perturbations, on the n = 4 RSAE seen in the experiment. However, this is not a validation study of the nonlinear gyrokinetic model and a comparison of simulation and measurement in beam particle transport, RSAE fluctuation amplitude, etc., as this requires details of the beam particle distribution not available to us. According to Fig. 2, such a validation study should include the n = 3, 5, 6 modes to account for possible interaction among the excited n > 0 RSAEs. But this is beyond the scope of the paper. The role of zonal structures is largely determined by how strong the instability is, insensitive to details of the beam particle distribution. We therefore continue to assume a Maxwellian beam distribution. Additionally, we assume the thermal density and temperature to have a flat profile, with n_e = 3 × 10¹⁹/m³ and T_e = T_i = 1 keV. This modification does not significantly affect the RSAE, but it eliminates the ion-temperature-gradient-driven (ITGs), in particular, the n = 8 ITG, which otherwise becomes dominant either linearly when the RSAE growth rate is small or nonlinearly when the RSAE has saturated. In such a situation, it is difficult to obtain the RSAE saturation amplitude. The n = 8 mode, in general, should be included, because it is also generated by the self-coupling of the n = 4 mode.

The effects of fluid nonlinearity have been studied by many authors.^{8,10,11,21,22} The hybrid model used in this work is based on the nonlinear gyrokinetic equation of Frieman and Chen.²³ Only nonlinear terms that are formally second order in the gyrokinetic ordering parameter are retained. Thus, the third-order nonlinear effect of ponderomotive force²² is not included in the present model. As typical in gyrokinetic simulations, we employ a quasi-neutrality condition that includes only the linear ion polarization density; thus, the effect of the ion polarization nonlinearity¹¹ is also not included. By including only the n = 0, 4, 8 modes, we study the fluid nonlinear effect that primarily arises from the RSAE instability coupling to itself, and not the coupling among different eigenmodes.^10,22

In the following, we scan over the n = 4 growth rate and compare the initial saturation amplitude of the n = 4 mode from two types of simulations. All nonlinear simulations are carried out with the hybrid model. Variation of the n = 4 growth rate is achieved by multiplying the beam density profile of Fig. 1 with a constant factor, adjusting the thermal ion density accordingly to maintain charge neutrality. In the first simulation, only the beam particle nonlinear effects are retained, the thermal ions and electrons being treated as linear. In the second simulation, all three species are nonlinear. n = 0, 8 are generated only in the fully nonlinear simulation. The initial saturation amplitudes, measured by the root-mean-square value of the radial magnetic perturbation at r/a = 0.45, are shown in Fig. 6. A typical evolution history of the three modes (n = 0, 4, 8) in a full nonlinear simulation is shown in Fig. 7 for ϕ and Fig. 8 for $A||$⁠. All simulations are collisionless and there is no beam particle source. The n = 4 mode amplitude first grows exponentially, reaches a maximum, and then decreases in time, due to background damping. We are only interested in this initial saturation, not the long term behaviour of the mode. The range of growth rate in Fig. 6 corresponds to the original beam density multiplied by a factor of 0.4 to 1.0. As can be seen, for growth rate γ/ω < 0.03, the results from the two nonlinear models agree well, and both approximately follow the trapping scaling δ B_r ∼ (γ/ω)², indicating that the beam nonlinear effects are dominant. The trapping scaling describes the steady state saturation of a driven mode when there is no background damping, as in the kinetic-MHD model.⁴ Here, the mode is not in a steady state after the initial saturation; nevertheless, the initial saturation amplitude also follows the trapping scaling.

As the growth rate increases to γ/ω > 0.03, the thermal nonlinear effects become important. With the original beam density profile (the cases with the largest growth rate in Fig. 6), the saturation amplitude is reduced to less than half of the amplitude with beam nonlinearity only. Comparing fully nonlinear simulations with and without the n = 8 mode shows that the n = 8 mode is not important. This is to be expected, as self-coupling of the n = 4 mode generates fluctuations at twice the n = 4 mode frequency, which mismatches the n = 8 mode and is damped. The zonal structures have zero frequency (after GAM oscillations are damped) and can strongly affect the driving mode. Figure 7 shows that the n = 0 mode grows with the n = 4 mode in the linear stage, with a growth rate twice that of the latter. This is typical of a force driven process.⁸ The various force generated zonal structures are shown in Fig. 9. The self-coupling of the n = 4 mode gives a nonzero surface-averaged radial particle flux in the linear stage at ω_At = 747 of Fig. 7. The radial E × B flow, $〈Ey×Bδne〉$⁠, is shown in Fig. 9(a), and the magnetic fluttering flow $〈δBrn0u||e〉$ is shown in Fig. 9(b). All the other radial structures in Fig. 9 are for the nonlinear stage at the time ω_At = 1178 in Fig. 7, where the n = 4 mode amplitude is near its peak. A n = 4 eigenfunction with a characteristic wave number k_r will lead to a zonal structure with a characteristic wave number of 2k_r. Thus, fine scale zonal structures naturally appear for RSAEs because of the localized mode structure. This can be readily seen in Figs. 9(a) and 9(b), as the radial electron flows are entirely determined by the linear eigenmode structure. All the other zonal structures in the linear stage are also thus determined, with a real frequency of zero, and grow at twice the n = 4 growth rate. The beam particle density perturbation in Fig. 9(f) indicates that beam particles are redistributed outward. The electron zonal density, once nonlinearly generated, tends to stay unless further modified by nonlinear sources. The main ion zonal density [Fig. 9(e)] has a magnitude close to that of the electrons [Fig. 9(g)], as a result of the ion neoclassical shielding effect.¹³ The zonal electric potential of Fig. 9(c) induces a sheared zonal flow with the shearing rate²⁴ $γE=(r/q)∂(qvE×B/r)/∂r, vE=Er/B$⁠, and this is shown in Fig. 9(h). The relative amplitude of the thermal species density perturbation is comparable to that of the beam ions but has a finer radial structure which can cause significant modification of the n = 4 mode structure and leads to damping. Spontaneous growth of the zonal structures is not present in Figs. 7 and. 8. This is very different from high-n drift waves. Figure 10 shows the evolution of the n = 30 ITG mode and the accompanying zonal flows. This ITG simulation uses the kinetic electron model. In the ITG instability, the zonal structures are also force generated initially, but at large ITG amplitude the zonal structure grows spontaneously as a secondary instability, which eventually suppresses the ITG growth via nonlinear E × B shearing.

For two of the growth rates in Fig. 6, one corresponding to the original beam density profile and the other corresponding to original profile multiplied by 0.6, we have carried out simulations with only the thermal ions and electrons treated as nonlinear, while the beam ions treated as linear. Not surprisingly, this leads to a much larger saturation amplitude in the weak instability case (about a factor of 4 larger), and only a small increase from the fully nonlinear result in the strong instability case. One can say that for near marginal instabilities, saturation is caused by the beam particle nonlinear effect, whereas for a strong instability, the background nonlinear effects determine the saturation amplitude.

Zonal structures exert two effects on the n = 4 mode. The first effect comes from the zonal flow and zonal field. Zonal flow can cause shearing of the n = 4 mode, similar to zonal flow shearing of drift waves. The zonal magnetic field leads to local modification of the q-profile. The second effect comes from zonal density perturbations, or more generally the n = 0 component of the perturbed distribution functions of the thermal ions and electrons. A zero-frequency zonal component in δf acts as a modification to the equilibrium density and temperature profile. In a particle-in-cell (PIC) simulation, the nonlinearly generated n = 0 component in $δf$ is always retained. It appears that the importance of zonal flow and zonal field can be ascertained by setting the $〈ϕ〉=〈A||〉=0$⁠. But this will also zero out the collisionless damping of the zonal flow¹³ and zonal density and will lead to the results that are difficult to interpret. Instead, we run the simulation with $〈ϕ〉$ and $〈A||〉$ retained, but their effects excluded in the particle motion. That is, we neglect the zonal flow and zonal field contribution to the nonlinear terms on the LHS of Eq. (5). The other nonlinear zonal structure effect on n = 4, described by

is still retained. Thus, nonlinear modification to the equilibrium density profiles is retained. Removing zonal flow shearing turns out to have little effect. For instance, for the case with the largest growth rate in Fig. 6, the fully nonlinear saturation amplitude is $δBr/B0=5.4×10−4$⁠, and the result without the zonal flow shearing is $δBr/B0=5.5×10−4$⁠. This proves that the effect of zonal structures on the n = 4 mode comes primarily from the n = 0 component in δf, and the n = 0 component of the perturbed density and parallel flow for electrons. This conclusion is somewhat surprising, since Fig. 9(h) indicates that the peak value of the zonal flow shearing rate is actually larger than the n = 4 linear growth rate. We have carried out linear simulations in which the zonal potential of Fig. 9(c), multiplied by a constant scaling factor, is applied as a fixed external potential. With a scaling factor of 0.25, 0.5, and 1.0, the n = 4 growth rate is reduced by 1.3%, 4.7%, and 11.7%, respectively. The reduction to the growth rate is relatively small. With a scaling factor of 3 (corresponding to a peak shearing rate about equal to the mode real frequency), the growth rate is reduced to 52% of the original value. These linear results are consistent with the small effect of zonal flow shearing in nonlinear simulations. It appears that, for Alfven waves, the zonal flow shearing rate needs to be comparable to the mode real frequency to have a significant effect on the mode growth rate.

A zero-frequency density or δf perturbation is equivalent to a local modification of the plasma equilibrium. Such a perturbation naturally has strong radial dependence with a characteristic radial wavenumber twice that of the driving instability and can readily cause modifications of the driving mode structure. Figure 11 shows the evolution of the m = 12 component of $ϕ$ in the fully nonlinear run, for the case of the largest growth rate in Fig. 7, at $ωAt=750$ (black), 1120 (green), and 1300 (blue). The change in the mode structure is quite large at later times. By comparison, Fig. 12 shows the same m = 12 component of $ϕ$ at the same times, for the case with beam nonlinearity only. There is no visible change to the mode structure, even though the mode amplitude in the nonlinear stage is larger. We can conclude that, for strong instabilities where zonal structures are important, it is the n = 0 perturbation in the distribution function of the thermal ions, and the n = 0 density and current perturbation in the case of fluid electrons that cause nonlinear modification of the mode structure and enhance damping. Zonal flow shearing is not important.

In conclusion, we have carried out nonlinear simulations of a single RSAE with the self-generated n = 0 fluctuations. The initial saturation amplitude is found to follow the trapping scaling for weak instabilities. For stronger instabilities with $γ/ω>0.03$⁠, the effects of n = 0 perturbations, including zonal structures, are important. The nonlinearly generated n = 0 distribution for thermal ions and the n = 0 density and current perturbation for electrons cause modifications of the RSAE structure. Spontaneous zonal flow generation is not observed, and zonal flow shearing is found to have little effect on the saturation amplitude. These conclusions are obtained for a specific DIII-D RSAE experiment. Whether it can be extended to other types of Alfven eigenmodes, or even RSAE in other experiments remain of course to be investigated. In particular, the gyrokinetic simulations of this paper use the δf-method, therefore will be difficult to apply to situations with very large fluctuation amplitude, $δB/B∼10−2$⁠, as in some kinetic-MHD simulations.⁸ 

The main conclusions are reached by comparing the two types of simulations presented in Sec. IV, namely, the beam nonlinearity only simulation and the all species nonlinear simulation. Here, we explain why the beam nonlinearity only simulation is chosen for comparison. Related to this procedure is the meaning of “single-n simulation.” The meaning varies with the model. An Eulerian simulation, in which only the n = 0 mode and an unstable n > 0 mode are kept in both the fields and the distribution, is very different from a particle-in-cell simulation in which only the fields are filtered to include the same n = 0 and n > 0 modes, while all components are retained in the distribution. This is because, in general, no filtering technique is applied to the particle weights in particle simulation. Even if a single toroidal component is retained in the fields, all harmonics of that toroidal component, including n = 0, are generated nonlinearly, and this is important for a proper description of particle trapping in the wave. In the present hybrid model, when the thermal species are nonlinear, both the thermal ion zonal density and the electron zonal density are generated. The electron zonal density, once generated, is little affected by the self-consistent zonal field perturbations. It can be artificially removed from δn_e. The ion zonal density cannot be removed and is subject to the collisionless damping incurred by the self-consistent zonal flow. Thus, even if the n = 0 fields, ϕ and $A||$⁠, and the n = 0 electron density and current are set to zero, the simulation is not a “single-n” simulation. The n = 0 component in δf_i can still be convected by the RSAE wave fields and leads to mode structure evolution. In fact, this is what we tried in the first place for a “single-n simulation,” and the thermal ion nonlinear effect is found to be strong, making it difficult to observe the trapping scaling. The reason is now clear. When the n = 0 field is set to zero, collisionless damping of the ion zonal density is also removed, therefore incurring larger modification to the mode structure. The thermal species nonlinear effects can only be removed by treating the thermal species as linear, so that no n = 0 density perturbations are generated.

This work was supported by the U.S. Department of Energy's SciDAC project “Center for Nonlinear Simulation of Energetic Particles in Burning Plasmas” (DE-FG02-08ER54987) and project “Gyrokinetic Turbulence Simulations of Turbulent Transport” (DE-FG02-08ER54954). This research used resources of the National Energy Research Scientific Computing Center, which was supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.