Main

On the road to a magnetic fusion power plant, a good confinement of a burning plasma, at temperatures >108 K and particle densities >1020 m−3, is mandatory. In tokamaks, these conditions are typically achieved in the high-confinement regime (H-mode)1. The H-mode has been selected as the most efficient and promising operational regime for the International Thermonuclear Experimental Reactor (ITER), which is planned to start operation within the next decade and should prove the scientific and technical feasibility of magnetic confinement fusion. H-modes are, however, characterized by steep edge pressure gradients that are accompanied by periodic edge plasma relaxations in the form of edge-localized modes (ELMs)2 that, if allowed to grow unabated, can severely damage the first wall of ITER and future fusion power plants3. Extensive efforts are dedicated towards the development of ELM control techniques and ELM-free high-performance regimes3,4,5,6,7,8,9,10. Linear magnetohydrodynamic (MHD) stability analyses indicate that ELMs are coupled peeling–ballooning modes driven by the plasma edge current and pressure gradients11,12,13,14. Nonlinear MHD models successfully reproduce some relevant experimental observations such as ELM cycles15, filamentary structure16,17,18 and collapse of the edge pressure profile19.

However, a detailed understanding of the ELM behaviour and consequences in a burning plasma with a notable fraction of energetic (suprathermal) ions is missing. Energetic ions are an essential source of momentum and energy that must be kept well confined until they slow down to the bulk plasma through Coulomb collisions20,21,22. Nevertheless, energetic ions are prone to show a rich variety of wave–particle interactions due to their large velocities and long mean free paths that can lead to an efficient exchange of energy and momentum with a broad spectrum of MHD fluctuations23,24,25,26,27. In fact, ELM-induced fast-ion losses are routinely observed in several tokamaks28 with scintillator-based fast-ion loss detectors (FILDs)29. Recent FILD observations of fast-ion acceleration during ELM crashes in the Axially Symmetric Divertor Experiment) Upgrade (AUG) tokamak highlight the strong interplay between the energetic particle population at the plasma edge and the electromagnetic perturbation developed during an ELM crash30. Moreover, many of the natural no-ELM or small-ELM regimes, such as the Enhanced D-Alpha high-confinement mode31 or type-II ELM regime (more recently dubbed quasi-continuous exhaust regime)32,33, feature plasmas with high collisionalities and, thus, negligible energetic particle contents. Extended Data Fig. 1 shows the plasma density perturbation during an ELM crash modelled with the MEGA code34 together with a typical fast-ion orbit (in white) at the edge of the AUG tokamak. Eventually, its orbit ends on the wall and deposits its energy at the plasma facing component. A good understanding of the interaction between energetic ions and ELMs is, thus, of paramount importance to optimize ELM control techniques while minimizing fast-ion losses.

The three-dimensional nonlinear hybrid kinetic–MHD code MEGA has been used to self-consistently model an ELM crash in the presence of a realistic energetic ion distribution. These simulations provide evidence of the key role that energetic ion kinetic effects play in the spatio-temporal structure of ELMs. The wave–particle resonances are responsible for the net energy exchange between fast ions and ELMs. Comparisons of experimental measurements of the magnetic perturbations and fast-ion losses induced by large type-I ELMs are in qualitative agreement with the simulations. The results presented here provide a better understanding of the physics underlying an ELM crash in the presence of energetic particles and contribute to developing robust ELM control techniques and/or ELM-free regimes for future burning plasmas with fusion-born alpha particles.

Experiments show that ELM and energetic particle dynamics are linked

The experiment analysed here has been performed with toroidal magnetic field Bt = 2.5 T, plasma current Ip = 0.8 MA, low collisionality (\({\nu }_{{\rm{i}}}^{* } < 0.5\) inside the pedestal top, where \({\nu }_{{\rm{i}}}^{* }\) is the ion collisionality, defined as the effective collision frequency normalized to the trapped particle bounce frequency) and 5 MW of neutral beam injection (NBI) as the main heating and fast-particle source.

