Simulation of fast-ion-driven Alfvén eigenmodes on the Experimental Advanced Superconducting Tokamak

Youjun Hu, 2 Y. Todo, Youbin Pei, Guoqiang Li, Jinping Qian, Nong Xiang, Deng Zhou, Qilong Ren, Juan Huang, and Liqing Xu Institute of Plasma Physics, Chinese Academy of Sciences, Hefei, Anhui 230031, China Center for Magnetic Fusion Theory, Chinese Academy of Sciences, Hefei, Anhui 230031, China National Institute for Fusion Science, Toki, Gifu 509-5292, Japan Abstract Kinetic-MHD hybrid simulations are carried out to investigate possible fast-ion-driven modes on the Experimental Advanced Superconducting Tokamak (EAST). Three typical kinds of fast-ion-driven modes, namely Toroidicity-induced Alfvén Eigenmodes (TAEs), Reversed Shear Alfvén Eigenmodes (RSAEs), and Energetic-Particle continuum Modes (EPMs), are observed simultaneously in the simulations. The simulation results are compared with the results of an ideal MHD eigenvalue code, which shows agreement with respect to the mode frequency, dominant poloidal mode numbers, and radial location. However, the modes in the hybrid simulations take a twisted structure on the poloidal plane, which is different from the results of the ideal MHD eigenvalue code. The twist is due to the radial phase variation of the eigenfunction, which may be attributed to the non-perturbative kinetic effects of the fast ions. By varying the stored energy of fast ions to change the fast ion drive in the simulations, it is demonstrated that the twist (i.e., the radial phase variation) is positively correlated with the fast ion drive.

0][11][12][13][14][15][16][17] The Experimental Advanced Superconducting Tokamak (EAST) 18 has recently been upgraded to include four deuterium neutral beam lines.Figure 1 shows the beam injection geometry on EAST.It is expected that, with the NBI generated fast ions, various AEs will be routinely observed on EAST.In the present work, kinetic-MHD hybrid simulations using MEGA code [19][20][21] are carried out to investigate possible fast-ion-driven Alfv en eigenmodes on EAST.The simulations use an equilibrium with a flat profile of safety factor.Anisotropic slowing down distributions are used to model the distribution of the fast ions from the neutral beam injection.Perturbations of multiple toroidal mode numbers are included in the simulations.Plenty of modes in the Alfv en frequency range are found in the simulation.The frequency and radial width of the modes observed in the simulations are compared with the Alfv en continua to categorize the modes, which indicates the modes include the toroidicity-induced Alfv en eigenmodes (TAEs), 1,3,4,19,22 reversed shear Alfv en eigenmodes (RSAEs), 7,8,[23][24][25][26][27][28][29][30][31] and energetic-particle continuum modes (EPMs). 5,20,32,33The simulation results are also compared with the results of an MHD eigenvalue code, 34 which shows agreement with respect to the frequency, dominant poloidal mode number, and radial location of the modes.
6][37] The radial phase variation makes the twodimensional mode structure on the poloidal plane take a twisted structure.9][40] In the present hybrid simulations, two-dimensional mode twist on the poloidal plane is observed, which is different from the results of the ideal MHD eigenvalue code.By varying the stored energy of fast ions to change the fast ions drive in the simulations, the dependence of the radial phase variation on the energetic particles drive is examined, which indicates that the radial phase variation is positively correlated with the fast ions drive.
The remainder of this paper is organized as follows.Section II reviews the model used in the simulation code MEGA.Section III gives the equilibrium used in the simulation.The simulation results are given in Sec.IV.A brief summary is given in Sec.V.

