Multi-phase hybrid simulation of energetic particle driven magnetohydrodynamic instabilities in tokamak plasmas

Magnetohydrodynamic (MHD) instabilities driven by energetic particles in tokamak plasmas and the energetic particle distribution formed with the instabilities, neutral beam injection, and collisions are investigated with hybrid simulations for energetic particles and an MHD fluid. The multi-phase simulation, which is a combination of classical simulation and hybrid simulation, is applied to examine the distribution formation process in the collisional slowing-down time scale of energetic ions for various beam deposition power ( P NBI ) and slowing-down time ( τ s ). The physical parameters other than P NBI and τ s are similar to those of a Tokamak Fusion Test Reactor (TFTR) experiment (Wong et al 1991 Phys. Rev. Lett. 66 1874). For P NBI = 10 MW and τ s = 100 ms, which is similar to the TFTR experiment, the bursts of toroidal Alfvén eigenmodes take place with a time interval 2 ms, which is close to that observed in the experiment. The maximum radial velocity amplitude (vr) of the dominant TAE at the bursts in the simulation is v r / v A ∼ 3 × 10 − 3 where vA is the Alfvén velocity at the plasma center. For P NBI = 5 MW and τ s = 20 ms, the amplitude of the dominant TAE is kept at a constant level v r / v A ∼ 4 × 10 − 4 . The intermittency of TAE rises with increasing P NBI and increasing τ s (= decreasing collision frequency). With increasing volume-averaged classical energetic ion pressure, which is well proportional to P NBI τ s , the energetic ion confinement degrades monotonically due to the transport by the instabilities. The volume-averaged energetic ion pressure depends only on the volume-averaged classical energetic ion pressure, not independently on P NBI or τ s . The energetic ion pressure profile resiliency, where the increase in energetic ion pressure profile is saturated, is found for the cases with the highest P NBI τ s where the TAE bursts take place.


