<figure><img src="image_27.png" /><figcaption>Figure 26. Alfvén continuum for $n = 12$ in ITER obtained by the ideal MHD code NOVA. The gaps in the continuum correspond to the BAE, TAE, EAE and BAAE modes. The low frequency BAE gap is due to Alfvén continuum upshift by the coupling of the Alfvénic and acoustic waves. The BAE gap hosts the GAM solution. This coupling also leads to the formation of a new gap near the rational magnetic surface called the BAAE gap shown in the insert. Dashed curves in the insert show the coupling of the Alfvénic and acoustic waves. Dashed horizontal lines correspond to the global eigenmode solutions. Reproduced from [30].</figcaption></figure>

analytically [458] which was later reproduced in Ref. [464]. The BAAE gap structure was recently reexamined [465, 466] where the coupling with the neighboring acoustic harmonics was ignored [467].

All branches are modified in the kinetic description compared to the MHD limit described so far, and new, purely kinetic branches arise. In particular, the low-frequency part can be strongly influenced by diamagnetic effects and kinetic effects of the thermal ion transit or bounce drift motion with characteristic frequencies $\omega_{*i}$ and $\omega_{ti,\theta}$. Here,

$$ \omega_{*i} = \omega_{*ni} + \omega_{*Ti} = \left( \frac{T_i b}{ZeB} \times \left( \frac{\nabla n_i}{n_i} + \frac{\nabla T_i}{T_i} \right) \right) \cdot k \\ \approx \frac{T_i}{ZeB} \frac{1}{n_i} \frac{\partial n_i}{\partial r} (1 + \eta_i) k_{m,\theta} \quad (55) $$

where $k_{m,\theta} \approx -m/r$ is the poloidal wave number and $\eta_i = \frac{\partial}{\partial r} \log T_i / \frac{\partial}{\partial r} \log n_i$.

General dispersion relations including finite orbit width effects have been derived [454, 462, 468–470] and implemented in various codes [29, 471, 472]. The general fishbone-like dispersion relation (see also section 7.2), relying on the separation of scales between the singular and regular layers, includes additional expressions for the EP contribution to the EP-drive and the thresholds to EPM branches that exist for

all the modes discussed above [32, 418]. The EPM branches often emerge from the least damped part of the spectrum where the drive can overcome continuum damping.

As already mentioned, an overview over some of the properties of AEs are given in table 4. It should be noted, however, that these expressions describe the gap locations, whereas the equilibrium non-uniformities have to be considered for calculating the global mode frequencies, since they remove the degeneracy with the continua. Thus, in many cases of practical interest global analyses are required for quantitative comparisons. Useful approximate formulas including finite beta effects, ellipticity and triangularity for many AEs are given in [473] and in the Appendix of [474]. For BAAEs and the low-frequency Alfvén modes (LFAMs) with frequencies $\omega \sim \omega_{ti,\theta} \sim \omega_{*i}$, the kinetic dispersion relation has to be solved, including ellipticity [475] and trapped particle contributions [461, 462, 464, 469, 476, 477].

### 6.3. AE eigenfunction

The AE eigenfunction is important because its structure determines both the resonance condition (Eq. 11) and the amount of energy the particle
<table><caption>**Table 4.** Summary of basic properties of various AE branches. TAE: Toroidal AE, EAE: Ellipticity-induced AE, NAE: Non-circularity-induced AE, RSAE: Reversed-shear AE, GAE: Global AE, BAE: Beta-induced AE, BAAE: Beta-induced acoustic AE. LFAM: Low-frequency Alfvén modes. ◦ RSAEs only exist at frequencies larger than the BAE frequency. † For the BAAEs in the MHD limit, we quote a simplified expansion around the rational surface according to [458,461,466]. ‡ For LFAM please refer to further details in the text.</caption><thead><tr><th>Mode</th><th>Δ<sub>m</sub></th><th>q</th><th>ω/ω<sub>A</sub></th></tr></thead><tbody><tr><td>TAE</td><td>1</td><td>(m + 1/2)/n</td><td>1/2q</td></tr><tr><td>EAE</td><td>2</td><td>(m + 1)/n</td><td>1/q</td></tr><tr><td>NAE</td><td>3</td><td>(m + 3/2)/n</td><td>3/2q</td></tr><tr><td>RSAE<sup>◦</sup></td><td>0</td><td>q<sub>min</sub></td><td>m/q<sub>min</sub>(t) - n</td></tr><tr><td>GAE</td><td>0</td><td>m/n</td><td>min[ω<sub>A</sub>(r)]/ω<sub>A</sub></td></tr><tr><td>BAE</td><td>0</td><td>m/n</td><td>√β<sub>i</sub>(7/4 + T<sub>e</sub>/T<sub>i</sub>)</td></tr><tr><td>BAAE<sup>†</sup></td><td>1</td><td>~ m/n</td><td>(n - m/q)/√(1 + 2q<sup>2</sup>)</td></tr><tr><td>LFAM<sup>‡</sup></td><td>0</td><td>m/n</td><td>~ 0.5ω<sub>*i</sub>/ω<sub>A</sub></td></tr></tbody></table>

exchanges with the wave. By definition, a resonant particle is one that encounters the same wave phase on multiple orbits [37]. Consequently, the net linear energy transferred per orbit $\Delta\mathcal{E}$ is the integral of the instantaneous power transfer (given in Eq. 12) averaged over a complete orbit,

$$ \Delta\mathcal{E} = Ze \oint \delta\mathbf{E} \cdot d\mathbf{l}. \qquad (56) $$

Equation (56) implies that the eigenfunction must be known accurately to calculate the drive term $\gamma_{\text{drive}}$ correctly so, in order to determine the stability threshold correctly, both the eigenfunction and the resonant orbits need to be known accurately.

Fortunately, the linear mode structures of all AEs are well understood. There are three directions to consider: the toroidal mode number $n$, the poloidal structure that is (in general) a combination of different mode numbers $m$, and the radial structure. The most unstable toroidal mode number occurs when the eigenfunction and EP orbit are comparable in size. In the simplest theory of AEs driven by the spatial gradient, the EP drive is proportional to the EP diamagnetic drift frequency $\omega_{*f}$

$$ \omega_{*f} = \left( \frac{T_{\text{f,eff}} \mathbf{b}}{ZeB} \times \left( \frac{\nabla f_{\text{f}}}{f_{\text{f}}} \right) \right) \cdot \mathbf{k} \qquad (57) $$

with $T_{\text{f,eff}} = \int v^2 f_{\text{f}} dv$. If $f$ is separable as $f(E, \mu, P_{\phi}) = n(r)f(E, \mu)$, a simplified expression is

$$ \begin{aligned} \omega_{*f} &= \omega_{*nf} + \omega_{*Tf} = \left( \frac{T_{\text{f,eff}} \mathbf{b}}{Z_{\text{f}} e B} \times \left( \frac{\nabla n_{\text{f}}}{n_{\text{f}}} + \frac{\nabla T_{\text{f,eff}}}{T_{\text{f,eff}}} \right) \right) \cdot \mathbf{k} \\ &\approx \frac{T_{\text{f,eff}}}{ZeB} \frac{1}{n_{\text{f}}} \frac{\partial n_{\text{f}}}{\partial r} (1 + \eta_{\text{f}}) k_{m,\theta} \end{aligned} \qquad (58) $$

where $k_{m,\theta} \approx -m/r$ is the poloidal wave number and $\eta_{\text{f}} = \frac{\partial}{\partial r} \log T_{\text{f,eff}} / \frac{\partial}{\partial r} \log n_{\text{f}}$. The EP diamagnetic drift frequency is linearly proportional to the poloidal mode

number $m \simeq n \cdot q$ [478]. However, the radial extent of the eigenfunction tends to decrease with increasing toroidal mode number $n$ so, when the EP orbit size significantly exceeds the spatial extent of the mode, the energy exchanged with the wave $\Delta\mathcal{E}$ decreases, reducing the EP drive [479]. The most unstable toroidal mode number occurs when the spatial extent and orbit size are comparable, an expectation that is roughly consistent with experiment [480]. Because of the large size and magnetic field of ITER, this scaling predicts that the most unstable toroidal mode numbers $n$ in ITER will be many times larger than in most contemporary experiments; an ITPA benchmark study (see section 6.4.4) predicts that mode numbers between $n \simeq 20$ and $n \simeq 30$ will be most unstable.

Available measurements of the poloidal structure agree with the expectations summarized in Table 4. For example, soft x-ray measurements of the poloidal structure of a GAE are consistent with a single $m$; in contrast, nearly equal values of $m$ and $m+1$ describe a TAE [481]. Similarly, electron cyclotron emission images of an RSAE at the start of the frequency sweep are consistent with a single $m$ [482,483], as expected.

The predicted radial structure has also been compared with experiment. In the example shown in figure 18 in section 4, the radial amplitude profile measured along the midplane agrees well with ideal MHD theory for both an RSAE (figure 18(b)) and a TAE (figure 18(c)) [280]. Another example of excellent agreement between the TAE radial amplitude profile and ideal MHD theory appears in [484]. RSAEs with the same toroidal mode number can appear with different radial mode numbers; the lowest radial mode is approximately Gaussian in shape but higher radial mode numbers contain nodes [483,484].

Many authors have noted that the AE radial phase profile depends upon the radial location of EP drive and wave damping. For example, a thorough theoretical overview of the connection between localized sources and sinks of energy, energy and momentum transport, and mode structure appears in [485]. Experimental observations of radial phase variation are commonplace for both TAEs and RSAEs [482, 484, 486, 487]. A recent paper [488] conducted a comprehensive survey of thousands of RSAE and TAE mode structure measurements from DIII-D beam-heated discharges. Amidst extremely large variability in the radial phase profile, the phase profile for a typical RSAE is approximately flat, while a typical TAE often has a phase ramp at large major radius. The authors speculate that the typical TAE radial phase profile reflects power flow from smaller minor radius (where the EP gradient and hence power flow from the EPs is large) to edge regions (where wave damping predominates).
During the current ramp phase of the discharge, RSAEs often first appear when $q_{min}$ is a rational number, then the mode sweeps upward in frequency as the q-profile evolves [268]. Observations of RSAEs that sweep downward in frequency to a minimum at the time of the rational $q_{min}$ crossing are relatively rare; for example, < 1% of JET cases surveyed were downward sweeping [489]. Proposed theoretical explanations for the enhanced stability of downward-sweeping RSAEs are summarized in [490]. Experimentally, the radial eigenfunctions of downward-sweeping and upward-sweeping RSAEs are very similar [490].

Despite its importance for correct mode identification and stability calculations, measurements of mode polarization are rare. Recently, however, a method to infer the parallel component of the vector potential $\delta A_\|$ and the effective parallel electrostatic potential $\delta\phi_\| = \delta\phi - \omega\delta A_\|/ck_\|$ from measurements of electron temperature and density $\delta T_e$ and $\delta n_e$ was developed [491]. Application of the method to instabilities in DIII-D shows that RSAEs and TAEs have small values of $\delta\phi_\|$ (as expected for shear Alfvén wave polarization) but lower frequency modes often have appreciable electrostatic potentials.

It should be noted that the measured AE mode structure is already in the nonlinear phase. Nonlinear simulations with the gyrokinetic code GTC found that the AE eigenfunction shifts radially in time as the EP gradient evolves, as observed experimentally [492]. A computationally expensive nonlinear simulation with the MEGA code achieved good agreement with experimental eigenfunction measurements for both mode amplitude and radial amplitude and phase profiles [493]. (This MEGA simulation modelled the EPs kinetically and the bulk plasma with resistive MHD.)

### 6.4. Alfvén frequency range: TAEs, EAEs, and RSAEs

Here we consider AEs in the intermediate frequency range on the order of the Alfvén frequency (TAEs and EAEs) and RSAEs or ACs with frequencies typically below the TAE frequency. Neutral beam, ICRF-accelerated, and energetic electron populations have driven these AEs unstable in virtually all toroidal magnetic confinement configurations [494]. TAEs and EAEs driven by electrons during lower-hybrid heating have been observed on EAST [495]. Alpha-driven TAEs were observed on TFTR [294] and in the recent DT campaign at JET using magnetics, soft x-ray, interferometry and reflectometry diagnostics [40]. At JET, it was clear that alphas drove the TAEs since they were observed in an afterglow experiment 50 ms after NBIs had been turned off and fusion alphas were still being produced, as illustrated in figure 27. The

observed mode was found on the outboard midplane and was driven by both trapped and passing alphas originating from the plasma core.

#### 6.4.1. Perturbative and non-perturbative simulations.
EP-driven instabilities are often described by a perturbative approach, i.e. the waves are linear eigenmodes of the bulk plasma and their frequencies and spatial structures are determined by the bulk plasma only. However, to model AEs with shear Alfvén polarization accurately, the descriptions of the interaction of EPs require the more sophisticated global, non-perturbative approach [496]. The term non-perturbative means that the simpler perturbative approach cannot be applied since the interaction between the EPs and the modes is too strong, and so the mode properties are not only determined by the bulk plasma, but also by the EPs. Non-perturbative eigenmode features also arise when the coupling to kinetic Alfvén waves is considered.

For the assessment of the damping, usually the gap structure is analyzed first in order to find regions without continuum damping. It is important to realize that some of these modes are in between rational surfaces $q = m/n$, e.g. the TAE $q_{TAE} = (m + 1/2)/n$, whereas others reside at rational surfaces with very small $k_\|$. Both electron and ion Landau damping depend very sensitively on $k_\|$, since the ratio of $\omega/k_\|v_{th,i,e}$ enters the damping rate exponentially [497]. Close to $k_\| = 0$, ion Landau damping can only occur through sidebands, and electron Landau damping can be dominated by trapped electrons. For TAEs, the high-energy tail of the thermal Maxwellian ions can contribute to the damping of the two main harmonics, whereas electron Landau damping depends critically on the background non-uniformities, i.e., the distance to the next rational surface [444]. Often trapped electrons dominate the damping, and thus also trapped-electron collisional damping has to be considered [498, 499].

The coupling of global modes to the kinetic Alfvén wave (KAW) introduces short wavelength-scales on the order of the ion Larmor radius $\rho_{Li}$ into the system, constituting the microscale. Formally, this coupling arises due to the ion polarisation term and finite Larmor radius contributions in the coupled system of quasi-neutrality and the gyrokinetic momentum equation causing finite $E_\|$ [500] that is dissipated by electron Landau damping. Since the damping scales with $\sim (k_\perp \rho_{Li})^2$, the mode structure needs to be calculated non-perturbatively. A non-ideal parameter [445, 501] can be derived for a local estimate of the KAW coupling, but, for global modes with various (radial) contributions to the radiative damping, this can lead to difficulties when interpreting
<figure><img src="image_28.png" /><figcaption>**Figure 27.** Mode observations during the afterglow in JET DT shot #99946. From top to bottom: Fusion neutron rate, NBI power, Fourier spectra from interferometry, and SXR emission. The TAE appears in both spectra at 7.05 s with frequency 115 kHz (circled). Reproduced from [40].</figcaption></figure>

experimental measurements. This issue was pointed out in [502, 503], but the strongly coupled TAE-KAW wavefield found in this work could not be reproduced with modern gyrokinetic codes, except when electron Landau damping is artificially reduced or omitted [504, 505].

This multi-scale, non-perturbative nature of the problem and the requirements for low-frequency modes (section 6.2) motivated various studies using global, fully self-consistent, gyrokinetic simulations that can simultaneously resolve micro- and mesoscales to map out the linear stability boundaries quantitatively and to accurately predict the nonlinear saturation amplitude and transport level for contemporary and future fusion devices. We will deal with non-perturbative approaches in section 7.

### 6.4.2. AE antenna measurements and benchmarks.
A powerful way to probe the linear stability of AEs is to launch waves with an external antenna [506]. By sweeping the antenna frequency across the resonance and detecting the response with fast magnetic coils, the resonant frequency $\omega$, net damping rate $\gamma/\omega$, and toroidal mode number $n$ of the mode can be inferred. These are then compared with theory and numerical simulations, improving extrapolations to future burning plasma operations such as in ITER and DEMO. Measurements of $\gamma/\omega$ are particularly

valuable, as the mode drive must exceed mode damping for linear instability. Since the method is most sensitive to modes with radial eigenfunctions that are appreciable near the antenna, TAEs, EAEs, or globally extended GAEs are usually studied [507]; for similar reasons, low-$n$ modes are more readily detected than high-$n$ modes.

The Alfvén antenna technique was introduced at JET [508]. After around 15 years and more than 100,000 individual damping rate $\gamma/\omega$ measurements, the original saddle coil system used for the active excitation of low-frequency AEs in JET was replaced with a set of eight, small and localized in-vessel antennas capable of driving AEs with higher toroidal mode numbers. Recall that these higher mode numbers are relevant for ITER and future burning plasma experiments. Details on the design of these new antennas and on the first results obtained with this system can be found in [509–515]. The system can drive, detect and track individual AEs for a given range of toroidal mode numbers [516, 517].

These first sets of measurements were then followed by a full campaign, totaling more than 30000 $\gamma/\omega$ measurements for individual toroidal mode numbers, now ranging from the low-$n$ previously probed with the saddle coil system (still a majority, covering about 50% of the measurements) to modes up to $|n| = 10 - 12$. Details of these measurements,
which covered the analysis of the dependence of $\gamma/\omega$ for higher $n$ modes on the edge plasma shape and on the plasma isotopic composition can be found in [518-520]. One important result is that the damping rate tends to increase with increasing plasma isotope mass (figure 28). The first detailed comparisons with theory for these higher-$n$ modes can be found in [521-524], and the first use of these higher-$n$ modes for MHD spectroscopy of the plasma isotopic composition can be found in [525].

<figure><img src="image_29.png" /><figcaption>**Figure 28.** JET antenna measurements of the damping rate and the mode frequency for $n = 7$ TAEs as a function of the effective isotope mass for modes with frequencies near the bottom (left), center (middle), and top (right) of the TAE gap. The damping rate increases strongly with increasing mass in all cases. Reproduced from [520].</figcaption></figure>

A dedicated ITPA EP group effort aimed to understand the damping rates as measured by the JET TAE antenna. Perturbative, fluid codes and non-perturbative, kinetic codes were benchmarked against each other and validated against experimental antenna results from JET for an elongation scan (figure 29 [526]). In this study, all codes captured the experimental trend, showing that radiative damping is the dominant contribution (see also [527]) for an open TAE gap (no continuum damping). However, figure 30 demonstrates that the damping in general cannot be estimated by one set of local parameters, since the global mode structure can cause both core and edge KAW coupling.

Over the past few years, an upgraded JET antenna diagnostic [531-536] has been operated on nearly a thousand JET plasma discharges, measuring thousands of stable AE resonances. Analyses of data collected during the 2019-2020 JET campaign of mainly deuterium plasmas has been reported, including comparisons with gyrokinetic codes [537-

<figure><img src="image_30.png" /><figcaption>**Figure 29.** ITPA benchmark and validation study: JET elongation scan comparing the damping rates of the TAE antenna with various code results; see also subsection 6.4.2. All involved codes captured the experimental trend of an increasing damping rate for increasing elongation. Reproduced from [526].</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>**Figure 30.** ITPA benchmark and validation study: the radial absorption profile for a global $n = 3$ TAE based on JET discharge #77788, 10.157s as calculated by LEMAN [528,529] and LIGKA [530], normalized to the LEMAN result $\gamma/\omega = 0.95\%$ [526]; see also section 6.4.2. The damping in general cannot be estimated by one set of local parameters. Reproduced from [526].</figcaption></figure>

541]. Statistical analysis of these data shows that, as expected, the damping rate is correlated with parameters that are related to continuum damping and radiative damping. Since the edge safety factor determines if the TAE gap is open or closed at the edge, the importance of continuum damping is reflected in a statistically significant dependence on $q_{95}$ [542]. ($q_{95}$ is the safety factor at the surface the contains 95% of the poloidal flux.) The measured damping rate also correlates strongly with the non-ideal parameter associated with the finite Larmor radius effects that cause radiative damping (figure 31) [542].

In recent work, a marginally stable AE was tracked by the JET antenna diagnostic in real-time
<figure><img src="image_32.png" /><figcaption>**Figure 31.** Normalized damping rate as a function of the non-ideal parameter $\lambda$ that is associated with radiative damping for a large database of JET antenna measurements. The correlation with $\lambda$ is the strongest dependence in the database. Reproduced from [543].</figcaption></figure>

in a discharge with 25 MW of external heating [543]. Interestingly, the AE was deduced to be an edge-localized EAE with a lab-frame frequency aligned with the core TAE gap. AEs were also monitored from destabilization through stabilization in a JET EP experiment [128].

AE antenna studies have also been conducted on MAST and C-Mod [544]. In C-Mod, damping-rate measurements for inner wall limited and diverted cases over a range of toroidal mode numbers $4 \le |n| \le 9$ were compared with NOVA-K calculations that included continuum and radiative damping [545]. Good agreement between the measurements and theoretical calculations were obtained but the results depended sensitively on the assumed $q$-profile.

### 6.4.3. Code benchmarks.
The ITPA EP group has performed several code benchmark studies on various machines. A detailed study of damping mechanisms for AEs, that successfully compares eigenvalue and initial value gyrokinetic codes for an AUG case, can be found in [546]. It should be noted that a straightforward comparison of numerical results with the general nonlinear theory 7.2 can be carried out by adopting the mode structure decomposition approach [547] and the ballooning formalism [471,548].

For benchmarking the linear drive, also in the non-perturbative regime, a synthetic test case with aspect ratio 10 has been set up within the ITPA EP group [549]. 11 codes with a broad variation in the physical and numerical models participated. The eigenfunctions as well as mode frequencies match in a satisfactory way. The growth rates were within around 20%. However, they are found to depend

strongly on the complexity of the model (see figure 32). Further code verification is necessary to improve the match of the results by including further physical effects. Furthermore, it was found that the mode frequency and mode structure of TAEs can change in the presence of a considerable EP population, as also documented for a DIII-D case [550]. This work also found symmetry breaking effects of TAEs due to the EPs. Such effects were analysed in more detail for BAEs and RSAEs [551–554]: an asymmetry in the EP gradient with respect to the rational surface of the modes leads to finite $k_r$ and tear-drop shaped 2D mode structures in excellent agreement with measurements [351,482,483,550,555–557], which can lead to enhanced radial energy transport [558].