II. SIMULATION MODEL
MEGA is a numerical code calculating the interaction of thermal plasmas and EPs in toroidal geometries. 20In MEGA, the thermal plasmas are described by the nonlinear full MHD equations, while the EPs are described by the drift-kinetic equation.The MHD equations solved by MEGA consist of the mass continuum equation the momentum equation (with modification by the effects of EPs) the Faraday's law and the equation of state where q, u, and p are the mass density, fluid velocity, and pressure of the thermal plasma, respectively, B and E are the magnetic field and electric field, the electric field E is given by E ¼ Àu Â B þ gðj À j eq Þ, the current density j is given by j ¼ l À1 0 r Â B, the vorticity X is given by X ¼ r Â u; l 0 is the vacuum magnetic permeability, C is the adiabatic constant (C ¼ 5/3 in the simulations presented here), g is the electric resistivity, and n are the artificial viscosity and diffusion coefficients chosen to maintain numerical stability (these dissipation coefficients also play a physical role of enhancing the damping of the modes in the MHD simulation that does not include kinetic damping from the thermal plasma), the subscript "eq" represents the equilibrium variables, and j 0 h is the current density of EPs without the contribution of E Â B drift (the contribution of E Â B drift disappears due to quasi-neutrality 19 ).Note that the effects of EPs on thermal plasma enter through the j 0 h term in the momentum equation (2) of the thermal plasma.This scheme of coupling EPs with thermal plasma (usually called "current coupling") is valid when the density of EPs is low so that the inertia of EPs can be neglected in the momentum equation.
In MEGA, the MHD equations ( 1)-( 4) are discretized in space by using a fourth-order finite difference scheme in the FIG. 1. Top view of the neutral beam injection geometry on EAST tokamak.The four beams all lie on the midplane of the device with the tangential radii of the beams being R tan 1 ¼ R tan 3 ¼ 1:26 m and R tan 2 ¼ R tan 4 ¼ 0:73 m (the tangential radius of a beam line is the perpendicular distance of the beam line to the axisymmetric axis of the device).The maximum source power per beam line is 2 MW.The maximum injection energy of particles is 80 keV.FIG. 2. Flux surfaces shape of EAST discharge #48916 at 4.5 s.The black rectangle indicates the computational box on the poloidal plane used in the simulation.The magnetic field at the magnetic axis B u0 ¼ þ1:72 T and the toroidal plasma current I p/ ¼ À417 kA.FIG. 3. The radial profiles of the safety factor, plasma pressure, and electron number density for EAST discharge #48916 at 4.5 s.The profile of the safety factor is flat and has a weak negative shear in the region ffiffiffiffiffiffi W t p 0:4 (the safety factor reaches its minimum value q min ¼ 2:373 at ffiffiffiffiffiffi W t p ¼ 0:40).The value of the safety factor at the magnetic axis q 0 ¼ 2.438.The value of the safety factor at the flux surface enclosing 95% of the poloidal magnetic flux q 95 ¼ 4.274.The stored plasma energy of the equilibrium is 127 kJ.right-handed cylindrical coordinates ðR; /; ZÞ.The time discretization uses the fourth-order Runge-Kutta method.The drift-kinetic description is employed for the EPs.The guiding-center motion of EPs is governed by the following equations:

022505
where X is the location of the guiding-centers, v k is the parallel (to the magnetic field) velocity, v rB is the rB drift given by k Þ; l is the magnetic moment, m h , Z h e, and X h are the mass, electric charge, and cyclotron angular frequency of EPs, respectively, B ? and B ? k are defined by respectively, where b ¼ B/B.In MEGA, the orbit of the EP guiding-centers is followed by using the fourth-order Runge-Kutta method.Expressed in terms of the guiding-center drift, the current density j 0 h appearing in Eq. ( 2) is written where f is the guiding-center distribution function of EPs and the last term on the right-hand side is the magnetization current.Note again that the E Â B drift does not appear in Eq. ( 9) due to the quasi-neutrality. 19In MEGA, the evolution of the distribution function of EPs is simulated by using the df particle-in-cell method.