Figure 1a shows an ELM-synchronized spectrogram of the radial magnetic field fluctuations (∂Br/∂t) measured with a Mirnov coil at the tokamak first wall for a plasma discharge with fast ions. In the experiment, the pre-ELM crash phase is characterized by a broadband frequency up to ~150 kHz but dominated by low frequencies well below 50 kHz. The ELM crash itself (at t = 0 ms) appears with a substantially larger perturbation amplitude and broadband frequency up to ~300 kHz. Right after the ELM crash (t ≈ 1.8–3.5 ms), the magnetic perturbation relaxes in amplitude and the frequency drops to f < 50 kHz. The main features observed here in the ELM magnetic spectrograms are, in fact, ubiquitous to most tokamak H-mode plasmas35,36. Figure 1b,c shows the spectrogram and temporal evolution, respectively, of the fast-ion losses measured during the ELM crashes shown in Fig. 1a. The FILD spectrogram shows striking similarities to the magnetic spectrogram. Whereas during the pre-ELM crash phase, low frequencies up to ~100 kHz dominate the spectrogram, during the ELM crash, large-amplitude high-frequency (up to ~300 kHz) losses are observed.

Fig. 1: ELM crash in AUG shot #37700.
figure 1

a,b, ELM-synchronized spectrograms of measured radial magnetic field fluctuations (∂Br/∂t, labelled as magnetics) for a plasma discharge with fast ions (a) and fast-ion losses at the low field side (LFS) midplane measured by the FILD (b), both for t = 4.0–4.2 s. c, The temporal evolution of the measured ELM-induced fast-ion losses (grey) during the time window used for the spectrograms shown in a and b, with the losses observed during a single ELM at t = 4.102 s highlighted in black. d, Thermal particle losses measured by \({I}_{{\rm{div}}}\), that is, the thermo-currents flowing through the outer divertor.

Full size image

The amplitude of the high-frequency losses during the ELM crash decays to their nominal levels before the ELM crash in about 2–3 ms in the experiment presented here. Figure 1c shows the time traces of the fast-ion losses used to construct the FILD spectrogram (in grey) with that of a single ELM highlighted in black. Several fast-ion bursts (filaments) are commonly observed during each ELM. Finally, Fig. 1d shows the thermal particle losses measured by the thermo-currents flowing through the outer divertor during the ELMs. While the maximum thermal particle losses peak slightly after the ELM crash in the magnetic or FILD spectrogram, the thermal particle losses decay with the same time scale as the high-frequency fast-ion losses in the ELM-synchronized FILD spectrogram.

Simulated ELM crash in the presence of energetic particles