<figure><img src="image_33.png" /><figcaption>**Figure 32.** Comparison of growth rates computed with various codes including FLR effects for an $n = -6$ TAE. The dashed line from CAS3D-K is valid in the limit of zero orbit width (small energies) and is shown for comparison. The shaded grey area marks the ±20% margin around the mean value. Reproduced from [549].</figcaption></figure>

Extensive studies of RSAEs have been carried out on DIII-D [494]. RSAEs on spherical tokamaks have been found to be largely suppressed due to their relatively higher $\beta$ [559]. A linear RSAE benchmark and a comparison between simulations and experiments based on a DIII-D discharge (figure 33) [555, 560] was done with modern gyrokinetic and gyrokinetic-MHD codes. RSAE frequencies and mode structures for $n = 4, 5, 6$ of eight participating codes agree very well with the experiment as measured by ECEI. The RSAE frequencies differ by less than 5% whereas the growth rates differ by about 20%. It is found that pressure gradients of the thermal plasma make a significant contribution to the growth rates. The differences become larger when the RSAE-TAE transition is approached by varying $q$. Due to
the open SAWC gap, various TAEs closely spaced
in the frequency domain but with different radial
structures can co-exist. Moreover, using a calculated
EP profile that takes diffusion by multiple unstable
modes into account, a TAE with *n* = 6 is found
to be unstable at the outer edge, which is consistent
with the experimental observations. Finally, electron
temperature fluctuations and radial phase shifts from
simulations show no significant differences with the
experimental data for the strong *n* = 4 RSAE, but
significant differences for the weak *n* = 6 TAEs. Other
DIII-D stability studies have found: (1) that flattening
the EP spatial gradient with off-axis beams stabilizes
RSAEs, consistent with TAEFL stability calculations
[561], (2) that creation of a bump-on-tail feature in
the beam-ion velocity distribution has little impact
on TAE and RSAE stability, consistent with MEGA
calculations [562], and (3) that, owing to their longer
slowing-down time, deuterium beams drive RSAEs
more strongly than hydrogen beams, as predicted by
MEGA [563].

Gyrokinetic simulations can incorporate realistic toroidal geometry and comprehensive physics such as effects of equilibrium currents [564] and compressible magnetic perturbations [565,566]. As discussed above, these simulations find that accurate damping and growth rate calculations require non-perturbative, fully self-consistent simulations to calculate the true mode structure. The radial position of the AEs can move with the location of the strongest EP pressure gradients as confirmed by experimental data [550]. Further validation of gyrokinetic simulations of TAEs have been carried out on many tokamaks including KSTAR [567], HL-2A [568], AUG [546, 569, 570] and JET [540,571].

6.4.4. *ITER predictions.* Based on scenario simula-
tions using the ASTRA code [572] (very similar to sim-
ulations stored under #131018, 50 in the ITER IMAS
database), the ITPA EP group motivated a study con-
cerning TAE stability and related EP transport for
the standard ITER *Q* = 10, 15 MA scenario. Var-
ious groups contributed to this study [31, 573-581],
extending previous studies [582, 583]. It is encourag-
ing to observe, that all studies using nearly identical
profiles predict only weakly unstable TAEs with small
or negligible EP transport. As mentioned above, the
range of the most unstable mode numbers lies between
*n* = 20 – 30. However, it has to be noted that the
*n*-number spectrum of the most unstable linear modes
varies significantly, depending on the model and the de-
tails of the safety factor profile. This leads to different
results when carrying out sensitivity scans in order to
determine the thresholds for benign AE transport. For
example, density peaking effects not included in the

<figure><img src="image_34.png" /><figcaption>**Figure 33.** Comparison of the linear dispersion relation calculated with various codes for an RSAE in DIII-D shot 158243 at 805 ms. Top: real frequencies. Bottom: growth rates. The plot markers are diamond, star, and circle for the gyrokinetic, kinetic-MHD, and perturbative eigenvalue codes, respectively. Reproduced from [560].</figcaption></figure>

ASTRA simulations may lead to steeper alpha gradi-
ents. It has been shown in [581] that global gyrokinetic
models [584] and gyrokinetic-MHD models [530, 585]
agree reasonably well on linear and nonlinear features,
emphasizing the importance of the linearly stable spec-
trum in the nonlinear phases. It should be noted
that for this comparison neither collisions nor zonal
flows were included for benchmark reasons. Relaxing
these simplifications is an ongoing effort (see subsec-
tion 7.3) needed for an increasingly realistic descrip-
tion and prediction of a burning fusion plasma. Recent
developments try to mitigate the problem of exponen-
tial sensitivity of the mode stability with respect to the
background profiles by setting up integrated automatic
workflows that combine equilibrium evolution and lin-
ear stability, in order to correctly capture parametric
dependencies (e.g during ramp-up or during heating
power scans) [577, 586].

Apart from the 15 MA scenario, also various
studies of the 9 MA steady-state ITER scenario have
been carried out [573, 587, 588]. Similar mode number
ranges are found to be unstable (again with a strong
dependence on the background profile details), but
in this case RSAEs, RSAE/TAE hybrids, and BAEs
are expected to be the most unstable modes, rather
then TAEs. Significant EP transport was predicted for
nominal parameters [587], hinting to the necessity to
include this EP transport in the scenario modeling in
order to obtain consistent background and EP profiles
that may deviate considerably from neoclassical values.
Further, the drive due to off-axis neutral beam heating
is expected to be more significant than in the 15 MA scenario [588].

## 6.5. High-frequency modes: CAEs and GAEs

Significant progress has been achieved in understanding and modeling of modes with frequencies in the ion cyclotron frequency range (higher than $\omega_A$) since the publication of the Progress in the ITER Physics Basis [20]. This frequency band includes the compressional and the global Alfvén eigenmodes (CAEs and GAEs). These modes do not cause strong transport of EPs in radius but they do cause transport in energy. They are primarily connected with the ion cyclotron emission (ICE) phenomena. The physics of CAEs in particular and its relationship with ICE is discussed in two review papers [589, 590]. Here we present a brief introduction to GAEs and CAEs. These modes have been found in low aspect ratio tokamaks [591–594] and in the conventional aspect ratio tokamaks DIII-D [595] and AUG [191, 596, 597]. Experiments have studied the parametric stability boundaries for a wide variety of EP-driven modes in NSTX [598].

### 6.5.1. CAEs
CAE linear stability theory has potentially attractive prospects for understanding ICE as a potential EP diagnostics (see section 6.5.3). The main elements of CAE linear stability theory have been developed [599].

CAE theory can be presented heuristically on the basis of the fast Alfvén wave dispersion relation, $\omega_{CAE} = k_{\perp}v_A$ for $k_{\parallel} \ll k_{\perp}$, as expected for these modes in tokamaks [600]. The CAE dispersion relation can be used to obtain an eigenvalue equation for the dominant magnetic perturbation $\delta B_{\parallel}$, i.e. $\hat{k}_{\perp}^2 \delta B_{\parallel} = (\omega_{CAE}^2/v_A^2) \delta B_{\parallel}$ where $\hat{k}_{\perp}^2$ is treated as a differential operator [601]. To leading order in inverse aspect ratio $\epsilon = r/R$, this equation becomes

$$ \frac{1}{r} \frac{\partial}{\partial r} r \frac{\partial}{\partial r} \delta B_{\parallel} + \frac{1}{r^2} \frac{\partial^2}{\partial \theta^2} \delta B_{\parallel} = \frac{\omega_{CAE}^2}{v_A^2(0)} \frac{n(r)}{n_0} (1 + \epsilon \cos \theta)^2 \delta B_{\parallel}, \quad (59) $$

where $r$ is the minor radius coordinate, $\theta$ is the poloidal angle, $n(r)/n_0$ is the plasma density normalized to its value on the magnetic axis, and $v_A(0)$ is the Alfvén speed on axis. The coefficient of $\delta B_{\parallel}$ on the right hand side plays the role of an effective potential and has a minimum absolute value in the low field side midplane. Within the heuristic approach we follow here, CAEs are located in the vicinity of this minimum.

Radially localized CAE solutions can be found for potential wells which are narrow and deep in the radial direction and shallow and long in the poloidal direction. The CAE poloidal wavelength is then shorter than the radial wavelength. This justifies the

choice of the eikonal for the following poloidal mode structure [600]:

$$ \delta B_{\parallel}(r, \theta) = b(r, \theta) \exp[-i\omega t + im(\theta + \epsilon_0 \sin\theta) - in\phi], \quad (60) $$

where $\phi$ is the toroidal angle, $(m, n)$ are the poloidal and toroidal mode numbers, and the subscript 0 denotes the value taken at the minimum of the local effective potential well, $r = r_0$. Because the potential well is typically poloidally elongated in tokamaks, the wave number components of CAEs have the ordering

$$ k_{\theta} \gg k_{\phi} \gg k_{r} \quad (61) $$

as also found in CAE stability theory.

Because of the complexity of the above eikonal, the realistic CAE dispersion relation is difficult to analyze and compare with experiments. For example, the theoretical methods employed in [602, 603] are based on the assumption for the eikonal and so are approximate. Instead, a heuristic dispersion relation of CAEs [604] leads to the identification of characteristic lengthscales and corresponding “quantum” numbers in each of three relevant directions: the toroidal mode number, $n$, associated with the major radius $R$; radial mode number, $S$, associated with the radial width of the effective potential in the radial direction, and the poloidal mode number, $M$, associated with the plasma minor radius. The CAE mode frequency is then

$$ \omega_{CAE}^2 \approx v_A^2 \left( \frac{M^2}{r^2} + \frac{S^2}{L_r^2} + \frac{n^2}{R^2} \right), \quad (62) $$

which is consistent with the tokamak ordering (61).

CAE solutions obtained numerically [604] using the ideal MHD code NOVA agree with the dispersion relation (62) which is consistent with the eikonal (60). It was also found to be consistent with numerical frequency splitting for low $n = (0, 1)$ but not for higher $n$ due to the strong coupling of the dominant compressional Alfvénic polarization of CAEs and the shear Alfvénic polarization of KAW structures. This may mask the observed frequency splitting because, in experiments, the most unstable modes are excited first [605], and they likely have high toroidal mode numbers, $n \gg 1$ [606–608]. Numerical simulations show the complicated CAE dispersion relation which only qualitatively agrees with equation (62). This is due to finite $k_{\parallel}$ and $\omega/\omega_c$ effects as was pointed out recently [609].

One new finding in CAE theory is that the shift of the eigenmode frequencies depends on the sign of CAE poloidal mode number $m$. For a plasma cross section with ellipticity $\kappa$, the CAE eigenfrequencies were found to be asymmetric with respect to the sign of the poloidal mode number [610]:

$$ \omega_{CAE} = \frac{m}{r} \sqrt{\frac{1 + \kappa^{-2}}{2}} v_{A*} $$
$$ \times \left[ \frac{\sigma_m v_A (\ln \rho_m)'}{2\omega_{ci}(\kappa)} + \sqrt{1 + \left( \frac{v_A \rho_m'}{2\omega_{ci}(\kappa) \rho_m} \right)^2} \right]_{r_*} , \quad (63) $$

where * represents location of the eigenmode, ' the radial derivative, and $\sigma_m = m/|m|$. The location of the eigenmode is given by

$$ 2 + r (\ln \rho_m)' - \sigma_m \frac{v_A}{\omega_{ci}(\kappa)} (r (\ln \rho_m)')' \sqrt{1 - \frac{2 + r (\ln \rho_m)'}{(r (\ln \rho_m)')'}} = 0, \quad (64) $$

where $\rho_m$ is the plasma density, $\omega_{ci}$ is the thermal ion cyclotron frequency.

The expression for CAE growth rate was obtained in the tokamak approximation, which neglects $k_\parallel$ and assumes high frequencies, $\omega/\omega_{ci} \gg 1$ [590,611]:

$$ \frac{\gamma}{\omega_{ci}} \approx \frac{\omega^3 \sqrt{2ecB}}{\omega_p^2 \omega_{ci}^2 \sqrt{\pi \Delta r_0 R_0}} \times \sum_{l, \sigma} \int dP_{\phi} d\mathcal{E} d\mu I^2 \frac{E_1^2 \mu l J_l^2}{E_0^2 z^2} \left[ \frac{\partial}{\partial \mathcal{E}} + \frac{l\omega_{ci}}{\omega B} \frac{\partial}{\partial \mu} \right] \tilde{f}, \quad (65) $$

where $\omega_p$ is the plasma frequency, $E_1 = E_0 f(r, \theta)$ is the CAE structure in the poloidal cross section required for proper averaging of the local growth rate expression, $\Delta$ is its radial mode width, $I^2 = 8\pi/|\frac{d}{dt}(\omega_{ci} + \omega_D)|$ is the resonant factor accounting for phase variation of wave particle interaction near the resonance, $J_l$ is the Bessel function of order $l$ with the argument $z = k_\perp \rho_{Lf}$, and $\rho_{Lf}$ is the Larmor radius of the EPs. The resonance condition is then $\omega - l\omega_{ci}(r, \theta) - \omega_D(r, \theta) = 0$ which is to be evaluated along the EP drift orbit.

### 6.5.2. GAEs.
GAE modes were found theoretically, with $q = m/n$ and the frequencies just below the minimum of the ideal MHD shear Alfvén continuum [612], and experimentally in tokamaks when plasma was heated by ICRF at frequencies lower than the cyclotron frequency [613]. These modes, called global Alfvén eigenmodes (GAEs), are localized in radial direction near the continuum minimum point outside the plasma centre [614]. These modes are also called the conventional GAEs, when their frequency is just below the Alfvén continuum minimum point where the GAE mode structure is localized. The so-called non-conventional GAEs (NGAEs) are found in stellarators [615] and in helical plasmas. They were studied theoretically in an attempt to explain sudden drops of the plasma beta during the low-frequency instabilities in the W7-AS shearless stellarator due to subsequent electron heating of the plasma periphery. According to [615], NGAEs have frequencies above the maximum of the Alfvén continuum.

The GAE dispersion relation, $\omega_A = k_\parallel v_A$, offers an efficient way to identify the modes, such as at NSTX [607], especially if the GAE frequency signal peaks overlap with each other, and their poloidal and toroidal numbers are known. Indeed, GAE modes with different $(m, n)$ mode number pairs not only have different frequencies but they also have distinguishable temporal frequency patterns which was confirmed by direct comparison with the Alfvén wave dispersion relation [616].

Estimates of the GAE growth rate were made from the time-dependence of the frequency chirps and compared with more traditional estimates based on the growth and decay rates of mode bursts [593, 617]. These estimates were found to agree well with numerical estimates by the HYM code. The frequency and mode number dependence of GAEs on the toroidal field and the EP distribution function has been documented in NSTX and NSTX-U [617].

In DIII-D, sub-cyclotron modes were originally identified as CAEs [595] but more recent analysis persuasively argues they were actually GAEs [618]. Recent experiments compared GAE theory and experiment at the stability threshold [608] and confirmed that the toroidal mode number of the unstable GAEs is $n \simeq 20$ [185].

A more accurate analytic evaluation of CAE and GAE growth rate drive in the sub-cyclotron frequency range was done recently with applications to spherical tokamak devices [619]. One important observation is that the addition of a small number of tangential NBI ions ($\sim 7\%$) stabilizes the counter-propagating GAE instability in NSTX-U [620], as observed experimentally [617,621].

These modes observed in NSTX can exhibit EPM features [622], as was found numerically, a new and surprising result in GAE stability theory. EPMs will be treated in sections 6.7 and 7. The frequency of the most unstable GAEs changes significantly with beam parameters such as injection energy and injection geometry. It was demonstrated that such changes depend on the Alfvén velocity and are in qualitative agreement with the Doppler-shifted cyclotron resonances driving the modes. This unexpected result holds for counter-propagating GAEs, which are routinely excited in NSTX, and for high frequency co-propagating GAEs, which had not been previously studied. Large changes in frequency without clear corresponding changes in the mode structure are signatures of an EPM, referred to as EPM GAEs [622].

### 6.5.3. ICE-based diagnostics for burning plasma devices.
As discussed in [589], EP-driven ICE has been detected in many fusion devices and is commonly
identified with CAEs: compressional Alfvén waves can be excited by non-Maxwellian populations of energetic ions through cyclotron resonances [623]. New observations of ICE have been reported from JT-60U [187,188], LHD [199,200,624], AUG [189–192], KSTAR [193,194], NSTX and NSTX-U [195], DIII-D [185,196], and EAST [197, 198]. EP-driven ICE in tokamak plasmas has been studied using either ICRF antennas [625] or dedicated RF probes [626]. Nevertheless, ICE detection is in principle possible using any technique for measuring electric, magnetic or density fluctuations with a sampling rate in the ion cyclotron range. Density fluctuations can be observed using microwave reflectometry, which is planned for ITER [627] and has been used in NSTX to determine the structure of eigenmodes with frequencies up to several hundred kHz [486] and at DIII-D with frequencies up to ~40 MHz [628].

High time-resolution measurements of the ICE bursts on MAST [591], NSTX and NSTX-U have revealed that frequency-chirping can occur during the short bursts of ICE. The growth rate determined from the frequency chirping was in good agreement with the direct measurement of the edge magnetic fluctuation amplitude growth [629].

Simulations are crucial to gain understanding of the ICE phenomenon. Fully nonlinear simulations of ICE in the ion cyclotron range and in the presence of energetic ion populations approximating those in tokamaks and stellarators have been carried out using the fully kinetic particle-in-cell (PIC) code EPOCH [630] and also using the kinetic-MHD code PROMETHEUS++ in which bulk and energetic ions are kinetic whereas electrons are treated as a charge-neutralising fluid [631]. A limitation of the latter approach is that the waves cannot exchange energy with the electrons. However, PIC simulations typically show that in the nonlinearly saturated state only a small fraction of the energy lost from the energetic ion population is transferred to electrons [632] while a substantial fraction is transferred to bulk ions [630, 632]. For this reason the excitation of ICE is one of several possible means of channeling the energy of fusion alphas directly into bulk ions rather than via electrons, thereby providing the possibility of sustaining the hot ion condition required for the thermonuclear fusion burn more effectively [633]. This idea for direct ion heating is sometimes called “alpha channeling”, which we will discuss in section 10.

A key result from the kinetic-MHD simulations of ICE reported in [631] is that the nonlinearly-saturated ICE intensity (as measured by the square of the perturbed magnetic field component parallel to the equilibrium field $\delta B_\|)$ was found to be proportional to the energetic ion concentration. This scaling

reproduces a linear dependence of ICE intensity on neutron flux over six orders of magnitude observed in JET [625], which indicates that ICE in future burning plasma devices could be used to obtain diagnostic information on the fusion alpha distribution.

Another conclusion of the PIC and kinetic-MHD studies of ICE is that nonlinear wave-wave interaction plays an essential role in the excitation of emission at lower ion cyclotron harmonics which are often found to be linearly stable. The three-wave interactions underlying these nonlinear excitation processes in the simulations can be studied using bicoherence analysis, as discussed for example in [630].

Most likely, ICE will be observed in ITER through the detection of magnetic fluctuations, as in the majority of contemporary experiments. Furthermore, ITER will have an ICRF heating system [141], and it has been demonstrated on JET that ICRF systems can be used to detect ICE either passively (during periods in which ICRF heating is not used [625]) or actively (during ICRF heating). In the latter case, detection was made possible in JET through the use of a sub-harmonic arc detection system [632]. The RF signals were sampled in the ICRF antenna transmission lines, which acted as resonators. The frequency response of these resonators depended on the antenna configuration and the length of the transmission line matching elements. As a result, the measured ICE spectra were strongly filtered, limiting the information that could be obtained. More complete spectra have been obtained using ICRF antennas in receiver mode (i.e. passively) [625] and it would be worth considering this option in the DT operation phase of ITER.

## 6.6. Low frequency modes: BAEs, BAAEs, LFAMs, EGAMs

### 6.6.1. BAEs.
In recent years, much attention was given to low-frequency modes in the BAE and BAAE regime. BAEs have been widely observed in present tokamaks. The excitation of BAEs is affected by both the energetic and thermal ion populations as well as coupling between the Alfvénic and ballooning mode branches [454, 634]. The BAEs can be divided into the three groups i-BAEs, e-BAEs and m-BAEs, depending on different energy sources. The letters ‘i’, ‘e’ and ‘m’ refer to mode instability drive by energetic ions, by energetic electrons and by magnetic islands, respectively.

The i-BAEs were first confirmed on DIII-D operated near the beta limit [453] and proved to be deleterious, causing large losses of EPs at levels similar to TAEs. The mode frequencies are approximately half the TAE frequency and display long-living or quickly-chirping behavior, depending on the plasma scenario. The i-BAEs are always localized in the core
<figure><img src="image_35.png" /><figcaption>**Figure 34.** ECEI spectrogram showing the bursting BAE modes associated with the *q* = 2 surface on AUG. (*b*) These modes are also observed by the FILD. Figures (*c*)–(*f*) show the 2D mode structures of the selected modes α–δ, respectively (in pairs of left the amplitude *A* and right *A*cos(Φ)). Modes α and β are 4/2 modes, modes γ and δ are 2/1 modes. The modes are observed to move radially outward with time. Reproduced from [483].</figcaption></figure>

