Gyrokinetic simulations of reverse shear Alfven eigenmodes in DIII-D plasmas

A gyrokinetic ion/mass-less ﬂuid electron hybrid model as implemented in the GEM code [Y. Phys. , 837 (2007)] is used to study the reverse shear Alfv (cid:2) en eigenmodes (RSAE) observed in DIII-D, discharge #142111. This is a well diagnosed case with measurement of the core-localized RSAE mode structures and the mode frequency, which can be used to compare with simulations. Simulations reproduce many features of the observation, including the mode frequency up-sweeping in time and the sweeping range. A new algorithmic feature is added to the GEM code for this study. Instead of the gyrokinetic Poisson equation itself, its time derivative, or the vorticity equation, is solved to obtain the electric potential. This permits a numerical scheme that ensures the E (cid:2) B convection of the equilibrium density proﬁles of each species cancel each other in the absence of any ﬁnite-Larmor-radius effects. These nonlinear simulations generally result in an electron temperature ﬂuctuation level that is comparable to measurements, and a mode frequency spectrum broader than the experimental spectrum. The spectral width from simulations can be reduced if less steep beam density proﬁles are used, but then the experimental ﬂuctuation level can be reproduced only if a collision rate above the classical level is assumed. V C 2013 American Institute of Physics . [http://dx.doi.org/10.1063/1.4775776]


I. INTRODUCTION
The reverse shear Alfv en eigenmode (RSAE), also called the Alfv en cascade (AC), is a class of shear Alfv en waves that have been observed in various tokamaks. 1-3 RSAE resides in a radial location where the safety factor is at a minimum, typically well inside the plasma core, and detailed measurements of the mode fluctuations are usually difficult. It is only recently that detailed measurements of the mode structure on DIII-D became available from electron cyclotron emission (ECE) and electron cyclotron emission imaging (ECEI) as first reported by Refs. 4 and 5. More experimental details are later reported on the mode frequency as it sweeps in time, the details of the electron temperature fluctuations inside the ECEI viewing range and the fluctuation amplitude. 6 This DIII-D plasma provides an ideal case for comparing simulation results with direct measurement. Such comparison has been attempted with simulation codes based on various theoretical models, such as the kinetic-MHD model, 7 the full gyrokinetic model, 8 and hybrid models where ions are gyrokinetic but electrons are fluid-like.
In this paper, we present simulations of the observed RSAEs using the gyrokinetic ion/massless fluid electrons hybrid model as has been implemented in the df particle-incell (PIC) code GEM. 9 Both linear and nonlinear simulation results are reported. The most salient feature of RSAE is the frequency sweeping in time as q min , the safety factor at its minimum, decreases in time as a result of current diffusion.
Reproducing this sweeping is the first step toward comparing simulation with the experiment. As a result of using the conventional df method, the magnetic equilibrium is held fixed throughout any particular simulation. The RSAE frequency sweeping in time in the experiments is reproduced in simulations by running many simulations, each with a different q-profile. Since in the experiment RSAEs are observed to persist over a long time ($20 ms), much longer than the inverse growth time 1=c (c is the linear growth rate), the mode structure and amplitude observed in the experiments must be assumed to correspond to the nonlinearly saturated state. If the perturbative effect of the beam particles on the mode structure is small, as is frequently assumed to be the case for the toroidicity-induced gap modes, 10 then it is reasonable to compare the linear mode structure from linear simulations with the experimental data. However, RSAE is usually considered to be non-perturbative, i.e., the mode structure is strongly influenced by the beam ion distribution, and strong modifications of the mode structure due to the nonlinearly modified beam distribution can happen. It is, therefore, essential to compare the observed mode features with the results of nonlinear simulations.
Extensive linear simulations of RSAEs using various models have been reported recently, including simulations with the gyrofluid code TAEFL, 6,11 the fully kinetic continuum code GYRO, 12 and the particle-in-cell code GTC with a second order perturbative approach for electrons. 13 The present study will, therefore, focus on nonlinear simulations. It should be borne in mind that the present study is only an initial attempt at directly comparing the hybrid model simulations with the experimental measurements. While some of the observed features of RSAEs, such as the up-sweeping of the mode frequency in time and the saturated electron temperature fluctuation level, are confirmed by simulations, there are significant disagreements between the simulation and observation that are not well understood. All simulations presented here are for the n ¼ 3 mode and the n ¼ 4 mode (n is the toroidal mode number) as are observed at t $ 725 ms in DIII-D discharge #142111. 6,11,14,15 We started with the code as was presented in detail previously, 9,16 but in time learned that, for the plasmas considered here, a significant change concerning the field equation for the electric potential u is needed. This change is motivated by consideration of the effects of the E Â B motion on generating electric charge, which is needed in the gyrokinetic Poisson equation (the quasineutrality condition). If the finite-Larmor-radius effect is neglected then all species undergo identical E Â B motions, leading to no net charge separation because of quasineutrality in the equilibrium. Charge separation due to E Â B convection of the equilibrium density profiles occurs only as a result of FLR effect. It is desirable for the numerical algorithm to take advantage of this fact, i.e., to ensure proper cancellation of the E Â B motion such that only the FLR effect is retained. This is particularly important for cases where the modes of interest have a low toroidal mode number and the beam ion density is high, as in the DIII-D discharge studied here. It is impossible, however, to make use of the exact cancellation between E Â B motions in conventional gyrokinetic simulations, where u is obtained from the gyrokinetic Poisson equation, with the charge density of each species calculated from the distribution of that species. On the other hand, if we take the time derivative of the Poisson equation to obtain the vorticity equation for @/=@t, the exact cancellation between the linear E Â B motions can be made explicit. This is found to greatly enhance the robustness of nonlinear simulations, allowing us to carry out long-term simulations (typically lasting a physical time of 5 ms) to obtain a good steady state, from which a timeaveraged fluctuation amplitude and the wave frequency spectrum can be calculated. These simulation results are then compared with ECEI data. The conclusion we draw from this comparison study is: (1) The mode frequency, in general, agrees with the experimental data to within 10%; (2) with the current experimental input for the beam density profile and the classical collision rate of beam particles, the fluctuation level (in terms of the electron temperature perturbation) in the simulations is consistent with ECEI measurement, but the mode spectrum is more broad than ECEI data; and (3) agreement with ECEI data in both amplitude and spectral width can be achieved if a more flat beam density profile is used and a somewhat anomalous collisional mechanism can be assumed. This paper is organized as follows. In Sec. II, the new algorithm for obtaining the electric potential is explained. Linear results are presented in Sec. III, and nonlinear results are presented in Sec. IV. Summary and discussions are given in Sec. V.