To understand the experimental observations, the nonlinear interaction between energetic ions and ELMs has been modelled using the hybrid kinetic–MHD code MEGA with realistic plasma equilibrium, thermal particle profiles and off-axis NBI energetic particle distributions, as seen in Extended Data Fig. 2. The simulation setup is described in Methods. The temporal evolution of the n = 10 (n denotes the toroidal mode number) mode energy (kinetic + magnetic)37 for a MHD single-n (including only n = 0 and 10 modes) simulation, that is, without fast ions, and three hybrid single-n simulations with different NBI injection energies (denoted by Ebirth hereafter) is shown in Fig. 2a. As shown in the inset, the linear growth rate of the mode decreases when Ebirth is increased, which delays the ELM crash. During the nonlinear phase of the hybrid kinetic–MHD simulations, the ballooning mode extends to inner regions (radial coordinate, ρpol < 0.9) and the mode energy takes values notably larger compared with the MHD case. In Fig. 2a, the dashed line represents the energy of the n = 10 mode (En=10) for ρpol > 0.9. We have also performed multi-n simulations (including n = 0, ⋯, 10 modes) to account for the nonlinear coupling of the modes. Figure 2b,c shows the time evolution of En, n = 1, ⋯, 10, for the MHD and hybrid simulations, respectively. The hybrid multi-n simulation has been performed for Ebirth = 60 keV. In the multi-n MHD simulation, n = 9 and n = 10 are the most unstable modes. However, in the hybrid simulation, although n = 10 begins as the most unstable harmonic, n = 8 becomes the most unstable mode in the nonlinear stage of the simulation with an energy almost five times larger than in the MHD simulation, indicating that the most unstable mode number is shifted to a lower one in the presence of fast-ion kinetic effects. Figure 2b,c also shows 0.5 × ∑nEn. In both the single- and multi-n MHD simulations, En=10 and ∑nEn take similar values, but in the hybrid simulations with Ebirth = 60 keV, the peak value of En=10 is larger than that of ∑nEn. Notice that the mode energy obtained in the multi-n simulations is divided by 10 because the mode energy is the volume integral of the mode energy density37 and, as explained in Methods, the volume enclosed by the simulation domain in the single-n simulations is ten times smaller than in the multi-n simulations. Note that the dominant mode numbers simulated here are in the range of those observed during ELMs in AUG plasmas38. Figure 2d shows how the thermal plasma pressure gradient, fast ions, viscosity and nonlinear mode coupling contribute to En=8 temporal evolution in the hybrid multi-n simulation with Ebirth = 60 keV. Whereas the ELM is, as expected, triggered by the steep edge thermal plasma pressure gradient, the abrupt crash is dominated by a large energy exchange between the mode and fast ions. The nonlinear mode–mode coupling contribution to the drive and damping of the n = 8 mode is included in Fig. 2d. This can be stabilizing or destabilizing depending on the relative phase of the different fluctuations’ terms. The nonlinear energy transfer plays a role in the evolution of high-n modes, as its contribution is comparable to those coming from the thermal plasma pressure gradient, power exchange with fast ions and viscosity.

Fig. 2: Simulated ELM crash.
figure 2

a, The time evolution of n = 10 toroidal mode energy (En=10) in single-n simulations (n denotes the toroidal mode number). The solid line corresponds to integration in the whole simulation domain and the dashed line to integration from the radial coordinate ρpol > 0.9. The case without fast ions (blue solid) is shown for reference. Inset: the linear growth rate (γl, normalized to the Alfvén frequency at the magnetic axis, ωA) versus NBI energy (Ebirth). b,c, The time evolution of each mode energy (En) and of the sum of the mode energies (∑nEn, multiplied by 0.5) in multi-n MHD (b) and hybrid (for Ebirth = 60 keV) (c) simulations. d, Different contributions to En=8 temporal evolution in the hybrid multi-n simulation shown in c. Ph (blue) denotes the energy transfer between n = 8 mode and energetic particles. The thermal plasma pressure gradient (red) is the volume integral of u ⋅ ∇p (dot product between plasma velocity and thermal plasma pressure gradient). ‘Viscous term’ is the dissipation brought by viscosity (black), and ‘NL mode coupling’ corresponds to the energy transfer between n = 8 and the rest of the toroidal modes (grey).

Full size image

The spatial structure of the thermal plasma pressure perturbation is shown in a poloidal plane in Fig. 3 at a time with the maximum perturbation amplitude, for a MHD simulation (Fig. 3a) and a hybrid simulation with Ebirth = 60 keV (Fig. 3b), both for the multi-n simulations. In the presence of energetic ion kinetic effects, the perturbation structure is radially extended to inner regions and strongly sheared. The radially extended ballooning structure causes a considerable fast-ion redistribution over a large plasma radius. This is observed in Fig. 3c, which shows the fast-ion pressure profile on a poloidal plane at the same timepoint as shown in Fig. 3b. The poloidal structure of the fast-ion pressure resembles the ballooning structure of the ELM with fingers in the energetic particle pressure extending into the far scrape-off layer.