Introduction
Nuclear fusion is a safe and environmentally friendly energy source in the next generation. Fusion reactors based on magnetically confined plasmas will harness the fusion reaction of deuterium and tritium in high temperature plasmas. Alpha particles born from the deuterium and tritium reaction with birth energy 3.5 MeV are expected to heat the plasma to maintain the high temperature. This type of high temperature plasma that is self-sustained by the fusion reaction is called 'burning plasma'. Magnetohydrodynamics (MHD) is a theoretical framework of plasmas that combines the fluid equations, the induction equation for the magnetic field, and the Ohm's law for the electric field. MHD explains well the macroscopic behaviors of fusion plasmas. Three types of wave exist in MHD, namely, slow and fast magnetosonic waves, and shear Alfvén wave. Shear Alfvén wave is an incompressible wave whose restoring force is the magnetic tension, while the slow and fast magnetosonic waves are compressional waves with the variations of the magnetic pressure and the plasma pressure. Alfvén eigenmodes (AEs) are MHD oscillations of the nonuniform magnetically confined plasmas. The spatial profile of AE peaks at the extremum of the Alfvén wave continuous frequency spectrum. AEs are classified on the basis of the extremum. For example, global Alfvén eigenmode [1] and reversed shear Alfvén eigenmode [2] are located at the extrema of the safety factor which is an index of the torsion of the magnetic field, and toroidal Alfvén eigenmode (TAE) [3,4] and beta-induced Alfvén eigenmode [5] appear at the gaps of the Alfvén continuous spectra, which are created by toroidicity and finite plasma pressure, respectively. The speed of alpha particle with energy 3.5 MeV exceeds the phase velocities of Alfvén wave and slow magnetosonic wave. The alpha particles can resonate with AEs in the collisional slowing-down process, and may destabilize and amplify the AEs. The alpha particle transport by the amplified AEs flattens the alpha particle spatial profile and leads to alpha particle losses. This will reduce the alpha particle heating efficiency and deteriorate the fusion reactor performance. Alpha particle driven AEs are one of the major concerns of burning plasmas. The interactions between energetic particles and AEs have been extensively studied using energetic ions generated by the neutral beam injection (NBI) and ion-cyclotron-range-of-frequency (ICRF) wave heating in tokamak and stellarator/heliotron plasmas [6][7][8][9][10].
Hybrid simulations of energetic particles and MHD are useful tools to analyze both the AEs and the energetic-particle transport [11][12][13][14][15][16][17]. In most of the hybrid simulations, the initial energetic particle distribution is assumed, and the evolution of AEs and energetic particles is simulated as an initial value problem for a short time of the growth and saturation of the AEs, which is typically ∼1 ms. However, we should notice that the energetic particle distributions in the experiments are formed in a longer time scale of the collisional slowingdown time for energetic particles interacting with the AEs. Then, comprehensive simulations that include the energetic particle source (birth of alphas, NBI, ICRF), collisions of energetic particles with the bulk plasma particles, energetic-particle losses, and the interactions between energetic particles and AEs are needed for the understanding of the energetic particle distribution in the experiments and the prediction of future experiments.
The author and his collaborators have tried to develop such comprehensive simulations of energetic particles and AEs with the energetic ion source, collisions, and losses [18][19][20][21][22][23][24][25][26]. We performed the first numerical demonstration of TAE bursts with parameters similar to a TFTR experiment [27] and reproduced many of the experimental characteristics [20]. However, the saturation amplitude of the magnetic field fluctuation normalized by the toroidal field was d~´-B B 2 10 2 , which is higher by one order of magnitude than the value d~-B B 10 3 inferred from the experimental plasma displacement measurements [20,28]. In the simulation of [20], the only nonlinearity retained was the nonlinearity in the energetic-particle orbits, while the nonlinear MHD effects were neglected. The nonlinear MHD effects on the evolution of a TAE was carefully investigated, and it was found that the energy transfer from the TAE to the nonlinearly generated zonal modes with toroidal mode number n=0 and higher harmonics reduces the saturation level of the TAE [29,30]. A hybrid simulation with both the MHD nonlinearity and the energetic ion source, collisions, and losses was constructed using the df particle simulation method for energetic particles to reduce the numerical noise [22]. The TAE bursts were reproduced with the df simulation, but the energetic ion velocity distribution was purely parallel to the magnetic field and the pitch-angle scattering was neglected because of the limitation of the df simulation model developed in [22].
The multi-phase simulation, which is a combination of classical simulation and hybrid simulation for energetic particles interacting with an MHD fluid, was developed in order to investigate a energetic ion distribution formation process with beam injection, collisions, losses, and transport due to the AEs with the MHD nonlinearity [24,25]. We use the MEGA code [15] for both the classical and hybrid simulations. We run alternately the classical simulation without MHD perturbations and the hybrid simulation with MHD perturbations. The multi-phase simulation is a comprehensive simulation, which deals with both the AEs and the energetic ion transport as self-consistently and realistically as possible, yet attainable on a tractable timescale. It was demonstrated with the multi-phase simulation of DIII-D discharge #142111 [31,32] that the energetic ion spatial profile is significantly flattened due to the interaction with the multiple AEs and that the energetic ion pressure profile is in agreement with that of the experiment with the root-mean-square of the deviations same as the error bar [25]. The predicted temperature fluctuation profiles of = n 3, 4, and 5 modes were quantitatively compared with ECE measurements, and it was found that the fluctuation profiles as well as phase profiles are in very good agreement with the measurements. Additionally, the saturated amplitudes are within a factor of 2 of those measured. The nonlinear MHD effects [22,29,30,33,34] that prevent the AE amplitude from increasing to a large amplitude observed in a reduced simulation [20] are included in the hybrid simulation. The energetic ion profile stiffness that was observed in the DIII-D experiments [35] was also investigated with the multi-phase simulation, and it was demonstrated that the resonance overlap of multiple eigenmodes [36][37][38] accounts for the profile stiffness with the sudden increase in energetic ion transport with increasing beam power [26].
In the present paper, we apply the multi-phase simulation to tokamak plasmas for various beam deposition power (P NBI ) and slowing-down time (t s ) with the other parameters similar to the TFTR experiment [27]. For the same beam deposition power ( = P 10 NBI MW) and the slowing-down time (t = 100 s ms) as those used in [20,22], we see that the TAE bursts take place with time intervals 2 ms, which is close to that in the experiment, and the maximum amplitude~´- . We examine the volume-averaged energetic ion pressure and the energetic ion profile for various P NBI and t s . We see the monotonic degradation of the energetic ion confinement and the profile resiliency, where the energetic ion pressure profiles are saturated with increasing volume-averaged classical energetic ion pressure which is well proportional to t P NBI s .