III. EQUILIBRIUM
The EAST equilibrium used here was reconstructed by the EFIT code 41 by using the constraints from experimental diagnostics in EAST discharge #48916 at 4.5 s. Figure 2 plots the flux surfaces of the equilibrium and the computational box on the poloidal plane used in the simulation.The computational box is chosen to enclose the flux surface with ffiffiffiffiffi ffi is the square root of the normalized toroidal magnetic flux.
The profiles of the safety factor q, plasma pressure p 0 , and electron number density n e0 of the equilibrium are plotted in Fig. 3.The profile of the safety factor has a weak negative shear in the region ffiffiffiffiffi ffi W t p 0:4.The electron number density n e0 is used here to determine the mass density of the thermal deuterium plasma through the approximate relation q 0 % n e0 m i , where m i is the mass of deuteron.
In this work, the equilibrium distribution of the fast ions from the deuterium neutral beam injection is modeled by the anisotropic slowing down distributions, which takes the following form: where C is a constant, which is chosen to achieve desired stored energy of fast ions; w p is the normalized poloidal flux; w scale is a quantity characterizing the radial gradient; v is the velocity of fast ions; v b is the injection velocity of the neutral beam; D v is a small velocity (compared with v b ), which is used to set the cutoff width near v b ; k is the normalized magnetic moment defined by k ¼ lB 0 /e, where B 0 is the strength of the equilibrium magnetic field at the magnetic axis and e is the kinetic energy of fast ions; k 0 and D k characterize the peak location and the width of the distribution over the pitch angle; and v crit is the critical velocity for the collisional friction of fast ions with thermal electrons and ions being equal, which is given by 42 where m e and v te are the mass and thermal velocity of the electrons.In this work, the beam velocity is chosen v b ¼ 2.35 Â 10 6 m/s, which is the velocity of a deuteron with kinetic energy of 58 keV; v crit ¼ 1.89 Â 10 6 m/s, which corresponds to the critical velocity given by Eq. ( 11) evaluated with the electron temperature T e ¼ 2 keV.The small cutoff width near the beam velocity is chosen as D v ¼ 0.21 Â 10 6 m/s.The central pitch angle variable k 0 is chosen as k 0 ¼ 0.5 with the expansion width D k chosen as D k ¼ 0.3, which represents a reasonable distribution over the pitch angle based on the beam injection geometry on EAST.The parameter w scale , characterizing the radial gradient of the fast ion pressure, is chosen as w scale ¼ 0.4.
The pressure p 0 given in Fig. 3 is the total pressure, which should include the contributions from both the thermal plasma and fast ions.However, in the present work, for simplicity, the pressure p 0 given in Fig. 3 is assumed to be the pressure of the thermal plasma.

IV. SIMULATION RESULTS
In the simulations, all the MHD perturbations are set to be zero at the computational boundary.Initial perturbation within the computational boundary is set to be random magnetic perturbation of order dB/B 0 % 10 À10 (via the curl of a random vector field, which ensures the divergence-free nature of the magnetic perturbation).The numbers of grid points for cylindrical coordinates ðR; /; ZÞ are (128, 64, 128).The number of particle markers loaded is 2 20 , which approximately corresponds to one marker per spacial grid.To reduce numerical noise, the MHD perturbations are filtered to keep only harmonics of low toroidal mode numbers.In this work, harmonics out of the range À8 n 8 are filtered out every 1000 time steps, where n is the toroidal mode number.The time step dt is chosen to be dt ¼ 0.1/X h0 ¼ 1.2 Â 10 À9 s, where X h0 is the cyclotron angular frequency of the fast ions at the magnetic axis.This small time step (compared with the cyclotron period of fast ions) is chosen to meet the Courant condition of the full MHD system (for the present case, the Courant number dtV max =dx % 0:9 < 1, where V max is the maximum of the Alfv en speed in the computational region and dx is one of the spacial grid intervals).The ratio of the time step dt to the Alfv en time t A is dt/t A ¼ 2.65 Â 10 À3 , where t A is defined by t A ¼ R axis /V A0 with R axis being the major radius of the magnetic axis and V A0 being the Alfv en speed at the magnetic axis.In MEGA, the fluid velocity u is normalized by V A0 .In the simulations, the artificial viscosity and diffusion coefficients are chosen as ¼ n ¼ 1:0 Â 10 À7 V A0 R axis .The electric resistivity is chosen as g ¼ 1:0 Â 10 À7 l 0 V A0 R axis .
Although MEGA code uses cylindrical coordinates ðR; /; ZÞ in advancing the MHD equations and the orbit of energetic particles, it uses flux coordinates ðw; h; /Þ in analyzing the simulation results, where w is a magnetic surface label (in this article, w is chosen as w ¼ ffiffiffiffiffi ffi ), / is the usual toroidal angle, and h is chosen to make magnetic field lines straight on ðh; /Þ plane (the transformation Jacobian of the coordinates is proportional to h(w)R 2 , where h(w) is a flux function).The zero point of h coordinate is chosen at the low-field side of the midplane and the positive direction of h is chosen in the counterclockwise direction when viewed in / direction.In ðw; h; /Þ coordinates system, the radial component of the fluid velocity u w can be expanded in two-dimensional Fourier series over h and /, which is written where u ðcÞ wmn ðw; tÞ and u ðsÞ wmn ðw; tÞ are cosine and sine components of the expansion coefficients, respectively.Note that, in this expansion, the poloidal mode number m is always positive, while the toroidal mode number n can be negative.
Figure 4 plots the two-dimensional mode structure observed in the simulation at time slice t ¼ 0.25 ms, which shows that the n ¼ À8 mode has the largest amplitude.This mode localizes in a narrow radial region and takes a ballooning structure with a large dominant poloidal mode number (m ¼ 24), as is shown in Figure 4(h).Furthermore, after examining the temporal evolution of the n ¼ À8 mode, we found it is a purely growing mode with zero frequency.Considering this, this mode is identified as an MHD ballooning mode, which will not be discussed further in this paper.Similarly, the n ¼ À7 mode is also an MHD ballooning mode.Excluding these two modes, the dominant modes among the modes given in Fig. 4 are the n ¼ À2, n ¼ À3, and n ¼ À4 modes, which will be analyzed in turn next.