II. OBTAINING U FROM THE VORTICITY EQUATION
We start with the gyrokinetic Poisson equation for the perturbed electric potential u Àn p ðx; yÞ ¼ dn i þ q b dn b À dn e : (1) Here, n p ¼ À P k ?
qn 0 T i ð1 À C 0 ðbÞÞ/ k ? expðik x x þ k y yÞ is the ion polarization density. The right-hand-side (RHS) contains the perturbed density for the main ions, the beam ions, and the electrons, all with a charge jqj ¼ e. The equilibrium density profiles also satisfy the charge-neutrality condition, n i ðrÞ þ n b ðrÞ ¼ n e ðrÞ. Equation (1) has been solved in all gyrokinetic simulations for u, including our previous studies of Toroidal Alfv en Eigenmodes (TAEs). 9,16 Equation (1) is also used in the early phase of this study, during which the following observations are made: (1) The thermal species kinetic equations for df contain the x Ã term, which is the instability drive for drift waves, but usually not important for energetic-particle-driven Alfv en waves. In fact, these terms are sometimes not included in simulations of energetic-particles-driven modes, in order to minimize the effect of thermal species drive and focus on the energetic particle physics. However, for the present DIII-D case, we find that doing so significantly changes the n ¼ 3 mode frequency and growth rate; (2) The E Â B motion is compressible in a toroidal field, which leads to a compressibility term in the electron continuity equation, @dn e =@t ¼ Á Á Á À n e r Á E Â b=B. This term is found to be unimportant in simulations of high-n modes, but important in the present case.
Both observations point to the need to carefully treat the E Â B motion. It is easiest to see this when we examine the vorticity equation. Taking the time derivative of Eq. (1) leads to where the summation is over the ion species and only linear terms directly related to the E Â B motion are shown. Here, hEi is the gyro-averaged electric field and the d-function indicates that deposition of the charge to the spatial grids is to be done for a charged ring. Without these FLR effects, the terms shown on the RHS of Eq. (2) cancel each other exactly if the equilibrium distributions for the ion species satisfy the charge neutrality condition. However, numerically such cancellation is not ensured. In particle-in-cell simulations, the marker particles are initially loaded according to the equilibrium distribution, but subsequently evolved in time without any procedure to enforce their charge neutrality. Deviation from the initially loaded distribution can arise from at least two sources. First, a local Maxwellian is not an exact equilibrium distribution for the thermal ions, because of the finite-orbit-width effect. Second, the marker distribution is subject to nonlinear modification if the fluctuation amplitude of the perturbed electromagnetic field is high. The finiteorbit-width effect is more significant for the energetic particles. If the energetic particle density is not very low compared with the thermal species density, the deviation of the energetic particle distribution from the assumed form can cause large artificial net charge. One can also see why this problem is especially severe for low-n modes. The ion polarization density is approximately given by n p ¼ Àn i r 2 ? / in the long wavelength limit, i.e., scales as $k 2 ? , whereas the terms on the RHS (including numerical errors!) scale as $k ? because E ? ¼ r/. The ratio between any numerical error on the RHS and the LHS scales as 1=k ? , i.e., large for low-n modes.
To remove the residual contribution from these linear E Â B terms in the absence of FLR effects, we can simply subtract the same ion terms from the RHS and drop the electron term Numerically, Eq. (3) can be implemented by looping through the ion markers twice in the deposition routine, treating the ions as gyrokinetic and drift-kinetic in turn, then take the difference between the two results. Equation (3) can be solved readily using the same Poisson solver as that used in the split-weight scheme. 17 Unlike in the split-weight scheme where the vorticity equation is solved in addition to the quasi-neutrality condition (Eq. (1)), here we only need to solve the vorticity equation, and / will be obtained by integrating @/=@t in time. Apart from this modification in the way / is obtained, the rest of the computational procedure for the hybrid model is the same as that described previously. 9,16 In summary, the gyrokinetic ion/fluid electron hybrid model used here neglects the perturbed parallel magnetic field perturbation. The effect of the equilibrium current as an instability drive has been implemented but not included in this work. The vector potential A jj is obtained by integrating the gauge equation @A jj =@t ¼ ÀE jj À r jj /, where E jj is obtained from the Ohm's law. The electron continuity equation provides the electron terms in the vorticity equation. The field quantities / and A jj are zero at the radial boundaries, while a particle moving across the boundaries re-enters the simulated volume at the point where its equilibrium trajectory intersects the boundary surface, with its weight unchanged. In linear simulations, the difference between Eqs.
(2) and (3) is often negligible, but the peak in the frequency spectrum of a nonlinear simulation using Eq. (2) can differ from that using Eq. (3) by 20%. Moreover, for cases with strong instability simulations using Eq. (2) often display nonlinear numerical instabilities leading to unphysical fluctuation amplitudes, whereas simulations using Eq. (3) do not exhibit such behavior.