Fig. 3: Two-dimensional poloidal cross-sections.
figure 3

a,b, Bulk plasma pressure perturbation (δp) in a multi-n MHD simulation at t = 0.16 ms (a) and a multi-n hybrid simulation with energetic particles with NBI energy Ebirth = 60 keV at t = 0.16 ms (b). The white line shows the trajectory of a characteristic resonant energetic particle. c, Fast-ion pressure profile (pEP) for Ebirth = 60 keV at t = 0.16 ms. d,e, Zooms of a (d) and b (e) at the plasma edge. In the figures, R and z correspond to the coordinates of the cylindrical coordinate system.

Full size image

The impact of the fast-ion kinetic effects on the spatial frequency pattern of the ELM crash is investigated in Fig. 4, which shows the spectrograms of n = 8 radial velocity in the multi-n MHD and hybrid simulations. Whereas in the MHD simulation the mode frequency takes rather low values (<20 kHz), in hybrid simulations, and during the nonlinear phase, the mode takes frequencies up to 200 kHz. In the presence of energetic ions, the mode becomes radially broad, as shown in Fig. 4b. The boost in frequency observed in the hybrid simulations is due to the imbalance of forces that appear in the MHD momentum equation. Our analysis reveals that the imbalance between the MHD and energetic particle forces, −∇p + j × B and \(-{{\bf{j}}}_{{\rm{h}}}^{{\prime} }\times {\bf{B}}\), respectively, is the main reason for the boost of frequency observed in the hybrid kinetic–MHD simulations. Likewise, the ELM radial broadening in the presence of fast ions is dominated by their drive (Supplementary Fig. 1).

Fig. 4: Spatial frequency ELM profile.
figure 4

a,b, Spectrograms of the radial velocity (ur) of n = 8 mode in the MHD (a) and hybrid (b) multi-n simulations in the radial coordinate (ρpol) and frequency (f) plane. The hybrid simulation shown in b considers the NBI energy Ebirth = 60 keV and βEP = 1.0%, where βEP is the ratio between fast-ion pressure peak and magnetic pressure.

Full size image

Interaction mechanism between ELMs and energetic particles

To understand the interaction between ELMs and fast ions, Fig. 5 shows the power transfer between the ELM and the energetic particles for Ebirth = 60 keV at t = 0.16 ms, when En=8 is largest in the hybrid multi-n simulation. The power transfer is given by39

$${P}_{{\rm{h}}}\left({\bf{x}},{\bf{v}}\right)=\sum _{l}{w}_{l}\frac{{\mathrm{d}}{E}_{l}}{{\mathrm{d}}t}S({\bf{x}}-{{\bf{x}}}_{l},{\bf{v}}-{{\bf{v}}}_{l}),$$
(1)

where wl = Vlδfl is the weight of the lth particle, Vl is the phase-space volume, δfl is the fast-ion distribution perturbation, \(\frac{{\mathrm{d}}{E}_{l}}{{\mathrm{d}}t}\) is the time derivative of the kinetic energy of the lth particle and S is a shape function that allows to select particles with x and v coordinates in the phase space. Since Ph < 0 at t = 0.16 ms, the ELM is extracting energy from the fast ions.

Fig. 5: Wave–particle interaction mechanism.
figure 5

Power transfer (Ph) for the hybrid simulation with NBI energy Ebirth = 60 keV at t = 0.16 ms, when En=8 is largest in the hybrid multi-n simulation. The dashed white and green lines denote the resonance conditions with \(n=10\left[p=\left(-9,-13\right)\right]\) and \(n=8\left[p=\left(-7,-10\right)\right]\), respectively, for μ = 1.6 × 10−15 J T−1 (μ is the magnetic moment of the particles). The black contours represent the fast-ion distribution. E and Pϕ denote the particle energy and toroidal canonical momentum coordinates, respectively.