A. Identification of the fast-ion-driven modes
The radial structures of the various poloidal harmonics of the n ¼ À2 modes are plotted in Fig. 5, which shows that the m ¼ 5 harmonics is dominant and the mode amplitude reaches its peak at the radial location ffiffiffiffiffi ffi W t p ¼ 0:2.A general harmonic in the Fourier expansion (12) can be further written as with the amplitude A given by A ¼ ju wmn j ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi q and the phase angle a given by a ¼ atanðu ðsÞ wmn ; u ðcÞ wmn Þ, where atan() is an arc tangent function with two arguments that can determine the correct quadrant of the results.Using this, the instantaneous angular frequency of a harmonic is obtained by calculating the temporal change rate of the phase angle a, i.e., x ¼ da/dt.The growth rate is given by c ¼ dA/dt.The frequency and growth rate of the n ¼ À2 mode can be obtained by calculating the frequency and growth rate of the dominant poloidal harmonic at the radial location where the amplitude of the mode is maximal.Figure 6 plots the temporal evolution of the m ¼ 5 poloidal harmonic at the radial location ffiffiffiffiffi ffi W t p ¼ 0:2.The results show that the frequency of the mode x/2p ¼ 69.0 kHz and the linear growth rate c ¼ 2.71 Â 10 4 s À1 .
Using the formula in Eq. ( 13) and considering the sign of the mode numbers (m, n), the sign of the frequency, and the defined positive directions of the poloidal angle h and toroidal angle /, we can determine the toroidal and poloidal propagation directions of the n ¼ À2 mode.The mode propagates poloidally in þ ĥ direction and toroidally in À / direction.In terms of the physical quantities, the mode propagates poloidally in the diamagnetic drift direction of the ions and toroidally in the co-current direction, which is consistent with the general rules for the propagation direction of Alfv en modes excited by fast ions. 43o identify which kind of fast-ions driven mode the n ¼ À2 mode discussed above belongs to, we plot the frequency and the half height full width (HHFW) of the mode on the graphic of the Alfv en continua, as shown in Fig. 7, which indicates that the mode lies in the TAE gap.Also plotted in Fig. 7(a) are the m ¼ 4 and m ¼ 5 Alfv en continua in the cylindrical limit, which shows that the two continua are well separated from each other, indicating there will be only weak coupling between these two harmonics in the corresponding toroidal geometry.Furthermore, as shown in Fig. 5, the mode is dominated by the m ¼ 5 harmonic with the m ¼ 4 harmonic being much smaller, thus excluding the possibility of being a TAE mode.Also note that the radial location of the mode is near the location where the safety factor reaches its minimal value ( ffiffiffiffiffi ffi W t p ¼ 0:4).Considering these characteristics of the mode, it is reasonable to identify the mode as a RSAE.
Similar analysis can be applied to the n ¼ À3 mode.The results are plotted in Figs. 8 and 9. Figure 8 shows the dominant poloidal mode number is m ¼ 7 and all other poloidal harmonics are negligible.The frequency of the n ¼ À3 mode is x/2p ¼ 69.0 kHz, which happens to be identical to the   Considering these characteristics of the mode, it is reasonable to identify the mode as an EPM, instead of a gap mode.Similar analysis can also be applied to the n ¼ À4 mode.The results are plotted in Figs. 10 and 11. Figure 10 shows that there are two dominant harmonics for this case, namely, m ¼ 9 and m ¼ 10.The frequency of the mode is x/2p ¼ 82.8 kHz. Figure 11 plots the frequency and the HHFW of the dominant m ¼ 10 harmonic on the graphic of the MHD continua, which indicates the mode lies in the TAE gap formed due to the coupling of the m ¼ 9 and m ¼ 10 harmonics.Considering this and that the two dominant harmonics m ¼ 9 and m ¼ 10 are comparable in amplitude, it is reasonable to identify the mode as a TAE.