Simulation model
2.1. Hybrid simulation model for energetic particles and MHD We use the MEGA code [15], in which the bulk plasma is described by the nonlinear MHD equations and the energetic ions are simulated with the gyrokinetic particle-in-cell (PIC) method. Several hybrid simulation models have been constructed [11][12][13][14][15][16][17] to study the evolution of AEs destabilized by energetic particles. An extended MHD model given in [39] has been implemented together with the equilibrium toroidal flow in MEGA [25,26]. In this paper, we use the standard MHD equations with the energetic ion effects: where m 0 is the vacuum magnetic permeability, g = 5 3is the adiabatic constant, ν, n n and χ are artificial viscosity and diffusion coefficients chosen to maintain numerical stability. In this work, the dissipation coefficients ν, n n , χ, and h m 0 are assumed to be equal to each other. The dissipation terms play a physical role to enhance the damping of AEs in the MHD simulation that does not include kinetic damping such as radiative damping [40] and thermal ion Landau damping. In this paper, we use one value of the coefficients,´-1 10 6 normalized by v R A 0 where v A is the Alfvén velocity at the plasma center, and R 0 is the major radius at the geometrical center of the simulation domain. The subscript 'eq' represents the equilibrium variables. The MHD momentum equation (equation (2)) includes the energetic ion contribution in the energetic ion current density ¢ j h that consists of the contributions from parallel velocity, magnetic curvature and gradient drifts, and magnetization current. TheÉ B drift disappears in ¢ j h due to the quasi-neutrality [15]. We see that the electromagnetic field is given by the standard MHD description. This model is accurate under the condition that the energetic ion density is much less than the bulk plasma density. The MHD equations are solved using a fourth order (in both space and time) finite difference scheme.
The energetic ions are simulated using the full-f PIC method, and a guiding-center approximation [41], where we employ the gyrokinetic approach to account for finite Larmor radius (FLR) effects. The electromagnetic fluctuations are averaged over the energetic ion gyro orbit for the energetic ion dynamics. It was demonstrated that the MEGA code with the full-f PIC method can be applied to energetic particle modes (EPMs) in JT-60U, although the numerical noise level is higher in the full-f PIC simulation than in the df PIC simulation [23,42]. The equations of motion for each computational particle are solved using a fourth-order Runge-Kutta method. MEGA code participated in the code benchmark of the Energetic Particle Physics Topical Group of the International Tokamak Physics Activity [43]. Good agreements were found in the spatial profile, frequency, and growth rate of a TAE among the seven codes without the energetic ion FLR effects and among the six codes with the FLR effects.

Equilibrium, beam injection, collisions and losses
We employ the tokamak equilibrium used in [22]. Cylindrical coordinates j ( ) R z , , are used in the simulations. For the purpose of the data analysis, magnetic flux coordinates j J ( ) r, , were constructed for the MHD equilibrium where r is the radial coordinate with r=0 at the plasma center and r=a at the plasma edge, and ϑ is the poloidal angle. The physical parameters are similar to the TFTR experiment [27]; a=0.75 m, , and the beam injection energy is 110 keV. Both the bulk and beam ions are deuterium. The bulk ion density is2.8 10 19 m −3 , and the bulk plasma beta value, which is the ratio of bulk plasma pressure to the magnetic pressure at the plasma center, is 1%. The beam injection velocity corresponds to = v v 1.1 b A . The safety factor profile is assumed to be the same as in [20,22] , and l D = 0.3. The numbers of grid points are (128, 128, 128) for j ( ) R z , , coordinates, respectively, and the number of computational particles is4.2 10 6 . It was confirmed in a reduced simulation of bursting evolution of five AEs with toroidal mode number n=1-5 that 2 million particles are sufficient for numerical convergence in burst interval, modulation depth of the stored fast ion energy at each burst, and saturation level of the stored fast ion energy [20]. Good convergence was also found with the similar number of particles per cell for the AEs and the energetic particle redistributions in an ITER plasma [44]. The reason why the relatively small number of particles works well may be that the phase space of energetic particles can be efficiently covered by the computational particles when we apply PIC only to to the energetic particles and not to thermal particles. Another reason for the hybrid simulation with full MHD is that the time step width is limited short by the Courant condition for fast magnetosonic wave. This may reduce the effective numerical noise for AEs because the energetic particle distribution is computed more frequently.
Collisions of energetic ions with thermal electrons and ions [45] are implemented in the MEGA code. The slowing-down is included in the time integration of total particle velocity v by where n s is the drag rate (inverse of the slowing-down time), and v c is the critical velocity above which the energetic ion collisions with electrons dominate the slowing-down process. The critical velocity is assumed to be = v v 0.65 c A . Pitch-angle scattering and energy diffusion are taken into account at the end of each time step using a Monte Carlo procedure [46], where a particle's pitch λ and total velocity v is altered according to the relations where n d is the pitch-angle scattering rate, m h is the energetic ion mass, Dt is the time step width, T e , T i are respectively electron and ion temperature, and±denotes a randomly chosen sign with equal probability for plus and minus. However, the energy diffusion given by equation (9) is turned off in this paper, and the pitchangle scattering rate is set to be n n = v v 2 d s c 3 3 in order to make the simulation condition same as the previous works [20,22].