Full size image

The nature of the interaction between ELMs and fast ions shown here has been investigated by computing the possible wave–particle resonances for the energetic particle phase-space volumes with the largest power exchange from Fig. 5. The resonance condition between the energetic particle drift orbits and a wave is given by40,41

$${\omega }_{n}-n{\omega }_{\varphi }-p{\omega }_{\theta }\approx 0,$$
(2)

where ωn is the frequency of the mode at the corresponding timepoint, ωφ and ωθ are the toroidal and poloidal particle orbital frequencies, respectively, and p is the bounce harmonic. The dashed lines shown in Fig. 5 represent the resonance conditions for the main bounce harmonics with n = 10 (white) and n = 8 (green) and for a given magnetic moment μ. The alignment of the power exchange structures with wave–particle resonances present in the plasma indicates that the wave–particle interaction observed here is likely to be mostly resonant. In fact, the orbital frequencies, ωφ and ωθ, for the energetic particles at the edge of the AUG plasma lie, mostly, within the range of frequencies observed in experimental magnetic and FILD spectrograms during ELM crashes. The trajectory of a typical resonant (passing) orbit with the ELM studied here is overplotted in white in Fig. 3b. The resulting power transfer observed in Fig. 5 is probably affected by a resonance overlap between the different modes considered in the simulation, given the closeness of the resonances.

The hybrid multi-n simulations presented here reproduce some outstanding experimental observations such as the ELM spatio-temporal structure, the filamentary and frequency pattern of the measured fast-ion losses and the temporal evolution of the measured thermal particle losses during the entire ELM crash. Figure 6 shows the synthetic signals of a magnetic pick-up coil, a FILD system and a thermal particle loss detector developed in MEGA for comparisons with the experimental observations shown in Fig. 1. The synthetic magnetic spectrogram shown in Fig. 6a qualitatively reproduces the overall frequency pattern shown in Fig. 1a with a large frequency band up to ~200–250 kHz dominating the crash at t ≈ 0.16 ms and with rather lower frequencies during the pre- and post-ELM crash. Similarly, the FILD spectrogram shown in Fig. 6b is dominated by a large frequency band up to ~200–250 kHz that decays in amplitude and frequency right after the crash to frequencies below 100 kHz. The temporal evolution of the synthetic FILD signal depicted in Fig. 6c shows several fast-ion filaments during the ELM crash, similar to the experiment (Fig. 1c). In contrast to the filamentary structures observed in the temporal evolution of the fast-ion losses, the synthetic signal of the thermal ion losses at the outer divertor shows a single well-defined peak extended over the entire ELM crash as in the experimental observation (Fig. 6d). The synthetic thermal ion losses are obtained by integrating the bulk plasma density over a small surface near the outer divertor, as the thermo-currents shown in Fig. 1 are measured at the outer divertor. Comparing the divertor current signal in Fig. 1d with the synthetic thermal ion losses in Fig. 6d, we note that the MHD model is accurate enough to describe the evolution of the divertor current. Note that the simulated ELM duration is one order of magnitude smaller than that observed experimentally. This is in accordance with the resistivity value used for these simulations, η = 20 × ηSpitzer (Supplementary Fig. 2 and ref. 42).

Fig. 6: Synthetic diagnostics for the multi-n hybrid kinetic–MHD simulation.
figure 6

a,b, Spectrograms of radial magnetic field fluctuations (∂Br/∂t) (a) and fast-ion losses (b) measured at the low field side (LFS) midplane measured by a synthetic FILD. c, The temporal evolution of the measured ELM-induced fast-ion losses. d, The temporal evolution of thermal particle losses to the outer divertor measured by a synthetic thermal ion loss detector (TILD). In a–d, the NBI energy is Ebirth = 60 keV.

Full size image

Significance of the interaction between energetic particles and ELMs in ITER

