Magnetohydrodynamic Vlasov simulation of the toroidal Alfven eigenmode

A new simulation method has been developed to investigate the excitation and saturation processes of toroidal Alfven eigenmodes (TAE modes). The background plasma is described by a magnetohydrodynamic (MHD) fluid model, while the kinetic evolution of energetic alpha particles is followed by the drift kinetic equation. The magnetic tluctuation of n=2 mode develops and saturates at the level of 1.8X 10m3 of the equilibrium field when the initial beta of alpha particles is 2% at the magnetic axis. after saturation, the TAE mode amplitude shows an oscillatory behavior with a frequency corresponding to the bounce frequency of the alpha particles trapped by the TAE mode. The decrease of the power transfer rate from the alpha particles to the TAE mode, which is due to the trapped particle effect of a finite-amplitude wave, causes the saturation. From the linear growth rate the saturation Ievel can be estimated. 0 1995 American Institute of Physics.


I. INTRODUCTION
The toroidal Alfven eigenmode (TAE mode)' has recently become the focus of attention for fusion physicists, since it can be excited resonantly with alpha particles of 3.52 MeV which are produced from deuterium-tritium reactions. TAE modes are destabilized when the alpha particle distribution has a density gradient. Fu and Van Dam' studied the linear stability of the TAE mode with alpha particles in terms of a linearized drift kinetic equation, and found that there are two conditions for the TAE mode to be unstable. The first condition requires that the scale length of the alpha particle density gradient be sufficiently small. The second condition requires that the alpha particle destabilization effects overcome the damping effects. Though the electron Landau damping is investigated in Ref. 2, there are many other important damping mechanisms such as ion Landau damping,3 continuum damping,4-6 collisional damping, and radiative damping.7 One of the unresolved, but significant problems of the TAE mode is the saturation mechanism and the saturation level. Sigmar et al.* analyzed alpha particle losses using a Hamiltonian guiding-center code for a given linear TAE mode, and found that for the amplitude of B,IB,>lO-", a substantial fraction of alpha particles can be lost in one slowing down time. This indicates that a precise knowledge of the saturation level and the saturation mechanism is crucial for ignited tokamak plasmas. Breizman et aZ.9 discussed the saturation level in the context of the nonlinear Landau damping which is studied by O'Neil." Wu and White" studied the same physics model as Ref. 9, simulating the nonlinear alpha particle dynamics and the evolution of a linear TAE mode employing the linear dispersion relation obtained by the linear analysis based on the alpha particle distribution. They found that modification of the particle distribution leads to mode saturation.
Computer simulation is a powerful tool to elucidate the saturation mechanism of the TAE mode. Spong et al.'2*'3 carried out linear and nonlinear simulations in which the back-ground plasma is described by a reduced magnetohydrodynamic (MHD) fluid model and the alpha particles by a gyrofluid model. They argue that the saturation occurs due to the nonlinearly enhanced continuum damping and the nonlinearly generated EXB flows. In their work, alpha particles are represented only by two components, namely, the density and the parallel flow velocity. Therefore, their method could not analyze the kinetic effects of the finite-amplitude TAE mode on the alpha particles which are the energy source for the mode. Park et al.14 developed another simulation technique in which they make use of an MHD fluid model and super particles as the energetic particles. Their method seems to suffer numerical noises which unavoidably arise because of a limited number of super particles, though it contains kinetic characters of alpha particles.
We have developed a new simulation method that enables us to investigate excitation and saturation of toroidal Alfvdn eigenmodes. In this method, the background plasma is described as a full-MHD lluid, while the kinetic evolution of the energetic alpha particles is followed by the drift kinetic equation. Both the MHD and drift kinetic equations are solved by a finite difference method. This new method can deal with the kinetic characters of alpha particles with nonlinear MHD waves free from numerical noises of particle discreteness. As is usual in nonlinear dynamics, high frequency modes are generated in the drift kinetic equation coupled with the MHD equations. As a matter of fact, the finite difference method eliminates extremely short scales which are comparable to the grid sizes and spoils the collisionless property of the drift kinetic equation in such scales. However, it can simulate the larger scales soundly.
In the present paper, we study the excitation and saturation process of the TAE mode using the new simulation method mentioned above. We focus particularly on the n =2 TAE mode and its nonlinear evolution including the generation of the n =0 mode. The plasma model and the simulation method are presented in Sec. II. In Sec. III we present the simulation results, comparing with the linear theory, and dis-cuss the saturation mechanism. We analyze the power transfer rate and show that the decrease of the rate causes the saturation of the instability. We show that the saturation is caused by the trapped particle effect of the finite-amplitude wave. Finally, a conclusion is given in Sec. IV.