Multi-phase simulation
We would like to investigate a energetic ion distribution formation process with beam injection, collisions, losses, and transport due to the AEs. A complicating factor is that the time scale of the classical processes without MHD perturbations is the slowing-down time, which is roughly ∼100 ms, and longer by four orders of magnitude than the typical oscillation period of AEs ∼0.01 ms. The time step width is limited by the Courant condition for fast magnetosonic waves in the hybrid simulation. On the other hand, in the classical simulation where the MHD part of the simulation is turned off, the time step width can be taken to be greater by one order of magnitude than in the hybrid simulation. To deal with this efficiently, a multi-phase simulation, where the classical simulation and the hybrid simulation are run alternately, was constructed [24]. In the classical phase of the simulation, the energetic ion distribution is built up with the beam injection and collisions. In the subsequent hybrid phase, the built-up energetic ion distribution destabilizes AEs leading to the relaxation of the distribution. We should note that the classical processes, beam injection and collisions, take place also in the hybrid phase. We repeat this combination of the classical and hybrid simulations until the stored energetic ion energy is saturated. After which, we continuously run the hybrid simulation to investigate the energetic ion distribution and AEs. In our previous work [24], multi-phase simulations with different classical phase duration (4 ms and 9 ms) were performed, and a good agreement was found between the two runs.

Simulation results
We have investigated ten cases for different beam deposition power (P NBI ) and slowing-down time (t s ). The multi-phase simulations are run with alternating classical phase for 2 ms and hybrid phase for 0.5 ms until stored energetic ion energy is saturated, after which the hybrid simulation is run continuously. The time step widths for the classical phase and the hybrid phase are 20 ns and 2 ns, respectively. Computational particles are injected at a constant rate over a beam injection period = t 100 inj ms. Table 1 summarizes the beam deposition power, the slowing-down time, and the starting time (t s ) and the ending time (t e ) of the last hybrid simulation for all the cases. We present the time evolution of stored energetic ion energy and MHD kinetic energy, the frequency spectra of temperature fluctuations, the spatial profiles of AEs, and the amplitude evolution of the dominant AEs for four cases A, D, F, and I.    The hybrid simulation is continuously run after t=47 ms. We see in figure 4(a) that the stored energetic ion energy is saturated, and in figure 4(b) that the bursty amplification of the MHD kinetic energy takes place repetitively after t=50 ms. Figure 5 shows the frequency spectra of the temperature fluctuations at   with the single dominant poloidal harmonic m=2 and the rapid variation of the sine part at   r a 0.5 0.6. This indicates the rapid phase variation of the spatial profile, which corresponds to the continuum damping [48,49]. All the other modes are TAEs.   figure 7. The hybrid simulation is continuously run after t=20 ms. We see in figure 7(a) that the stored energetic ion energy is saturated, and in figure 7(b) that the bursty amplification of MHD kinetic energy takes place repetitively after t=20 ms. Figure 8 shows the frequency spectra of the temperature fluctuations at   [20,22]. The time evolutions of stored energetic ion energy and MHD kinetic energy are shown in figure 10. We see in figure 10(a) that the stored energetic ion energy is saturated, and in figure 10(b) that the bursts take place recurrently with time intervals about 2 ms, which is close to that observed in the TFTR experiment [27]. The neutron rate from the beam and bulk deuterium reaction is evaluated before and after the burst att 26 ms using the fusion cross-section [50]. The drop in neutron rate is 9%, which is also close to the measurement at the TFTR experiment. Figure 11 shows    Figure 13 shows the time evolution of the dominant radial velocity fluctuation harmonics for cases A, D, F, and I at each peak location. We see in figure 13(a) that the mode amplitude stays around a steady level´v v 4 10 r A 4 . For the longer slowing-down time (t = 100 s ms) in case D, we see in figure 13(b) that the evolution becomes bursty. The classical stored energetic ion energy is well proportional to the product of the beam deposition power and the slowing-down time ( t = P NBI s ). The slowing-down time in case D makes the the classical stored energetic ion energy five times larger in case D than in case A. The higher energetic ion pressure may lead to the bursty evolution of AEs. We see a bursty evolution also for case F shown in figure 13(c). We see repetitive bursts in figure 13(d) for case I, where the classical stored energetic ion energy is the largest. The bursts take place with a time interval ∼2 ms. The maximum amplitude at each burst is~´v v 3 10 r A 3 . The time interval between the busts is close to that observed at the TFTR experiment [27], and the maximum amplitude of radial velocity fluctuation is close to the value that was inferred for the TFTR experiment in [20].