To discuss the implications of these results for ITER, a preliminary estimation of an efficient resonant interaction between ELMs and fast ions has been performed for representative energetic particle distributions and expected ELM parameters43. To this end, the Larmor radius ρL of fusion-born alpha particles and NBI ions has been compared with the expected edge ballooning mode wave number \({k}_{\perp }\approx \frac{nq}{r}\) (ref. 44), with q being the safety factor and r the radial coordinate. If ρLk⊥ < 1, a substantial interaction between the mode and fast ions is expected26. For the standard H-mode of ITER, it is observed that ρLk⊥ < 1 (with toroidal mode numbers n = 3, ⋯, 10) for both NBI ions (Ebirth = 1 MeV) and alpha particles. The mode numbers examined here are in the range of those expected during ELMs in ITER if they are not mitigated45. For the same particles and mode numbers, the drift orbit resonance condition is satisfied for a broad range of bounce harmonics p. This indicates that a considerable, and even stronger than in present machines, interaction between ELMs and energetic particles is expected in ITER. Whereas the majority of the NBI ions in ITER are born near the separatrix, thus suggesting a strong interaction with ELMs, most fusion-born alpha particles are born on-axis, suggesting a weak interaction with edge perturbations. A strong on-axis peaked alpha particle distribution is expected to trigger a broad spectrum of energetic particle-driven instabilities that will lead to a large alpha particle transport over the entire plasma radius46,47,48,49,50. Such a transported alpha particle population might interact with unabated ELMs in ITER.

Methods

Nonlinear hybrid kinetic–MHD model MEGA and simulation setup

The hybrid kinetic–MHD simulations of ELMs presented in this manuscript have been performed with the code MEGA34. MEGA calculates the interaction between MHD instabilities and energetic particles self-consistently by adopting a drift-kinetic description51,52,53 that couples the nonlinear MHD equations, which evaluate bulk plasma dynamics, with a δf method that describes the energetic particle population. In the hybrid simulations, finite Larmor radius effects of energetic particles are included by using the standard technique of four-point gyroaveraging around a particle’s guiding centre54.

The equilibrium profiles that are used as initial conditions for the bulk plasma and geometry are taken from AUG shot #33616 at 7.2 s, just before the ELM crash (that is, pre-ELM)35. Extended Data Fig. 2 shows the profiles of the thermal electron density (ne), thermal electron temperature (Te) and the safety factor (q) as a function of the radial coordinate ρpol. Here, \({\rho }_{{\rm{pol}}}=\sqrt{{\varPsi }_{{\rm{N}}}}\) with ΨN being the normalized poloidal flux with respect to the last closed magnetic flux surface. The pressure-constrained pre-ELM equilibrium has been reconstructed with the CLISTE code55. Extended Data Fig. 2a shows the experimental uncertainties of ne, Te and q. The error bars of q at the core (ρpol < 0.3) are not well defined because of the limitation of the available diagnostics. However, such error bars at the core do not affect the observations at the edge of the plasma reported in this paper. In the MEGA simulations presented here, the resistivity is given by \(\eta (T)={\eta }_{0}{\left(T/{T}_{0}\right)}^{-3/2}\), where T0 = Te,0 + Ti,0 = 6.6 keV is the temperature and η0 = 10−7 Ωm ≈ 20 × ηSpitzer is the resistivity, both at the magnetic axis. The viscosity follows the same profile, leading to a constant magnetic Prandtl number, \({\Pr }_{{\rm{m}}}=\frac{\nu }{\eta /{\mu }_{0}}=10\). The perpendicular thermal diffusivity is inferred from charge exchange measurements with \({\chi }_{{\perp }_{0}}=3.6\,{{\rm{m}}}^{2}\,{{\rm{s}}}^{-1}\) and a radial profile that mimics an edge transport barrier56,57. The parallel thermal diffusivity is given by the ratio \({\chi }_{{\parallel }_{0}}/{\chi }_{{\perp }_{0}}=1{0}^{5}\) at the magnetic axis, as used in previous ELM simulations18, and with a radial dependency given by \({\chi }_{\parallel }={\chi }_{{\parallel }_{0}}\times {\left(T/{T}_{0}\right)}^{5/2}\) (ref. 58). Notice that both the MHD and hybrid kinetic–MHD simulations shown in this work consider the same resistivity, parallel and perpendicular thermal diffusivity profiles to highlight the impact of fast-ion kinetic effects on ELMs. Moreover, the impact of parallel gradients on the fast-ion kinetic effects on ELMs reported in this work is negligible because, in the MEGA code34, the effect of energetic particles on the thermal plasma is accounted for in the MHD momentum equation of the thermal plasma through the energetic particle force term, \(-{{\bf{j}}}_{{\rm{h}}}^{{\prime} }\times {\bf{B}}\), which acts in the direction perpendicular to the magnetic field. In addition, the fast-ion distribution considered in the simulations is toroidally axisymmetric.