III. LINEAR SIMULATION RESULTS
The equilibrium density profiles at the time $725 ms for DIII-D discharge #142111 are shown in Fig. 2. 11 The beam density at r=a ¼ 0:3 is a significant fraction of the thermal ion density. The beam distribution in velocity is assumed to be an isotropic, slowing-down distribution ÞÞdv is a normalization constant so that n b is the local beam density, obtained from classical fast ion transport calculation. There is no direct measurement on the velocity distribution. Since the neutral beam injection is tangential to the magnetic field, one expects an anisotropic velocity distribution, especially in the vicinity of the injection speed. In linear simulations, we have tested one anisotropic distribution function model with a Gaussian factor in the pitch angle variable, 18 and an isotropic Maxwellian distribution. In all nonlinear simulations presented in this paper, the isotropic distribution of Eq. (4) is used for simplicity. In the simulations, we characteristically vary the beam distribution function, including the radial profile, the critical velocity v I , and the cut-off velocity v b to show the effect of beam distribution. As a consequence, care is needed when comparing simulation results that depend on the details of the distribution with the experimental data, as, for instance, when we try to interpret the nonlinear results in Sec. IV. The linear growth rate is most sensitive to the beam distribution, but one expects the spatial mode structure and real frequency to be less sensitive.
Collisions between the beam particles and the thermal species determine the critical velocity v I , in Eq. (4), which is given by using the experimental parameters for the main ions. The critical velocity is the velocity where the fast-ion drag on thermal electrons is equal to that on thermal ions. It is not sensitive to the impurity concentration as long as the impurities have the same charge-to-mass ratio as the primary gas, which is normally the case for plasmas of deuterium, carbon, and oxygen (the DIII-D case). Once v I is given, linear results depend little on the beam particle collision frequency. Beam particle collisions in the form of drag and pitch-angle scattering are important in determining the saturated mode amplitude, and the collision frequency (or equivalently the slowing-down time) must be specified in nonlinear simulations. The actual form of the beam particle collision operator is given elsewhere. 16 For most of the simulations in this paper, a slowingdown time of 50 ms is used, which is estimated using the electron temperature and density in the core. The simulation domain is chosen to be 0:1 < r=a < 0:9. For linear simulations in this section, a total of 1 048 576 particles are used per ion species. The time step is X c ᭝t ¼ 10, where X c is the proton gyro-frequency. Both the thermal ions and the beam ions are deuterium. The grid setting is ðN x ; N y ; N z Þ ¼ ð128; 32; 32Þ in the three dimensions of the field-aligned coordinates. 17 We remind the reader that GEM is an initial value code, therefore only the most unstable mode can be clearly seen in linear simulations. The frequency of a subdominant mode can be traced in the parameter space only with an eigenmode tool. 12 If a flat beam density profile is used the plasma is stable. Ideally, the frequency spectrum of a stable plasma obtained in the simulation should contain discrete peaks that correspond to discrete eigenmodes, but the spectrum is frequently found to be dominated by high frequency (>200 kHz) fluctuations, yielding no clear results in the frequency range of interest, 20 kHz < f < 100 kHz.
The safety factor profile has a minimum q min near r=a ¼ 0:33, which decreases in the time window 700 ms < t < 770 ms (Ref. 11) according to q min ¼ À0:00337ðt À 700Þ þ3:33. RSAEs are characterized by the frequency upsweeping as typically seen in the experiments, or equivalently, rapid increase in frequency as q min decreases. The experimental spectrogram of Fig. 1 shows up-sweeping n ¼ 3 and n ¼ 4 RSAEs in the time window 700 ms < t < 750 ms. The mode frequencies of these two modes obtained from linear simulations are plotted in Figs. 3 and 4. Fig. 3 is intended to show the strong effects of the beam distribution function on the n ¼ 3 mode frequency, and Fig. 4 compares the simulation results with the experimental data for both the n ¼ 3 and the n ¼ 4 modes, and with the theoretical Alfv en continuum at the q min surface. A Doppler shift of 7:5 kHz for n ¼ 3 and 10 kHz for n ¼ 4 is added to all the simulation data, corresponding to a plasma rotation frequency of 2:5 kHz at the mode location. The experimental mode frequencies of the n ¼ 3 and n ¼ 4 RSAEs, calculated from the ECEI data, are also plotted. These curves represent the centers of best-fit Gaussians to the n ¼ 3 and n ¼ 4 peaks in the spectrogram, respectively, at each time point. The spectrogram itself is the sum of the spectral power density for all ECEI channels, with each timepoint representing a 1024-point fft window.
For the n ¼ 3 mode, simulations are first run with the beam density profile of Fig. 2. The mode frequency increases as q min decreases. This mode continues to be the most unstable mode until the minimum safety factor decrease to q min ¼ 3:17, at which point a jump in the frequency (of the most unstable mode) is seen. Examination of the corresponding eigenmodes at q min < 3:17 indicates that the most unstable mode becomes a TAE near the edge, with frequency nearly constant as q min decreases. These TAE data are not shown in Fig. 3. It is seen that, using the original beam density profile, the n ¼ 3 mode frequency sweeps from $60 kHz at t $ 700 ms to $70 kHz at t ¼ 740 ms, with a sweeping range of less than 10 kHz. This range is too small compared with the experimental sweeping range of about 40 kHz for the n ¼ 3 mode. More rapid sweeping can be produced by modifying the beam density profile. All the other simulation FIG. 2. The experimental safety factor, beam density, electron density, and electron and ion temperature profiles.
Phys. Plasmas 20, 012109 (2013) results of n ¼ 3 in Fig. 3 are obtained with the following analytic beam density profile In Eq. (6) s ¼ r=a, s 0 is the location where n b has maximum gradient, characterized by the local density scale length L h , and ᭝s specifies the radial width in which the gradient is localized. The advantage of using Eq. (6) is that one can choose s 0 to focus on mode activities at a given location and choose L h to match the local instability drive of the experimental profile. Modes far away from the s 0 surface can be avoided by choosing a small ᭝s. The frequencies from simulation with modified beam density profiles in Fig. 3 are obtained with a=L h ¼ 1:86, chosen to match the local radial dependence of the beam density profile of Fig. 2 at r=a $ 0:4. We have carried out linear simulations for both the n ¼ 3 and n ¼ 4 mode with other variations of the beam density profiles, e.g., by multiplying the original profile with a factor of e À2t (t is the square root of the normalized toroidal flux) or setting ᭝s ¼ 0:3 if Eq. (6) is used. We have also used a model beam particle distribution that is anisotropic in velocity. 18 All simulations with modified beam density profiles display frequency up-sweeping at a rate comparable to but smaller than the experimental observations. Fig. 3 shows the effects on the mode frequency of the beam particle density profile and the distribution function in velocity. The cutoff velocity of v b =v A ¼ 0:38 (v A is the Alfv en velocity at q min ) corresponds to the injection energy of 80 KeV. At this cut-off velocity, the simulated mode frequency is 73.5 kHz for q min ¼ 3:24 (including the Doppler shift effect), compared with the experimental value of 88 kHz. If the cut-off velocity is increased to v b =v A ¼ 0:54 and the total beam density reduced by half (with a moderately reduced total beam pressure), the mode frequency from simulation becomes 78 kHz. Also shown are the results from using a Maxwellian beam distribution with an effective beam temperature of 21 KeV, which are close to the results obtained with the slowing-down distribution with V b =V A ¼ 0:54. Fig. 3 confirms the non-perturbative nature of the RSAE mode 19 and highlights the need for more accurate beam distribution in code validation. The simulation data in Fig. 4 are obtained with the analytical beam density profile of Eq. (6) with a=L h ¼ 1:86, s 0 ¼ 0:33, ᭝s ¼ 0:2. The beam density at q min is half of the beam density of Fig. 2, and the cut-off velocity is v b =v A ¼ 0:54. The q min value is mapped into the discharge time according to q min ðtÞ ¼ À0:00337ðt À 700Þ þ 3:33 and the frequency is now plotted as a function of time. For reference, the theoretical Alfv en continuum given by the formula 20 is also plotted, with m ¼ 10 for n ¼ 3 and m ¼ 13 for n ¼ 4.
The RSAE mode frequency is expected to deviate slightly from the continuum due to the effects of toroidal geometry, 3 beam ions, 19 and plasma pressure gradient. 20 As can be seen, this theoretical continuum predicts the frequency sweeping rate for both modes with remarkable accuracy. In comparison, the sweeping rates from simulation are smaller. Fig. 4 shows that while the gyrokinetic ion/fluid electron hybrid model used here captures important features of the RSAEs, significant difference still exists between the simulated frequency and the experimental data. Fig. 3 suggests that some of the difference might be attributed to inaccuracy in the beam distribution used in simulations. We also note that the theoretical continuum result of Eq. (7) is obtained with the ideal MHD, in which the compressional component of the perturbed magnetic field, dB jj , is included. This effect is neglected in the current simulation model. It remains to be seen whether the inclusion of dB jj can account for the disagreement in the mode frequency between simulations and observation. The mode is characterized by a dominant poloidal harmonic, the m ¼ 10 component for the n ¼ 3 mode and the m ¼ 13 component for the n ¼ 4 mode, during most of the time. For n ¼ 3, the m ¼ 9 harmonic becomes significant only for q min < 3:18, and becomes comparable in magnitude to the m ¼ 10 mode when q min % 3:13. The frequency sweeping slows down as q min decreases near this value, while the m ¼ 9 harmonic grows in importance, suggesting a transition from RSAE (characterized by a single poloidal harmonic) to TAE (characterized by coupling between neighboring harmonics). 20 The n ¼ 3 mode structures in the poloidal plane and the radial dependences of the dominant poloidal harmonics for q min ¼ 3:33, 3.28, 3.23, and 3.18 are shown in Fig. 5. The plotted quantity is the potential /, which in the core region should have the same spatial structure as the electron temperature fluctuation dT e (see Sec. IV).
The poloidal mode structures of Fig. 5 (left column) show noticeable poloidal shearing as the radius increases, which has been identified as the fast ion induced shearing. 6 The comparison between the simulation and experiment on the mode shearing direction has been somewhat controversial recently. The shearing direction as represented by the q min ¼ 3:33 case (top left in Fig. 5) is robustly observed in the experiments. In this case, the toroidal magnetic field points out of the page (i.e., in the toroidal angle f direction if the coordinate system ðr; h; fÞ is right-handed). The mode shearing direction is clockwise as the radius increases, or in the ion diamagnetic direction, which is what has been seen in the experiment. However, Fig. 5 shows that this feature is not robust in simulations. In general, the shearing direction can change not only with q min but also with the beam density profile. The only solid conclusions we can draw from simulations are that, first, the shearing direction is determined by the direction of the toroidal field. Reversing the toroidal field produces an eigenmode with reversed shearing. Second, the shearing direction is independent of the plasma current direction relative to the toroidal field. In GEM, the plasma current direction is indicated by the sign of the safety factor, with a positive sign for a plasma current parallel to the toroidal field. Changing the sign of q does not change the mode structure. These observations are in agreement with the linear simulation results using the GTC code. 13 We also note that these simulation results suggest that the mode shearing is not determined by coupling of the RSAE to the kinetic Alfv en wave, as the latter is included in the hybrid model through the gyrokinetic ions.
It must be noted that as q min or n b ðrÞ is varied in simulations all other equilibrium quantities are fixed. Given the sensitivity of the mode characteristics to these equilibrium parameters, one expects that in the experiment the qÀprofile might evolve in time in a way that cannot be captured by elevating or lowering the entire q-profile, as is done here when q min is varied. The beam distribution can also nonlinearly evolve in response to the mode. We speculate that the observed mode shearing as a robust feature might be the result of the simultaneous evolution of the magnetic equilibrium and the self-consistently evolved beam distribution.