Amplitude evolution
We have investigated the time-averaged amplitude and the standard deviation of amplitude for the dominant n=2 radial velocity fluctuation harmonics for all the cases. Figure 14(a) shows time-averaged amplitude versus beam deposition power P NBI . We see in the figure a linearly proportional increase of the timeaveraged amplitude to the beam deposition power, which is similar to the observations on DIII-D [35]. Figure 14(b) shows the standard deviation of amplitude versus volume-averaged classical energetic ion beta, which is the ratio of the energetic ion pressure to the magnetic pressure at the plasma center. The standard deviation increases with the volume-averaged classical energetic ion beta for low energetic ion beta, but it is . The ratio of the standard deviation to the time-averaged amplitude is plotted versus slowing-down time in figure 14(c). This ratio is an index of intermittency of the AEs. For the cases with P NBI =5 MW represented by squares in figure 14(c), we see the intermittency rises with increasing slowingdown time (=decreasing collision frequency). For the cases with t s =20 ms, we see a trend that the intermittency is enhanced with increasing beam deposition power (=increasing drive to AEs). Two types of evolution, a steady state or pulsations were predicted depending on the relationships between the source, the  mode damping rate, and the collisions [51]. The steady state was predicted for larger collision frequency and lower drive to AEs. The simulation results shown in figure 14(c) are qualitatively consistent with the theoretical prediction. It has been predicted that the effect of micro-turbulence-induced radial diffusion of energetic particles is similar to and greater than that of the pitch-angle scattering on the AE evolution [52]. Quantitative comparisons between simulations and recent experiments will promote a better understanding of the evolution of AEs interacting with energetic particles.