The initial fast-ion distribution is modelled as an off-axis slowing-down NBI energetic ion distribution. The spatial dependence of the fast-ion distribution is given by the fast-ion pressure profile as

$${p}_{{\rm{EP}}}={\beta }_{{\rm{EP}}}\frac{{B}_{0}^{2}}{2{\mu }_{0}}\exp \left[-{\left(\frac{{\varPsi }_{{\rm{N}}}-{\varPsi }_{{{\rm{N}}}_{0}}}{{\sigma }_{{\varPsi }_{{\rm{N}}}}}\right)}^{2}\right],$$
(3)

and the slowing-down NBI component is given by

$$\begin{array}{rcl}I(v,\varLambda )&=&\frac{1}{{v}^{3}+{v}_{{\rm{crit}}}^{3}}\frac{1}{2}{\rm{erfc}}\left(\frac{v-{v}_{{\rm{birth}}}}{{\varDelta }_{v}}\right)\\ &&\exp \left[-{\left(\frac{\varLambda -{\varLambda }_{0}}{{\varDelta }_{\varLambda }}\right)}^{2}\right].\end{array}$$
(4)

In equations (3) and (4), βEP = 0.01 is the ratio between fast-ion pressure peak and magnetic pressure, μ0 = 4π × 10−7 H m−1 is the vacuum magnetic permeability, B0 = 2.5 T is the magnetic field at the magnetic axis, ΨN is the normalized poloidal flux, \({\varPsi }_{{{\rm{N}}}_{0}}=0.55\) is the centre of the off-axis, \({\sigma }_{{\Psi }_{{\rm{N}}}}=0.3\) is the spatial width, vcrit is the critical velocity for the collisional friction of fast ions with thermal electrons and ions being equal59, vbirth is the NBI injection velocity, erfc denotes the complementary error function, Δv = 0.05vA is the distribution width in velocity space, vA = 4.4 × 106 ms−1 is the Alfvén speed at the magnetic axis, Λ0 = 0.5 is the pitch angle for the distribution peak and ΔΛ = 0.2 is the distribution width. Extended Data Fig. 2a also shows the energetic particle pressure (pEP) profile. The initial fast-ion distribution in energy (E) and normalized pitch-angle (Λ = μB0/E) space is shown in Extended Data Fig. 2b. A realistic fast-ion distribution in position and velocity space has been implemented using the analytical formula shown in this section. For the sake of clarity, to isolate the interaction between fast ions and ELMs, and to save computing resources in the simulations presented here, the fast-ion density in the plasma core has been artificially reduced and maximized at the edge. Single-n hybrid simulations with a realistic fast-ion distribution in the region of interest \(\left({\rho }_{{\rm{pol}}} > 0.5\right)\) as calculated with the TRANSP code60,61 have been carried out to check that this simplified distribution does not affect the findings presented in this Article.

