Introduction to the interaction between energetic particles and Alfven eigenmodes in toroidal plasmas

This article is a tutorial review of the interaction between energetic particles and Alfvén eigenmodes (AEs) which is one of the important research issues for fusion burning plasmas. The destabilization mechanism of AEs is a kind of inverse Landau damping through the resonant interaction with energetic particles. The important properties of the AE instability, such as resonance condition, conserved variable during the interaction, and particle trapping by the AE, are explained. The time evolution of AEs is classified into various types, steady state, frequency splitting, frequency chirping, and recurrent bursts. Berk and Breizman presented both a one-dimensional weakly nonlinear theory for marginal stability and a reduced simulation model that qualitatively explain the various types of time evolution. Berk–Breizman’s theory and reduced simulation model are introduced, and their limitations and the future works are discussed in this article. In addition, energetic particle transport by AEs is illustrated with surface-of-section plots. The particle trapping by the AE creates phase space islands and leads to the local flattening of the energetic particle spatial profile. The resonance overlap of multiple AEs and the overlap of higher-order resonances of a single AE lead to the emergence of stochasticity in phase space and the global transport of energetic particles.


Introduction
Alpha particles born from the deuterium-tritium (D-T) reaction with energy 3.5 MeV are expected to heat the fusion burning plasmas to maintain the high temperature that is needed for the D-T reaction. Magnetohydrodynamics (MHD) is a theoretical framework of plasmas that combines the fluid equations, the induction 1 Page 2 of 33 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, shear Alfvén wave, slow magnetosonic wave, and fast magnetosonic 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 variations of the magnetic pressure and the plasma pressure.
Alfvén eigenmodes (AEs) are magnetohydrodynamic oscillations of the nonuniform magnetically confined plasmas. The spatial profile of AE peaks at the extremum of the shear Alfvén wave continuous frequency spectrum. AEs are classified on the basis of the extremum. For example, global Alfvén eigenmode (GAE) (Appert et al. 1982) and reversed shear Alfvén eigenmode (RSAE) (Fukuyama and Akutsu 2002) 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) (Cheng et al. 1985;Cheng and Chance 1986) and beta-induced Alfvén eigenmode (BAE) (Heidbrink et al. 1993) 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 shear Alfvén wave and magnetosonic waves. 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 (Fasoli et al. 2007;Heidbrink 2008;Toi et al. 2011;Sharapov et al. 2013;Gorelenkov et al. 2014). When the energetic particle drive is so strong that the drive overcomes the continuum damping (Rosenbluth et al. 1992;Zonca and Chen 1992), energetic particle mode (EPM) can be destabilized (Chen 1994). The spatial peak of EPM is located on the shear Alfvén continuous spectrum.
In this article, we give a tutorial review of the interaction between energetic particles and AEs. The basics of the interaction is explained in Sect. 2. Section 3 is devoted to the Berk-Breizman's theory. Berk and Breizman presented a one-dimensional weakly nonlinear theory for marginal bump-on-tail instabilities. It has been demonstrated that the theory can be applied to the AE instabilities. Berk and Breizman also constructed a reduced simulation model for the nonlinear evolution of the bump-on-tail instability. The reduced simulation model can be applied to AE instabilities. An elementary derivation of the reduced simulation model and the simulation results are presented in section 4. The AEs show various types of time evolution, steady state, frequency splitting, frequency chirping, and recurrent bursts. Berk-Breizman's theory and reduced simulation model qualitatively explain the various types of time evolution of AEs. Section 5 is devoted to the energetic particle transport by AEs. Energetic particle transport is illustrated with surface-of-section plots. The resonance overlap of multiple AEs and the overlap of higher-order resonances leads to the emergence of stochasticity in phase space and the global transport of energetic particles. The limitations of the reduced simulation model and the future works are discussed in Sect. 6.

Interaction with Alfvén eigenmodes
The parallel electric field to the magnetic field is zero in the ideal MHD. The energy transfer between charged particles and Alfvén eigenmodes (AEs) take place through the magnetic (grad-B and curvature) drifts and the perpendicular electric field. For the interaction with most of the AEs, the magnetic moment is an adiabatic invariant because the AE frequency is sufficiently lower than the ion Larmor frequency. The energetic particle transport is brought about by the E × B drift ( E ) and the perturbative magnetic field. For parallel particle velocity v ∥ , the transport velocity by the perturbative magnetic field ( ) is given by v ∥ ∕B . The ratio of E × B velocity to the Alfvén velocity ( v A ) is comparable to the ratio of the perturbative magnetic field to the equilibrium field for shear Alfvén wave ( v E ∕v A ∼ B∕B ). Since the parallel particle velocity is comparable to the Alfvén velocity v ∥ ∼ v A for the resonant particles, the transport by the E × B drift is comparable to that by the magnetic per- . Both the effects of the electric field and the perturbative magnetic field should be considered for the energetic particle transport by AEs.