Energetic ion profile resiliency
Volume-averaged energetic ion beta b á ñ h in the simulation results is plotted versus volume-averaged classical energetic ion beta b á ñ h classical in figure 15(a). The line in the figure represents b á ñ h = b á ñ h classical . We see that the reduction of b á ñ h from b á ñ h classical increases monotonically for higher b á ñ h classical , which indicates the monotonic degradation of energetic ion confinement with increasing b á ñ h classical . We see a reduction of energetic ion beta by . It is interesting to note that the data points look located on one curve, while the data points represent the results for various beam deposition power P NBI and slowing-down time t s . This indicates that the most important parameter that determines the energetic ion confinement is b á ñ h classical rather than the independent parameters, P NBI or t s . Volume-averaged classical energetic ion beta b á ñ h classical is well proportional to the product of the two parameters t P NBI s . The ratio of the neutron rate to the classical neutron rate is also plotted figure 15(b). We see the monotonic degradation in the neutron rate with increasing b á ñ h classical . The neutron rate is about 35% of the classical neutron rate for b á ñ~0.015 h classical with P NBI =10 MW and t s =100 ms. This is close to the results of [20] where the stored energetic ion energy is 40% of the classical simulation result. The energetic ion beta profiles are compared in figure 16 for the cases D, H, and I with the highest b á ñ h classical . We see that the energetic ion beta profiles are very similar to each other, and are saturated with the increase of b á ñ h classical by factor of 2 from case D to case I. This saturation of the energetic ion pressure profile with increasing beam power is called profile resiliency [35], and the energetic ion profile is regarded as 'stiff' when the energetic ion transport suddenly increases for the beam power or the energetic ion profile gradient above critical values. In our simulations based on a DIII-D equilibrium [26], it was demonstrated that the profile stiffness is brought about by resonance overlap of multiple AEs. The profile resiliency was not observed in the simulations presented in [26]. It was discussed that the profile resiliency would take place if the phase space was fully covered by the resonance overlap from the plasma center to the edge. We would like to point out that the recurrent AE bursts occur with a constant time interval for case I as is shown in figure 13. This type of AE bursts was not observed in the simulations presented in [26]. The phase space may be fully covered by the resonance overlap when the AE bursts take place as was demonstrated with the simulations of AE bursts in [20].

Summary
We performed the multi-phase hybrid simulations for energetic particles and an MHD fluid in order to investigate MHD instabilities driven by energetic particles in tokamak plasmas and the energetic particle distribution formed with the instabilities, neutral beam injection, and collisions. The multi-phase simulation, which is a combination of classical simulation and hybrid simulation, was applied to examine the distribution formation process in the collisional slowing-down time scale of energetic ions for various beam deposition power (P NBI ) and slowing-down time (t s ). The physical parameters other than P NBI and t s are similar to those of a TFTR experiment [27]. For P NBI =10 MW and t s =100 ms, which is similar to the TFTR experiment, the TAE bursts take place with a time interval 2 ms, which is close to that observed in the experiment. The maximum radial velocity amplitude of the dominant TAE at the bursts in the simulation is~´v v 3 10 r A 3 . For P NBI =5 MW and t s =20 ms, the amplitude of the dominant TAE is kept at a constant level´v v 4 10 r A 4 . It was found that the intermittency of TAE rises with increasing P NBI and increasing t s (=decreasing collision frequency). With increasing volume-averaged classical energetic ion pressure, which is well proportional to t P NBI s , the energetic ion confinement degrades monotonically due to the transport by the instabilities. The volume-averaged energetic ion pressure depends only on the volume-averaged classical is well proportional to the product of the two parameters t P NBI s . We found the energetic ion pressure profile resiliency, where the increase in energetic ion pressure profile is saturated, for the cases with the highest t P NBI s where the TAE bursts take place.
Since the employed MHD model lacks a complete physical description that includes kinetic damping of AEs, we employ artificial dissipation coefficients to control the damping. In this paper, we used one value of the coefficients´-1 10 6 . We investigated the dependence on dissipation coefficients in our previous works [22,24]. We found that the AE evolution transitions from a steady state to recurrent bursts with increasing dissipation coefficients. The time interval between the AE bursts increases with increasing dissipation coefficients. The energetic ion pressure and stored energy are higher for the larger diffusion coefficients which lead to higher damping rate of the AEs. We discussed in [25] that we may need a more sophisticated simulation model that handles kinetic damping for a perfect reproduction of all the AEs observed in the experiment, on the other hand, the energetic ion profile may not be sensitive to the individual damping rate of each mode in multimode systems.
The multi-phase hybrid simulation presented in this paper has been applied to DIII-D [24][25][26] and JT-60U [53]. For the DIII-D experiments, the energetic ion pressure profile and the electron temperature fluctuations were reproduced with the simulations, and the energetic ion profile stiffness was also demonstrated. For the JT-60 experiments, the repetitive frequency chirping of EPM has been successfully simulated. The multi-phase hybrid simulation is now being applied to the LHD experiments where AE bursts were observed [54]. The simulations validated on the basis of the experiments will enable us to predict the behaviors of AEs and the distributions of energetic particles in burning plasmas.