A. Obtaining the electron temperature perturbation in simulations
The mode structure and fluctuation amplitude obtained with ECEI are for the electron temperature dT e , which can be obtained in the massless fluid electron model by the isothermal condition 21 along the perturbed field line ðB 0 þ dB ? Þ Á rðT 0 þ dT e Þ ¼ 0; (8) or linearized to give B 0 Á rdT e ¼ ÀdB ? Á rT 0 : Equation (9) can be combined with the ideal MHD relation dB ¼ r Â ðn Â B 0 Þ and dn ? =dt ¼ E Â b=B 0 to yield an evolution equation for dT e @dT e @t ¼ À Equation (9) has been incorporated in the Ohm's law equation for E jj and dT e does not appear in the dynamical equations. 16 Equation (10) is evolved in the simulations for the sole purpose of comparison with the experimental data. It is preferred to directly inverting the magnetic differential equation, Eq. (9), because the problem of solvability does not arise. However, if in the nonlinear stage a stationary perturbation with x ¼ 0 is present, Eq. (10) will produce a temperature perturbation that grows secularly with time. To remove such a slow growing component, a damping term, Àc Te dT e , is added to the RHS of Eq. (10), with c Te ( x RSAE . Since Eq. (10) is completely passive in the simulation, such a damping term has negligible effect on the dT e in the frequency range of interest. B. Saturation amplitude of the temperature perturbation Parameters, such as grid sizes and the time step, for the nonlinear simulations are identical to that used in linear simulations, except that the number of particles is increased to 4 194 304 per ion species. The beam distribution function is the same as that of case (2) in Fig. 3, i.e., using the analytical beam density profile of Eq. (6), with the local beam density magnitude and gradient equal to the profile in Fig. 2, and a cut-off velocity of v b =v A ¼ 0:38. It is widely believed that the crucial nonlinear effects that determine the single-n mode saturation level is the energetic particle nonlinearity, 10 with the thermal plasma nonlinear effect negligible. However, in the present model with multiple species, the fluctuation potential / is calculated with the quasi-neutrality condition (specifically, from the vorticity equation derived from the quasi-neutrality condition), neglecting the nonlinear effects of the thermal species while keeping that of the beam particles is inconsistent, due to the same consideration that leads to Eq. (3). The nonlinear E Â B convection from all species should nearly cancel each other if FLR effects are small. In a multi-species model with the charge density of each species separately calculated, keeping the E Â B nonlinearity of one species while neglecting that of other species leads to numerical charge accumulation. We, therefore, include the nonlinear effects (the E Â B nonlinearity and the magnetic fluttering nonlinearity, but not the so-called parallel nonlinearity, $E jj @df =@v jj ) of all species in the following simulations, but at the end of the section give an example of simulation with only the beam nonlinearity to show the difference.
In the following simulations, the analytic beam density profile, Eq. (6), is used. With s 0 ¼ 0:4, ᭝s ¼ 0:2, and a=L h ¼ 1:865, this analytic beam density profile near r=a ¼ 0:4 is the same as that of Fig. 2, but flat outside this region. It is found that with the equilibrium profiles of Fig. 2, a TAE-like mode near the edge of the simulation box, r=a $ 0:8, can become unstable, but the growth rate of this mode is affected when the boundary is moved slightly. The analytic beam profile is chosen to avoid this edge mode and allows us to focus on the RSAE in the core region.
All nonlinear simulations are run for 100 000 time steps, corresponding to a physical time of 5 ms. At each time step, the root-mean-square (rms) value of the temperature perturbation on each flux surface, rmsðdT e Þðr; tÞ, is calculated. This quantity is then averaged over the time window of 1:5 ms < t < 5 ms to give an estimate for the steady state fluctuation level. The temperature fluctuation levels (in keV) at four radial locations for various q min values are given in Table I for n ¼ 3 and Table II for n ¼ 4.
At the time of t ¼ 725 ms (q min ¼ 3:246), the relative rms fluctuation level is also obtained: dT e =T e ðr=a ¼ 0:4Þ ¼ 1:1% for n ¼ 3 and 0.75% for n ¼ 4. The time history of rmsðdT e Þðr; tÞ at r=a ¼ 0:4 for n ¼ 3 at the time of t ¼ 725ms (q min ¼ 3:246), with large fluctuations of the amplitude in time, is shown Fig. 6. Typically, the linear stage is followed by a large initial overshoot of the mode amplitude, followed by a steady state at significantly lower amplitude. This is a feature seen in all of the nonlinear simulations. The average magnitude of thermal ion weights follows the same temporal pattern, because thermal ions are not resonant with the mode and respond adiabatically, with the weight proportional to the mode amplitude. The beam ions, on the other hand, respond to the mode non-adiabatically, and their average weight magnitude stays roughly at the level immediately after the initial overshoot. One, therefore, is concerned whether the discrete particle noise due to the large magnitude in the beam particle weights can lead to large noise that overwhelms the low level of physical fluctuation at later times. As mentioned above, the number of particles per ion species in the nonlinear simulations is already increased to four times that used in linear simulations. Fig. 7 compares the results of three simulations, with a total of 4 194 304 (black), 8 388 608 (red), and 1 048 576 (green) particles per species, respectively, for the q min ¼ 3:23 case. The timeaveraged temperature fluctuation level is dT e ðr=a ¼ 0:4Þ ¼ 0:0085 KeV, 0:0094 KeV, and 0:01 KeV, respectively. Given the large fluctuation of dT e in time and the finite simulation time, we consider the result to be converged in the number of particles.
The beam density profile at the end of the simulation for the n ¼ 3, q min ¼ 3:246 case is shown in Fig. 8, in comparison with the initial profile. The RSAE mode causes a global relaxation of the profile. The beam density at the inner boundary (r=a ¼ 0:1) is reduced by $10% from the initial value. We found that if a fixed boundary condition, i.e., df ¼ 0 at the boundaries, is used for the beam particles, the steady state fluctuation level is much increased. This is because fixing df ¼ 0 at the boundary is equivalent to an unphysical source of beam particles at the inner boundary that prevents the profile from relaxing, and this leads to a larger saturation amplitude.
For near-marginal instabilities, the saturation amplitude most sensitively depends on the density scale length L h , but the scale length of a=L h ¼ 1:865 used in the above simulations is not near marginal (the instability threshold is about a=L h ¼ 0:6), and we found that the saturation amplitude is insensitive to a=L h for moderate variations in L h . For instance, for the n ¼ 3 mode and for q min ¼ 3:246, a value of a=L h ¼ 1:3 leads to a temperature fluctuation of dT e ðr=a ¼ 0:4Þ ¼ 0:0086 KeV, a moderate reduction from the value of 0.0137 at a=L h ¼ 1:865. The critical velocity v I in the slowing-down distribution can have a large effect on the saturation amplitude. Again for the n ¼ 3, q min ¼ 3:246 case, reducing the critical velocity from v I =v Te ¼ 0:07 to v I =v Te ¼ 0:035 reduces dT e ðr=a ¼ 0:4Þ from 0.0137 to 0.0069. The effect of reducing v I is the decrease in the number of beam particles at high energy, given the total beam particle density fixed. Since resonant particles are at high energy, reducing v I leads to weaker instabilities. Figs. 9 and 10 show the frequency spectra of the n ¼ 3 and n ¼ 4 mode at r=a ¼ 0:4, respectively. The frequency spectrum is obtained by first decomposing / at the outer midplane (h ¼ 0) in the yÀdirection, then Fourier transform the k y ¼ 2p=L y component in time. This causes the asymmetry in the frequency spectrum, as the mode propagates in the ion diamagnetic direction, x=k y < 0.