Resonance condition
Alfvén eigenmodes (AEs) are destabilized by energetic particles through resonant interactions. The resonance condition is given by where and n are the AE frequency and the toroidal mode number, and l is an arbitrary integer. The poloidal and toroidal orbit frequencies are defined by = 2 ∕T and by = ∕T , respectively, where T and are the time for each particle to complete one round in the poloidal angle and the toroidal angle which the particle passes in T , respectively. With these definitions, we can find that the resonance condition Eq. (1) is equivalent to which indicates that when a resonant particle passes one round in the poloidal angle, the phase of the AE at the location of the particle should change by a multiple of 2 (Todo and Sato 1998). A generalized resonance condition for stellarator-heliotron plasmas is given in Kolesnichenko et al. (2002).
The resonance condition Eqs. (1) and (2) can be generalized to the higher-order resonance where the phase of the AE is the same when the particle passes K (K > 1) times around the poloidal angle: The higher-order resonance is also called nonlinear resonance or fractional resonance. The higher-order resonance is negligible for the linear stability analysis, but may transfer substantial energy between finite-amplitude waves and particles.

Conservation of E'
In axisymmetric time-independent fields, energy E and toroidal momentum P = e h + m h Rv are conserved along the particle orbit. Here, is the poloidal magnetic flux, and e h and m h are charge and mass of the particle. In the presence of a wave with angular frequency and toroidal mode number n, neither energy nor toroidal momentum is conserved whereas magnetic moment is an adiabatic invariant if << h ≡ e h B∕m h is satisfied. However, a combination of energy and toroidal momentum E � = E − P ∕n is conserved . The reason is explained as follows. The energy and toroidal momentum evolution with the equilibrium field Hamiltonian H 0 and the wave field Hamiltonian H 1 are given by because H 0 t = H 0 = 0 holds. Supposing that the wave amplitude is constant, the wave field Hamiltonian can be written as H 1 =Ĥ 1 (R, z)e in −i t in cylindrical coordinates (R, , z) . The energy and momentum evolution is given by Then, is satisfied. For the wave-particle interaction in tokamak plasmas, the conservation of E ′ relates the variations in particle energy ( E ) and poloidal flux ( ) As the variation in poloidal flux corresponds to the radial transport, the radial transport is related to the energy transfer between the wave and the particle. It is qualitatively suggested that the primary role of the waves with low and high n such as the ion temperature gradient (ITG) modes is the radial transport rather than the energy transfer while the waves with high , such as the ion cyclotron range of frequency (ICRF) and electron cyclotron heating (ECH) waves, work for heating rather than for transport. One may ask the question why AEs can be destabilized by the resonant interaction with energetic particles although the slowing-down distribution (and also the Maxwell distribution) has negative gradient in energy f E < 0 . In spatially uniform plasmas, the distribution with f E < 0 leads to the Landau damping of the wave. However, in toroidal plasmas, the spatial gradient of the energetic particle distribution destabilizes the AEs. When we take the energy derivative of the distribution function, we should keep E � = E − P ∕n constant. Then, the energy derivative is 1 Page 6 of 33 direction [ sgn( f ∕ r) = −1 ], and the plasma current is + direction [ sgn(B ) = −1 ], the triple product =1 means the AE propagates to + direction.
When both the poloidal and the toroidal mode numbers m and n are given for the AE, the poloidal propagation direction is determined from the toroidal propagation direction through the sign of m/n. For low frequency AEs such as TAE and RSAE, the parallel wave number is low |k ∥ | = |(mB ∕r + nB ∕R)∕B| << |n∕R| , and the signs of mB and nB are opposite to each other. Then, sgn(m∕ ) = −sgn(n∕ ) ⋅ sgn(B ∕B ) gives the poloidal propagation direction. When sgn(e h ) ⋅ sgn( f ∕ r) ⋅ sgn(B ) = 1 , the AE propagates to − direction. The result is often summarized as a simple conclusion: "AEs destabilized by energetic particle spatial gradient propagate to the energetic particle diamagnetic drift direction." Figure 1 shows fast ion distribution functions for co-going and counter-going particles to the plasma current in a tokamak plasma. The major and minor radii are 1.8 and 0.6 m, respectively. The magnetic field strength is 2T at the plasma center, and the safety factor is 1.1 at the center and 3.0 at the edge. The fast ion species is deuterium and the energy distribution is a slowing down distribution with the critical energy 30 keV. The maximum energy of the distribution is 80 keV, and the spatial distribution is exponential in P with the scale length 0.4 0 where 0 is the poloidal magnetic flux at the plasma center, and = 0 at the plasma edge. The magnetic moment is constant and chosen so that |v ∥ ∕v| = 0.7 for E = 80 keV at the plasma center. The white lines in the figure represent E � = const. for a wave with n = 4 and frequency 70 kHz. The distribution functions along the E � = const. lines are plotted versus E in panel (c).We see in the figure that d f ∕ d E < 0 along the E � = const. lines. This demonstrates the destabilization mechanism of AEs through the spatial gradient of energetic particles.
In Fig. 1a, b, white lines for E � = const. are drawn from the plasma center and E = 80 keV to the plasma edge. We see in Fig. 1c that the energies at the plasma edge are 60 and 69 keV for co-going and counter-going particles, respectively. This indicates that the transport of counter-going particle from the plasma center to the edge needs less energy transfer than that of co-going particle. This property is further enhanced by the loss region that appears in Fig. 1b. This region represents trapped particles whose orbit width is large enough to reach the plasma edge. Counter-going particles with energy lower than about 70 keV at the plasma center and the same magnetic moment are lost by entering the loss region across the passingtrapped boundary. We can say that counter-going particles are easier to be lost than co-going particles.
Alfvén eigenmodes with toroidal mode number n = 0 are not destabilized by energetic particle spatial gradient. However, when the energetic particle distribution function is not isotropic in velocity space and depends on pitch angle variable = B 0 ∕E , the energy derivative is given by The second term on the right hand side can lead to destabilization of n = 0 modes such as the energetic-particle driven geodesic acoustic mode [EGAM (Fu 2008)].