region. The mode structures have been measured by ECEI diagnostic at AUG, as shown in figure 34. The 2D imaging suggests that the mode numbers are either *m/n* = 4/2 or 2/1. The modes move radially outward, following the outward moving *q* = 2 surface. Such a localized mode structure can also be observed by the microwave reflectometer on HL-2A [635]. Linear excitation threshold analysis of BAEs based on observations at Tore Supra indicates that ion Landau damping is important for the mode in the acoustic frequency range. The BAEs will be driven unstable only when EP-drive exceeds ion Landau damping [636]. A major study of BAE stability in DIII-D found that the observed mode frequencies are usually close to analytic estimates of the BAE accumulation point, that the modes occur in bursts with rapid frequency chirping, and that BAEs are most likely to be unstable when the poloidal beta exceeds 0.5 [637]. A DIII-D study of isotope dependence found that, with deuterium NBIs, BAEs are at least as unstable in mixed hydrogen-deuterium plasmas as in deuterium plasmas; however, with hydrogen NBIs, BAEs are stable [638]. In HL-2A, modes were more readily excited in low-density plasmas [635] but the density dependence is relatively weak in DIII-D.

The e-BAEs were first observed on HL-2A [639]. They can be excited in both Ohmic and ECRF-heated plasmas with low densities (figure 35). Since these modes appear in plasmas with low ion temperatures, the mode frequencies are always lower than for i-BAEs and range from 10 to 30 kHz. The e-BAEs are driven unstable by barely circulating or deeply trapped particles. The mode characteristics are affected not

<figure><img src="image_36.png" /><figcaption>**Figure 35.** An e-BAE observed on HL-2A. (a)–(c) Counts of HXR photons in different energy ranges, (d) central line-averaged density, (e) magnetic probe signal, (f) spectrogram of magnetic probe signal. Reproduced from [639].</figcaption></figure>

only by the population of the energetic electrons, but also by their energies and pitch angles. Theoretical analysis further indicates that the frequency and growth rate are sensitive to the energetic electron temperature, and there exists a maximum growth rate [640]. Moreover, the growth rates, the mode width and its radial asymmetry can be affected by the energetic electron density [472].

The BAEs occurring in the presence of magnetic islands (m-BAE) were first reported on FTU [641] and subsequently found on HL-2A [642] and EAST [643]. The m-BAE is a mode formed when counter-propagating Alfvén waves form a standing-wave structure within a magnetic island [644]. The m-
BAEs can be driven only when the magnetic islands width exceeds a threshold, and then they appear in pairs. Two m-BAEs usually propagate in different diamagnetic drift directions, and their mode numbers are $m/n = 2/1$ and $-2/-1$ for modes propagating in ion or electron diamagnetic drift directions.

The BAE dispersion relation and the related global modes as obtained by kinetic models [26, 451, 454, 469, 557, 646–650] could be clearly linked to experimental observations [469, 476, 483, 637, 650–653]. Alternative interpretations as Alfvén-acoustic modes [466] seem less likely, given that a clear relation to rational surfaces is observed experimentally [483, 637, 653], and that all branches with predominantly acoustic or mixed Alfvén-acoustic polarisation are generally very strongly ion Landau-damped, as demonstrated in direct MHD-gyrokinetic comparisons [654].

### 6.6.2. BAAEs and LFAMs.
As for the BAE frequency range, it has been recognised [29, 655] that kinetic theory is needed to describe modes in the BAAE frequency regime, because the mode frequencies become comparable to diamagnetic frequencies, and ion Landau damping is very effective. Modes with frequencies below the BAE frequency are observed on DIII-D and NSTX [645, 655] and have been initially dubbed BAAEs. However, experiments at ASDEX-Upgrade [469, 476, 483, 652] and DIII-D [645] show a clear connection of the observed instabilities to the evolution of rational surfaces, which contradicts the interpretation as Alfvén-acoustic gap modes. Accordingly, these modes were dubbed low-frequency Alfvén modes (LFAMs). Figure 36 shows a typical example of LFAMs on DIII-D [645]. LFAMs appear in ascending patterns in plasmas with EPs and high electron temperature but modest beta. The mode frequencies are in the diamagnetic frequency range, and the toroidal mode numbers are in the range $n = 3–12$. The mode occurrence is correlated with rational values of the safety factor $q$. The gyrokinetic dispersion relation [454, 469, 472, 656] was successfully applied to interpret the experimental frequency and mode number patterns of LFAMs and their connection to the kinetic ballooning branch [29, 464, 476, 656]. Extensions including the trapped particle response [462] are needed for accurate quantitative analyses. The excitation mechanisms were further investigated [464, 477, 656–658], demonstrating that EPs are not necessary to destabilise LFAMs (as reported in DIII-D experiments [645]), i.e. the drive can be provided by the thermal background pressure gradients ($\sim \omega_{*i}$). The modes exhibit a predominantly Alfvénic polarisation, are favoured by high $T_e$ and can have a reactive character. LFAMs are often observed at low

$\beta$ [645], which can be understood from the coupling of Alfvénic, diamagnetic and acoustic branches [29, 462, 464]. The same type of analysis also demonstrates that for increasing $\beta$ the BAE branch is favoured over the LFAM branch which leads to the conclusion that no strong EP transport is expected from LFAMs under reactor-relevant conditions. Modes with similar properties, but at higher thermal ion beta and $\omega_{*i}$ than in DIII-D, have been recently reported at JET [571]. Additionally, modes between the BAE and TAE frequencies at JET have been interpreted as global perturbations within higher-order geodesic Alfvén-acoustic gaps [463].

Due to different interpretations as BAAEs or LFAMs in the literature, simulated mode properties need to be carefully analysed when comparing to the experiment. However, for most typical experimental conditions, the BAAE branch is always more stable than the LFAM/KBM branch [656]. Recent global gyrokinetic simulations [658] find that the LFAM mode structure and many of its parametric dependencies are consistent with the theoretical analyses and the experimental observations: the linear growth rate increases rapidly with increasing electron temperature; the LFAM can be excited without EPs and has a frequency in the BAAE gap; an antenna scan confirms that it is not the conventional BAAE. Instead, the LFAM is an interchange-like electromagnetic mode excited by non-resonant drive of pressure gradients. Trapped electrons and equilibrium current have modest effects on the LFAM [658]. Recent polarisation measurements on DIII-D show that LFAMs exhibit a fast changing mixed Alfvén-acoustic polarisation during their lifetime, indicating that LFAMs indeed consist of a mainly Alfvénic component and electrostatic sidebands. While the location of the rational surface evolves, the measurements pick up a rapidly varying mix of electrostatic and electromagnetic components [491].

In the BAE/BAAE frequency range, there are two other instabilities, i.e. kinetic ballooning modes (KBM - typically at high $n$) due to finite diamagnetic drifts [659] and Alfvénic ion temperature gradient (AITG) modes due to kinetic thermal ion compressibility and wave-particle resonance [660, 661]. The AITG mode, a new branch connecting KBM and BAE modes, was first experimentally identified in Ohmic plasmas on HL-2A [662], see figure 37. These modes appear in plasmas with peaked density profiles and weak magnetic shear, which indicates that corresponding instabilities are excited by pressure gradients. The time trace of the fluctuation spectrogram is either a frequency staircase with different modes excited at different times, or multiple modes may simultaneously coexist. AITG-like modes that trigger disruptions are
<figure><img src="image_37.png" /><figcaption>**Figure 36.** LFAMs during discharges #178631 (left) and #132710 (right): (a) ECE spectrum from channels near $q_{min}$ on the reference discharge, each mode is labeled by the (m, n) value (b) $q_{min}$ from the equilibrium reconstructions (line) and fitted rational q values m/n (symbols) vs time. (c) Inferred frequency in the plasma frame $f_{lab} - nf_{rot}$ vs fitted toroidal mode number. Reproduced from [645].</figcaption></figure>

also observed in high ion temperature (> 10 keV) plasmas [663].

<figure><img src="image_38.png" /><figcaption>**Figure 37.** AITG modes with different mode numbers in an Ohmic plasma at HL-2A. The 2D patterns are spectrograms of ECEI (left) and SXR (right) signals. Left: with the frequency staircase. Right: with the multi-mode coexistence. Reproduced from [662].</figcaption></figure>

### 6.6.3. EGAMs.
EGAMs are similar to the standard GAMs [456], and as for standard GAMs, their existence is based on the geodesic curvature of the magnetic field lines. The main difference with respect to the standard GAMs is that they are driven by the EP population and therefore the EGAM frequency and radial structure depend on the EP distribution. EGAMs have been described analytically [443, 664–675] and have been observed on DIII-D [676–678], LHD [679–681], HL-2A [682, 683] and AUG [684, 685]. Similar modes were also observed on JET as global geodesic acoustic modes (GGAMs) [686, 687]. EGAMs have also been

studied in a series of numerical gyrokinetic simulations using GYSELA [667, 688–690], ORB5 [668, 691, 692] GENE [693] and GTC [694], as well as in kinetic-MHD simulations using MEGA [351, 550].

EGAMs are axisymmetric oscillations ($n = 0$) of the electrostatic potential, density, pressure and magnetic perturbations. The latter have been measured, for example, at LHD [681], HL-2A [682], DIII-D [678] and AUG [684]. Their location, mode structure and non-linear mode structure evolution has been measured by SXR e.g. in HL-2A [682] and AUG [684]. The electrostatic potential is dominated by a zonal structure ($m = 0$), the density and pressure exhibit a poloidally up-down asymmetric structure (with $m = 1$), as illustrated in figure 38, and the magnetic perturbation is characterized by $m = 2$ structure. The expected $m = 1$ structure of the density was imaged in [695], the potential in [681], and the radially broad $E_r$ structure was measured in [678].

EGAMs on LHD have been found to contribute to nonlinear destabilization of a subcritical mode (figure 39). Moreover, the EGAMs may couple with turbulence and significantly degrade the plasma confinement [688, 689]. Theoretical analysis reveals that EGAMs can be regarded as a potential energy channel to transfer the fusion-born alpha energy to the thermonuclear plasma, referred to as GAM channeling [666]. Such effects will be discussed in section 9. EGAMs are in some cases accompanied by strong bursting and frequency chirping. Also, significant density perturbations and large drops in neutron emission can be observed in the presence of EGAMs
<figure><img src="image_39.png" /><figcaption>**Figure 38.** Calculated $m = 1$ poloidal structure of a DIII-D EGAM. Reproduced from [676].</figcaption></figure>

[676, 677], which suggests EP losses, as corroborated by modeling [677, 690].

<figure><img src="image_40.png" /><figcaption>**Figure 39.** EGAM observed on LHD: spectrograms of (a) magnetic field fluctuations and (d) electric potential fluctuations; waveforms extracted by bandpass filters with a passband of 25 – 55 kHz [(b), (e)] and 58 – 92 kHz [(c), (f)], respectively. The bold curves show the envelopes. Reproduced from [696].</figcaption></figure>

Since EGAMs are axisymmetric, they are linearly driven by positive gradients in velocity space. The EGAM linear dispersion relation $D(\Omega) = 0$ has been derived in several works making various assumptions for the equilibrium distribution function of EPs. For instance, in a kinetic-MHD approach the thermal plasma is modelled as a fluid, and the response of EPs is computed using the perturbed distribution function obtained from the drift kinetic equation [443, 665].

A kinetic approach for both the thermal plasma and the EPs [667–669] can explain the existence of two types of EGAMs, as illustrated in figure 40, where the

zeros of the dispersion relation are represented in the complex plane. Solving the kinetic dispersion relation of the standard GAM (without EPs) results in a branch of highly Landau damped modes ($Im(\omega) < 0$) and in a mode marginally stable ($Im(\omega) \approx 0$), which is the so-called GAM. When EPs are introduced, the distribution of zeros in the complex plane is modified. Depending on the conditions of the plasma (safety factor, energy and mass of EPs) the driven mode can originate either from the standard GAM (top panel of figure 40) or from the Landau branch (bottom panel of figure 40). In the figure, the thick black solid lines with arrows represent the evolution of the excited mode when the density of EPs is increased. The dashed lines represent the evolution of a secondary mode getting closer to marginal stability. These two types of EGAM were found in gyrokinetic simulations using the ORB5 [668] and the GYSELA [690] codes, which were in very good quantitative agreement. Moreover, using a kinetic approach, other branches can be found depending on the resonance between characteristic EP frequencies and the mode, e.g. due to the resonance between EPs and the magnetic drift frequency [671]

$$ \omega_d = k_{\text{GAM}} (R\omega_{\text{ci}})^{-1} (v_{\parallel}^2 + v_{\perp}^2/2) \sin\theta, \quad (66) $$

with $k_{\text{GAM}}$ the GAM radial wavenumber, $R$ the major radius, $\omega_{\text{ci}}$ the ion cyclotron frequency and $\theta$ the poloidal angle.

Finally, in a fluid approach [670], the EPs are modelled by a distribution with a narrow width in energy, leading to the excitation of a reactive branch. Since the frequency is determined by EP kinetic effects, the EGAM is qualitatively different from the pure MHD mode GGAM observed in JET [664, 687]. The radial structure is determined by the EP drift orbit width [443, 665]. The effect of the elongation of the plasma on the linear excitation of EGAMs has been studied in a joint investigation with GENE and ORB5 [693].

The linear excitation of EGAMs has been extensively studied in numerical simulations in different frameworks. Depending on the framework, different strategies can be employed to model EPs. In full-$F$ global gyrokinetic simulations like GYSELA [697] or ORB5 [698], an initial value problem can be solved. In that case, EPs are initialized using a distribution function which depends on the motion invariants [667, 668, 690, 691]. Alternatively, an external source can be used to force the distribution function to be as close as possible to a distribution required to excite EGAMs [688, 689], which is useful to excite EGAMs in the presence of turbulence. On the other hand, one can use $\delta F$ codes such as GENE [699], where an equilibrium distribution function is imposed and does not evolve.
<figure><img src="image_41.png" /><figcaption>**Figure 40.** Evolution of the zeros of the dispersion relation (represented by thick black solid lines) in the complex plane for an increasing EP density. The top panel corresponds to the case where the EGAM is excited from the standard GAM and the bottom panel corresponds to the case where the EGAM is excited from the Landau branch. Reproduced from [669].</figcaption></figure>

## 6.7. Linear stability of EPMs

The frequencies and eigenfunctions of the EP-driven AEs discussed so far are determined by the bulk plasma, almost unaffected by the presence of EPs, but the EPs usually provide the drive through gradients in the distribution function. These EP-driven AEs are usually described by a perturbative approach. New types of modes, the energetic particles modes (EPMs) appear when the EP pressure is similar to the bulk plasma pressure, such that the EPs significantly influence the dielectric response of the plasma. The EPs then affect the very existence of the mode and not just the mode growth rate. To describe EPMs, for example fishbones, we need the more sophisticated non-perturbative approach [496]. Due to the larger EP pressure, the EP drive is often strong enough to overcome continuum damping, and so EPMs can exist at frequencies in the Alfvén continuum, and

their frequencies and drive depends strongly on EP distribution. The frequencies of EPMs are those maximizing the power transfer of the interaction between the EPs and the wave and are often given by the characteristic EP drift orbit frequencies discussed in section 2.

One of the technical difficulties, associated with the EPMs requiring a non-perturbative approach, is that their spatial structure and frequency can change significantly and rapidly in the nonlinear phase of the instability. Rapidly changing mode frequency on short timescales relative to the timescales for changes to equilibrium parameters, often on the order of milliseconds, is referred to as chirping.

Chirping events are routinely observed [700, 701], and their salient features have been reproduced qualitatively using theoretical models [702–706]. EPMs can be driven by runaway electrons [707] as well as fast ions. A quantitative understanding of experiments requires detailed numerical modeling which is challenging but should be a feasible task to be discussed in section 7.

Experimental studies of the stability of EPMs range from parametric studies of the presence of chirping modes (e.g., [700]), direct measurements of mode growth rates based on growth and decay rates of bursting modes (e.g., [629]), and indirect measures of mode growth rates (e.g., [594]) based on theoretical frequency-chirping rates [702, 703]. Finally, theoretical predictions of methods to reduce the size of mode bursts and frequency chirps by scattering resonant particles have been tested experimentally [708–711].

Estimates of growth rates based on the frequency chirping of EP-driven instabilities and their related EPM branches have been made for TAEs [712–714], GAEs [594], and ICE [629]. The correlation of mode amplitude and chirping rate has been studied on START, MAST and NSTX [706, 712, 715]. The parameter dependence of chirping TAEs has been characterized on NSTX [712]. The TAE instability is excited by the free energy in the EP distribution present in gradients in radius or energy. Flattening the EP distribution, e.g., by off-axis NBI, is often considered a method to suppress TAEs. However, experiments have found that far off-axis co-current NBI, that is when creating a hollow EP profile, can destabilize counter-propagating TAEs [716, 717]. More generally, it is shown in [718, 719] that especially in situations with large orbit width, the mode drive (equation 45) needs to be evaluated carefully, as even non-inverted EP profiles may lead to the destabilization of counter-propagation TAEs.

Similarly, on AUG, off-axis NBI heating without any central NBI heating leads to a hollow EP profile with large ratios of $\beta_{NBI}/\beta_{th} \sim 1$ and $\epsilon_{NBI}/T_{i,e} \sim 100$.
In this regime, nonlinear coupling signatures between co- and counter-propagating AEs and their corresponding EPM branches have been observed in the ramp-up and flat-top phases, as well as the coupling of zonal modes (EGAMs) and AEs/EPMs [684, 685, 716, 717]. Linear and nonlinear modeling and code validation efforts are ongoing [570, 692, 720, 721].

## 6.8. Impact of anisotropy and toroidal flows

The introduction of auxiliary heating in tokamak plasmas can introduce toroidal and poloidal flows and pressure anisotropy as well as change the magnetic configuration. In this subsection, we discuss how these effects can change the linear MHD stability of the plasma. Since NBI and ICRF heating generate highly anisotropic EP distribution function with a kinetic EP pressure comparable to the bulk plasma pressure, the total plasma pressure can become anisotropic. NBI heating in MAST can cause anisotropy up to $p_{\perp}/p_{\parallel} = 1.7$ [723], and ICRF heating in JET can cause anisotropy up to $p_{\perp}/p_{\parallel} = 2.5$ [724], where $p_{\parallel}$ and $p_{\perp}$ are the parallel and perpendicular pressures. At AUG, the anisotropy was inferred to be $p_{\perp}/p_{\parallel} = 1.2$ in an NBI heated plasma based on velocity-space tomography [725]. However, most MHD equilibrium codes do not account for anisotropy and solve the usual static isotropic Grad-Shafranov equation. MHD equilibria with anisotropy and flows are described by the generalised Grad-Shafranov equation, written as (e.g. [726] and references therein)