Comparison with ECE Imaging
ECEI data in the time window between t ¼ 0.7245 and t ¼ 0.7255 are analyzed to yield the experimental mode characteristics. Quantifying the spectral width of a sweeping FIG. 7. Convergence test with respect to particle numbers. mode with a sweeping rate comparable to the spectral width is nontrivial. In this study, we have simply taken a singletimepoint (t ¼ 725 ms) slice through a spectrogram with an appropriately chosen fft window length. The length of the fft windows in this case was empirically chosen to be 1024 points long (about 1 ms) to minimize the spectral width of the resulting mode peaks. In practice, the spectrogram could be constructed using longer fft windows (resulting in an artificial broadening in the time domain) or shorter fft windows (resulting in an artificial broadening in the frequency domain). By choosing the optimum balance between these two effects, while still acknowledging some unavoidable broadening due to sweeping, the resulting mode widths can, thus, be interpreted as an upper-bound.
The ECEI data thus analyzed indicate the following mode properties: for the n ¼ 4 RSAE, the spectral peak is located at f 0 ¼ 66:9 kHz, with a width of fwhm $1 kHz; The root-mean-square amplitude for the electron temperature perturbation is 1.3%. For the n ¼ 3 RSAE, the peak is at f 0 ¼ 80:7 kHz, with a width of fwhm $1 kHz and root-meansquare amplitude of 1.1%. The rms amplitude is obtained by averaging over the entire ECEI array of pixels (20 horizontal rows of pixels, 8 in each row). The horizontal row of ECEI pixels covers a radial range of approximately ᭝r=a ¼ 0:21 at the outer mid-plane. The ECEI fluctuation levels show weak dependence on the radial location, which suggests that the ECEI window is well aligned with the RSAE mode location, in our simulation roughly in the region 0:2 < r=a < 0:5. Here, we compare the ECEI data with simulation results at r=a ¼ 0:4, which is about the same as r=a ¼ 0:3 in Tables I  and II. We note that direct comparison between measurement and simulation at an exact physical location should be avoided, as this will be meaningful only if the magnetic equilibrium used in the simulation, in particular the Miller approximation of the surface shapes, has very high accuracy.
We can now summarize the nonlinear simulation results in comparison with ECEI data. The fluctuation level of both modes, dT e =T e ðr=a ¼ 0:4Þ ¼ 1:1% for n ¼ 3 and 0.75% for n ¼ 4, is in reasonable agreement with the ECEI data analyses. The spectral width of each mode, however, is significantly larger than ECEI data ($1 kHz for each mode), as can be seen from Figs. 9 and 10. We note that the mode frequency sweeping rate is about 1 kHz/ms, so that on a fft window of about 1 ms the uncertainty of any frequency estimate is about 1 kHz. However, allowing an uncertainty of $1 kHz in the estimated experimental mode width will not change the conclusion. We have also verified that small variations (e.g., within 20%) of the beam density gradient scale length and the velocity dependence (in terms of the critical velocity) will not change the conclusion either.