Particle trapping
The most important nonlinear process of the (inverse) Landau damping is "particle trapping" where the resonant particles are trapped by the finite-amplitude wave.
After the particle is trapped by the wave, the particle executes a bounce motion in the wave potential well and does not transfer net energy to the wave in a time longer than the bounce period. This leads to the saturation of the (inverse) Landau electrostatic wave. The importance of particle trapping for the saturation of AE was predicted and the bounce frequency for Alfvén eigenmode in tokamak plasma is given by , . It is well known that the phase space with a finite amplitude wave is divided into the "passing (open)" and "trapped (closed)" regions with the separatrix. It is easy to draw a figure of the phase space structure for one-dimensional problems. However, also for the waves and the particles in toroidal plasmas, we can analyze the phase space structure utilizing the surface-of-section plots where only one eigenmode is taken into account and the amplitude and the frequency of the wave are kept constant. We choose particles which have a constant value of E � = E − P ∕n , because E ′ is conserved in the interaction with a constant-amplitude wave with frequency and toroidal mode number n. Then with a single mode we have a conserved variable and we can make surface of section plots that are easily interpretable. Otherwise, with more than one wave we would obtain phase oscillations about Kolmogorov-Arnold-Moser (KAM) surfaces that ruin the simplicity of the output so that we would have difficulty resolving KAM boundaries. Figure 2 shows an example of surface of section plots for a toroidal Alfvén eigenmode (TAE) in a tokamak plasma . The horizontal axis is the toroidal angle multiplied by the toroidal mode number of the TAE in the co-moving frame with the TAE. The vertical axis is the major radius in the weak field side. In the surface of section plots we print the major radius R/a and phase, n − t , of a particle each time the poloidal angle of the particle reaches the midplane = 0 or . The particles plotted in the figure have the same E ′ and magnetic moment. We see an island structure in the figure that indicates that the particles are trapped by the TAE.

Saturation of the energetic particle driven instabilities
When an energetic particle driven instability grows to a finite amplitude, it is expected that the instability is saturated by the particle trapping. The saturation level will be determined by a condition b ∼ L where b is the bounce frequency for Fig. 2 Surface of section plots for a toroidal Alfvén eigenmode in a tokamak plasma . Reproduced from  with the permission of AIP Publishing the finite-amplitude wave and L is the linear growth rate of the instability. Since the bounce frequency is in proportion to the square root of the wave amplitude, the saturation level is in proportion to the square of the linear growth rate. The saturation due to the particle trapping was demonstrated by the earlier works on nonlinear simulation of Alfvén eigenmode (Wu and White 1994;Fu and Park 1995;Todo et al. 1995;Briguglio et al. 1995). The quadratic scaling of the Alfvén eigenmode saturation level on the linear growth rate is also investigated in detail (Briguglio et al. 2017;Slaby et al. 2018). The quadratic scaling means that larger phase space islands are created for larger growth rate. However, the AE spatial profile limits the radial width of the islands. This leads to a weaker scaling of the saturation level for larger growth rate.

Berk-Breizman's theory
In fusion plasmas, fast ions are created by the neutral beam injection, the ion-cyclotron-range-of-frequency heating, and the fusion reaction for alpha particles. The fast ions slow down and are scattered in velocity space due to the collisions with the bulk electrons and ions. The fast ion losses also take place when the particle enters the loss cone and/or are transported by the electromagnetic perturbations to the plasma edge. Then, the fast ion distribution is formed with source and sink. Since the unstable distribution is gradually formed with source and sink in experiment, we need to consider marginally unstable states. No unstable distribution appears suddenly. We should consider energetic particle distributions close to the marginal stability together with source and sink.
Berk, Breizman, and their collaborators established a weakly nonlinear theory for the marginally unstable energetic particle driven waves (Berk et al. 1996). The theory applies to the 1-dimensional electrostatic problem, but it also applies to general energetic particle driven instabilities including Alfvén eigenmodes in toroidal plasmas . Let us consider an electrostatic wave with the frequency and the wave number k. Suppose an electron distribution function shown in Fig. 3. The gradient of of the distribution function is positive at v = ∕k . This leads to the growth of the wave through the resonant interaction. This instability is called the bump-on-tail instability.
We follow Berk et al. (1996) to derive the cubic equation that describes the time evolution of the marginally unstable energetic particle driven waves. The Vlasov equation of the distribution function F(x, v, t) for a one-dimensional electrostatic problem is given by where Ê (t) and are the electric field amplitude and the phase of the wave. The source is represented by S(v), and the sink or dissipation is modeled with the Krook collision frequency . The distribution function F is expanded in a Fourier series, The time evolution of Ê (t) is given by where d is the intrinsic damping rate of the wave. In Eq. (19), it is taken into account that the total wave energy is double the electric field energy, and the other half is the electron kinetic energy for the electrostatic wave. We introduce the bounce frequency B (t) = √ ekÊ(t)∕m , and assume that the electric field amplitude is so small that B ≪ L , d , , (1∕t) is satisfied. The linear growth rate L is defined later. The perturbed distribution functions are assumed to be so small that (15) yields the equations for the time evolution of f 0 and f 1 . We used a formula for delta function, A factor of 1/2 has been multiplied for the derivation of Eq. (24) since t � − t takes only negative value in the integration. Then, Eqs. (19) and (24) give the equation for the time evolution of the electric field amplitude without the damping term Next, we consider the contribution from f * 1 to f 0 , The 3rd order f 1 is generated from f 0 , We integrate the 3rd order f 1 given by Eq. (30) in velocity space, The other contribution to f 0 from the 1st order f 1 is accompanied by a delta function (t � − t + t 2 − t 1 ) . This term does not generate any net contribution to the 3rd order f 1 because t � − t + t 2 − t 1 = 0 is satisfied only when t � = t and t 2 = t 1 . Then, . (29) Eqs. (19), (24), and (32) give the equation for the time evolution of the electric field amplitude up to the 3rd order of the wave amplitude, (33) is rewritten with a normalized parameter ̂= ∕( L − d ) (Berk et al. 1996),