$$ \nabla \cdot \frac{(1-\Delta)\nabla\psi}{R^2} = -\frac{F(\psi)F'(\psi)}{(1-\Delta)R^2} - \mu_0\rho \\ \times \left[ k_B T_{\parallel}^{\prime}(\psi) + H^{\prime}(\psi) + R^2\Omega(\psi)\Omega^{\prime}(\psi) - \left(\frac{\partial W}{\partial \psi}\right)_{\rho,B} \right] \quad (67) $$

where five constraints are used: $F(\psi)$, $\Omega(\psi)$, $H(\psi)$, $T_{\parallel}(\psi)$, and $\Theta(\psi)$. Here, the anisotropy is reflected in the flux function $\Theta(\psi) = B(1/T_{\parallel} - 1/T_{\perp})$, where $T_{\parallel}$ and $T_{\perp}$ the parallel and perpendicular temperatures of a plasma with a bi-Maxwellian velocity distribution function [727].

Several studies have examined the impact of anisotropy and flows on the plasma equilibrium (for an analytic approach, see [728]). For MAST equilibria, the $J_{\parallel}$ contribution can reach 20% of the total current for $p_{\parallel}/p_{\perp} = 1.5$ [722]. This is significant as it results in a 10% change in the current profile and consequently the q-profile, which influences the plasma stability. Figure 41 shows the differences in the equilibrium reconstruction using the EFIT TENSOR code for MAST.

EFIT TENSOR [729] solves the generalised Grad-Shafranov equation, equation (67) using magnetic and kinetic constraints. The unique solution to the generalised Grad-Shafranov equation is determined by the plasma boundary contour $C_A(R, Z)$ and the profiles of the flux functions $\{F(\psi), \Omega(\psi), H(\psi), T_{\parallel}(\psi), \Theta(\psi)\}$. A recent study examined the impact of toroidal flow and anisotropy for the same magnetic configuration (q-profile) and plasma stored energy [726] using the fixed boundary solver HELENA+ATF [722]. As the anisotropy is varied for a constant q-profile, the most striking observation of the equilibrium profiles are that contours of constant $p_{\perp}, p_{\parallel}$ and $\rho$ shift outboard of the magnetic axis for $T_{\perp} > T_{\parallel}$ and inboard of the magnetic axis for $T_{\perp} < T_{\parallel}$. Figure 42 shows this shift as a function of $\Theta = 1 - T_{\parallel}/T_{\perp}$.

A recent study [730] extends pressure anisotropy and toroidal flow equilibrium models to include the EP orbit width. An equilibrium model with q-solver constraint enabled a wider systematic study of the parametric dependencies for the same q-profile. The addition of finite orbit width effects reduces the shift in magnetic axis, reducing the change in the TAE gap frequency from the centre of the gap, reducing continuum damping and thereby destabilising the TAE.

The impact of pressure anisotropy can be investigated using kinetic energy principles [731] [732]. Chew, Goldberger and Low (CGL) [733] introduced the now widespread form of the pressure tensor and derived the double-adiabatic CGL closure. The corresponding energy principle was later derived and studied [731, 734, 735]. CGL ignores the heat flow when the mode frequency is comparable to or smaller than the particle streaming frequency, especially in the vicinity of marginal stability boundary [736]. Alternative fluid closures without this drawback are, e.g., the double polytropic closure, a higher-order-momentum closure, and, recently, the single adiabatic closure [729] which produces the same result as MHD for an isotropic equilibrium.

MHD linear stability was computed [726] using MISHKA-AD [732], which uses the CGL closure and the single adiabatic closure [729]. As the q-profile remains unchanged, the shear Alfvén continuum, the gap modes and the gap frequencies remain largely unchanged. However, the compressional continuum changes significantly. Figure 43 illustrates this change by plotting the continuum for $T_{\parallel}/T_{\perp} = 0.8, 1.2$.

The equilibrium and stability of $n = 1$ TAE modes in MAST were computed using EFIT-TENSOR and MISHKA-AD [737], and here substantial changes in the equilibrium compared to the isotropic EFIT were found, and consequently the stability was significantly affected. The anisotropic solution had a q-profile with
<figure><img src="image_42.png" /><figcaption>**Figure 41.** (a) Pressures on the mid-plane in the anisotropic reconstruction (two pressures, shown with the solid and dashed-dotted lines) and in the isotropic reconstruction ($p^*$, shown with the dotted line) for MAST discharge #18696 at 290ms. (b) The reconstructed $J_\phi$ profile and $q$-profile on the mid-plane. Reproduced from [722].</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>**Figure 42.** Calculation of the displacement of peaks in density $\rho$, parallel pressure $p_\parallel$ and perpendicular pressure $p_\perp$ in MAST from the magnetic axis as a function of $\Theta = 1 - T_\parallel/T_\perp$. Reproduced from [726].</figcaption></figure>

<caption>**Figure 43.** Compressional Alfvén continuum for $n = 1$ with $T_\parallel/T_\perp = 0.8$ (blue) and $T_\parallel/T_\perp = 1.2$ (red). Reproduced from [726].</caption>

reversed shear, the $n = 1$ shear Alfvén continuum gap opening near the core, and the core safety factor lowering, resulting in a broader TAE mode. Mode drive was computed the wave-particle interaction code HAGIS [585]. The resonant regions of phase space were significantly modified between the isotropic and anisotropic cases. The anisotropic $l = 1$ bounce resonance shifted radially inward for given particle energy relative to the isotropic case. The linear growth rate in the anisotropic case was 35% larger than in the isotropic case (see figure 44), and the saturation amplitude in the anisotropic case was 18% smaller than in the isotropic case. The linear growth rate was larger in the anisotropic case because the radial gradient of the distribution function at resonance was larger. The saturation amplitude in the anisotropic case was slightly smaller, despite the larger linear growth rate. This may be caused by the differences in the safety factor and magnetic shear affecting the

nonlinear bounce frequency.

Recently, the impact of anisotropy on the stability of infernal modes, driven when $q \approx 1$ over a wide region, was studied [738]. Infernal instabilities are instabilities in plasmas with weak shear, such that bands of unstable $n$-values can form, even when standard ballooning theory predicts stability. The guiding centre plasma motion was described by fluid equations, whereas the motion parallel to the magnetic field was described by a collisionless kinematic equation. In this model, the plasma anisotropy enters through a modification to the magnetic well, yielding better stability for tangential injection. A stability criterion for the linear stability of plasma equilibria with incompressible flow parallel to the magnetic field, constant mass density and constant $\sigma_d = \mu_0(P_\parallel - P_\perp)/B^2$ can be found in [739].

In all cases studied, the change in equilibrium affects the change in plasma stability. In ITER,
<figure><img src="image_44.png" /><figcaption>Figure 44. Wave amplitude $\delta B/B_0$ for the isotropic (red solid line) and anisotropic (blue dashed line) cases. Reproduced from [737].</figcaption></figure>

strong flows or anisotropy are not anticipated to have a significant effect on the equilibrium reconstruction, and uncertainty in the experimental inputs to the equilibrium reconstruction will likely predominate. Nevertheless, such effects may play a role in other burning plasma devices.

## 7. Nonlinear mode evolution, theory and simulation

This section describes progress made in the field of nonlinear EP physics since the topical review in the Progress of the ITER Physics Basis [20]. The build-up of the EP population in fusion plasmas is typically slow compared to the growth times of EP-driven instabilities. This scale separation suggests that we need nonlinear studies of unstable waves in the near-threshold regimes. We need to characterize the long-time behavior of the waves and resonant particles in the presence of particle sources and sinks. In some observations of EP-driven instabilities and EPMs, the unstable modes grow to a level at which they enhance transport and cause anomalous losses of the EPs. In other observations, the losses are small, but the modes exhibit a complex nonlinear behavior, such as the generation of sidebands, quasi-periodic bursts, or changing mode frequencies. The fishbone instability exhibits such features, which we will describe in the following subsection 7.1. This section presents a first-principles physics basis for modeling these phenomena as well as reduced models.

As explained in section 2, the resonance condition for interaction of EPs with AEs in a tokamak is

$$ \omega - n\omega_{\phi}(\mathcal{E}, \mu, P_{\phi}) - l\omega_{\theta}(\mathcal{E}, \mu, P_{\phi}) = 0. \quad (68) $$

Here we state the arguments of the toroidal and poloidal frequencies $\omega_{\phi}(\mathcal{E}, \mu, P_{\phi})$ and $\omega_{\theta}(\mathcal{E}, \mu, P_{\phi})$ of

the unperturbed ion motion explicitly to highlight that equation (68) selects a surface in the 3D phase space $(\mathcal{E}, \mu, P_{\phi})$. A single resonance in $(\mathcal{E}, \mu, P_{\phi})$ thus covers a relatively thin region of phase space, and it causes transport only across this thin surface, which can flatten the distribution function only locally. But when the phase space is covered by overlapping resonances, there can be a significant transport of particles via stochastic diffusion [740].

The nonlinear simulation of EP-driven modes remains a significant challenge for future devices. In these systems it is expected that a dominant component of alpha particle transport will be driven by interactions with various instabilities rather than classical orbit losses. There are several issues that make modeling the fusion reactor regime difficult. These include the large timescale separation between the instabilities (on the order of microseconds) and the alpha slowing-down times (on the order of a second). Next, disparate spatial scales are involved ranging from the ion gyroradius, (on the order of 0.01 meter) to the device size, (on the order of several meters). Finally, there are no existing experiments that operate in the reactor parameter regime, so there is no possibility for model validation until DT operations in ITER or elsewhere.

These challenges have resulted in a variety of approaches, ranging from global full-f gyrokinetic, fluid-kinetic, gyrofluid, quasilinear, as well as various semi-analytic approaches. Some of the physics issues that must be addressed include: mode-coupling (wave-wave) nonlinearities vs. wave-particle nonlinearities, global vs. local (flux tube), interaction of alpha particles with both Alfvén instabilities and core microturbulence, inclusion of multiple fast ion species (beams, ICRF, alphas), roles of zonal flows and currents, perturbative quasilinear vs. non-perturbative nonlinear, collisionless vs. collisional, absence vs. presence of external sources and sinks, and linear critical gradient profiles vs. nonlinearly flattened profiles. Currently, there is no single approach that includes all known relevant effects and is sufficiently computationally efficient and robust to run for even a small fraction of an alpha slowing down time. All of them involve some degree of simplification, but each tries to address at least some of the above challenges in ways that should be complementary to the others.

In this section some of the nonlinear simulation models currently in use will be reviewed. We will start with a description of fishbones, which is a prominent example of an instability requiring a nonlinear description. We will then progress from gyrokinetic theory and simulations, which is the most fundamental description, to reduced descriptions requiring an increasing amount of modeling. We will review nonlinear gy-
rokinetic theory of wave-particle and wave-wave inter-
action, discuss progress in nonlinear gyrokinetic sim-
ulation, nonlinear kinetic-MHD simulations, gyrofluid
model simulations, and nonlinear simulations with re-
duced models. Development is continuing, and it is
expected that increasingly comprehensive and sophisti-
cated simulations will be available by the time of ITER
DT.

7.1. *Fishbones*

Historically, the so-called fishbone was the first
EPM observed experimentally. It was discovered in
experiments with NBI heating on the PDX tokamak
[741] and then found in most tokamaks with EP
populations (either ions or electrons). Figure 45
shows how this perturbation appears on a Mirnov
coil measurement; the similarity to a fish skeleton
for multiple bursts resulted in the name “fishbones”.
Such bursts can enhance EP losses [22, 741, 742],
as also observed for fusion alphas in the recent
DT experiments and JET [47], and deteriorate the
confinement of the thermal plasma via triggering kink
modes and neoclassical tearing modes (see section 5)
[592, 743–747].

<figure><img src="image_45.png" /><figcaption>**Figure 45.** Oscillations of the perturbed magnetic field time derivative during a fishbone burst (JET discharge #16341). The gradual lengthening of the period associated with frequency chirping is visible to the eye. Reproduced from [748].</figcaption></figure>

Fishbones are an EPM that has the mode
structure of an (m/n) = (1/1) internal kink. They are
usually observed in a repetitive bursting cycle when the
central safety factor q(0) has fallen below unity when,
in the absence of EPs, the MHD internal kink growth
rate γMHD drives instability. As discussed in section 5.2
and [372], because the EP precession frequency is

large compared to the MHD mode frequency, EPs
on trapped-particle orbits stabilize the internal kink
that triggers the sawtooth crash. However, when the
EP population becomes too large, the EPs destabilize
a new branch of the internal kink, the fishbone
(figure 46a). In present-day tokamaks with so-called
"hybrid" scenarios of improved performance [750],
fishbones often appear when q(0) ≳ 1.

The fishbone bursts involve oscillations with
a time-dependent frequency, as seen in figure 45.
A qualitative physics explanation of the downward
frequency chirping for fishbones has been given in
[28]. There are two kinds of internal-kink fishbone
instability. The first are the so-called “precessional”
fishbones [751] appearing when the mode frequency
ω resonates with the bounce-averaged precession
frequency of the trapped EPs, ω<sub>φ,tr</sub>, and it is much
greater than the thermal ion diamagnetic frequency,
ω<sub>*i</sub>. The precessing EPs then destabilize the (1/1)
mode in the Alfvén continuum. At the onset of
the fishbone pulse, the drive from the EPs is almost
balanced by the continuum damping coming from two
Alfvén resonances near the q = 1 surface: ω =
±k<sub>‖</sub>(r<sub>A</sub>)v<sub>A</sub>(r<sub>A</sub>) (see section 6). To overcome strong
continuum damping, the precessional fishbones require
relatively high values of the EP pressure. The
fluid nonlinearity at the Alfvén resonances tends to
reduce the continuum damping. This is the dominant
nonlinear effect at the onset of the fishbone pulse. It
leads to an explosive growth of the pulse so that the
precessional fishbones follow a hard excitation scenario
[752].

The second kind of fishbones are those with ω ≈ ω*<sub>i</sub> [753, 754]. This mode lies in a low-frequency gap in the shear Alfvén continuum and consequently is not damped by continuum damping. These fishbones are one of the two oscillatory kink modes stabilized by FLR effects. These modes are unstable within the framework of ideal MHD [755].

The explosive growth of the fishbone oscillations
and the significant change of the oscillation frequency
call for a systematic nonlinear description of fishbones
with a self-consistent treatment of kinetic and MHD
nonlinearities [752]. This is a challenging technical
issue for numerical modeling. One of the difficulties
here comes from the need to incorporate an accurate
description of the narrow phase-space resonances into
global MHD simulations. For linear problems, this
difficulty is only a moderate obstacle since the resonant
response of the system is insensitive to the width of
the resonance and can be treated in terms of Landau
damping. In contrast, nonlinear problems typically
need much better resolution to calculate the resonant
response appropriately, which is very demanding for
the existing global codes. Several attempts have
<figure><img src="image_46.png" /><figcaption>**Figure 48.** (a) Magnetic probe and (b) neutron signals during an off-axis fishbone in DIII-D. As with the classic fishbone of figure 45, the gradual lengthening of the period associated with frequency chirping is evident to the eye but, unlike the classic fishbone, the magnetic waveform becomes increasingly distorted by higher frequency harmonics as the burst evolves. The coincident drop in neutron rate is caused by ejection of trapped EPs. Reproduced from [763].</figcaption></figure>

with mode amplitude for both [763, 765]. The non-ambipolar EP loss causes a coincident, rapid change in radial electric field [763].

MEGA simulations that treat the EPs kinetically and the bulk plasma with resistive MHD have successfully reproduced an off-axis fishbone burst [766] and the nonlinear distortion of the waveform caused by coupling to higher harmonics [767]. In both theory [768] and experiment [769, 770], transport of EPs by off-axis fishbones can trigger ELMs.

## 7.2. Nonlinear gyrokinetic theory

Nonlinear gyrokinetics [771–773] provides the most comprehensive framework for addressing EP physics in reactor relevant fusion plasmas. It accounts for the resonant excitations of instabilities, the ensuing transport processes as well as the broad range of spatio-temporal scales associated with their nonlinear dynamics. Employing gyrokinetic theory is not only necessary for the correct assessment of wave-particle resonant interactions, which provide crucial driving and damping mechanisms, but is also mandatory for dealing with the short scales that are spontaneously produced by the spatial phase mixing of the SAW continuous spectrum [774]. KAWs [500, 775], excited by mode conversion at resonances with the SAW continuum, enhance nonlinear wave-wave couplings and modify the spectral features of the fluctuation spectrum qualitatively and quantitatively [776–778].

Nonlinear gyrokinetics is the foundation of the unified and self-contained theoretical framework describing these physics, as discussed in a recent review [32]. The self-consistent theory treats the fluctuation spectrum evolution and EP transport on the same footing. The first aspect is dealt with by the so-called general fishbone-like dispersion relation [634, 779],

which consists of the weak formulation of nonlinear gyrokinetic quasi-neutrality and vorticity equations. This can be seen as a gyrokinetic energy principle, valid in a wide frequency interval, ranging from the low-frequency MHD up to the Alfvén wave frequencies. It has been systematically verified numerically and validated experimentally [32, 634, 779]. In the high-n (toroidal mode number) limit, the general fishbone-like dispersion relation has the form

$$i\Lambda(\omega, r) = \delta\bar{W}_f(\omega, r, k_r) + \delta\bar{W}_k(\omega, r, k_r), \quad (69)$$

where $\Lambda(\omega, r)$ is the generalized inertia, depending on the mode frequency and the radial coordinate. $\delta\bar{W}_f(\omega, r, k_r)$ and $\delta\bar{W}_k(\omega, r, k_r)$ are the fluid and kinetic potential energies, which also depend on the radial wave vector. Equation (69) can be cast as a nonlinear Schrödinger-like equation, describing the amplitude evolution of short wavelength Alfvénic fluctuations [780]. The EP transport description is based on the derivation of the renormalized EP response, i.e., the EP distribution function in the presence of a finite fluctuation level [781–783]. The evolution equation for the renormalized EP response self-consistently accounts for EP transport in phase space and is cast in the form of a Dyson-like equation [781]. Summarizing, a key point is that there are two “routes” to nonlinear dynamics of EP-driven Alfvénic fluctuations in magnetized fusion plasmas [32]; i.e., nonlinear wave-wave and wave-particle interactions [777]. Each of these routes is discussed in the following two sub-sections.

### 7.2.1. Nonlinear wave-wave interactions.
Nonlinear wave-wave interactions have been first addressed by Alfvén [784], who demonstrated that, in uniform, incompressible ideal MHD plasmas, SAWs can exist in the pure “Alfvénic state”, independently of their amplitude, due to the cancellation of the Reynolds and Maxwell stresses [785, 786]. The pure “Alfvénic state” is closely linked, in realistic nonuniform plasmas, with the existence of the SAW continuum spectrum (section 6). As a consequence, the equilibrium magnetic field geometry and the plasma nonuniformity play crucial roles in the nonlinear dynamic evolution of the system [32] along with previously mentioned short-scale kinetic effects, plasma compressibility as well as deviation of the Alfvén eigenmode frequency from the local continuum [777, 787, 788]. Therefore, the nonlinear gyrokinetic approach is crucial to qualitatively and quantitatively assess the nonlinear wave-wave coupling process. In fact, it may also be applied to the polarization nonlinearity at long wavelengths as well as to Reynolds and Maxwell stresses in the short wavelength kinetic regime [776, 777]. Among the various nonlinear wave-wave interactions, excitation of zonal field structures (ZFs)
[780, 789] and frequency cascading via ion induced scattering [790, 791] are two channels expected to significantly influence the SAW instability nonlinear dynamics in fusion plasmas. These physics have been analyzed taking TAEs [792] as an example. That general approach, developed for TAEs, can be applied to other SAW instabilities based on their corresponding linear properties [793–795].

ZFs are toroidally and predominantly poloidally symmetric variations of scalar and parallel vector potentials and are connected with the radial corrugations of the equilibrium profiles such as zonal flows and fields/currents (ZF & ZC) [796, 797]; in other words, zonal flow and fields/currents are low-frequency $n=0$ and predominantly $m=0$ structures characterized by appreciable radial wavenumber $k_r$. ZFs may regulate fluctuations by radial mode structure shearing and scattering the fluctuations into the short-wavelength stable domain, as typically occurs for drift wave turbulence [797]. In reactor relevant fusion plasmas, EPs dominate the local power balance and, thus, play a unique role as mediators of cross-scale couplings [798] as they may linearly and nonlinearly excite ZFs, acting, thereby, as generators of nonlinear equilibria, or zonal states (ZS), that generally evolve on the same timescale as the underlying fluctuations [781–783, 798].

Nonlinear excitation of zero frequency ZFs by TAEs was investigated in [789], which found that both electrostatic zonal flows and electromagnetic zonal current can be excited. The former results from noncancellation of Reynolds and Maxwell stresses due to toroidicity, while the latter is related to dynamo effects. The condition for the spontaneous excitation of zonal field structures is more easily satisfied when the neoclassical polarizability enhancement by trapped particles is properly accounted for [799], causing the zonal current to be preferentially excited with a much smaller threshold on the TAE amplitude being $|\delta B_r/B_0| \sim \mathcal{O}(10^{-4})$, compatible with the fluctuation amplitude in tokamak experiments [800].

Due to the typically weak ballooning feature of SAW instabilities, zonal field structures excited by TAEs have a fine radial scale structure [793, 801], in addition to the usual well-known mesoscale structures. This additional fine radial structure may lead to enhanced nonlinear coupling and a corresponding much lower TAE saturation level. It was further demonstrated that the nonlinear contribution of resonant EPs to plasma compressibility may dominate the thermal plasma contribution to the Reynolds and Maxwell stresses. In that case, zonal flows may be forced driven [793, 802] rather than spontaneously excited by finite amplitude TAEs, with the zonal flow growth rate being twice of the instantaneous TAE growth rate [802], as observed in kinetic-MHD [803]

as well as PIC simulations [804]. A thorough discussion of zonal field structures generation by TAEs is given in [801], with emphasis on the various underlying physics, e.g., zonal flow vs. zonal current generation, the roles of fine- vs. mesoscale structures, as well as spontaneous decay vs. forced driven processes. Further important applications of the same theoretical framework to, e.g., zonal flow generation by BAEs [700, 793, 805] as well as finite frequency GAM [456] generation by TAEs [806–808] have also been investigated. Simulations find that zonal flows and currents play an important role in the saturation of AEs are discussed further in section 7.3.1.

Another important channel of wave-wave coupling is TAE frequency cascading via thermal-ion induced scattering, also called nonlinear ion Landau damping [809]. In reactor-relevant realistic geometry, there exists $\mathcal{O}(n^2q)$ TAEs, with $n \gtrsim \mathcal{O}(10)$ [810] (section 6). Thus, many TAEs with overlapping radial structures and slightly different frequencies co-exist. The TAE spectral cascading was first investigated using a drift kinetic theory [790], and then generalized using a nonlinear gyrokinetic approach [791]. In a single process, a test TAE couples with a counter-propagating fluctuation of the TAE spectrum with slightly lower frequency and generates a low frequency ion-sound-wave quasi-mode, with ion Compton scattering and shielded-ion scattering contributing to the process on the same footing. The nonlinear evolution of the TAE spectrum can then be cast as wave-kinetic equation by summing all the strongly interacting background TAEs in the continuum limit. This describes the downward spectral energy transfer and TAE saturation due to enhanced continuum damping. The nonlinearly saturated TAE spectrum can be derived from the wave-kinetic equation, which yields the estimated TAE saturation amplitude as $|\delta B_r/B_0| \lesssim \mathcal{O}(10^{-4})$; i.e., typically one order of magnitude smaller than that in the MHD limit [790] due to enhanced nonlinear coupling in the kinetic regime [777]. The resulting transport of circulating EPs can be derived from quasi-linear transport theory [811], and be estimated in the range of $\sim 1 - 10 \text{ m}^2/\text{s}$ for typical reactor parameters [791]. The resulting bulk ion heating rate can also be derived [812]. Specific applications to ITER to quantitatively assess the impact of wave-wave interactions on AE saturation and ensuing EP transport will require dedicated numerical simulations accounting for all the various nonlinear processes mentioned above, since they all compete on the same footing and estimating their effect independently would produce an incorrect prediction.

### 7.2.2. Nonlinear wave-particle interactions.
Nonlinear wave-particle interactions between Alfvénic fluctu-
ations and EPs are an essential element for the assessment of EP transport in fusion plasmas. The theoretical framework discussed in [32] suggests that burning plasmas in ITER will require lifting the transport description to phase space in order to capture the complex nonlinear behaviors due to the many interacting degrees of freedom [781, 798]. A new approach has been recently proposed within this framework [782, 783] to effectively compute the self-consistent evolution of plasma equilibria and fluctuation spectra on the energy confinement time scale. This approach is consistent with the general nonlinear gyrokinetic description and is verified and validated in a number of proof-of-principle cases. First results are briefly discussed in the following, showing that they provide a practical tool to make feasible predictions in realistic ITER plasma conditions [813, 814] (see also section 6 for recent results on ITER simulations). A hierarchy of reduced models with verifiable fidelity can be constructed which should lead to applications of increasingly higher sophistication and reliability before ITER operations, using numerical simulation tools that are integrated within the IMAS framework.

Wave-particle interactions must be described in phase space, where resonant processes dominate due to the relatively small fluctuation levels ($|\delta B_r/B| \lesssim 10^{-4}$). For the proper definition of nonlinear EP equilibrium evolving on the spatio-temporal mesoscales, the concept of phase-space zonal structures (PSZS) has been introduced [32, 781, 783, 798], that is, the part of the distribution function remaining undamped by collisionless processes on the characteristic non-linear timescale. The PSZS, accounting for the “renormalized” (nonlinear/evolving) equilibrium distribution function, depend only on the invariants of motion, $\bar{F}_z(0, \mathcal{E}, \mu, P_\phi; t)$ Introducing the magnetic-drift/banana center pull-back operator $e^{-iQ_z}$, with $Q_z = g(\psi)(v_\parallel/\Omega)k_z/(d\psi/dr)$ and $k_z \equiv -i\partial_r$ the zonal field structures radial wave-number, the evolution equation for the phase-space zonal structures is given by [780, 782, 783]

$$ \partial_t e^{iQ_z \bar{F}_0} = - \frac{e^{iQ_z} F(\psi)}{B_0} \partial_t \langle \delta A_{\parallel g} \rangle_z \frac{\partial}{\partial \bar{\psi}} \bar{F}_0 + \frac{e^{iQ_z} [C_g + S]}{\tau_b} - \frac{1}{\tau_b} \frac{\partial}{\partial \psi} \left[ \tau_b e^{iQ_z} \delta \dot{\psi} \delta F \right] - \frac{1}{\tau_b} \frac{\partial}{\partial \mathcal{E}} \left[ \tau_b e^{iQ_z} \delta \dot{\mathcal{E}} \delta F \right]. \quad (70) $$

Here, $\tau_b$ denotes bounce/transit time$^+$, $[\dots]$ denotes

bounce/transit averaging and, thus, $e^{iQ_z} [\dots]$ denotes orbit averaging [780]. The first term on the right hand side represents the nonlinear equilibrium change associated with the zonal field structures due to the gyro-averaged parallel vector potential, $\langle \delta A_{\parallel g} \rangle_z$. $\delta$-quantities stand for fluctuations, while $C_g$ and $S$ on the right hand side represent gyro-center collision and source terms [815–818]. Equation (70) can generate all spatio-temporal scales, from micro- via meso- to macroscales. For the sake of convenience and numerical implementation, it is often useful to separate slow meso- and macroscale variations, which define the mean equilibrium evolution, from the fast scales that describe the zonal state deviations about it, as discussed above in section 7.2.1 [780, 782, 783]. In addition, it serves as fundamental equation defining EP transport in the phase space. It provides the basis for any reduced approach that seeks to describe EP dynamics over transport timescales (for a practical implementation see e.g. [813, 814]). In these works, various simplifying assumptions can be adopted in the quantitative evaluation of the phase space fluxes on the RHS. The crucial element is that all these models can be readily verified/falsified within a unified theoretical framework [32, 783] by means of nonlinear GK codes that incorporate specialized diagnostics, as demonstrated by the work of [819]. Meanwhile, verification/falsification is based on comparisons of fluctuation spectra as well as the corresponding phase space fluxes, which provide a credible assessment of the predictive capability of self-consistent evolution of fluctuation spectra and equilibrium profiles.

The approach based on equation (70) and the self-consistent solution of the Alfvénic fluctuation spectrum allow describing the non-perturbative dynamics of EP nonlinear equilibria that generally evolve on the same timescale $\tau_{NL}$ of the underlying fluctuations; i.e., $\gamma\tau_{NL} \sim \omega_B\tau_{NL} \sim \mathcal{O}(1)$, with $\gamma$ the instantaneous (nonlinear) growth rate and $\omega_B$ the wave-particle trapping frequency. By construction, thus, this approach can handle both weak and strong drive regimes and naturally recovers the approaches that assume $\omega_B\tau_{NL} \gg 1$ near marginal stability, as shown in [32, 781, 820].

Equation (70) reduces to the quasilinear diffusion equation for a broad spectrum of overlapping resonances [821, 822]. More generally, for a narrow spectrum of quasi-coherent fluctuations, it describes a broader class of transport processes [823, 824] and can be cast into the form of a Dyson-like equation, which accounts for convective transport and avalanches in the EP phase space [32, 781]. A particular application of practical interest of the Dyson-like equation, coupled with the nonlinear Schrödinger-like equation for

section 2; we name it $\tau_b$ here to comply with the cited references.

<footnote>¶ The COMs are here defined per unit mass, i.e. $\mathcal{E} = v^2/2$, $\mu \simeq v_\perp^2/2B_0$ the leading order expression of the magnetic moment and $P_\phi \simeq (e/c)(g(\psi)(v_\parallel/\Omega) - \psi) \equiv -(e/c)\bar{\psi}$ the leading order toroidal canonical angular momentum, where $g(\psi) = RB_t$ is the covariant component of the toroidal magnetic field and $\Omega = eB/(mc)$ is the cyclotron frequency.</footnote>

<footnote>+ The bounce/transit time $\tau_b$ is the same as $\tau_\theta$ defined in</footnote>
the self-consistent description of the fluctuation amplitude (section 7.2), is the nonlinear dynamics of the EPMs [421], (section 6) which is shown to yield convective amplification of the EPM wave packet as a soliton leading to an EP avalanche. This case was discussed earlier in [20], while the EPM chirping rate, $\dot{\omega} \sim \omega_B^2$, and the radial structure of the EPM wave packet are analyzed in [32, 781]. Extensive verifications of the linear scaling of chirping rate with the mode frequency have been recently reported in [825, 826]. In particular, it can be shown that the self-similar shape of the EPM envelope, $A(r,t) = U(\xi)e^{\int^t \gamma(\tau)d\tau}$, with $\xi = k_{NL}(r - r_0 - \int^t v_g d\tau)$, $k_{NL}$ the EPM nonlinear wave number, $r_0$ the location of the linear instability, and $v_g$ the EPM group velocity, obeys the nonlinear equation [32]

$$ \partial_{\xi}^{2} U = \lambda U - 2iU|U|^{2}, \qquad (71) $$

with $\lambda = -\sqrt{2}/3 + i(4/3)$, which is a particular case of the complex Ginzburg-Landau equation [32]. In this case, the linear scaling of the theoretical prediction of frequency chirping [32, 781], which can be cast as

$$ \dot{\omega} \approx \pm \frac{1}{2} \omega_B^2 \qquad (72) $$

$$ \omega_B^2 = \left| en \frac{\partial \omega_{\text{res}}}{\partial P_{\phi}} + ew \frac{\partial \omega_{\text{res}}}{\partial E} \right| \left| \frac{\mathbf{v}_d \cdot \nabla \delta \phi}{\omega} \right| \qquad (73) $$

Here, $\omega_{\text{res}}$ is the wave particle resonance frequency and $\delta\phi$ is the EPM scalar potential fluctuation, while (...) indicates equilibrium orbit averaging as defined above [783]. The $\pm$ sign shows that both up- and downward chirping are possible, although the downward chirping is the most typical case due to equilibrium non uniformity [32]. Another successful application of the Dyson-like equation is the fishbone burst cycle [32], which is consistent with recent nonlinear kinetic-MHD simulation results [827, 828] and proves the secular loss mechanism conjectured by [829], reducing, in the proper limit, to the model equations originally introduced by Chen-White-Rosenbluth [751]. This fundamental understanding of the secular loss of EP to be expected in ITER and reactor relevant fusion plasmas imposes the necessity of identifying operation regimes where these phenomena are mitigated in suitably designed operation scenarios and/or, possibly, even controlled by means of the interplay of ZFs and PSZS [759]. In this case, the fishbone downward frequency chirping can be obtained by the same expression given above, with a proper extension to take the effect of zonal flow decorrelation properly into account [830].

The Dyson-like equation approach to self-consistently compute the EP response is not restricted to the strong EP drive case, but it properly recovers the weak EP drive in the relevant limit [32,781]. One of such limiting cases is investigated in [831] and applied

to the fishbone chirping where it was interpreted as a Doppler shift, which is self-generated by the fluctuation itself. In the light of equation 70, this highlights the importance of self-interactions in the plasma nonlinear dynamic evolution [32, 781, 783, 798].

### 7.3. Nonlinear gyrokinetic simulations of EP-driven instabilities and EPMs

Effects of collective EP instabilities (AE/EPM) on the EP confinement depend ultimately on the self-consistent nonlinear evolution of the fluctuation spectrum. More precisely, this self-consistent evolution will be reflected by the properties of the zonal state (ZS), introduced above, which consists of the ZFs and corresponding PSZS, which depend critically on the nonlinear phase-space dynamics of EPs as well as on the nonlinear mode-mode couplings among the multiple EP-driven modes, which are mesoscale (defined as fast ion gyroradius $\rho_{Lf}$) with intermediate toroidal mode numbers, typically $n \sim (10, 40)$ in ITER plasmas [31, 575]. Both nonlinear effects, in turn, depend on the global features of wave-particle resonances and mode structures whose description requires to accurately account for kinetic effects of thermal particles; e.g., the existence of kinetic Alfvén waves. Furthermore, wave-wave interaction and resonance overlap in EP phase space will induce cross-scale coupling between AE/EPM and macroscopic MHD modes such as the $n = 1$ fishbone instability. Finally also, the coupling between EP-driven turbulence and the ubiquitous drift-Alfvén wave turbulence driven by thermal particles, which are micro-scale (defined as thermal ion gyroradius $\rho_{Li}$) with high toroidal mode numbers, typically $n \sim (100, 200)$ in ITER [581] needs to be retained. A unified simulation model treating consistently macro-, meso- and microturbulence is needed to explore new physics frontiers associated with the complex dynamics of cross-scale couplings (see section 7.2) [32].

For fully self-consistent simulations, we highlight that they therefore must incorporate three physics elements: kinetic effects of thermal particles, nonlinear interactions of many mesoscale modes, and cross-scale couplings of macro-meso-microturbulence. The large dynamical ranges of spatial-temporal processes further require global simulation codes that are efficient in utilizing massively parallel computers at the exascale. Global gyrokinetic simulation [773, 832] is a suitable approach. Since the publication of the EP chapter in the Progress in the ITER Physics Basis [20], several gyrokinetic codes with comprehensive physics and realistic geometry have been developed, verified, and partially validated [555, 560, 584, 699, 833–836]. In particular, for PIC codes these advances were facilitated by the development and implementation of
advanced noise control [837] and 'pull-back' schemes [833,838-840]. These nonlinear gyrokinetic simulations are providing more complete physics insights on nonlinear dynamics of EP instabilities regarding their saturation by ZF, regulation of EP-driven turbulence by microturbulence, and fast frequency chirping. Improved understanding from gyrokinetic simulations provides physics foundation for reduced EP transport models. For example, effects of zonal flows have been incorporated in the critical gradient model (CGM) and resonance broadening quasilinear (RBQ) models [574] and effects of microturbulence have been implemented in the RBQ model [841].

### 7.3.1. Saturation of EP instabilities by wave-particle and wave-wave nonlinearities.
Since EP nonlinear dynamics depends linearly on the wave amplitude while nonlinear wave-wave coupling depends on the wave intensity, it is generally believed that wave-particle nonlinearity dominates over wave-wave nonlinearity near marginality where the wave amplitude is low. Consequently, the saturation mechanism for the EP instability near marginality has been attributed to the nonlinear wave-particle trapping [842], a 1D model [702] that has been successfully utilized to explain many simulation results and experimental observations near marginality. In this classic view, an EP instability saturates when the nonlinear trapping frequency equals the linear growth rate. In the presence of multiple resonances, overlap of phase space islands leads to EP diffusion and associated flattening of the EP distribution function at resonance, which diminishes the instability drive. Another mechanism of wave damping is the resonance broadening which arises from EP scattering by an ensemble of non-interacting waves.

Far away from marginality, wave-wave coupling could become the dominant nonlinear process. The transition from wave-particle to wave-wave nonlinearity needs to be quantitatively determined by gyrokinetic simulations. Gyrokinetic simulations with the GEM code [844] show that when the RSAE growth rate $γ > 0.03ω_r$ (here $ω_r$ is real frequency), wave-wave nonlinearity becomes important. In the wave-dominated weak turbulence, unstable modes can saturate through energy transfer to damped modes in a three-wave coupling process. A nonlinearly-generated ZF often plays an important role in the saturation of the instability in the toroidal geometry. The ZF is associated with the perturbed distribution function (called zonal structure) [781] averaged over the flux-surface, i.e., the $n = m = 0$ component of the perturbed density and flows that generate ZF such as zonal flows and zonal currents. Zonal flows have been found to often dominate the nonlinear saturation of the toroidal drift wave instability [845,846]. On the other hand, nonlin-

ear toroidal coupling is sub-dominant, but can lead to an inverse cascade to the lower toroidal mode number $n$ at a longer timescale [847].

ZFs and PSZS have been shown to be generated by, and in turn, suppress AEs in global simulations using gyrokinetic codes GTC [557, 843, 848], GEM [844], ORB5 [650], and kinetic-MHD codes TAEFL [849], FAR3D [850], and MEGA [803]. Effects of zonal flows are typically stronger than that of the zonal currents, but zonal structures have also been reported to dominate the AE saturation [844]. Local gyrokinetic GYRO simulations [851] found that strong zonal flows generated by microturbulence are needed to saturate RSAE. Global gyrokinetic GTC simulations using realistic plasma profiles and geometry of DIII-D experiments [657] found that zonal flows dominate the nonlinear RSAE saturation even in the presence of multiple unstable modes ($n = 3-9$) [852]. These results suggest that effects of ZFs and PSZS on EP instabilities are universal and could play an important role in EP transport.

### 7.3.2. Regulation of mode saturation by microturbulence.
An outstanding issue in global simulations of mesoscale EP instabilities is the absence of a steady state EP-driven turbulence and related transport. The nonlinear dynamic of a huge initial burst followed by a quickly diminished AE amplitude and EP transport is a common phenomenon in global gyrokinetic [804,844] and kinetic-MHD [803,849,853] simulations, where quasi-steady state EP transport can only be sustained by artificially large dissipation due to, e.g., resistivity or scattering by Coulomb collisions. However, collisions are found to have negligible effects on EP turbulence when using realistic experimental parameters. This issue has recently been addressed in global gyrokinetic simulations coupling mesoscale AE turbulence and microturbulence, where microturbulence is found to regulate EP-driven modes, resulting in a quasi-steady state EP transport [843,848]. In this work, larger EP transport has been found for stronger microturbulence even though the microturbulence directly drives little EP transport due to gyro-averaging effects, as expected by conventional wisdom [854-856]. When background microturbulence is artificially suppressed in the simulation, the RSAE amplitude and EP transport are much higher than experimental levels at nonlinear saturation, but quickly diminish to very low levels after the saturation (figure 49a). In contrast, in simulations coupling micro-mesoscales, the RSAE amplitude and EP transport decrease drastically at the initial saturation but later increase to the experimental levels in the quasi-steady state with bursty dynamics due to regulation by thermal ion temperature gradient (ITG) microturbulence (figure 49b). The RSAE
<figure><img src="image_47.png" /><figcaption>**Figure 49.** Time history of the zonal flow shearing rate $\omega_E/\gamma_4$, electron density perturbation $\delta n_e/n_e$ [%], and effective EP diffusivity $D_f$ [m$^2$/s] from the simulation of an RSAE only (panel a) and from the simulation of an ITG-RSAE (panel b). Panel c shows radial profiles of the electron temperature perturbation $\delta T_e$ from GTC simulations and ECE measurements in DIII-D. In panel d, the thermal ion temperature gradient is increased by 30%. Reproduced from [843].</figcaption></figure>

amplitude in the quasi-steady state mediated by ITG-RSAE turbulence agree well with experimental measurements (figure 49c). When the thermal ion temperature gradient is increased, owing to the stronger microturbulence, the quasi-steady state EP transport is larger (figure 49d).

Pressure gradients of thermal particles excite various drift-wave instabilities, leading to ubiquitous microturbulence responsible for turbulent transport of thermal plasmas. Despite the separation in the spatial and temporal scales, there can be strong cross-scale coupling between AEs and microturbulence. Zonal flows can be nonlinearly generated by AEs, and in turn, suppress both AEs and microturbulence. Microturbulence can dampen the zonal flows and zonal structures generated by the AEs. EP scattering by the microturbulence [857] can affect the phase space dynamics in nonlinear AE-EP interactions [858]. Drift waves and AEs can also nonlinearly interact through wave-wave coupling. Furthermore, mesoscale EP turbulence can also drive significant thermal transport as shown in recent gyrokinetic ORB5 simulations that find large electron heat fluxes driven by BAEs excited by EPs [804].

7.3.3. *Fast frequency chirping in gyrokinetic simulations.* Increased EP transport by AEs has been correlated with a fast frequency oscillation (chirping) with a sub-millisecond period that has been observed in many experiments. An analytic model for the chirping based on the one-dimensional (1D) nonlinear wave-particle interaction near marginal stability has been constructed, and single [859] and repetitive [860] bursts of chirping have been observed in kinetic MHD simulations with sources and sinks. Global gyrokinetic GTC simulations [861] of BAEs found fast and repetitive frequency chirping without sources and sinks. The wave frequency exhibits a fast, repetitive and mostly downward chirping with a sub-millisecond period and a 90° phase shift from the amplitude oscillation (figure 50). The frequency chirping is induced by the evolution of coherent structures in the EP phase space. The dynamics of the coherent structures is controlled by the competition between the phase-space island formation due to the nonlinear particle trapping and the island destruction due to the free streaming process, an intrinsically 2D dynamics in toroidal and radial directions [861, 862].
<figure><img src="image_48.png" /><figcaption>**Figure 50.** Time evolution of (a) BAE amplitude $|e\delta\phi/T_i|$ (red) and dominant frequency $\omega$ (black), and (b) frequency power spectrum. The y-axis on the left is $\omega/\omega_0$. The unit of the power intensity in panel (b) is arbitrary. Reproduced from [861].</figcaption></figure>

### 7.4. Nonlinear kinetic-MHD simulations of EP-driven modes

Kinetic-MHD simulations for EPs interacting with an MHD fluid have now become irreplaceable for understanding and predicting EP behavior [355, 828, 863–872]. In kinetic-MHD models, the bulk plasma is described as an MHD fluid, and a particle simulation method is applied to EPs. The MHD fluid and the EPs are coupled through EP pressure or EP current in the MHD momentum balance equation. Both the MHD nonlinearity and the nonlinear EP dynamics are included in the simulations. Nonlinear dynamics of EPs interacting with MHD waves is also studied in reduced simulations [585, 840, 873] and in gyrokinetic simulations [650, 855, 874, 875]. Gyro-Landau closure models such as FAR3d [850] and TGLF-EP [876] constitute a further type of kinetic-MHD hybrid model and are described in the next section. The reduced simulations are computationally more efficient, but they can be sensitive to simplifying assumptions and are less detailed in predictions. The gyrokinetic codes, on the other hand, offer more complete description of the bulk plasma but they are computationally more demanding than the kinetic-MHD simulations. One of the most capable kinetic-MHD codes is the MEGA code. It uses a multi-phase simulation method of the classical collisional processes and EP interactions with an MHD fluid. This approach makes it possible to simulate a population of EPs taking into account the injection of a neutral beam, collisional drag, pitch-angle scattering and energy diffusion, beam losses and transport processes due to MHD waves with the effects of MHD nonlinearity and EP finite Larmor radius [756]. We present some examples of the MEGA

modeling. In DIII-D experiments with neutral beam injection, there was a significant flattening of the EP profile during AE activity as shown in figure 51 [800].

<figure><img src="image_49.png" /><figcaption>**Figure 51.** EP pressure profiles and FIDA density profiles versus normalized minor radius at two different times in DIII-D. The dashed lines are the classical pressure profile predicted by TRANSP. Reproduced from [800].</figcaption></figure>

In these experiments, there was a rich spectrum of TAEs and RSAEs during the current ramp-up phase with reversed magnetic shear [482, 484, 800, 877, 878]. The ORBIT code calculations of the EP for the measured mode amplitudes [879, 880] have shown that the profiles of NBI-produced EPs flatten, consistent with the measurements. These calculations also reveal that resonance overlap of multiple AEs is essential for EP transport. The multi-phase MEGA simulation of DIII-D discharge #142111 with a self-consistent calculation of the mode amplitudes has demonstrated that the simulated flattened EP pressure profile matches with the experimental profile within the error bars (see figure 52) [493].

<figure><img src="image_50.png" /><figcaption>**Figure 52.** Comparison of an EP pressure profile for DIII-D discharge #142111: a multi-phase MEGA simulation (circles), a classical simulation (triangles), and experiment (squares), also showing an error bar. Reproduced from [493].</figcaption></figure>

The simulated temperature fluctuation profiles
were quantitatively compared with the ECE measure-
ments, and they, as well as phase profiles, showed good
agreement with the measurements. Additionally, the
saturated AE amplitudes were δB/B ~ 10⁻⁴, which
is within a factor of 2 of those measured. [879, 880].
Experiments in DIII-D also show that EP transport
suddenly becomes “stiff” above a critical threshold in
the presence of many overlapping small-amplitude AEs
(see figure 53) [881–883].

<figure><img src="image_51.png" /><figcaption>**Figure 53.** Time-averaged divergence of modulated flux, i.e., transport, inferred from (a) solid-state neutral particle analyzer (SSNPA) signal and (b) neutron emission for the first half (triangles pointing up) and second half (triangles pointing down) of the modulation period at DIII-D. Error bars are the standard error of the time average over the half period. The onset of transport corresponds to the theoretical level of stochasticity [solid squares in (a)]. Reproduced from [881]</figcaption></figure>

The EP profile is regarded as “stiff” when the
increase of the EP profile gradient above a critical
value is lower than the proportional increase to the
beam power. The threshold is phase-space dependent
and occurs when particle orbits become stochastic due
to resonances with AEs. The multi-phase MEGA
simulations predict the EP pressure profile and the EP
transport flux for different beam power levels [884].
They show stiffness and a monotonic degradation of
EP confinement with increasing beam power. The
confinement degradation and profile stiffness are due
to the presence of multiple AEs when the EP pressure
gradients exceeds a critical value. The critical pressure
gradient and the corresponding beam power depend
on radial location. The EP pressure gradient stays
moderately above the critical value, and the profiles of
the EP pressure and EP transport flux spread radially

outward from the inner region, where the beam is
injected. The resonance regions in the EP phase
space have been analyzed for the eigenmodes observed
in the MEGA simulations. Figure 54 shows the
particle trajectories in the phase space of normalized
major radius and energy with the AEs present in the
simulations for DIII-D experiments for different beam
deposition power levels [884]. With increasing beam
power, the resonance regions broaden due to the larger
amplitude of eigenmodes leading to the overlapping.

This analysis verifies that the overlap of multiple
resonances is the underlying reason for the sudden
increase in EP transport with increasing beam power.
The multi-phase MEGA simulations [885] of the TAE
bursts and EP losses in TFTR have reproduced the
observed synchronized bursts of multiple TAEs and
the corresponding variations of the stored beam energy
with regular time intervals close to the experimental
value, as shown in figure 55.

The mode saturation amplitude is now lower than
in the earlier reduced simulations with a linearized
description of the bulk plasma [886]. The inclusion
of nonlinear MHD effects [787, 803, 849] in MEGA
prevents the AE amplitude from excessive growth.

In experiments with super-Alfvénic beam ions,
repetitive bursting events consisting of multiple AEs
or EPMs are observed. In JT-60U, negative-ion
beams drove a series of repetitive bursts dubbed “fast
frequency sweeping modes” that culminated in an
“abrupt large event” (ALE) [887]. Similarly, in NSTX,
a sequence of bursting AEs often culminated in a
larger burst with multiple harmonic content dubbed
an “avalanche” [888, 889]. The avalanche events
caused substantial drops in neutron [888, 890] and
FIDA signals [889, 891]. Many features of the JT-
60U ALE [860, 892, 893] and NSTX avalanche [894]
are reproduced by simulations, including the entire
repetitive burst cycle in simulations with MEGA [860].

The EP-driven AEs and EP transport in ITER
plasmas have been studied numerically in [31,573,578,
580, 876, 895]. It was found that the EP drive of
AEs can overcome the thermal ion Landau damping
only in the outer half of the plasma, and AE-induced
redistribution is not expected to affect the fusion burn
in the ITER 15 MA baseline scenario [31]. For the
steady-state scenario with 9 MA plasma current, BAEs
with low toroidal mode number were found to dominate
in the nonlinear phase although many TAEs with n ~
15 are more unstable in the linear MEGA simulations
[573].

7.5. *Nonlinear kinetic-MHD simulations of dynamics and frequency chirping of single-n Alfvén modes*

The nonlinear dynamics of a single toroidal Alfvén
mode and, in particular, the saturation mechanisms
<figure><img src="image_52.png" /><figcaption>**Figure 54.** Particle trajectories in the phase space of normalized major radius $\hat{R} = (R - R_{\text{axis}})/(R_{\text{edge}} - R_{\text{axis}})$ and energy $\mathcal{E}$ (keV) for beam deposition power (a) 1.56 MW, (b) 3.13 MW, (c) 6.25 MW, and (d) 15.6 MW. Particle orbits are followed with the electromagnetic perturbations of a single TAE with fixed amplitude and frequency, and $\hat{R}$ and $\mathcal{E}$ are recorded when the particle passes the mid-plane from bottom to top. Only the particles trapped by the TAE are plotted in the figure. The eigenmodes are represented by colors: n = 1 (blue), n = 2 (purple), n = 3 (green), n = 4 (orange), and n = 5 (red). Reproduced from [884].</figcaption></figure>

<figure><img src="image_53.png" /><figcaption>**Figure 55.** (a) Radial profile evolution of EP energy transport flux in the radially outward direction, and amplitude evolution of radial MHD velocity fluctuation for (b) the dominant n = 2 (red) and n = 3 (blue) AEs and (c) the other AEs with n = 1–3. The unit of color bar is MW and the beam injection power is 10 MW. Reproduced from [885].</figcaption></figure>

have been investigated by particle-in-cell simulations performed with the gyrokinetic-MHD code XH-MGC [866, 896, 897] and the so-called Hamiltonian-mapping technique [898], both for constant and chirping frequencies. This technique samples the regions of phase space that provide the most important contribution to the growth or damping of modes by means of a population of test particles evolving in the electromagnetic fields obtained from self-consistent particle simulations. Wave-trapped particles form an island structure, around the resonance radius, in the 2D space $(\Theta, r)$, with $\Theta$ being the wave phase seen by the particle and $r$ its radial coordinate (figure 56). This yields a density flattening delimited by a sufficiently steep negative density gradient. In general, mode saturation occurs as the flattening region extends over the whole radial region where the mode-particle power transfer can take place [870, 898]. A variety of phase space diagnostics can be used to characterize the trapping-detrapping process by multiple resonances, see e.g. [825, 826, 899, 900].

In the constant-frequency case, two regimes can be distinguished. In the first regime, the power-transfer region is mainly limited by the need to satisfy the resonance condition $|\omega - \omega_{\text{res}}(r)| \lesssim \gamma$, where $\omega_{\text{res}}(r)$ is the resonance frequency at radius $r$, and $\gamma$ is the linear growth rate of the mode. This is
<figure><img src="image_54.png" /><figcaption>**Figure 56.** Hamiltonian mapping of test-particle markers in the $(\Theta, r)$ plane in the nonlinear stage. Each marker is coloured according to the birth $r$ value of the particle. The formation of an island-like structure around the resonance radius, including the trajectories of particles trapped in the potential well of the wave, is clearly seen. Reproduced from [825].</figcaption></figure>

the case occurring for low values of the growth rate, sharp profiles of the resonance frequency or extended mode structure. The second regime, occurring in the opposite limits, sees the power-transfer region limited by the finite radial width of the mode structure. The two different regimes correspond to different scalings of the saturation mode amplitude with the linear growth rate: the former (*resonance detuning* regime) exhibits a quadratic scaling; the latter (*radial decoupling* regime), a linear scaling; this is shown in figure 57 [870].

In the case of chirping frequency, the dynamics, in a first nonlinear stage, is analogous to that analysed in the constant-frequency case. If the frequency is allowed to change, however, the resonance region may move such that the coverage of the power-transfer region by the density-flattening region is delayed and the mode can extract more power from the resonant particles. In turn, the island and, then, the density-flattening region reconstitute around the new resonance radius, counteracting this retarding effect. This process goes on until a further change in frequency becomes unable to produce a significant displacement of the resonance region. To further grow, the mode has to tap a different resonant structure, possibly making use of additional frequency variations. The phenomenology described above has been observed in simulations referring to a variety of physical situations, ranging from AE modes driven unstable by EPs (both fusion alphas [898] and EPs from auxiliary heating [825,870,901,902]) as well as energetic electrons (as in the case of electron fishbones [827]).

<figure><img src="image_55.png" /><figcaption>**Figure 57.** Scaling of the saturation amplitude of the scalar potential (defined as the radial peak, in arbitrary units, of the dominant poloidal harmonic) versus the linear growth rate $\gamma$ for XHMGC simulations. The reference quadratic and linear $\gamma$ scalings associated with the resonance-detuning and radial-decoupling regimes are also shown. Reproduced from [870].</figcaption></figure>

## 7.6. Nonlinear gyrofluid simulations of EP-driven instabilities

Gyrofluid models are a further step in the evolution of the MHD-kinetic hybrid paradigm for EP instability analysis. This approach is motivated by the fact that all gyrokinetic models must communicate kinetic information to field equations through some form of moments integral (e.g., charge/current density, pressure, etc.). Since these low order moment integrals will inherently average away some of the fine structure of the kinetic distribution function, this motivated replacing the kinetic component with a hierarchy of moment equations. The advantages of this approach are: a significant dimensionality reduction from 5D/6D to 3D, leading to fewer degrees of freedom and high computational efficiency; no discrete particle noise; no issues with time-evolving uneven particle phase space densities (cavitation); and a mathematical structure that is compatible with direct eigenmode solvers. While moments hierarchy or MHD models have been developed and applied for many years in plasma physics, they were not capable addressing the phase-mixing and instability feedback effects associated with wave-particle resonant phenomena, such as Landau damping/growth until the development of gyro-Landau closure methods. These new forms of moments closures, developed originally by Hammett and Perkins [903] demonstrated that Landau resonance phenomena could be included in fluid-like equations. This approach has been adapted to the analysis of EP instabilities through models such as TAEFL [865], FAR3d [904], MAS [905] and TGLF-EP [906, 907], which use closure relations optimized for Alfvénic instabilities. FLR stabilization effects are also included. FAR3d and TAEFL were based on coupling a set of EP gyrofluid equations to an existing
<figure><img src="image_56.png" /><figcaption>**Figure 58.** (a) Time evolution of poloidal magnetic field component at $\rho/\rho_{\text{edge}} = 0.2$, $\theta = 0$, $\zeta = 0$, (b) spectrogram of magnetic field signal in (a), with the dominant mode indicated at $\simeq 100$ kHz. Reproduced from [850].</figcaption></figure>

global reduced MHD model. MAS [905] is also a global Landau closure model that has been tested in the linear regime for drift waves (ITG and KBM), internal kinks, and Alfvén instabilities (KAW, TAE and RSAE). TGLF-EP was developed starting from the flux-tube model TGLF [908] that incorporated gyro fluid closures for the modeling of drift waves.

The global EP gyrofluid model, FAR3d has been used to model a variety of EP instability observations in both tokamaks and stellarators [909]. It includes both linear and nonlinear options. The nonlinear saturated state is achieved through the inclusion of convective and $J \times B$ nonlinearities [850, 910, 911]. This model has been used for the long-time scale simulation of nonlinearly saturated AE modes observed in the DIII-D experiment, and has shown reasonable agreement. An example of the simulated poloidal magnetic field fluctuations and their associated spectrogram is shown in figure 58 [850]. In the EP density evolution equation the convective nonlinearity drives flattening of the fast ion profile in the region where the EP instability amplitude is strongest; EP transport fluxes can be derived from these profile changes. 2D zonal flow and current generation effects also are present and play a role in regulating the nonlinear state. Mode coupling effects are well-resolved in this model and indicate that for the case of simulations with multiple toroidal modes, linear growth rates can be strongly altered by nonlinear mode coupling effects.

The TGLF-EP model [906, 907] provides a rapid method to evaluate local EP instability properties and has mostly been used for linear stability evaluations. TGLF-EP goes higher in the moments hierarchy than FAR3d and includes trapped particle and drift

effects in the resonance conditions, in addition to the usual passing particle resonance. This model has been particularly useful in the development of critical gradient models for EP transport. Local variations in the EP density gradient drive are made until a local marginal stability state is achieved. By repeating this process over the full range of flux surfaces, a marginal stability EP density profile can be reconstructed and EP transport rates inferred from the profile flattening effects.

### 7.7. Nonlinear simulations of EP-driven modes with reduced models

Apart from the cutting-edge numerical models aimed at rigorous first principle coverage of all known essential physics effects, there has been significant progress in less rigorous but fast modeling that involves some phenomenological assumptions and free input parameters. Less accurate quantitative interpretation of wave-particle phenomena can be expected from such an approach with respect to first principle codes, especially for predictive simulations that are not constrained by experimental data. Nonetheless, these simplified models allow inexpensive scans across multiple discharges for a quick assessment of the relative role of instabilities in EP transport and losses. Such reduced models distill information from theory and first-principles codes and implement efficient numerical methods to cover EP transport in long-timescale integrated simulations [912–914].

The degree of model reduction determines which physics aspects are retained and to which degree of fidelity. The most common simplifications are assumptions about the instability spectrum and the transport mechanisms, a simplified representation of the EP population, and simple ad-hoc transport coefficients. When a Monte Carlo approach is used to simulate the EP evolution, such as in the NUBEAM module of TRANSP [915, 916], EP transport coefficients can be represented by matrices that describe the impact of modes on EP orbits as kicks in EP variables during the simulation. This EP kick model [917, 918] has been implemented in NUBEAM. The kick matrices are defined for the EP constants of motion (COM) ($\mathcal{E}$, $\mu$, $P_\phi$) [919]), see section 2. For each ($\mathcal{E}$, $\mu$, $P_\phi$) region (or bin), a 2D matrix represents the conditional probability $p(\Delta\mathcal{E}, \Delta P_\phi | \mathcal{E}, \mu, P_\phi)$ of correlated kicks in $\mathcal{E}$ and $P_\phi$ resulting from EP interaction with instabilities. The kick matrix comes from particle following codes such as ORBIT [918, 919]. For Alfvénic modes, MHD codes such as NOVA [920] provide a radial mode structure to use in ORBIT. Analytic representations can also be used to represent kink modes, fishbones, and NTMs. Kick amplitudes in NUBEAM are time-dependent according
to a user-supplied waveform. If the mode damping rate is known, the model can adjust the amplitude accordingly [918]. The kick model has been applied to several experiments to study EP transport by Alfvénic modes [232,882,918,921], NTMs [232,317,369,922] and internal kinks [392, 396, 399], including scenarios with multiple instabilities present simultaneously [335, 341]. By retaining full phase-space details and realistic EP distribution functions from NUBEAM, the model is suitable for detailed comparisons between simulation results and EP diagnostics [232, 317, 369, 392, 396] via synthetic diagnostic codes such as FIDASIM [923,924]. Application of the kick model to ITER requires more work to meet the remaining challenges. Work is in progress to extend and validate the model to multi-species plasmas (D, T, or fusion products) [46, 925, 926], leveraging new data from the recent JET D-T campaign [50, 51]. Also, depending on the granularity of $(\mathcal{E}, \mu, P_\phi)$ bins, matrices can be memory-consuming, and the model requires amendments to cover multiple modes.

For multiple Alfvénic modes, in scenarios with either isolated or overlapping resonances, a resonance-broadened quasi-linear theory, originally proposed in [927], provides a way to reduce the dimensionality of the problem while still capturing the essential dynamics in the EP diffusive regime [928]. The Resonance Broadening Quasi-linear model RBQ [841, 929] computes the diffusion coefficients to relax the distribution function of EPs while simultaneously evolving the amplitudes of multiple modes. The model, originally developed to account for diffusion along $P_\phi$ only [841, 929], has recently been upgraded to resolve the $(\mathcal{E}, P_\phi)$ correlation (equation 13) for resonant interactions, $\Delta P_\phi / \Delta \mathcal{E} = n/\omega$ [930]. RBQ has motivated recent basic analytic work on the formulation of quasilinear theory in a broad sense. For instance, sufficiently close to marginal stability, collisions regulate the wave evolution, phase memory is poorly retained and the EP dynamics becomes time-local [931, 932]. This understanding led to the formulation of a self-consistent quasi-linear theory [933, 934] that recovers the exact saturation level of the more complex nonlinear theory [935], while being considerably less computationally demanding. An essential ingredient of the theory, the shape of the resonance (window) function, previously employed with an arbitrary ad-hoc shape, emerges spontaneously in the self-consistent derivation [933, 934]. This self-consistent window function is used in RBQ [858], yielding better agreement when compared with the nonlinear Vlasov code BOT [936]. Verification of the width of the wave-particle interaction has been performed using the guiding center following code ORBIT [937]. The model can be used to find

the relaxed EP profile with self-consistent mode amplitudes, based on AE structure and spectrum from MHD codes such as NOVA [920] or LIGKA [530], and the EP drive calculated for an analytic slowing down distribution. RBQ has been used to assess NBI ion transport by AEs on DIII-D [938]. At present, RBQ communicates with TRANSP/NUBEAM via the same 5D transport matrix formalism as the kick model. The RBQ diffusion coefficients provide Gaussian probabilities for $(\Delta\mathcal{E}, \Delta P_\phi)$. Work is in progress to provide a direct interface and reduce the memory usage with respect to the kick model. Further heuristic elements lead to critical gradient models such as TGLF-EP [876, 906, 939] and CGM [940], which seek a solution for the EP radial flux such that the EP drive balances the mode damping at each radial location. The radial diffusion coefficient computed from the EP flux can be used as an input for TRANSP.

As mentioned above in subsection 7.4, the TGLF-EP reduced model [906, 907] provides an additional method for evaluating nonlinear Alfvén instability driven EP transport using an approach known as a critical gradient model. TGLF-EP uses the TGLF [908] local flux tube gyro-Landau closure model to evaluate AE instability growth rates on each flux surface in the ballooning limit and find the marginal stability point (or point where the growth rate is below some minimum level) by varying the EP gradient drive. Under the assumption of stiff transport (i.e., rapidly increasing transport as AE modes destabilize), a steady-state solution of the EP continuity equation (including beam/fusion source, slowing-down sink, and diffusive AE flux) gives the marginally stable EP density profile. The underlying stiff-transport assumption has been checked [851, 907] against simulations using the gyrokinetic code GYRO [941]. The calculated effective diffusivity can be coupled into integrated simulations. This model has been confirmed with EP profiles measured on DIII-D from FIDA measurements and neutron deficit data [563, 881, 942]. It has also been applied to cases from ITER [876].

Finally, reduced modeling can also lead to simulations of the EP distribution function fast enough to allow real-time control of discharges. The RABBIT code is such a model allowing very fast computations of the NBI distribution functions neutral-beam populations based [943–945].

## 7.8. Strengths and weaknesses of nonlinear simulation approaches

Here we attempt to summarize typical strengths and weaknesses of various nonlinear simulation approaches treated in the previous subsections: gyrokinetic simulations, kinetic-MHD simulations, gyrofluid simulations,
tions, and ad hoc reduced model simulations, e.g. the kick model implemented in TRANSP. For each simulation approach, several codes are available with developed or developing user communities. Naturally, it is hardly possible to provide fully correct lists encompassing all codes of their category. Particular codes might have overcome some of the typical weaknesses of their category, or not fully exploit typical strengths. Nevertheless, we can provide a starting point for further investigation. In addition to the typical strengths and weaknesses pointed out below, numerical aspects should be considered, such as particle sampling noise or the need to introduce hyperviscosity or time step restrictions due to numerical stability. Since these numerical issues depend on the implementation, we will not list them here.

### Gyrokinetic simulations - strengths:

* They currently provide the most complete EP physics model available. Highest fidelity.

* They can take into account the detailed phase-space structure of the distribution function and finite orbit width effects.

* Efficient CPU/GPU parallelization is achieved with particle and domain decompositions.

* They can analyze coupled EP instability/core microturbulence effects.

* Both linear and nonlinear options available.

### Gyrokinetic simulations - weaknesses:

* Gyro-motion is averaged out.

* Typically only a subset of *E*- and *B*-field components are included. For a detailed recent discussion in the framework of gyrokinetic field theory, see [946].

* Heavy computational requirements for nonlinear full-*f* runs (significant CPUs/GPUs + lengthy run times).

* Methods have not been developed to insure evenly distributed phase space densities through the simulation - may limit simulation time.

### Kinetic-MHD simulations - strengths

* Relative to gyrokinetic, this is the next best level of completeness in terms of the physical model.

* Guiding-center orbit width effects are taken into account, and most models can utilize a range of EP distribution functions.

* The full set of *E*- and *B*-field-field components are typically included.

* Can be run in hybrid mode (i.e., alternate between nonlinear instability and classical collisional physics to achieve long simulation times).

* Efficient CPU/GPU parallelization is achieved with particle decomposition.

* Both linear and nonlinear options are available.

### Kinetic-MHD simulations - weaknesses

* Gyro motion is averaged out.

* Heavy computational requirements for nonlinear runs (significant CPUs/GPUs + lengthy run times).

* The thermal plasma is generally treated in the fluid approximation so coupled EP instability/core microturbulence effects cannot be modeled.

* Some damping mechanisms of EP-driven modes such as radiative damping and thermal ion Landau damping are missing while continuum damping is retained.

* $E_\parallel$-effects are not treated consistently, as the electron dynamics is missing.

* Methods not used to insure evenly distributed phase space densities through the simulation - may limit simulation time.

### Gyrofluid simulations - strengths

* Fast execution speed.

* These are continuum models. Relatively accurate treatment of mode-coupling effects is achieved – important for zonal flow effects.

* Both initial value and eigenmode options available for linear stability analysis.

* Includes approximations for FLR effects and finite orbit width effects.

* Both linear and nonlinear options are available.

* CPU/GPU parallelization via domain decomposition - nonlinear runs can be lengthy, but do not require as many computational nodes as gyrokinetic or kinetic-MHD.

### Gyrofluid simulations - weaknesses

* Limited flexibility with respect to EP distribution function model (mostly near Maxwellian models are used, but may be extended by including higher order moments).

* Full field components are not included, most models are based on reduced MHD.

* Some damping mechanisms of EP driven modes such as radiative damping and thermal ion Landau damping are included via perturbative approximations while continuum damping is retained.

* $E_\parallel$-effects are not treated consistently, as the electron dynamics is missing.
* Nonlinear wave-particle trapping which plays an essential role in the saturation process of EP driven modes is not present per se, but modeled by convective cell trapping of fast ion density.

* Developed from gyrokinetic models that average over gyromotion.

### Ad hoc reduced model simulations - strengths

* Large user community, e.g. the kick model implemented in TRANSP.

* Fast execution time.

* They provide answers for transport rates and profile flattening that can readily be incorporated into integrated models.

### Ad hoc reduced model simulations - weaknesses

* Lowest physics fidelity

* Models are based on many simplifying assumptions that may not apply to all regimes, in addition to gyroaveraging.

* Ranges of applicability not always well known.

* Do not include mode coupling effects, zonal flows, and currents.

* Reduced models generally give asymptotic states (e.g., critical gradient profiles) based on linear theory with not much information on the time intervals to get to such states.

Besides these ad hoc reduced models, recent studies have suggested approaches for developing reduced descriptions based on first principle theory derivations with different levels of approximations that are justified with assumptions, which can be checked at each level of the hierarchical approach (cf. Sec. 7.2 and [32]). The feasibility of the corresponding workflow has recently been demonstrated [813].

## 7.9. Application of nonlinear simulation models to ITER

The parameters of ITER introduce new regimes for nonlinear Alfvén turbulence with unstable modes extending to higher toroidal mode numbers (i.e., smaller $\rho_{Li}/a$ allows instability at higher wavenumbers before FLR stabilization becomes active). This results in many interacting modes and closely spaced Alfvén gap locations. As described above, there are a variety of models now under development to address the nonlinear physics of ITER burning plasma scenarios. A recent collaboration of nonlinear EP models focused on simulating these effects as part of an ITPA joint simulation activity [587]. Both an ITER pre-fusion baseline scenario case (#101006) and a steady-state case (#131041) with a reversed shear $q$-profile [947] were considered. The pre-fusion baseline case indicated

that a weak fishbone instability was present, but that Alfvén modes were stable; only weak neutral beam ion redistribution was predicted. The steady-state case (in the DT phase) contained both anisotropic neutral beam ions and isotropic alpha particles. Due to the stronger drive and reversed $q$-profile, it showed stronger Alfvénic growth rates than for cases with monotonic $q$-profiles. The global nonlinear models applied to this case included GTC, MEGA, FAR3d, M3D-C1 [948], [587]. Also, the reduced RBQ and TGLF-EP models [949] were applied in this study. The most unstable modes were BAEs and RSAEs and were located near the $q_{min}$ surface. All models [587] indicated high levels of transport (50 to 70 m²/s) both for the beams and alphas. The simulations were limited to relatively short time intervals (~ 0.25 ms). In contrast, NBI-driven cases from DIII-D with similar $q$-profiles only resulted in a few unstable AEs. If this level of transport extrapolates to longer timescales, it would be difficult to sustain a steady-state burning plasma. In particular, the maintenance of reversed shear $q$-profiles will be important to achieve non-inductive (bootstrap current aligned) steady-state operation. Future directions for this type of modeling should include longer timescale integrated simulations that include cross-scale coupling between MHD, AEs, microturbulence and the self-consistent effects of alpha and beam source functions.

## 7.10. Soft and hard nonlinearities of near-threshold instabilities

The advent of burning plasmas might suggest an increase focus on alpha-driven instabilities in burning plasmas. However, a large part of the underlying physics is generic and understandable at a very basic level within properly idealized models of wave-particle interaction. Even though the EP pressure can be comparable to the bulk plasma pressure, the EP density is usually much smaller than the bulk particle density. EPs hence interact mainly with the bulk plasma rather than among each other. EPs have, by definition, a non-Maxwellian distribution which has free energy available to excite waves via wave-particle resonances. Depending on the problem of interest, the free energy can come either from inverted velocity-space gradients, anisotropy, or from spatial gradients of the EP distribution function (section 6).

EP instabilities have typically much shorter intrinsic timescales than the particle or energy confinement time. Consequently, linear stability analysis of a given initial configuration does not cover how the instability saturates in the presence of particle sources and the slow collisional relaxation processes. To address saturation in general, the theory must be nonlinear, which is still a challenge. However, given the difference in timescales and given that the system
is weakly driven, the instability should occur in a near-threshold regime, and then the nonlinear stage is often tractable analytically.

There are two basic sources of nonlinearity: the nonlinearity of the EP motion in the field of the excited waves and the nonlinearity of the bulk plasma. In each of these two cases, it is critical to determine whether the near-threshold instability exhibits a soft or a hard nonlinear regime, i.e., whether it can be stabilized by weak nonlinearity or not. First, the soft or hard regime require different theoretical and numerical tools to describe the system, and second, there are immediate experimental implications. In particular, the fishbone instability [741] generally requires treatment of wave-particle resonances and bulk plasma nonlinearities. For AEs, the nonlinear interaction of EPs with linear waves dominates the problem. In what follows, we first present a theory of single resonance phenomena, and then discuss particle transport in the presence of multiple modes.

### 7.11. Weakly nonlinear near-threshold dynamics

The interplay of the key ingredients in the problem of wave-particle interaction can be understood within a simple electrostatic bump-on-tail model that exhibits the characteristic nonlinear scenarios. Within this model, the bulk plasma is represented by cold electrons, and we assume that there are sources and sinks that create an unstable energetic electron tail. The tail provides an instability drive $\gamma_{drive}$ due to a positive gradient of the velocity distribution function $F$ (see figure 59).

<figure><img src="image_57.png" /><figcaption>Figure 59. Unstable bump-on-tail distribution function with a positive gradient at the resonant velocity.</figcaption></figure>

The unstable mode in this case is a plasma wave, and its eigenfrequency, $\omega_{pe}$, is the cold-electron plasma frequency. The cold-electron collision frequency $\nu_{cold}$ provides a linear damping rate $\gamma_{damp} = \nu_{cold}/2$ of the mode that determines the the minimum slope of $F_0(u)$ needed to excite a mode, i.e. the instability threshold.

The spectrum of EP-driven modes in a tokamak is generally discrete due to periodicity in the toroidal and poloidal directions and the radial boundary conditions (section 6). To take that into account in the bump-on-tail model, we consider a single electrostatic mode with a given wavelength $\lambda$ and wavenumber $k = 2\pi/\lambda$. The electric field of the mode can then be written as

$$E = \frac{1}{2} [\delta E(t) \exp(ikx - i\omega_{pe}t) + \text{c.c.}], \quad (74)$$

where $\delta E(t)$ is a slowly varying complex amplitude and c.c. denotes the complex conjugate. The basic equations for the bump-on-tail problem are the kinetic equation for the energetic electrons and the wave evolution equation

$$\begin{aligned} \frac{\partial F}{\partial t} + u \frac{\partial F}{\partial x} + \frac{e}{2m} [\delta E(t) \exp(ikx - i\omega_{pe}t) + \text{c.c.}] \frac{\partial F}{\partial u} \\ = \left[ \frac{\nu^3}{k^2} \frac{\partial^2}{\partial u^2} + \frac{\alpha^2}{k} \frac{\partial}{\partial u} - \beta \right] (F - F_0) \end{aligned} \quad (75)$$

$$\frac{\partial \delta E}{\partial t} = -4 \frac{\omega}{k} \pi e \int \langle F \exp(-ikx + i\omega_{pe}t) \rangle du - \gamma_{damp} \delta E, \quad (76)$$

where $\langle ... \rangle$ denotes averaging over wavelength. The right-hand side of the kinetic equation describes three different collision models for the resonant tail electrons: velocity-space diffusion, drag, and Krook-type collisions [950]. The characteristic rates for these collisions are characterized by the quantities $\nu$ (velocity-space diffusion), $\alpha$ (collisional drag), and $\beta$ (Krook model collisions), respectively (note that these specific meanings for $\nu$, $\alpha$, $\beta$ are specific for this subsection and do not apply in other parts of this paper). The collision operator also includes the source and sink terms written as $(\frac{\nu^3}{k^2} \frac{\partial^2}{\partial u^2} + \frac{\alpha^2}{k} \frac{\partial}{\partial u} - \beta) F_0$, which sets up an equilibrium distribution function $F_0$ in the absence of the wave field. The appropriate collision operator for the problem is determined by what collisional process is dominant at the wave-particle resonance in phase space. For EPs in a tokamak, Coulomb collisions can be described as a combination of pitch angle scattering and electron drag (section 2) [78, 951]. The former can be represented by a diffusive operator, while the latter introduces a slowing-down operator to the kinetic equation. The near-threshold regime of wave excitation makes it possible to expand the perturbed distribution function $F$ in powers of the wave amplitude $\delta E$ and solve the kinetic equation iteratively. The actual expansion parameter is $\omega_B t$, where

$$\omega_B \equiv (k|\delta E|e/m)^{1/2} \quad (77)$$
is the bounce frequency of the resonant particles trapped in the wave, and $t$ is the time interval of interest. The first term in the power series for $F$ gives the linear instability drive $\gamma_{\text{drive}}\delta E$ in the wave equation (76). The difference between $\gamma_{\text{drive}}$ and $\gamma_{\text{damp}}$ is small in the near-threshold limit, which allows the lowest order nonlinear correction to compete with this difference. It follows from the expansion procedure that the nonlinear correction to the wave growth rate scales as $\gamma_{\text{drive}}(\omega_{\text{B}}t)^4$ whereas the linear growth rate itself is $\gamma_{\text{drive}} - \gamma_{\text{damp}} \ll \gamma_{\text{drive}}$. Consequently, the lowest order nonlinearity becomes important when $(\omega_{\text{B}}t)^4 \approx (\gamma_{\text{drive}} - \gamma_{\text{damp}})/\gamma_{\text{drive}} \ll 1$. At this level, the next-order nonlinear term, $\gamma_{\text{drive}}(\omega_{\text{B}}t)^8$, is still negligible. Thus, the inequality $(\gamma_{\text{drive}} - \gamma_{\text{damp}})/\gamma_{\text{drive}} \le (\omega_{\text{B}}t)^4 \ll 1$ defines a window in which the dynamics are already nonlinear, but the nonlinearity can still be treated perturbatively. The ensuing relation between the perturbed distribution function and the wave field involves a sequence of time integrations. Once this relation is used in equation (76), we obtain a cubic integro-differential equation for the wave amplitude, which can be written in the following dimensionless form:

$$ \frac{dA}{d\tau} = A(\tau) - \frac{1}{2} \int_{0}^{\tau/2} dz z^2 A(\tau - z) \times \int_{0}^{\tau-2z} dx \exp \left[ -\hat{\nu}^3 z^2 \left( \frac{2z}{3} + x \right) - \hat{\beta} \left( 2z + x \right) - i\hat{\alpha}^2 z (z+x) \right] \times A(\tau - z - x) A^*(\tau - 2z - x). \quad (78) $$

The dimensionless amplitude $A$ in this equation is defined as

$$ A = (k\delta E e/m) \left( \frac{\gamma_{\text{drive}}}{\gamma_{\text{damp}}} - 1 \right)^{-1/2} (\gamma_{\text{drive}} - \gamma_{\text{damp}})^{-2}, \quad (79) $$

the dimensionless time is $\tau = (\gamma_{\text{drive}} - \gamma_{\text{damp}})t$ and the normalized relaxation rates are defined as $\hat{\nu} \equiv \nu/(\gamma_{\text{drive}} - \gamma_{\text{damp}})$, $\hat{\alpha} \equiv \alpha/(\gamma_{\text{drive}} - \gamma_{\text{damp}})$, and $\hat{\beta} \equiv \beta/(\gamma_{\text{drive}} - \gamma_{\text{damp}})$. The cubic nonlinear equation (78) was originally derived in [935, 952] for the diffusive and Krook-type collisions, and it has been generalized in [953] to include the effect of drag. According to equation (78), the initial linear instability can develop into a soft or hard nonlinear regime. In the soft case, the amplitude $A$ saturates at a finite level. In the hard case, the solution for the amplitude $A$ rapidly grows in a finite time. Without drag ($\hat{\alpha} = 0$), a saturated solution of equation (78) is

$$ |A|^2 = 2 \left[ \int_0^\infty \frac{z^2 dz}{\hat{\beta} + \hat{\nu}^3 z^2} \exp \left( -2\hat{\nu}^3 z^3 / 3 - 2\hat{\beta}z \right) \right]^{-1} \quad (80) $$

for $\tau \to \infty$. If the annihilation rate $\hat{\beta}$ and/or diffusion rate $\hat{\nu}$ are above some threshold, the amplitude approaches that solution. However, if $\hat{\beta}$ and $\hat{\nu}$ are below this threshold, the steady saturated solution is unstable. The solution then has a characteristic “pitchfork splitting”, a periodic limit-cycle behaviour. An experimental example of pitchfork splitting driven by ICRF-accelerated tail ions appears in figure 60 [954].

<figure><img src="image_58.png" /><figcaption>Figure 60. Nonlinear splitting of AEs in JET [954, 955]. Reproduced from [954].</figcaption></figure>

A further reduction in relaxation rates leads to period doubling bifurcations, resulting in a chaotic mode amplitude evolution and explosive growth of the mode [954–956]. The cubic nonlinear term in equation (78) destabilizes the mode for pure drag (where $\hat{\beta} = \hat{\nu} = 0$), so the mode is then explosive. Then there are no saturated solutions for equation (78), and the mode amplitude grows beyond the applicability range of equation (78). When both drag and diffusion are present, steady saturated solutions are only prohibited when the integral in equation (78) has a negative real part for $\tau \to \infty$. This occurs when $\hat{\nu}/\hat{\alpha} < 1.043$ (dashed line in figure 61). However, some of the formal steady solutions for $\hat{\nu}/\hat{\alpha} > 1.043$ have been shown to be unstable [953]. The solid line in figure 61 represents the stability boundary between stable and unstable steady solutions. Recent work shows that, in contrast to diffusion, drag fundamentally shifts the resonance condition in the quasilinear regime, producing shifts
and splitting of the resonance lines that alter the saturated amplitude [934].

<figure><img src="image_59.png" /><figcaption>**Figure 61.** The boundaries in parameter space that give stable, unstable and no steady state solutions to equation (78). The unstable solution lies in between the solid and dashed lines. Reproduced from [953].</figcaption></figure>

The explosive growth described by equation (78) allows the growing wave to reach the level of resonant particle trapping by the wave. One might, therefore, expect that particle trapping will flatten the distribution function near the resonance and eliminate the instability drive, after which the wave will decay quickly because of the background damping. However, the solution of equations (75) and (76) reveals a remarkably different behaviour. It exhibits formation of long-living coherent structures with time-dependent frequencies, as discussed next.

### 7.12. Spontaneous frequency sweeping, phase-space holes and clumps

The tendency for the mode frequency to change in the strongly nonlinear regime is already seen in the explosive solution of the reduced cubic equation (78). The explosive solution is in fact oscillatory, and the period of oscillations in the wave amplitude A shortens as the solution approaches the singularity. These oscillations indicate that the wave tends to split into the upshifted and downshifted sidebands. The fully nonlinear set of equations (75) and (76) prevents the mode from growing indefinitely. However, the trend for frequency sweeping continues, as found numerically in [702, 703] and shown in figure 62.

The time-dependent wave spectrum (lower plot in figure 62) correlates perfectly with the evolution of the spatially averaged particle distribution function (upper plot). The latter exhibits an upward moving depletion (hole) and a downward moving protrusion (clump). The hole and clump contain particles trapped in the upshifted and downshifted waves, respectively. Unlike the fast explosive onset of holes and clumps, the evolution shown in figure 62 takes place over

<figure><img src="image_60.png" /><figcaption>**Figure 62.** Evolution of the particle distribution function and wave spectrum after explosive formation of holes and clumps in the hard nonlinear regime of near-threshold instability. Reproduced from [703].</figcaption></figure>

many bounce periods of the particles that are trapped in the wave, so that these particles respond to the wave field adiabatically. Conservation of the adiabatic invariant preserves the values of the particle distribution function within the trapped particle areas of phase space. The reasons for frequency sweeping shown in figure 62 are the presence of dissipation in the bulk plasma and the availability of free energy in the EP distribution.

<figure><img src="image_61.png" /><figcaption>**Figure 63.** Adiabatic motion of hole and clump releases kinetic energy of the EPs. The increased slope of the distribution function in the wake (thick line) facilitates formation of subsequent holes and clumps. Reproduced from [936].</figcaption></figure>

A schematic snapshot of the distribution function in figure 63 shows that the particle kinetic energy decreases when the hole and clump move away from the original resonance with the constant values of the
distribution function at the bottom of the hole and the top of the clump. This energy release balances the dissipation in the background plasma to allow the wave to last over hundreds of linear damping times. As shown in [702], each hole and clump represent a nonlinearly saturated wave with

$$ \omega_B = (16/3\pi^2) \gamma_{\text{drive}}. \qquad (81) $$

Calculation of the power release via frequency sweeping relies on the knowledge that the distribution function of passing particles is smooth around holes and clumps. This calculation gives the following square root time-dependence for the frequency shift in the absence of EP collisions [702]:

$$ \delta\omega = (16/3\pi^2) \gamma_{\text{drive}} \sqrt{2\gamma_{\text{damp}}t/3}. \qquad (82) $$

The results shown in figure 62 are consistent with this relation. More recent simulations with improved computational accuracy [936] reveal that holes and clumps are produced continuously in the collisionless case. This can be understood by noting that the slope of the distribution function at the original resonance steepens somewhat after the hole and clump move away from it. There is, therefore, a tendency for a recurrent instability. The presence of drag and velocity space diffusion adds interesting new features to the behavior of holes and clumps [936]. The drag alone breaks the symmetry of the sweeping pattern, as demonstrated in figure 64.

<figure><img src="image_62.png" /><figcaption>Figure 64. Effect of drag on holes and clumps. The plots show an asymmetric frequency spectrum and two snapshots of the particle distribution in phase space. Reproduced from [936].</figcaption></figure>

The source term in the drag collision operator acts to enhance a phase space hole and weaken, or even suppress, a phase space clump. Also, the combined effect of drag and velocity space diffusion can produce the repetitive pattern of hook-shaped frequency chirping events shown in figure 65.

<figure><img src="image_63.png" /><figcaption>Figure 65. Hooked frequency spectrum of holes and clumps represents interplay of drag and diffusive collisions. Reproduced from [936].</figcaption></figure>

<figure><img src="image_64.png" /><figcaption>Figure 66. Typical examples of non-perturbative modes on MAST (a) and (b) and START (c). Reproduced from [715].</figcaption></figure>
### 7.13. Long-range sweeping

The initial theory for phase space holes and clumps was limited to the case of small frequency deviations from the bulk plasma eigenfrequency [702, 703]. However, there are multiple experimental observations of frequency sweeping events in which the change in frequency is comparable to the frequency itself [594, 709, 715]. Figure 66 shows examples of that. A non-perturbative theoretical formalism [957] is needed to interpret such a long-range frequency sweeping event.

Since the EPs typically have a much lower density than the bulk plasma, at first sight it seems unlikely that they can change the Alfvén eigenmode frequency significantly. However, a small but coherent group of EPs can generate a detectable signal at a different frequency than the bulk plasma eigenfrequency. Consider, for example, a modulated beam in a plasma. The initial modulation occurs spontaneously because of instability at the plasma eigenmode frequency.

However, as the coherent structure evolves due to dissipation, the trapped EPs gradually decelerate while maintaining their coherency. Consequently, the resulting frequency undergoes a significant shift from the initial frequency. This can be described as a nonlinear Bernstein-Greene-Kruskal (BGK) mode [958]. A solution of this type was obtained for a one-dimensional bump-on-tail model with immobile ions. The perturbed electrostatic potential $\phi$ takes the form

$$ \phi \equiv -\frac{1}{|e|} U(x - s(t), t). \quad (83) $$

The electron potential energy $U$ is a periodic function of its first argument $(x - s(t))$ and a slowly varying function of the second argument $t$. Additionally, the wave phase velocity $\dot{s} \equiv ds(t)/dt$ is a slowly varying function of time, characterized by a sweeping rate $\ddot{s}$.

For small deviations of $\dot{s}$ from $\dot{s}_0$ (early phase of frequency sweeping), a sinusoidal mode with constant amplitude occurs at the beginning of frequency sweeping [702, 957]. For large deviations of $\dot{s}$ from $\dot{s}_0$, the amplitude and the mode structure change significantly. Figure 67 illustrates the separatrix between the passing and trapped particles which changes its shape.

While the separatrix shrinks, some of the originally trapped EPs now become passing particles. The EPs that remain trapped decelerate to lower velocities and supply their energy to the wave. The power extracted from the EP population balances the power dissipated in the bulk plasma, and thus determines the rate of sweeping needed to compensate for collisional dissipation of the BGK-mode. This power balance condition reproduces the square root scaling of frequency sweeping during the initial phase of the sweep [702]. Later on, the mode phase velocity $\dot{s}$

<figure><img src="image_65.png" /><figcaption>**Figure 67.** Evolution of the phase-space bucket during sweeping event. The plot shows the initial separatrix (upper shaded area) and the shrunken separatrix at half of the initial mode phase velocity (lower shaded area). Reproduced from [28].</figcaption></figure>

deviates gradually from this square root scaling. This process is interpreted as a transformation of the initial plasma wave into a nonlinear EPM. EPMs might be generated by Alfvén wave instabilities in an analogous way [421, 496, 959]. Such a nonlinear modification of the mode structure appears to be essential to generate EPMs from AEs.

The analysis of the 1D electrostatic bump-on-tail problem leads to an analogous approach for understanding frequency-sweeping phenomena in tokamaks. In experiments, these events can be linked to the resonant excitation of TAEs. For a linear mode, the resonance condition is

$$ \omega - n\omega_{\phi}(P_{\phi}, P_{\theta}, P_{\psi}) - l\omega_{\theta}(P_{\phi}, P_{\theta}, P_{\psi}) = 0, \quad (84) $$

where $\omega$ is the mode frequency, $\omega_{\phi}(P_{\phi}, P_{\theta}, P_{\psi})$ and $\omega_{\theta}(P_{\phi}, P_{\theta}, P_{\psi})$ are the toroidal and poloidal transit frequencies, and $n$ and $l$ are integers. The canonical action-angle variable pairs $(P_{\phi}, \phi)$, $(P_{\theta}, \theta)$, and $(P_{\psi}, \psi)$ describe the integrable unperturbed motion. However, $(P_{\psi}, \psi)$ describes the Larmor radius gyration which is much faster than the wave period, so that gyration does not resonate with shear Alfvénic perturbations. For an isolated linear resonance, the perturbed particle Hamiltonian is a sinusoidal function of $\omega t - n\phi - l\theta$. The transition to the nonlinear case now generalizes the Hamiltonian to

$$ H = H_0 + U \left( \int_0^t \omega(\tau) \, d\tau - n\phi - l\theta, t \right), \quad (85) $$

where the function $U$ to be determined numerically is a not necessarily sinusoidal but still periodic function of its first argument. Note that $P_{\psi}$ and $P = lP_{\phi} - nP_{\theta}$ are constants of motion and that since the function $U$ evolves slowly, it should preserve an adiabatic invariant
for trapped particles. These three conservation laws suggest a relationship between the trapped particle distributions at any two locations of the resonance as illustrated in figure 68.

<figure><img src="image_66.png" /><figcaption>**Figure 68.** Transport of resonant particles during frequency sweeping. The shaded areas are snapshots of the moving resonant region in the momentum space. The shades of grey mark different values of the particle distribution function. The trapped resonant particles form a locally flat distribution across the resonance and preserve the value of their distribution function when the resonance carries them along the dotted lines. Reproduced from [28].</figcaption></figure>

The macroscopic quantities, like the perturbed EP pressure, now become known functions of the potential energy profile U and the unperturbed distribution. To determine the wave profile U, we can now solve a set of linear MHD equations for the bulk plasma response with an analytic nonlinear source from the EPs. The frequency sweeping rate can then be obtained from the power balance condition.

### 7.14. Multiple modes and global transport of EPs

A single unstable Alfvén mode tends to be benign in terms of global losses of EPs. The fundamental reason for this is that the resonances associated with a single mode occupy only a small fraction of the particle phase space. Many modes are usually needed to achieve resonance overlap in a large part of phase space and thereby produce global diffusion. When more than one unstable mode is present in the system, these modes are basically independent if their wave-particle resonances do not overlap. Then the waves should flatten the particle distribution function locally near each resonance. However, when the resonances overlap, the particle motion becomes stochastic, which allows individual particles to diffuse in phase space over many resonances and flatten the distribution function over a larger phase space area. This stochastic motion gives rise to both the velocity space transport and real

space transport. The relative importance of the two is different for different instabilities but the mathematical description of the resulting transport is essentially the same from the technical standpoint. Particle diffusion over a set of overlapping resonances can be described by quasilinear theory. The constants of motion ($\mathcal{E}$, $P_\phi$, $\mu$) of unperturbed particle orbits are no longer guaranteed to be constant if there are overlapping resonances. However, the particle magnetic moment $\mu$ remains almost constant for wave-particle interaction with AEs since the wave frequency is much smaller than the cyclotron frequency. In addition, the particle energy remains almost constant if the wave frequency is much smaller than the other two terms in the resonance condition

$$ \omega - n\omega_{\phi}(\mathcal{E}, P_{\phi}, \mu) - l\omega_{\theta}(\mathcal{E}, P_{\phi}, \mu) = 0, \quad (86) $$

which consequently almost balance each other. The waves then affect mostly the toroidal canonical angular momentum $P_\phi$ and hence the radial position of the particle orbit described by the flux surface in the poloidal cross-section of the tokamak. The width of the resonance (86) in $\delta P_\phi$ can roughly be estimated from

$$ \delta P_{\phi} \frac{\partial}{\partial P_{\phi}} \left( \omega - n\omega_{\phi}(\mathcal{E}, P_{\phi}, \mu) - l\omega_{\theta}(\mathcal{E}, P_{\phi}, \mu) \right) = \omega_{b}. \quad (87) $$

Here, $\omega_b$ is the nonlinear bounce frequency for a resonant particle in the wave field. Whereas $\omega_B$ in subsection 7.11 represents an electrostatic mode, $\omega_b$ represents the perturbed magnetic field of an electromagnetic mode. Both $\omega_B$ and $\omega_b$ are proportional to the square root of the mode amplitude.

The 1D quasilinear diffusion equation in $P_\phi$ is

$$ \frac{\partial f}{\partial t} - \frac{\partial}{\partial P_{\phi}} D \frac{\partial f}{\partial P_{\phi}} = -\nu(f - f_0). \quad (88) $$

$D$ is the diffusion coefficient, which is proportional to the wave intensity and is given by

$$ \frac{\partial D}{\partial t} = 2 \left( a \frac{\partial f}{\partial P_{\phi}} - \gamma_{\text{damp}} \right) D. \quad (89) $$

The damping rate $\gamma_{\text{damp}}$ sets an instability threshold and is assumed to be smaller than $\gamma_0$. $f_0$ is a classical equilibrium distribution function, and the Krook-type relaxation $-\nu(f - f_0)$ drives the distribution towards this equilibrium distribution. The gradient of this equilibrium distribution in $P_\phi$ can drive the waves unstable with a linear growth rate $\gamma_0 = a \frac{\partial f}{\partial P_\phi}$. The factor $a$ depends on the modes that resonate with the particles for a given value of $P_\phi$.

If there are no waves, the distribution function $f$ will converge to $f = f_0 \sim \frac{\gamma_0}{a} P_\phi$. But if waves are present, the steady-state solution of (88) and (89) gives $f = \frac{\gamma_{\text{damp}}}{a} P_\phi < f_0$, i.e. smaller than $f_0$. The diffusion coefficient can then be estimated as

$$ D \approx \frac{a}{\gamma_{\text{damp}}} \nu f_0 P_{\phi} \approx \nu \frac{\gamma_0}{\gamma_{\text{damp}}} P_{\phi}^2. \quad (90) $$
To obtain the resonance overlap constraint on the steady regime, we consider a set of barely overlapping resonances for which the correlation time is the inverse of the bounce frequency $1/\omega_b$, and the diffusion coefficient is

$$D \approx \omega_b (\delta P_\phi)^2 \\ \approx \frac{\omega_b^3}{\left\{ \frac{\partial}{\partial P_\phi} [\omega - n\omega_\phi(\mathcal{E}, P_\phi, \mu) - l\omega_\theta(\mathcal{E}, P_\phi, \mu)] \right\}^2}, \quad (91)$$

which can be written as

$$D \approx (\delta P_{\phi})^3 \frac{\partial}{\partial P_{\phi}} [\omega - n\omega_{\phi}(\mathcal{E}, P_{\phi}, \mu) - l\omega_{\theta}(\mathcal{E}, P_{\phi}, \mu)]. \quad (92)$$

If there are $N$ resonant modes in the full range of $P_\phi$, then the overlap condition becomes

$$D > \left(\frac{P_{\phi}}{N}\right)^3 \frac{\partial}{\partial P_{\phi}} [\omega - n\omega_{\phi}(\mathcal{E}, P_{\phi}, \mu) - l\omega_{\theta}(\mathcal{E}, P_{\phi}, \mu)]. \quad (93)$$

Equations (90) and (93) suggest that the resonance overlap condition requires a source strength at least given by

$$\nu \frac{\gamma_0}{\gamma_{\text{damp}}} > \frac{1}{N^3} P_{\phi} \frac{\partial}{\partial P_{\phi}} [\omega - n\omega_{\phi}(\mathcal{E}, P_{\phi}, \mu) - l\omega_{\theta}(\mathcal{E}, P_{\phi}, \mu)] \quad (94)$$

If the source strength is weaker than that, the global transport completely terminates, or it becomes intermittent if the individual modes can reach the overlap condition somewhere during their nonlinear evolution.

If there is no resonance overlap, the EPs are not transported from one resonance to the next, i.e., they do not move across the Kolmogorov-Arnold-Moser (KAM) surfaces between the resonances. If such transport barriers can be created on purpose, especially at the plasma edge, they may lead to good global confinement of EPs even if there is a local instability in the core.

We may also obtain an overlap condition for individually saturated modes, which follows immediately from equation (87) and the estimate for the nonlinear bounce frequency obtained in [960],

$$\omega_b \approx \gamma_0 \left( 1 + \frac{\nu}{\gamma_{\text{damp}}} \right). \quad (95)$$

This overlap condition for individually saturated modes is

$$\gamma_0 \left( 1 + \frac{\nu}{\gamma_{\text{damp}}} \right) > \frac{1}{N} P_{\phi} \frac{\partial}{\partial P_{\phi}} [\omega - n\omega_{\phi}(\mathcal{E}, P_{\phi}, \mu) - l\omega_{\theta}(\mathcal{E}, P_{\phi}, \mu)], \quad (96)$$

which is much more restrictive than equation (94). The substantial difference between equations (94) and (96) for large $N$ arises from the enhanced energy release per mode for overlapping resonances compared to isolated modes. Consequently, even if neighboring modes are linearly stable, an overlap of two closely spaced resonances could couple the neighbouring modes and trigger an avalanche-type relaxation event.

<figure><img src="image_67.png" /><figcaption>Figure 69. Cartoon illustrating intermittent quasilinear relaxation due to resonance overlap in the multi-mode regime. The resonances broaden and eventually overlap as the modes grow above the instability threshold. The EP population drops to subcritical values and then is replenished by the source to a metastable level that exceeds the linear threshold level due to separation of neighboring linear resonances. Reproduced from [28].</figcaption></figure>

During such an event, rapid quasilinear diffusion can decrease the EP density to a subcritical value below the linear instability threshold. This process is illustrated in figure 69. Subsequently, the waves will decay within a linear damping time, and the system will be quiescent until a particle source makes the EP population unstable again, triggering the next avalanche.

In intermittent diffusion, the bursts of different modes are synchronized due to the triggering effect. Due to these bursts, the EP population is close to the marginally stable level, in the steady as well as the intermittent quasilinear regimes. The primary difference lies in the time behavior of the turbulence level. The profile stiffness resulting from this marginal stability condition causes the turbulence level to self-adjust, maintaining a similar profile regardless of the strength of the particle source. The time-averaged transport coefficients are given by the injection rate and gradients of the marginally stable profile. However, the turbulence responsible for enhanced transport is caused by the small deviations away from marginal stability. Predicting the turbulence level based on a few macroscopic parameters such as the EP pressure or density gradient becomes challenging but might turn out to be unnecessary. Instead, a better route for theory-experiment comparison might be to
examine marginal stability constraints and resonance overlap criteria. This requires numerical tools for linear stability assessments for realistic magnetic configurations and plasma parameters.

## 8. 3D effects on EP confinement and losses

3D effects are unavoidably present in all toroidal magnetic fusion devices due to a variety of sources: the discreteness of toroidal field (TF) coils, 3D ferritic inserts, externally applied 3D fields for plasma control, internal MHD perturbation developed by the plasma itself or, in ITER, test blanket modules (TBMs). In an idealized, perfectly axisymmetric tokamak, the EPs move along their drift orbits which from closed loops in the poloidal plane and guarantee optimal confinement (section 2). In actual tokamaks, internal and external perturbations break the perfect axisymmetry and deteriorate the confinement. Energetic ions are particularly vulnerable to enhanced cross-field transport from symmetry breaking due to their large velocities and, consequently, wide orbits and low collisionalities. The 3D effects on energetic ion confinement have to be understood and mitigated to ensure good performance. Transport in 3D fields is further discussed in chapter 2 of the volume [14].

### 8.1. Introduction to 3D perturbations

ITER has 18 toroidal field (TF) coils which leads to a non-axisymmetric geometry. The resulting TF ripple in the magnetic field is larger than 1% in the plasma periphery in standard full-field scenarios. Such a magnetic ripple causes stochastic diffusion since the drift orbits no longer close upon themselves in the poloidal plane. For deeply trapped ions, it even leads to direct losses due to $\nabla B$-drift, leading to localized power loads on plasma-facing components that need to be avoided. Ferritic steel inserts placed between the TF coils and the plasma reduce the peak TF field at the plasma edge near the TF coil. The optimised ferritic inserts were developed and demonstrated to reduce the TF ripple by a factor of 4 at JT-60U [961].

However, experiments at JFT-2M have shown that ferritic inserts must be carefully designed to avoid enhanced transport from higher harmonics of the ripple [961, 964]. Additional geometrical constraints, such as the toroidal distribution of ports, may hamper the ferritic insert design and optimisation, possibly leading to localised enhanced losses.

Since the unmitigated TF ripple in ITER is unacceptably strong, ITER is equipped with ferritic inserts at each coil as illustrated in figure 70. In the 15 MA baseline scenario, the optimized ferritic inserts reduce the TF ripple substantially, see figure 71: the TF ripple at the outboard midplane separatrix is

mitigated to 0.3% except for near the NBI ports where it is mitigated to 0.6%.

ITER will be the first tokamak to study tritium breeding, a critical technology for the production of tritium fuel to sustain operation of a fusion power plant. It is planned to test a range of design concepts for test blanket modules (TBMs) in the equatorial port plugs. The TBMs will be fabricated using ferritic steels and will consequently affect the edge magnetic field. Two pairs of TBMs, at different toroidal locations, are available in ITER for testing tritium breeding, producing localised perturbations of the magnetic field. Like the TF ripple, these symmetry-breaking fields could deteriorate the confinement of energetic ions. It is also expected that structural components of future reactor systems will predominantly utilize martensitic stainless steels, which are ferritic and may induce symmetry-breaking effects.

Since ITER is foreseen to operate in H-mode, ELMs are expected and must be mitigated. External coils are widely used in present devices to suppress ELMs by means of externally applied resonant magnetic perturbations (RMPs). In most present tokamaks, RMPs are generated by two or three toroidal rows of window-frame coils. The current flowing through the coils is modulated toroidally, and the relative toroidal phase between the waveforms can be adjusted to vary its effect on the plasma. The relative toroidal phase between toroidal rows is thus used to modify the poloidal mode spectrum of the perturbation; this has been observed to affect the plasma stability in several experiments such as MAST [965], AUG [966] and DIII-D [967]. Further information on the interaction of ELMs and RMPs are found in chapters 3 and 5 of this volume [15, 17].

Furthermore, an efficient coupling between internal MHD fluctuations and externally applied 3D fields via particle transport, finite orbit width effects and mode couplings can lead to synergistic effects on EP transport with unexpected consequences. Recent advances in diagnostics and modeling techniques together with collaborative efforts between the tokamak and stellarator communities have led to significant progress in understanding the effects that symmetry breaking 3D effects have on tokamak plasma stability and related EP transport and losses [968]. In the following, the current state of experimental observations and modeling will be presented.

### 8.2. Experiments

8.2.1. Toroidal field coil ripple. Extensive experimental efforts carried out at JT-60U [969], JET [373,970], TFTR [24,971–973], JFT-2M [961, 974, 975] and Tore Supra [976, 977] have led to a solid understanding, allowing us to keep ripple losses at acceptable levels in
<figure><img src="image_68.png" /><figcaption>**Figure 70.** Illustration of spatial location of symmetry-breaking elements in ITER: finite TF coils (white), ferritic inserts (light grey), NBI ports (no color) and test blanket modules (light grey). Reproduced from [962].</figcaption></figure>

<figure><img src="image_69.png" /><figcaption>**Figure 71.** The toroidal field strength at the outer midplane at ITER, near the separatrix, as a function of toroidal angle showing the cases without (blue) and with (red) ferritic inserts. The ferritic inserts are found to be very effective in mitigating the TF ripple down to 0.3%, but around the neutral beam ports, $\phi \approx 60^\circ - 80^\circ$, the mitigation is incomplete at around 0.6%. At about $\phi \approx 30^\circ, 310^\circ$ and $350^\circ$, the influence of the three pairs of TBMs is visible. Reproduced from [963]. Later, the number of pairs of TBMs has been reduced to two.</figcaption></figure>

present tokamaks as well as to make robust predictions for future devices. For a given device, ripple losses are numerically predictable but highly dependent on the EP population and the *q*-profiles. The TF ripple can cause rapid losses of deeply trapped particles at the edge. These losses are usually lower for low-*q*, monotonic *q*-profiles than for high-*q*, reversed *q*-profiles, such as those used in “advanced tokamak” scenarios [24, 971–973, 978, 979]. Optimised ferritic inserts have successfully reduced the measured ripple losses at JT-60U [961] and EAST [980].

### 8.2.2. Test blanket modules.
The effect of TBM-like perturbations was studied using error fields produced by an external coil system at DIII-D [982]. The features of the ITER TBM perturbation were mimicked as closely as possible, but the magnitude of the perturbation was significantly larger. Heating

<figure><img src="image_70.png" /><figcaption>**Figure 72.** Infrared measurements of heat flux on DIII-D graphite tiles with TBM simulation coil turned off and on. Reproduced from [981].</figcaption></figure>

was observed on the protective tiles of the TBM mock-up surface when NBIs and the TBM fields were engaged (figure 72). The EP core confinement was not significantly affected. Different orbit-following codes predicted the formation of a hot spot on the TBM mock-up surface arising from beam ions deposited near the edge of the plasma [981, 983]. The codes are in good agreement on the total power deposited at the hot spot, predicting an increase in power with decreasing separation between the plasma edge and the TBM surface.

Simulations carried out for the TBM parameters in ITER full-field, standard scenarios indicate that the TBM-induced losses can be kept below 1% for alphas, below 2% for NBI ions, and below 3% for ICRF ions [963, 984]. However, in reduced magnetic field scenarios, alpha losses can be up to 10% and NBI losses can be up to 3%.

### 8.2.3. ELM mitigation coils.
Application of the RMP coils has been observed to cause significant EP losses at AUG [985, 986], DIII-D [985] and KSTAR [391]. Figure 73 shows the measured losses at DIII-D with modulated perturbation amplitude for a counter-
<figure><img src="image_71.png" /><figcaption>**Figure 73.** Midplane FILD signal showing modulation of EP loss correlated with rotating $n = 2$ field in DIII-D discharge #146121. Reproduced from [985].</figcaption></figure>

<figure><img src="image_72.png" /><figcaption>**Figure 74.** AUG discharge #28061. (a) Time trace of the heat load on the detector aperture. (b) 2D image of heat load due to ELM mitigation coils on FILD detector head. (c) 2D image of detector head temperature. Reproduced from [986].</figcaption></figure>

current and a co-current NBI. The measured power load due to the EP losses on a FILD detector head reaches up to 2 MW/m² as figure 74 shows [986]. However, the measured losses depend strongly on the existing EP distribution and on the poloidal and toroidal spectra of the applied perturbation. At AUG, striking differences in the velocity space of the EP losses have been measured by FILD systems for different NBI geometries (figure 75). Dedicated experiments and numerical simulations have revealed that the observed EP losses in the presence of externally applied RMPs are produced by an edge resonant transport layer with a high density of linear

and nonlinear resonances between the EPs and the externally applied static 3D fields [987]. In AUG [988] and DIII-D [989], the EP displacement during a single pass through the RMP-induced fields was measured using the “light-ion beam probe” [990, 991] technique. For comparison with theory, the plasma response was calculated by MARS-F and M3D-C1, respectively. In both studies, the EP displacement depends sensitively on the poloidal mode spectrum. In DIII-D, the displacement also depends sensitively on the normalized beta $\beta_N$, since this parameter impacts the plasma response to the perturbation.

<figure><img src="image_73.png" /><figcaption>**Figure 75.** AUG discharge #28061. Velocity-space of escaping ions measured by FILD1 with RMP coils off, (a), (c) and (e) and with RMP coils on, (b), (d) and (f). Reproduced from [986].</figcaption></figure>

8.2.4. *Synergistic effects between internal MHD fluctuations and external 3D fields.* Externally produced 3D perturbations predominantly affect only the edge EP population. However, if an efficient coupling to the edge exists, internal MHD fluctuations that otherwise would cause only a minor particle redistribution, could now cause significant losses. DIII-D experiments explored the synergy between EP transport caused by the simulated test blanket model of figure 72 and NTMs, AEs, sawteeth, and RMP fields. A definitive synergistic effect was observed at sawtooth crashes where, in the presence of the TBM, the localized heat flux at a burst increased from $0.36 \pm 0.27 \text{ MW-m}^{-2}$ to $2.6 \pm 0.5 \text{ MW-m}^{-2}$ (figure 76). Similarly, AUG experiments show that, under certain conditions, the edge resonant transport layer produced by RMPs can couple to NTM-induced EP losses to channel particles from the core of the plasma to the wall. NTM induced EP
losses were measured with FILD systems only when the RMP was active (figure 77) [986]. This effect must be taken into account when applying external RMPs in future devices.

Experiments at NSTX have shown that RMPs can also be used to control AEs modifying the resonant EP distribution. RMPs were observed to reduce the amplitude and frequency chirp of global Alfvén eigenmodes and to increase the bursting frequency [992]. More recently, experiments at AUG have shown that RMPs can be used to not only mitigate TAEs but also to excite them by modifying the gradients in the EP distribution that drive the modes unstable [993] (figure 78). To optimize the control capabilities of the externally applied 3D fields at certain plasma locations and to minimize the any losses that may be caused by the same symmetry-breaking fields, the location of the geometrical resonances, the targeted fast-ion phase-space volume, and the transport properties can be varied by means of the applied 3D fields spectrum and initial phase [987,994].

<figure><img src="image_74.png" /><figcaption>**Figure 76.** Measured heat flux near the TBM mockup coil (dashed rectangle) in a sawtoothing DIII-D discharge before (left) and after (right) application of the TBM field. Reproduced from [995].</figcaption></figure>

### 8.3. Modeling of 3D equilibria, stability and EP confinement

In the following, we will review the modeling tools available for 3D magnetic field effects on equilibrium, EP confinement, and stability. The efforts in modeling 3D effects have been increasing with the advent of large stellarator experiments. The same tools developed for stellarators can often be used to also model tokamak equilibria, EP confinement and stability, accounting for the 3D perturbation fields due e.g. TF coils, ferritic inserts, TBMs and ELM control coils.

#### 8.3.1. Modeling of 3D equilibria.
In non-axisymmetric devices, 3D equilibria with closed, nested flux surfaces are not guaranteed to exist. Magnetic islands can appear where the rotational transform, $\iota$, is a rational number, which can deteriorate the confinement. The

<figure><img src="image_75.png" /><figcaption>**Figure 77.** AUG discharge #28061. Spectrogram of (a) magnetic pick-up coil and (b) EP losses measured by FILD1. Reproduced from [441].</figcaption></figure>

VMEC code [345, 996] solves the equilibrium problem in general geometry by minimizing the plasma energy, assuming closed, nested flux surfaces. VMEC has proven useful for the computation of 3D tokamak equilibria with applied RMPs [997–1009] as well as 3D equilibria in ITER [1010]. In plasmas with a large EP population, anisotropic 3D equilibria can be calculated using energy principles [1011–1013]. The 3D equilibrium models HINT [1014–1017], PIES [1018–1020], and SIESTA [1021–1023] do not need the assumption of closed, nested flux surfaces, allowing magnetic islands to form. 3D equilibrium effects can also be addressed using linear and nonlinear time-dependent MHD models, which are useful to study time-dependent effects such as time-varying RMP fields [993, 1024, 1025].

The simplest approach to generate the 3D magnetic fields in tokamaks for the purpose of studying EP confinement (see next subsection) is to calculate the 3D perturbations from vacuum solutions and superimpose these on a 2D Grad-Shafranov equilibrium [984]. At the next level, VMEC can provide a self-consistent, finite $\beta$ plasma response to the external fields [1010, 1026]. For time-dependent external 3D fields, such as from ELM coils, the plasma response becomes dynamic, which has been computed in [985, 1027].