C. Collisional effects
The role of energetic particle collisions with the background plasmas is well-known for the nonlinear saturation of a single-n, energetic particle driven instability. 10 The mode initially saturates due to resonant particle trapping in the wave field and flattening of the beam distribution in the vicinity of resonance. In the absence of collisions, non-resonant particles cannot become resonant and the mode amplitude will decay to zero because of finite background damping. Collisions, both pitch-angle scattering and velocity drag, can bring nonresonant particles into resonance with the wave, leading to saturated steady state with a finite amplitude. As the collision rate increases, the mode amplitude increases. These theoretical findings have been verified by simulations repeatedly for the toroidicity-induced Alfven eigenmodes. 22,23 All the previous results are obtained with a collision frequency of ¼ 1=s s ¼ 20=s, or =x A ¼ 3:5 Â 10 À5 . For n ¼ 4, q min ¼ 3:246 the effect of collisions is shown in Fig. 11. Similar to TAEs, the RSAE saturation level is also seen to increase with the collision rate.
Although collisions can have a large effect on the fluctuation amplitude, this is not through changing the linear growth rate, as the collision rate is much smaller than the growth rate. In linear simulations of RSAEs, the linear growth rate has been compared with the estimated experimental spectral

012109-9
Chen et al. Phys. Plasmas 20, 012109 (2013) width, 11,13 with the implicit assumption that the spectral width from a nonlinear simulation is correlated with the linear growth rate. In general, our nonlinear simulations bear out this correlation between the linear growth rate and the nonlinear spectral width, but as seen above the nonlinear spectral width is much larger than the experimental width. These considerations suggest the following possibilities which can explain the large discrepancy between the spectra of Figs. 9 and 10 and the ECEI measurement of the spectral width. The actual beam density profile at the time of t ¼ 725 ms differs from the profile of Fig. 2, either with lower beam density at the mode location or with weaker density scale length (or both), so that the mode has smaller growth rate, consequently the spectral peak becomes narrower. On the other hand, beam ion collisional effects are not limited to classical pitch-angle scattering and velocity space drag (these are the only collisional processes included in simulations). Small scale drift wave turbulence can act as an anomalous collisional mechanism for the beam ions, as demonstrated in a recent study. 24 Perhaps more importantly, an additional scattering mechanism can arise from the many simultaneously present Alfv en waves, as can be seen in Fig. 1, and these waves can lead to a multitude of wave-particle resonances and cause substantial additional stochastic transport in fast-ion phase space. 25 This can be demonstrated with a simulation of the n ¼ 4 mode, with the beam density scale length doubled (a=L h ¼ 0:93) and the collision rate increased to ¼ 80=s. The fluctuation level is dT e =T e ðr=a ¼ 0:4Þ ¼ 0:7%, similar to the previous value. However, the mode spectrum, shown in Fig. 12, is very different from Fig. 10. Not only does the spectral width compare more favorably with the ECEI estimate, but also the peak frequency at $60 kHz (including a Doppler shift of 10 kHz for n ¼ 4) is in better agreement with the measurements.