II. SIMULATION MODEL
In our model, the background plasma is described by the ideal MHD equations and the electric field is given by the MHD description which is a reasonable approximation under the condition that the alpha density is much less than the background plasma density: dP ,=-wpv,. (1) where h is the vacuum magnetic permeability, y is the adiabatic constant, and all other quantities are conventional. The energetic alpha particles are described by the drift kinetic equation constructed with the following equation of motion of an alpha particle under the guiding-center approximation with the EXB, grad-B, and curvature drifts, a $ .c=eavd.E+p -g B, E= )m&;+pB, ExB V,=F ' (10) vp$ i P a ";K---ihiE xb, 1 where UII is the velocity parallel to the magnetic field, p is the magnetic moment which is the adiabatic invariant, and K is the magnetic curvature vector. From Eqs. (6)-( 11) we can obtain the drift kinetic equation which describes the temporal evolution of the alpha distribution function in the phase space hq ~-4: &A x,q ,pEL)= -g [(y +v,+~~)f.l-~(afL (12) avll (13) To complete the equation system in a self-contained way, we should take account of the effects of the alpha particles on the background plasma. For this purpose, we invoke the Maxwell equation for the perpendicular component of the electric field: A E&F2 Erj where, in Eq. (14), the second and third terms are the polarization current and the diamagnetic current of the background plasma, respectively, and the fourth term is the current of alpha particles. The polarization current of alpha particles is negligible since the alpha density is much less than the background plasma density. Neglecting the displacement current and multiplying Eq. (14) by B, we obtain the following equation: We finally arrive at the magnetohydrodynamic momentum equation by approximating the left-hand side of Eq. (16)  where vE is rewritten as v. This relation is nothing but the MHD momentum equation in which the contribution of the alpha particle current is extracted from the total current. This relation is essentially the same as the model of Park et a&l4 except for -e,n,E in their force term, though their derivation is different from ours. We consider the alpha particle current as the sum of the parallel current, the magnetic drift current, and the diamagnetic current: ' .
We can easily show that the total energy is conserved in our model. The temporal evolutions of energy for MHD part and alpha particles are, respectively, given by FlG. 1. The q profile and the toroidal shear Alfven continuous spectrum of the two-mode-coupling model 'J for (n=Z, m=2) and (n=2, m=3) modes as functions of the minor radius. The (I value is an average, since the magnetic surface is not a concentric circle. The continuous spectrum is obtained assuming concentric circular magnetic flux surfaces, expanding the toroidicity effect to first order in the inverse aspect ratio u/R.

0%
We sum these equations using the Maxwell equation (dldt)B = -V X E and a vector identity, and obtain the temporal evolution of the total energy in the following conservative form:

04
In the remainder of the paper, we set the magnetic moments of alpha particles to zero. We solve these equations using a finite difference method of second-order accuracy both in time and in space.
The aspect ratio of the system is 3 and the poloidal cross section is rectangular. The cylindrical coordinate system (R,cp,z) is used. The simulation region is 2aGRG4a, --aSz sa where a is the minor radius. The simulation region in the rp direction is O=+Sr, since we focus on the n=2 TAE mode and its nonlinear evolution. We make use of (65, 20, 65, 60) grid points for (R,p,z,v~~) coordinates, respectively. The simulation region in the u 11 direction is --3 v A<~ 11<3 u A with the grid size of 0. Iv A.
The initial condition is an MHDRequilibrium where both the background plasma beta and the plasma beta of alpha particles are 2% at the magnetic axis. We obtain the initial condition using an iterative method. As an initial guess for the iteration, we numerically solved the Grad-Shafranov equation neglecting the alpha particle current, and set the parallel pressure distribution of alpha particles to be proportional to the MHD pressure with a Maxwellian distribution whose thermal velocity u a= (T,lm,) I" is equal to the Al-f&n velocity at the magnetic axis. We obtain an exact equilibrium by the iteration both for the MHD force balance and the distribution function of the alpha particles. The magnetic axis locates at R=Ro=3.16a, z=O.
(4 E'P We show, in Pig. 1, the initial 4 profile and the continuous spectrum of the two-mode-coupling mode11Y2 of (n=2, m =2) and (n =2, m =3) modes as functions of the minor radius. The q value is an average, since the magnetic surface is not a concentric circle. The continuous spectrum is obtained assuming concentric circular magnetic flux surfaces, expanding the toroidicity effect to first order in the inverse aspect ratio a/R. The q=(2m+ 1)/2n=S/4 surface is located at Y= 0.35a on average. The scale length of the alpha particle pressure gradient at the q= (2m+ 1)/2n surface, which is denoted as L, is 0.42a. The parallel Larmor radius of alpha particles, which is denoted as pm, is set to al 16, which leads to the drift frequency of alpha particles at the q=(2m+ 1)/2n surface, wP,vcY **.a--=----=1 670 2rL .
AT (25) where oA=uA/Ro and vA is hereafter the Alfven velocity at the magnetic axis. For a simulation run described in the following section, it took 100 hours of CPU time at NEC SX-3/24R, maximum computational performance of which is 12.8 GPLOPS.

A. Linear growth
We continued a simulation run up to t= 1294wi' . The TAE mode which consists mainly of (n ~2, m =2) and (n ~2, m=3) modes appeared. Shown in Fig. 2 are the toroidal electric field and the alpha particle distribution function at u([=-1.05u,, l.OSrJ,, and 2.05~~ on a poloidal cross section at t=423wi' from which the initial one is subtracted. According to the linear theory of TAE, the q value for coupling of (n =2, m =2) and (n =2, m =3) modes is q =.5/4. It is to be noted that its amplitude is peaked near the q=5/4 surface. For the distribution function, the m=3 component is dominant for u/i= -l.O5u, whereas the m = 2 component is conspicuous for VII= 1.05v,. At the q=5/4 surface the phase velocities of the m =2 and m =3 modes have the same absolute value ( =uA) with an opposite sign to each other. Thus when m =3 mode has a phase velocity of -vA , m =2 mode has that of uA. and they interact with alpha particles at each phase velocity.
In Fig. 3 we plot the temporal evolution of the (n=2, m =3) mode at r = 0.35~ where the mode structure is peaked. It is seen that the instability saturates at t=64Oo,' . The real frequency is 0.32~~ and the linear growth rate is 8.8X IO-"wA. After saturation the mode executes an amplitude oscillation.
Let us compare the simulation results with the analytical formula for the linear growth rate based on a localized analysis which is given by Fu and Van Dam2 FIG. 4. Temporal evolutions of (a) the TAE mode energy, and (bj the ratio of the power transfer rate (-j&-E} to the TAE mode energy, which is divided by a factor of 2 to relate directly to the growth rate. The decrease of this ratio leads to the saturation of the TAE instability.
This formula is applicable for alpha particles which have an isotropic Maxwellian distribution, For the distribution with only the parallel velocity component as considered here, the formula is reduced to This formula yields the growth rate of 1.10X 10-'wA for x=1/a WY: u= 1,67w,, and p,= 1.5% at the q=5f4 surface which coincides within a factor of 1.25, with the rate obtained from our simulation results.

B. Saturation level and mechanism
When the instability growth stops at t=640wi', the maximum value of the fluctuating magnetic fieId becomes 1.8X 1 Oe3 of the equilibrium field. This value is significantly large and thus a significant loss of alpha particles may be anticipated." Let us now consider the saturation mechanism. There are two possibilities for saturation mechanisms: (a) mode coupling effects in the MHD component which include generation of EXB flows and enhanced continuum damping;r3 and (b) effects of finite-amplitude waves on the alpha particles.g In order to identify the saturation mechanism we analyze the temporal evolution of the power transfer rate from alpha particles to the MHD component, namely, (-j,.E} (( > means volume integration). We show the temporal evolution of the TAE mode energy in Fig. 4(a) and the ratio of the power transfer rate to the TAE mode energy (divided by a factor of 2 to relate directly to the growth rate) in Fig. 4(b). The ratio fluctuates at the early stage of the simulation since the initial perturbation is not an eigenfunction but an arbitrary one, and it gradually converges to a constant value during the linear growth stage. At t =47Owi' the ratio begins to decrease, thus leading to saturation of the instability. It is evident that the decrease of the power transfer rate is the cause of the saturation.
After saturation both the TAE mode energy and the power transfer rate show oscillatory behaviors. Specifically, when the ratio reaches above zero at t= 1OlOoi' the mode energy begins to increase and when the ratio becomes negative at t= 118Owi' the mode energy begins to decrease. We have confirmed that the oscillation of the power transfer rate is caused by the evolution of the phase between j, and E.
It is' welI known that resonant particles in the nonlinear phase of Landau damping are trapped by the potential well and execute a bounce motion in it." The bounce frequency for the TAE mode is given by (see Ref. 15), where I'* is the minor radius of the resonance and B, is the magnetic field strength at the magnetic axis. The Fourier amplitude of the electrostatic potential #+,,J is defined by I 2~ d6' d%1,1= 0 ; qb,Jr,+Ab cos t))e "', Ap(u;+ ~u~)q(r,)lu,,fL (32) with a=e,BOlm, (alpha particle Larmor frequency). It is necessary to use a mode number 1 in addition to the poloidal mode number 112, since a particle orbit does not coincide with a magnetic surface.
The frequency of the amplitude oscillation of the TAE mode, as well as that of the power transfer rate in Fig. 4, was, is 1.9X 10-'(,jA. It should be remarked here that w, and o,, coincide with each other within a factor of 1.4. This remarkable coincidence indicates that the oscillation observed is caused by the bounce motion of the particles trapped by the TAE mode. Furthermore, the growth rate yL=8.8X lo-'w, also coincides fairly well with u+, and w,,$. This coincidence among oh, w,,, , and yL suggests that the instability saturates in the manner of the nonlinear Landau damping of a finite-amplitude wave which is as suggested by O'Neil." We can estimate the saturation level from the linear growth rate yL having the condition tib at saturation equal to hy, where h is a factor of --O(l) and for the present results h=3.0. Making use of Eq. (28) we obtain the saturation level, 4,,L,l+)2 g!!i. For the present case, A,*, which makes the major contribution to the instability, decides the saturation level of the entire TAE mode. In general, the saturation level of the entire TAE mode will depend on the mode structure and hb , but we can say, at least, that it will be in proportion to the square of the linear growth rate.

C. Distribution function at the final stage
In Fig. 5 the same quantities are plotted as in Fig. 2, but for t = 12940,' . Comparing it with Fig. 2, the distribution function shows a very different structure, while the toroidal electric field is still similar. The (n=O, m=O) quasilinear modes are generated through the nonlinear coupling of the n = 2 TAE mode and the n =2 modes of the distribution function, and they have become dominant at t = 1294~;' . The quasilinear modes are the manifestations of the transport of alpha particles by the TAE modes which smoothes the spatial gradient of alpha particle distribution. Though the amplitude of the quasilinear mode at t = 1294wi' is a few percent of the initial distribution at u11=0 at the magnetic axis, it affects the spatial gradient significantly. Figure 6 shows distributions of alpha particles at UI~= -1 .05uA and z=O as a function of R which are averaged in the cp direction at t =0 and f= 1294~2' , respectively. The spatial gradient of the distribution function is reduced to half near R = 2.8 a and 3.3 a.

iv. CONCLUSION
A new simulation method has been developed to investigate the excitation and saturation processes of the TAE modes. In this method, the background plasma is described by the full-MHD model, while the kinetic evolution of energetic alpha particles is represented by the drift kinetic equation. It is demonstrated that the n =2 TAE mode is excited, the linear growth rate of which is in good agreement with the linear theory of Fu and Van Dam.* After saturation the mode amplitude shows an oscillatory behavior with the bounce frequency of alpha particles trapped by the TAE mode. Moreover, the bounce frequency is in agreement with the linear growth rate. The saturation is caused by the decrease of the power transfer rate from alpha particles to the TAE mode. Thus, we conclude that the growth of the unstable mode is suppressed by the trapped particle effect of a finite-amplitude wave.
The saturation of the magnetic field fluctuation can reach to a significant level, e.g., 1.8X low3 of the equilibrium field intensity when the initial beta of alpha particles is 2% at the magnetic axis, which is supposed to lead to a non-negligible alpha particle loss in one slowing time.* We can estimate the saturation level for each 4m,c from the linear growth rate yL using the condition that @b at saturation equals X yL where X is a factor of -c2( I). The saturation level of the component which makes the major contribution to the instability will decide the saturation level of the entire TAE mode. It will be in proportion to the square of the linear growth rate.
The (n =O, m =O) quasilinear mode of the alpha particle distribution is generated through the nonlinear coupling of the n=2 TAE mode and the n =2 mode of alpha particle distribution. This quasilinear mode spatially flattens the distribution function, reducing the free energy source of the instability.
This new simulation technique will be useful to study other fast particle physics such as the fishbone and sawtooth oscillations.