Applications to Alfvén eigenmodes: steady amplitude and frequency splitting
Equation (34) depends on only one parameter, ̂ . For large ̂>̂c r = 4.38 , A( ) reaches the steady saturation. The steady saturation amplitude is A 0 = 2 √ 2̂2 which gives As ̂ decreases below ̂c r , the transition of the solution takes place from the steady saturation to periodic limit cycle solution, chaotic regime, and explosive growth. The explosive growth goes beyond the applicability range of the weakly nonlinear analysis and leads to the saturation levels that are expected from particle trapping The weakly nonlinear theory was extended with diffusion in velocity space (pitchangle scattering) . The extended theory was applied to Alfvén eigenmodes observed in JET (Fasoli et al. 1998;Heeter et al. 2000), and successfully demonstrated that the periodic limit cycle solution explains the "pitchfork splitting" of the Alfvén eigenmodes. The theory was also extended with drag in velocity space (slowing down) (Lilley et al. 2009). The scaling of steady saturation levels with collision frequency was confirmed for Alfvén eigenmodes with numerical simulations (Lang et al. 2010;Slaby et al. 2018). Micro-turbulence in fusion plasmas may work as an effective collision or dissipation, and may affect the time evolution of Alfvén eigenmodes. The effect of micro turbulence on the Alfvén eigenmode evolution has been theoretically predicted to be larger than that of collisions for present tokamaks and ITER (Lang and Fu 2011;Duarte et al. 2017).

Nonlinear evolution of single Alfvén eigenmode
In the weakly nonlinear theory discussed in the previous section, the wave amplitude was assumed to be so small that the expansion in the wave amplitude converges within the third order. We know that the explosive solution given by the weakly nonlinear solution grows to a large amplitude beyond the applicability range of the analysis. In this section, a nonlinear simulation model that can deal with the particle trapping by a finite-amplitude wave is presented. The simulation model is not restricted by the condition B ≪ L , d , , (1∕t) that is assumed for the weakly nonlinear theory presented in Sect. 3. On the other hand, the simulation uses a perturbative approach where the wave spatial profile is assumed fixed, while amplitude and phase of the wave and the energetic-particle nonlinear dynamics is followed self-consistently. The simulations based on the model clarify the saturation of the instability caused by the particle trapping and the spontaneous frequency chirping of the wave with the creation of the hole-clump structure in phase space.

Nonlinear reduced simulation model for the bump-on-tail problem
Since the destabilization mechanism of Alfvén eigenmodes (AEs) is the inverse Landau damping, the evolution of the bump-on-tail instability is illuminating for the understanding of the AE evolution. In the bump-on-tail instability, an electrostatic wave is destabilized by energetic particles that form a bump on the tail of the distribution function. If the velocity gradient of the distribution function is positive at the resonant velocity that is equal to the wave phase velocity, the wave grows through the inverse Landau damping. Berk and Breizman constructed a reduced simulation model to investigate the nonlinear evolution of the bump-on-tail instability . The derivation of the reduced simulation model is based on the Lagrangian formalism. For a one-dimensional electrostatic wave and energetic electrons, the reduced simulation model can be derived in a relatively simple way as follows. Let us start with the linearized momentum equation of electron fluid and the Maxwell equation, where m and −e are electron mass and charge, n 0 is the uniform equilibrium density, and u(x, t) and E(x, t) are the electron fluid velocity and the electric field, respectively. We assume that the fluid velocity is zero in the equilibrium and the temperature is neglected for simplicity. We know the following solution of Eqs. (36) and (37) where p is the electron plasma frequency, and Ê and û are real. When we consider energetic electrons, Eq. (37) is extended with the energetic electron current density j h (x, t), When Eqs. (36) and (42) are integrated explicitly in time, the time step width is limited by the Courant condition for the plasma oscillation. If we assume that the evolution of the electric field amplitude is sufficiently slower than the plasma oscillation period and the electrostatic wave is well approximated by the linear solutions given by Eqs. (38) and (39), we can eliminate the Courant condition for the plasma oscillation utilizing the linear solutions. Equations (36) and (42)  The second order term in cancels out with the first term on the right hand side. We retain the first order terms in , and define new variables Ê cos =Ê cos(− ) and Ê sin =Ê sin(− ) . Then, we have the following reduced equations.
− mn 0 2 1 2iû where we have introduced the wave damping rate d and L is the system length. The advantage of the reduced model is that the time scale is the growth and damping times, which are longer than the wave oscillation period. Then, the reduced model enables simulations in a longer time than the original equations that describe the wave oscillation. An important property of the reduced model is that it incorporates the phase ( ) evolution in addition to the amplitude evolution. Equation (44) indicates that the term with the time derivative of phase ( ̇ ) is the same order of as the term with the time derivative of amplitude. Then, the time derivative of phase (= frequency shift) should be retained together with the amplitude evolution. Another important property of the reduced model is that the energetic particle current j h ,cos and j h ,sin are divided by a factor of 2 on the right-hand sides of Eqs. (46) and (47). This factor appears through the mathematical derivation procedure, but physically indicates the equi-partition of the plasma oscillation energy to the electric field and the electron fluid motion.
The time evolution of the energetic particle distribution function f h (x, v, t) is described by the Vlasov equation, and the energetic particle current j h (x, t) is given by the integration of the distribution function f h (x, v, t) in velocity space where f 0 (x, v) = f (x, v, t = 0) and the term with is the Krook collision operator that appeared in the weakly nonlinear theory presented in Sect. 3. The Vlasov equation can be simulated with either Lagrangian approach (particles) or Eulerian approach. The particle simulation results are presented in the following subsections.