B. Mode twist
Figure 12 compares the two-dimensional mode structures calculated by MEGA and GTAW (an ideal MHD eigenvalue code 34 ), which shows agreement with respect to the radial location and dominant poloidal harmonics of the modes.However, the n ¼ À4 mode in the hybrid simulation takes a twisted structure on the poloidal plane, which is different from the results of the ideal MHD eigenvalue code.9][40] By varying the stored energy of fast ions to change the fast ions drive in the simulations, we examined the dependence of the radial phase difference of the n ¼ À4 TAE on the fast ions drive, which is shown in Fig. 13(a).The radial phase difference is defined to be the phase difference of the m ¼ 10 harmonic between the radial location ffiffiffiffiffi ffi W t p ¼ 0:173 and ffiffiffiffiffi ffi W t p ¼ 0:356, which is the radial range where the mode amplitude is significant, as is shown in Fig. 14. Figure 13(b) shows that the growth rate (indication of the fast ions drive) increases with the increasing of the stored energy of fast ions, while the frequency remains the same.Figure 13(a) shows that the phase difference increases with the increasing of the stored energy of fast ions, i.e., the mode twist is positively correlated with the fast ions drive.
V. SUMMARY Kinetic-MHD hybrid simulations are carried out to investigate possible fast-ion-driven modes on the EAST tokamak.Toroidicity-induced Alfv en eigenmodes, reversed shear Alfv en eigenmodes, and energetic-particle continuum modes are observed simultaneously in the simulations.The slow-sound approximation of the Alfv en continua proves to be useful in identifying the modes found in the simulations.The agreement between the hybrid simulations and linear eigenvalue analysis provides confidence in the simulation results.It is demonstrated numerically that the radial phase variation of the toroidicity-induced Alfv en eigenmodes is positively correlated with the fast ion drive in the hybrid simulations.The present work is limited to the linear properties of the modes and does not provide any analysis for the resonant interaction between the fast ions and the modes.This subject will be investigated in future by examining the phase space structure of the fast ion distribution function.

FIG. 4 .
FIG. 4. Contour of the toroidal harmonics of the radial velocity u w on the poloidal plane for different toroidal mode number n (the toroidal mode numbers are indicated on the corresponding figures) at time t ¼ 0.25 ms.The initial stored energy of the fast ions within the boundary flux surface is 35 kJ for this case.The back line on every figure indicates the last closed flux surface.

FIG. 5 .
FIG.5.The radial profiles of (a) amplitude ju wmn j, (b) cosine components u ðcÞ wmn , and (c) sine components u ðsÞ wmn of the n ¼ À2 mode.The time for plotting these figures is chosen to be at the moment when the sine component of the dominant m ¼ 5 harmonic reaches zero at the radial location where the amplitude of the mode is maximal (i.e., ffiffiffiffiffiffi W t p ¼ 0:2 for this case).The poloidal mode number range for this case is 0 m 64.

FIG. 7 . 45 FIG. 8 .
FIG. 7. (a) n ¼ À2Alfv en continua with the frequency and the HHFW of the n ¼ À2 mode plotted on.Also plotted on (a) are the m ¼ 4 and m ¼ 5 Alfv en continua in the cylindrical geometry limit.(b) Enlarged part of the Alfv en continua in the frequency range [60 kHz: 70 kHz] to show the tip of the continua near the zero magnetic shear point.The Alfv en continua is calculated in the slow-sound approximation.44,45

FIG. 9 .
FIG.9.n ¼ À3 Alfv en continua with the frequency and the HHFW of the n ¼ À3 mode plotted on.Also plotted are the m ¼ 7 and m ¼ 8 Alfv en continua in the cylindrical geometry limit.
frequency of the n ¼ À2 RSAE analyzed above.Figure9plots the frequency and the half height full width of the dominant m ¼ 7 harmonics on the graphic of the Alfv en continua, which indicates the mode intersect with the m ¼ 7 Alfv en continua.

FIG. 12 .
FIG. 12. Contour of the radial fluid velocity on the poloidal plane calculated by MEGA ((a) and (b)) andGTAW ((c) and (d)).The n ¼ À2 mode is an RSAE with frequency f ¼ 69.0 kHZ (67.9 kHz given by GTAW).The n ¼ À4 mode is a TAE with frequency f ¼ 82 kHz (83 kHz given by GTAW).The initial stored energy of fast ions used in the MEGA simulations is 36 kJ.
is maximal (i.e., ffiffiffiffiffiffi W t p ¼ 0:29 for this case).The poloidal mode number range for this case is 0 m 64.