D. The effects of thermal species nonlinearity in the hybrid model
Finally, we test the effect of thermal species nonlinearity in the hybrid model as used here. Fig. 13 shows the time history of dT e ðr=a ¼ 0:4Þ for two identical runs with n ¼ 4, a=L h ¼ 0:93, and ¼ 200=s, differing only in whether thermal species nonlinear terms are included. It is clearly seen that without the thermal species nonlinear terms the temperature perturbation becomes much larger. It is tempting to conclude, based on this observation, that the mode saturates primarily due to background nonlinearity. However, the term "background nonlinearity" is usually understood differently. In kinetic-MHD simulations where the thermal plasma is described by the ideal MHD equations, 7 the background nonlinear effect is usually understood to be the nonlinear terms in the ideal MHD equation, and it is usually true that such background nonlinear terms are not important in determining the Alfv en wave saturation amplitudes, at least for near marginal instabilities. Because of the difference between the MHD model and the hybrid model used here, neglecting the background nonlinearity in MHD is not equivalent to neglecting the thermal species nonlinear terms in the hybrid model. Therefore, the simulation results of Fig. 13 do not lead to the conclusion that the background nonlinearity, as the term is usually understood in kinetic-MHD simulations, is important. Based on the same consideration that leads to the modified vorticity equation in Sec. II, the simulation without the thermal species nonlinearity while keeping the beam particle nonlinearity should be viewed, for the multiple species model studied here, as numerically inconsistent.