In MEGA, the time evolution is calculated using the fourth-order explicit Runge–Kutta method, and the spatial derivatives are evaluated using a fourth-order finite difference scheme. The boundary conditions are periodic in the toroidal direction. At the boundaries of the simulation domain, the perfect conducting no-slip wall boundary conditions are applied. In the simulations reported in this Article, the grid resolution is NR × Nφ × Nz = 512 × 16 × 512 in the single-n and NR × Nφ × Nz = 512 × 160 × 512 in the multi-n simulations, respectively. The time step is Δt = 4.2 × 10−11 s, and the number of computational particles for the kinetic module is NCP = 1.8 × 106 in the single-n simulations; in multi-n simulations, NCP = 1.8 × 107. In the single-n simulations, the toroidal angle ranges from 0 to 2π/N, with N = 10. In the multi-n simulations, the toroidal angle ranges from 0 to 2π. Therefore, in the multi-n simulations, the volume enclosed by the simulation domain is ten times larger than in the single-n simulations. A Fourier filter is applied in the simulations and, therefore, only mode numbers n = 0 and 10 are studied in the single-n simulations, whereas all n = 0, ⋯, 10 modes are retained in the multi-n simulations. The digital filter (smoother + compensator) described in page 438 of ref. 62 is applied in the simulations to reduce the numerical noise.

FILD

The scintillator-based FILD acts as magnetic spectrogram and uses the magnetic field to scatter charged particles that then impact on different positions on a scintillator plate, depending on their energy and pitch angle29,63. The light emitted by the scintillator material is directed by means of an optical setup to the data acquisition systems. The light beam is divided with a beam splitter and guided to a photomultiplier array and to a charge-coupled device camera. The measurements from the photomultiplier array with a high temporal resolution (down to 1 μs) are used in this work to analyse the fast-ion losses during an ELM cycle. The FILD spectrogram shown in this Article has a temporal resolution of 512 μs and a frequency resolution about 30 Hz. The spatial resolution of FILD is not defined because it is placed at a fixed location.

Magnetic spectrograms

Magnetic pick-up coils are one of the basic diagnostics to track MHD instabilities by measuring magnetic field perturbations64,65. The measurement principle is based on Faraday’s law; a voltage is induced in a wire loop when immersed into a temporally varying magnetic field. The Fourier transformation of the coil signals provides magnetic frequency spectrograms that allow to identify the frequency of the MHD instability. At AUG, several magnetic pick-up coil arrays are installed at the low field side and at the high field side of the tokamak and allow instabilities with frequencies up to 2 MHz to be detected. The magnetic spectrogram shown in this work has a temporal resolution of 256 μs, and its frequency resolution is about 30 Hz. The spatial resolution of this diagnostic is not defined because it is placed at a fixed location.

ELM synchronization method

The measured data from magnetic pick-up coils and FILD measurements are synchronized with respect to the ELM onset times. To this end, the divertor shunt current measurements are first used to identify the onset times of individual ELM crashes within a given time window of a discharge. From these measurements, the ELM crashes become apparent as spikes in the signal because heat and particles flushed out by the ELMs increase the currents flowing onto the divertor. The start times of the ELM crashes are then determined by the onset of the rise of the shunt current signal. Once the ELM start times are identified, a correlation analysis of the shunt current rise profiles is performed. This is done to select only those ELMs with similar characteristics for synchronization. Here, only ELMs with a correlation factor above 0.7 with respect to the initial ELM crash in the time window are selected. In this way, small ELMs and other small plasma bursts are ignored. Once a set of the most similar ELMs has been identified, a new time base for the ELM-synchronized signal is created. This is done by subtracting the selected ELM onset times from the original time base of the corresponding signal and rearranging the data chronologically. In this way, timepoint 0 in the synchronized time base corresponds to the ELM onset time, with negative and positive timepoints corresponding to the times before and after the ELM onset, respectively.