Saturation due to particle trapping
Let us consider a 1-dimensional problem where the system size in x direction is L = 8 v t h ∕ , the wave number of the wave is k = 2 ∕L = ∕4v th and the resonance velocity with the wave is v res = 4v t h . The energetic particle distribution is chosen so that the linear growth rate is L ∕ = 2.42 × 10 −2 . The time evolution of the electric field amplitude A = √Ê 2 cos +Ê 2 sin is shown in Fig. 4 for various cases. Let us start with a special case for d = 0 and = 0 . The time evolution of the electric field amplitude is shown in Fig. 4a. We see in the figure linear growth, saturation, and amplitude oscillation after the saturation. The saturation level A = 0.024 gives the bounce frequency B = 0.077 ∼ 3.2 L . Figure 5 shows the distribution function for the different phases of the evolution. We see the vortex structure formation in phase space that indicates the particle trapping by the wave. This leads to the saturation of the instability. The amplitude oscillation is brought about by the bounce motion of the trapped particles by the electric field.

Hole-clump creation and frequency chirping
The most surprising discovery of the nonlinear reduced simulation was the spontaneous formation of the hole and clump structure in phase space and the consequent frequency chirping (Berk et al. , 1998(Berk et al. , 1999. Before this study, it was expected that the wave damps due to the intrinsic damping after the saturation of the instability was caused by the particle trapping. Figure 4b shows the result for d ∕ = 5 × 10 −3 and = 0 . We see in the figure the saturation of the instability and the damping of the wave amplitude to the noise level. This is what was expected before the discovery of the hole-clump creation. However, the story is different if the damping rate is close to the linear growth rate. Figure 4c shows the wave amplitude evolution for d ∕ = 2 × 10 −2 and = 0 . The damping rate ( d ) is close to the linear growth rate L ∕ = 2.42 × 10 −2 . We see in the figure that the wave amplitude does not monotonically damp, but the time-average amplitude is kept at a constant level. The time evolution of the wave frequency spectrum is shown for another case with d ∕ = 1 × 10 −2 and = 0 in Fig. 6. The frequency chirps both upward and downward. The frequency up-shift and down-shift are associated with the formation of the hole and clump structure in phase space. Figure 7 shows the distribution function for various moments. We see in the figure that the holes and the clumps, which are islands in phase space, are formed and shift to higher and lower velocities in time, respectively. The holes and the clumps are types of the Bernstein-Greene-Kruskal (BGK) solution that is consistent with the electric potential distribution. The velocity up-shift (down-shift) of the hole (clump) releases energy. The frequency chirping rate is determined so that the released energy balances with the energy dissipation due to the wave damping (Berk et al. 1999). The frequency chirping is in proportion to the square-root of time and given by, An exact solution and a theoretical framework for long-time behavior of the holes and clumps were developed in Breizman (2010), Nyqvist et al. (2012). A helpful explanation for the understanding of the frequency chirping is given in Breizman (2011). The spontaneous frequency chirping takes place for both the assumed constant mode damping and the self-consistent damping due to heavy stabilizing species. Comprehensive simulation studies have been devoted to the nonlinear evolution of the bump-on-tail problem including the hole-clump pair creation (Vann et al. 2005;Lesur et al. 2009;Lilley et al. 2010). The frequency chirping of Alfvénic modes has been observed in many tokamak plasmas (Shinohara et al. 2002;Fredrickson et al. 2006;Gryaznevich and Sharapov 2006).
The frequency chirping of an energetic particle mode and Alfvén eigenmodes in tokamak plasmas have been investigated by computer simulations, and the relation between the frequency chirping and the hole-clump pair creation has been discussed Pinches et al. 2004;Zhu et al. 2013. Asymmetry between the upward chirping and the downward chirping are observed in the experiments and the simulations. The asymmetry is attributed to the non-uniform distribution of the free energy source (Pinches et al. 2004; and to the interaction with the shear Alfvén continuous spectra. The time evolution of energetic particle driven geodesic acoustic modes (EGAMs) is often accompanied by the frequency chirping (Boswell et al. 2006;Nazikian et al. 2008;Ido et al. 2016). The frequency chirping and the sudden excitation of EGAMs observed in LHD experiments were successfully reproduced by kinetic MHD hybrid simulations (Wang et al. , 2018. Though the reduced simulation model may not be applied to the EGAMs for which energetic particles have the non-perturbative effects on the real frequency and the spatial profile, the creation of hole and clump was demonstrated for a chirping EGAM with the kinetic MHD hybrid simulation ).

Effect of collisions
Next, we introduce a finite Krook collision frequency . Figure 4d shows the result for ∕ = 1 × 10 −2 . The linear growth rate and the damping rate are L ∕ = 2.42 × 10 −2 and d ∕ = 2.4 × 10 −2 . The growth rate that was observed in the simulation result is ∕ = 5.43 × 10 −4 . The normalized collision frequency is ̂= ∕ = 18.4 . Here, we replaced L − d by observed in the simulation. We see the steady amplitude in the figure. The weakly nonlinear theory predicts the steady saturation, and the saturation level given by Eq. (35) is A = 1.7 × 10 −4 . The saturation level in the simulation result is A = 2.0 × 10 −4 , which is in good agreement with the theoretical prediction. The collision term in Eq. (50) of the reduced simulation model plays two roles, source and sink. The term f is a sink because it dissipates fine structures in phase space, while the term f 0 is a source because it restores the equilibrium distribution function. The steady state shown in Fig. 4d is established by a balance between the distribution flattening caused by the wave and the source and the sink provided by collisions.

Energetic particle transport by Alfvén eigenmodes
The interaction between energetic particles and Alfvén eigemodes (AEs) leads to energetic particle transport and losses. The particle trapping by a single AE flattens the energetic particle spatial profile. The flattening due to a single AE takes place in a localized region near the resonance. When multiple AEs exist and their amplitudes are sufficiently high, the resonance overlap may occur leading to the emergence of stochasticity and the global transport. Stochasticity may arise even for a single AE with large amplitude. The higher-order (nonlinear, fractional) resonance discussed in Sect. 2 leads to the emergence of stochasticity for a single AE. Energetic particles may be lost when they are transported to the plasma edge or the particle orbit transitions to another type of orbit with a large orbit width during the interaction with the AEs. Both the particle trapping and the stochasticity in phase space can lead to energetic particle losses. We discuss in this section the energetic particle transport by a single AE and multiple AEs, and energetic particle loss. Figure 8 shows surface-of-section plots for counter-passing fast ions in the strong field side of a tokamak plasma for a TAE with n = 3 . The plasma center and the edge are located at R∕a = 3.2 and R∕a = 2.2 , respectively. Figure 8a shows the plots for the field amplitude B∕B = 2 × 10 −3 . We see islands that are the trapped regions by the TAE. The energetic particle distribution is flattened in the islands. In addition to the islands, we see a stochastic region near the plasma edge. The emergence of stochastic region can be understood when we investigate another case with lower TAE amplitude. Figure 8b shows the surface-of-section plots for the amplitude B∕B = 8 × 10 −4 . Now the particle dynamics hardly have any region of stochasticity. Instead, we see the emergence of second-and fourth-order islands around R∕a = 2.27 and R∕a = 2.35 , respectively, in addition to the two first-order islands around R∕a = 2.22 and R∕a = 2.32 . With increasing field amplitude these islands overlap to eventually destroy the KAM surfaces and create the stochastic region that appears in Fig. 8a. The stochastic regions can be created by the overlap of higher-order islands of a single AE. This effect was first demonstrated by Refs. , . Fast ion losses by energetic particle driven geodesic acoustic modes through the higher-order resonance were observed in DIII-D experiments (Kramer et al. 2012).

Fig. 8
Surface of section plots for counter-passing fast ions in the strong field side for a tokamak plasma for a TAE with n = 3 when the field amplitude is fixed in time at: a B∕B = 2 × 10 −3 and b B∕B = 8 × 10 −4 . The plasma center is located at R∕a = 3.2 and the plasma edge at the strong field side is R∕a = 2.2 . Reproduced from  2888 with the permission of AIP Publishing

Resonance overlap among multiple Alfvén eigenmodes
When multiple AEs exist, each eigenmode has trapped regions (islands) in phase space. The trapped regions of the different AEs may overlap each other. When the trapped regions overlap, the particles inside the overlapped regions can move around the whole of the overlapped region and their motions become stochastic. The surface-of-section plots presented in the previous subsection use particles with the same E ′ value. Since E ′ depends on the AE frequency and the toroidal mode number, we should make the plots separately. The surface-of-section plots for different AEs consist of particles with different energy and toroidal canonical momentum. The resonance overlap of multiple AEs cannot be discussed in an exact sense if we focus on the particles with the same E ′ . We need to extend the analysis to twodimensional phase space (P , E).
A new analysis method to distinguish resonance overlap has been developed (Todo et al. 2016). In Todo et al. (2016), fast ion transport by multiple AEs in a DIII-D plasma was investigated for various beam deposition power with hybrid simulations for energetic particles interacting with an MHD fluid. The resonance regions in fast ion phase space was analyzed for a single eigenmode with fixed amplitude and frequency that are observed in the simulation results. Figure 9 shows the particle trajectories in the phase space of normalized major radius R = (R − R axis )∕(R edge − R axis ) and energy E [keV] with the AEs observed in the simulations for DIII-D experiments for beam deposition power (a) 1.56 MW, (b) 3.13 MW, (c) 6.25 MW, and (d) 15.6 MW (Todo et al. 2016). Here, P is replaced by R to clarify the spatial location. The AE amplitude is larger for higher beam deposition power. The particle orbits were followed with the electromagnetic perturbations of a single eigenmode for 2 ms, and the position in (R, E) space was recorded as each particle passed the mid-plane from bottom to top. Collisions are turned off to clarify the resonance regions. The particles are co-going to the plasma current with the same magnetic moment which gives v ∥ ∕v = 0.63 for E = 70 keV. The peak of the fast ion distribution is located at v ∥ ∕v = 0.63 . Only the particles trapped by the eigenmode are plotted in the figure. The time for this analysis ( = 2 ms ) is 1-30 times longer than the bounce period of the wave-particle trapping. The particles are initially located uniformly in (R, E) space with intervals 0.01 in horizontal axis and 1 keV in vertical axis. The particles plotted can be regarded as resonance regions in phase space. The eigenmodes are represented by colors: n = 1 (blue), n = 2 (purple), n = 3 (green), n = 4 (orange), and n = 5 (red).
In Figure 9a for P NBI = 1.56 MW, only the resonance region of the n = 1 mode (blue) overlaps the resonance regions of the other modes, but the overlapped regions are small in the phase space. For P NBI = 3.13 MW shown in Fig. 9b, the resonance regions of n = 2-5 modes broaden due to the larger amplitude of eigenmodes. However, we see that gaps exist between the resonance regions of n = 5 (red), n = 4 (orange), and n = 3 (green) for 0.4 <R < 0.6 . These gaps prevent the fast ion transport flux from increasing. For P NBI = 6.25 MW shown in Fig. 9c, the resonance regions overlap substantially for 0.4 <R < 0.6 . This leads to the enhancement of fast ion transport flux for 0.4 < r∕a < 0.6 . The overlapped region spreads outward up to R = 0.8 for P NBI = 15.6 MW shown in Fig. 9d, which spreads the fast ion transport flux profile up to r∕a = 0.8 . This analysis demonstrated that the enhanced fast ion transport observed in the simulations can be attributed to the resonance overlap among the multiple AEs.
The importance of the resonance overlap of the multiple AEs was first predicted in , , . In DIII-D experiments, anomalous flattening of the fast ion spatial profile was observed with a rich spectrum of toroidal Alfvén eigenmodes and reversed shear Alfvén eigenmodes in the current rampup phase ). The electron temperature fluctuations due to the AEs were observed with the electron cyclotron emission measurement. The electron temperature fluctuations are well reproduced by analyses based on MHD models (Van Zeeland et al. 2006Zeeland et al. , 2009Todo et al. 2015). The magnetic fluctuation amplitude of the AEs is of the order of B∕B ∼ O(10 −4 ) . Numerical analyses of the fast ion transport demonstrated that the multiple AEs observed in the experiment with the amplitude B∕B ∼ O(10 −4 ) bring about the significant fast ion transport and the profile flattening (White et al. 2010a, b). The significantly flattened fast ion pressure profile due to the multiple AEs was reproduced by the comprehensive kinetic-MHD hybrid simulations with neutral beam injection and collisions of fast ions (Todo et al. 2015). Experiments in the DIII-D tokamak show that fast-ion transport suddenly becomes stiff above a critical threshold of beam power in the presence of many overlapping smallamplitude AEs (Collins et al. 2016). The sudden increase in fast ion transport flux and its phase space dependence was also reproduced with the comprehensive hybrid simulations (Todo et al. 2016). In the hybrid simulation results, monotonic degradation of fast ion confinement and fast ion profile stiffness is found with increasing beam deposition power. The confinement degradation and profile stiffness are caused by the sudden increase in fast ion transport flux that is brought about by multiple AEs for fast ion pressure gradients above a critical value. The redistribution of fast ion profile with multiple AEs was also observed on JT-60U tokamak ) and the Large Helical Device (LHD) (Osakabe et al. 2006).

Fast ion loss
Fast ion losses induced by AEs have been measured in tokamak and stellarator/ heliotron devices with the scintillator-based detectors (Darrow et al. 1997;Weller et al. 2001;Isobe et al. 2006;Darrow et al. 2008;García-Muñoz et al. 2010;Van Zeeland et al. 2011;Ogawa et al. 2012). The lost particles measured by the detectors on the tokamaks are mainly trapped particles that were originally counter-passing particles to the plasma current, lost energy due to the interaction with the AEs, and became trapped particles. Since the orbit width of trapped fast ion is large in the present devices, the trapped fast ions may reach the fast ion loss detectors located outside the plasma. We discussed this particle loss mechanism for Fig. 1b in Sect. 2. For the LHD experiments, the lost fast ions measured with the detector were co-passing particles and transition particles (Ogawa et al. 2012).
Two types of fast ion loss have been observed. The fast ion loss rate of the first type is linearly proportional to the AE amplitude, while the second type has a quadratic dependence on the AE amplitude (García-Muñoz et al. 2010;Ogawa et al. 2012). The fast ion loss signal for the first type oscillates with the same frequency as that of the AE(García-Muñoz et al. 2010). In test particle simulations where the AE amplitude and the frequency are assumed to be constant, the quadratic dependence of fast ion loss rate on the AE amplitude was predicted when stochasticity emerges in phase space ). The first type of fast ion loss is related to convective transport and the second type is called diffusive transport. Both the linear and the quadratic dependences were reproduced by test particle simulations (Van Zeeland et al. 2011;Ogawa et al. 2012). Reduced simulations of multiple AEs (Schneller et al. 2013) and hybrid simulations for the bursting evolution of multiple AEs  reproduced the quadratic dependence of the fast ion loss rate on the AE amplitude. We can attribute the first type of fast ion loss to a single low-amplitude AE without stochasticity in phase space and the second type to a single large-amplitude AE with stochasticity or multiple AEs with resonance overlap which generates stochasticity.
However, we should notice that fast ion transport flux is based on the product of the perturbed distribution function ( f ) and the radial particle velocity perturbation ( E + v ∥ ). Neither the phase space analysis such as the surface-of-section plots 1 Page 26 of 33 nor the test particle simulations consider the effect of f whose evolution should be consistent with that of the AEs. The two types of fast ion response to an EPM were observed with the directional Langmuir probe located inside the plasma of the Compact Helical System (CHS) (Nagaoka et al. 2008). The first type of response oscillates with the same frequency as that of the EPM and is linearly proportional to the EPM amplitude, while the second type has zero frequency and is proportional to the square of the EPM amplitude. It was clearly demonstrated that the phase difference between the first type of response and the EPM generates a net transport of fast ions (Nagaoka et al. 2008). When we consider a single AE destabilized by fast ions, the perturbed distribution function is given by an equation similar to Eq. (16) where the linear and the quasi-linear responses of the fast ion distribution function are given by f 1 and f 0 , respectively. We know in Eqs. (16), (23), and (28) that f 1 oscillates with the same frequency as that of the AE and has linear dependence on the AE amplitude, and the frequency and the dependence on the AE amplitude of f 0 is 0 and quadratic, respectively. These are the general properties of the wave-particle interaction and the same as those observed with the fast ion loss detectors and the directional Langmuir probe, but are not related to the emergence of stochasticity. More works would be needed to conclude the relationship between the fast ion transport processes and the fast ion loss measurements.