V. DISCUSSION AND SUMMARY
The nonlinear simulations of Sec. IV all start from a white noise perturbed beam distribution (random initial

012109-10
Chen et al. Phys. Plasmas 20, 012109 (2013) values of particle weights) at sufficiently low amplitude. Typically, the instability grows out of this initial noise, goes through a linear stage with exponential growth, at the early nonlinear stage reaches a large peak in amplitude and then decreases rapidly to a steady state. In comparing the temperature fluctuation level of the steady state with the experimental measurement, the validity of a multi-scale analysis in time is tacitly assumed, namely, the observed mode activity in the experiment at any time corresponds to the saturated steady state in the simulation that uses the equilibrium profiles at this time. In order to have good statistics in obtaining the time-averaged fluctuation level, we typically run the simulations for 5 ms or longer, but on this time scale q min would have changed by an amount of $0:017 according to q min ðtÞ ¼ À0:00337ðt À 700Þ þ 3:33, a significant change for RSAEs. Moreover, when setting up the simulation for a different time, the entire q-profile of Fig. 2 is moved down by the amount of change in q min , with the shape of the profile unchanged. No attempt is made during the simulation to modify the evolution of the q-profile over the simulation time scale. The same assumption of the validity of such multi-scale analyses is made in global drift-wave turbulence simulations of anomalous transport of the core species. It is very difficult to verify the assumption directly, because to do so requires simulating the self-consistent evolution of both the waves and the background equilibrium. One can say that in the real plasma the RSAEs are always already in the nonlinear stage as q min decreases in time, and this nonlinear state evolves in response to the evolving q-profile. In the experiment, a robust shearing direction of the mode structure is seen, but the linear mode structure in the simulation varies with q min as shown in Sec. III. It is quite possible that the robust shearing of the mode is the result of selfconsistent evolution of the mode in the evolving magnetic equilibrium.
In summary, we have carried out linear and nonlinear simulations of the n ¼ 3 and n ¼ 4 RSAEs observed in DIII-D #142111, and compared the simulation results with ECEI measurement. The linear mode frequency and the frequency sweeping range are in reasonable agreement with the ECEI data. The saturated mode amplitudes also compare well with measurements, but the mode frequency spectral width obtained the original beam density profile is much broader than the ECEI data. These results represent an initial attempt at comparing the hybrid model with the direct experimental measurements. Nevertheless, this comparison exercise already identifies important areas for future research. The linear simulations indicate that there is still a disagreement between the simulation and measurement, in both the mode frequency and the frequency sweeping range, and this disagreement is unlikely due to uncertainties in the experimentally inferred beam distribution. This suggests the need of improving the simulation model, perhaps by including the parallel magnetic field perturbations. The nonlinear simulation results indicate that the beam density profile should be less steep (more close to stability threshold) than the estimated beam density profile, and a physics-based model for anomalous collisional processes should be included in nonlinear simulations.