Discussion and summary
We have explained the basics of the interaction between energetic particles and Alfvén eigenmodes (AEs), which is an important research issue for burning plasmas. Experimental, theoretical, and computational studies of the interaction have been extensively conducted. The various types of time evolution of AEs such as steady state, frequency splitting, frequency chirping, and recurrent bursts have been observed in the experiments in tokamak and stellarator/heliotron plasmas. Berk and Breizman presented both a one-dimensional weakly nonlinear theory for marginal stability and a reduced simulation model that qualitatively explain the various types of time evolution. The fast ion profile flattening and losses brought about by AEs have been also measured in the experiments. These experimental results promoted the development of theory and computer simulation, especially the extension with source (NBI, ICRF, and birth of alphas in future work) and sink (collisions and losses). When we turn our eyes on the academic aspects, the AE instability is a kind of inverse Landau damping that is the fundamental research subject of plasma physics. We can say that the studies of AEs have pioneered a new frontier in plasma physics on the nonlinear evolution of (inverse) Landau damping with source and sink.
Let us discuss the limitations of Berk-Breizman's reduced simulation model and the future works. In this article, we focused on the nonlinear wave-particle interaction and neglected the nonlinear "wave-wave" interaction. The author and his collaborators performed the first numerical demonstration of TAE bursts with parameters similar to a TFTR experiment (Wong et al. 1991) using the reduced simulation model presented in section 4, and reproduced many of the experimental characteristics . However, the saturation amplitude of the magnetic field fluctuation normalized by the toroidal field was B∕B ∼ 2 × 10 −2 , which is higher by one order of magnitude than the value B∕B ∼ 10 −3 inferred from the experimental plasma displacement measurements Durst et al. 1992). In the simulation of Ref. , 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 to the higher harmonics reduces the saturation level of the TAE (Todo et al. 2010. The zonal flow generation induced by AEs was first found in a gyrofluid simulation (Spong et al. 1994) and was found also in gyrokinetic simulations (Zhang and Lin 2013;Biancalani et al. 2016). Theory of the zonal flow generation by TAE has been also developed .
Another important missing subject in this article is the energetic particle modes (EPMs) for which energetic particles play an essential role on the mode frequency and the spatial profile. When the energetic particle distribution is modified by EPMs, the frequency and the spatial profile of the EPM are also modified. The reduced simulation model presented in Sect. 4 cannot be applied to EPMs. The reader is referred to a review article that discusses EPM in detail . Comprehensive kinetic MHD hybrid simulations have been performed for EPMs in JT-60U plasmas (Bierwage et al. 2014(Bierwage et al. , 2018. For the evolution of EPM, the interaction with the shear Alfvén continuous spectra is also an essential issue. The interaction with the shear Alfén continuous spectra affects also the time evolution of the gap modes such as TAE and RSAE when they have the continuum damping and/ or the chirping frequency hits the continuum. A theoretical model of the frequency chirping of AEs with the interaction with shear Alfvén continua has been developed (Wang and Berk 2012).
An interesting correlation between the amplitude and the chirping frequency of TAE was observed in NSTX experiments (Podestà et al. 2011). A similar correlation was observed for BAE in gyrokinetic simulations (Zhang et al. 2012). Even for AEs whose frequency is located inside the continuum gap, it has been pointed out that energetic particles have non-perturbative effects on the mode spatial profile and the real frequency Todo et al. 2005;). These modes cannot be simulated with the reduced model. Kinetic-MHD hybrid simulations (Park et al. 1992;Spong et al. 1992;Todo et al. 1995;Briguglio et al. 1995;Todo and Sato 1998;Wang et al. 2010;) and gyrokinetic simulations (Mishchenko et al. 2009;Lang et al. 2009;Zhang et al. 2010;Bass and Waltz 2010) will be powerful tools to reproduce and predict the nonlinear behaviors of energetic particles and Alfvénic modes beyond the limitations of the reduced simulation.