Gyrokinetic investigation of the nonlinear interaction of Alfvén instabilities and energetic particle-driven geodesic acoustic modes

The next generation of fusion relevant machines (ITER,^1,2 DEMO³) will be characterized by a large population of energetic particles (EPs). These are fusion products and charged particles generated by external power sources (like neutral beam injection, NBI). EPs are characterized by velocities much higher than the thermal velocity of the particles of the background plasma, $vEP≫vth,i,e$⁠, where the subscripts i and e refer, respectively, to the main ion population and to the electrons present in the plasma. The transport and confinement of the EPs in Tokamaks and the physics related to EPs are of primary importance in order to achieve self-heating plasma. In fact, EPs can drive unstable, via wave–particle interactions, symmetry breaking electromagnetic perturbations whose presence can redistribute the EPs population, expelling them out of the plasma before they can thermalize.⁴ This can consequently lead to a less effective heating. In addition, the violent migration of the EPs toward the walls caused by generated short but intense perturbations dubbed abrupt large-amplitude events⁵ can possibly damage the machine. The understanding of the intensity of the fields of the saturated induced instabilities together with the transport of the EPs represents a key topic to reach fusion.

In fusion devices, the EPs have a characteristic periodic motion (bounce frequency, transit frequency, etc.) of the same order of magnitude as that typical of the shear Alfvén waves (SAWs).⁶ When the drive provided by the EPs exceeds the damping of these plasma fluctuations, a broad spectrum of Alfvén waves (Alfvén modes, AM) can be excited.^7,8 These driven instabilities are classified into two types: Alfvén eigenmodes (AEs) and energetic particle modes (EPMs). The AMs are characterized by frequencies located inside the frequency gaps of the SAWs continuum spectrum created by the field geometry and by plasma nonuniformities.^9,10 The EPMs are nonnormal modes of the SAWs continuum spectrum, emerging as discrete fluctuations at the frequency where wave-EP power is maximized.¹¹ The linear and nonlinear dynamics of SAWs-driven instabilities has been reviewed in Ref. 4.

In addition, EPs can drive unstable modes with frequencies comparable with that of the geodesic acoustic modes (GAMs)¹² characterized by $(m,n)=(0,0)$ scalar potential and $(m,n)=(1,0)$ up-down antisymmetric density perturbation (being m and n, respectively, the poloidal and toroidal mode numbers). These driven modes are the energetic particle geodesic acoustic mode (EGAM), excited via free energy associated with velocity space gradients in the EPs distribution, as it has been shown analytically in Ref. 13. Studies on the nonlinear dynamics of EGAMs have been recently presented in Refs. 14 and 36.

The comprehension of the dynamics of the great zoology of modes presents in plasma, their interaction and the redistribution of EPs, is crucial to understand the properties of burning plasma. The study of the saturation levels of these instabilities is fundamental to be able to be predictive regarding the intensity of the fields that will be present in future reactors. This motivates the interest in the study of the nonlinear evolution of these instabilities. Traditionally, the two main saturation mechanisms of such instabilities have been recognized as (a) reduction of the drive due to the redistribution of the EP population¹⁵ and (b) the transfer of energy to other modes via mode–mode coupling mediated by the thermal plasma nonlinearities.¹⁶ More recently, a new mechanism has also been studied, consisting in mode–mode coupling mediated by the EP nonlinearity. This mechanism has been found to be responsible, for example, of the excitation of zonal structures (ZSs) by Alfvén modes, and it has been called forced-driven excitation.^17,18

In the present paper, the interaction between AMs and EGAMs is studied showing the results of numerical simulations obtained with the nonlinear, gyrokinetic, electromagnetic, and particle-in-cell (PIC) code ORB5.¹⁹ The dynamics that EGAM $(m,n)=(0,0)$ and dominant AM $(m,n)=(2,1)$ exhibit in numerical simulations where only one toroidal mode number is retained (⁠$n={i} , i=0,1$⁠) is compared with the dynamics observed in simulations where both toroidal mode numbers are kept. The nonlinearities have been maintained only in the EPs dynamics, while the background particle species (deuterium and electrons) follow their unperturbed trajectories. The EPs have a double-bump-on-tail distribution function in velocity space with radial density gradient. The background plasma species have Maxwellian distribution functions.

The paper is structured as follows. In Sec. II, the model of ORB5 is presented. In Sec. III, the equilibrium in use in the simulations is shown, together with the simulations details. In Sec. IV A, we investigate the coupling of AMs and EGAMs in a simplified configuration (circular flux surfaces) to simplify the physics and more easily compare with analytical theory. In Sec. IV B, the analytical interpretation is provided and the theoretical predictions compared with more simulation results. In Sec. IV C, the application to a more realistic configuration (ASDEX Upgrade experimental magnetic equilibrium) is described. Here, the results of numerical simulations obtained with a realistic scenario, the so-called NLED-AUG case,^20,21 are presented. Here, we will show that our simulations have similar features as those presented in the section where a simplified equilibrium was considered. This will lead us to retain that the theory exposed in Sec. IV B can be applied even when a realistic magnetic equilibrium is taken into account. This section represents an extension of the studies detailed in Refs. 22–24. There, the dynamics, respectively, of EGAMs and AMs has been individually investigated, retaining only one toroidal mode number. In this work, the main novelty is that, in order to study the interaction of the EGAM, with the dominant Alfvén mode (AM), $(m,n)=(2,1)$ (which is identified as an EPM), both the toroidal modes $n={0,1}$ are retained in the performed simulations.

ORB5¹⁹ is a nonlinear, global, electromagnetic, particle-in-cell (PIC) code, which solves the gyrokinetic Vlasov–Maxwell system of equations,^25,26 accounting for the presence of collisions and sources. The code uses a system of straight field-line coordinates, $(s,θ∗,φ)$⁠. The poloidal flux ψ, normalized at its value at the edge ψ₀, plays the role of radial coordinate (⁠$s=ψ/ψ0, 0≤s≤1$⁠). $φ$ is the toroidal angle, while the poloidal magnetic $θ*$ angle is defined as

being q(s) the safety factor profile, $θ′$ the geometric poloidal angle, and B the background magnetic field, linked to the magnetic potential $A0$ through the equation $B=∇×A0$⁠.

In ORB5, all the physical quantities are normalized to four reference parameters: the mass and charge of the main ion species (m_i and $qi=eZi$⁠, being e the elementary charge and Z_i the atomic number of the ith ion species), the values of the electron temperature at a radial location s₀ and of the magnetic field on-axis (⁠$Te(s0)$ and B₀, respectively). The derived units are obtained from these four parameters. For example, the time is given in units of the inverse of the ion cyclotron frequency $ωci=qiB0/(mic)$ (being c the speed of light in vacuum), and the velocities are multiple of the ion sound velocity $cs=Te(s0)/mi$⁠, and so on.

In ORB5, it is possible to consider both analytical equilibrium, comprising circular magnetic surfaces, and ideal-MHD equilibria. This is a solution of the Grad–Shafranov equation, calculated through the CHEASE code.²⁷ 

The species distribution function f_s is divided into a prescribed time-independent background distribution function $F0,s$ and into a time-dependent part $δfs$⁠. The subscript s refers to the particle species (i.e., $s=i,e,f$⁠, respectively, background ions, electrons, and fast particles). The time-dependent part of the distribution function is sampled with numerical particles, called markers, representing a portion of the phase space. The gyrokinetic Vlasov equation for the perturbed distribution function is as follows:

for the bulk ions, $F0(R,E)$ only, while for a bump-on-tail $F0(R,E,v∥)$⁠. Details can be found in Ref. 28. In Eq. (2), R is the gyrocenter trajectories, $v∥$ and $v⊥$ are the parallel, and perpendicular velocities respect to the background magnetic field. The equations of motions in mixed-variable formulation of the gyrocenter characteristics $(R,E,v∥,μ)$ are as follows:²⁹ 

where the terms proportional to $ϵδ$ are the nonlinear terms, corresponding to the perturbed equations of motion and $⟨…⟩$ is the gyroaverage. Equations (3)–(5) are formulated in mixed-variables formulation.³⁰ Here, the perturbed magnetic potential $A∥$ has been split into its symplectic and Hamiltonian parts: $A∥=A∥s+A∥h$⁠. In Eqs. (3)–(5), the modified vector potential is present $A*$⁠, the modified parallel magnetic field $B∥*$ and the unit vector in the direction of the magnetic field, are present

The equations of motion are coupled with the following field equations: the quasineutrality condition, the parallel Ampère's law, and the ideal Ohm's law

where

The integrals in Eqs. (8) and (9) represent, respectively, the mixed-variable gyrocenter density and the mixed-variable gyrocenter current. These integrals are calculated in the phase-space volume $dW=B∥*dv∥dμdα$ (being α the gyrophase). $ρ$ is the particle gyroradius, while $ρs=msTs/(qsB)$ is the thermal gyroradius. The equations are solved through the mixed-variable pullback algorithm, presented in Ref. 30, which has been able to mitigate the so-called cancelation problem.³¹ The typical modes of interest are mainly aligned with the magnetic field line $m≃nq(s)$⁠, so a filter is applied to the Fourier coefficients of the perturbed density and current.³² In this way, all the nonphysical modes introduced by charge and current deposition are filtered out. For each toroidal mode, $n∈[nmin,nmax]$⁠, and only for the poloidal modes, $m∈[−nq(s)±Δm]$⁠, where $Δm$ is the width of the retained poloidal modes.

In the present paper, we will study the nonlinear interaction of modes, analyzing simulations where only the EPs full dynamic is retained. This means that only in the EPs equations of motion, the nonlinear terms will be present, that is, the terms proportional to $ϵδ$ in Eqs. (3) and (4). Electrons finite Larmor radius effects are neglected, given $ρe→0$⁠.

The discharge No. 31213@0.84s has been chosen as a base case for linear and nonlinear EPs simulations. The uniqueness of this scenario, the so-called NLED-AUG case,²⁰ is due to the fact that it exhibits a neutral beam (NB)-induced fast-ion β comparable to that of the background plasma. In addition, the fast ions have energy 100 times larger than the thermal background. This unexplored corner of plasma parameters has been chosen to match the realistic ratios of plasma parameters that are going to be met in future fusion machines and to obtain a scenario where the transport of fast particles and the induced mode dynamics can be mainly attributed to the presence of the EPs being the effects of the background plasma minimized.²¹ This scenario is rich of nonlinear physics. A toroidal Alfvén eigenmode (TAE) burst is observed to trigger EGAMs, suggesting the coupling of these modes via the velocity space (EPs avalanches) and via mode–mode coupling processes. The great variety of nonlinear physics present here, together with the fact that the mode dynamics is mainly mediated by the EPs, makes of this an important scenario for the validation of theoretical tools and codes.

The safety factor profile has a reversed shear, with a minimum located at the radial position $s≃0.5$⁠. In Fig. 1, the q-profile of the NLED-AUG case is shown (blue curve). The background plasma temperature profile and the electron density profile of the NLED-AUG case are shown in Fig. 2.

In Table I, the values of some important constants used in the simulations are shown.

In Secs. IV A and IV C, we will show the results of simulations obtained with the radial electron density profile and temperature profile of the background species of the NLED-AUG case. The main ion species and the EPs are deuterium plasma. The equilibrium quasi neutrality is fulfilled by keeping constant the radial electron density profile and varying the EPs concentration, together with the deuterium concentration, satisfying: $ne=∑iZini$ (being n the density profile of the ith species). The bulk plasma species (electrons and deuterium) have Maxwellian distribution functions, while the EPs have a double bump-on-tail distribution function (as in Refs. 22, 23, 28, and 33) because an anisotropy in velocity is needed to drive unstable EGAMs¹⁴ 

The local maximum of the distribution function, in the velocity space, is located at $v∥=v∥,0$⁠, while tval is the width in velocity space of the two shifted Maxwellians.

In Sec. IV A, a circular magnetic equilibrium will be considered and will be characterized by a safety factor profile having a radial dependence close to that of the NLED-AUG case (see Fig. 1, orange dotted line). The EPs will have an on-axis radial density profile (see Fig. 2 on the right, the orange dotted curve). Through this approximation, we have that the EPs provide a contribution with constant sign to the growth rate of the driven modes, through the spatial derivative of their distribution function,^8,24,34 being

In Eq. (13), the poloidal flux ψ plays the role of radial coordinate, and f_i and q_i are, respectively, the distribution function and charge of the ith species of the plasma.

In Sec. IV C, we will adopt the realistic ASDEX Upgrade magnetic equilibrium. The EPs will have an off-axis radial density profile (see Fig. 2 on the right, the green curve), as modeled by TRANSP.³⁵ 

In this section, we focus on the basic physics of the interaction between AMs and zonal structures (ZSs) indicating, with this term, axisymmetric perturbations, in general, can be zero frequency zonal flows (ZFZFs), geodesic acoustic modes (GAM), or energetic particle-driven geodesic acoustic modes (EGAM). To this aim, we consider here a simplified configuration, where the flux surfaces are circular and concentric. This allows us to neglect the secondary correction due to the geometry.

For computational reasons, in Sec. IV A the electron mass has been then taken 500 times lighter than the ions mass $me=mi/500$⁠. In Table II, the values of some important parameters used in the simulations are shown. The EPs have an on-axis radial density profile (see Fig. 2 on the right, the orange dotted line).

In Fig. 3, we show the simulation results obtained in two reference regimes: at low (⁠$⟨nEP⟩/⟨ne⟩=0.0379$⁠) and high (⁠$⟨nEP⟩/⟨ne⟩=0.114$⁠) EPs concentration (being $⟨..⟩$ the volume average). For each regime, we present the mode dynamics observed in simulations where only a single toroidal mode was retained, that is, only $n={0}$ or $n={1}$⁠. These correspond, respectively, to the green and red curves in Fig. 3. We compare this evolution with the dynamics observed in simulations where both toroidal modes are present $n={0,1}$ (blue and orange curves in Fig. 3). In this work, the results of nonlinear simulations are discussed. We will refer to the time interval where the mode is observed to grow exponentially before it saturates, as exponential growth phase. For instance, considering the simulation with $n={0}$ at low EPs concentration, its exponential growth phase lies in the temporal range $t∈[10 000;38 000]ωci−1$ (Fig. 3 on the left, green curve). We discuss now the dynamics observed in the two reference regimes of EPs concentration (⁠$⟨nEP⟩/⟨ne⟩=0.0379$ and $⟨nEP⟩/⟨ne⟩=0.114$⁠).

The measured growth rates and frequencies will be provided in units of Alfvén frequencies measured on axis: $ωA0=1R0B04πρm,0$⁠, being $ρm,0$ the value of the background plasma density on axis: $ρm,0=mini+mene$⁠. The frequencies will be numerically calculated by Fourier analyzing the harmonic components of the scalar potential $δϕm,n∈ℂ$⁠. Positive/negative frequencies correspond to clockwise/anticlockwise rotations in the poloidal plane. Note that since ZS are dominant in $(m,n)=(0,0)$⁠, they will have symmetric frequency spectrum.

In the regime at low EPs concentration, we observe the ZS (green curve in Fig. 3 on the left) to be more unstable than the dominant AM (red curve in Fig. 3 on the left). The AM has an exponential growth phase $t∈[25 000;65 000]ωci−1$ that is longer than that of the ZS, $t∈[10 000;38 000]ωci−1$⁠. In Fig. 4, their mode structure is displayed at two chosen times. The first chosen time $t ∼31 000ωci−1$ is inside the exponential growth phase of both simulations. The second (⁠$t ∼38 000ωci−1$⁠) corresponds, for the ZS, to a time taken before the saturation level is reached (compare with the green curve in Fig. 3 on the left). For the AM, it represents a time taken in the middle of its exponential growth phase (compare with the red curve in Fig. 3 on the left). The measured growth rates and frequencies are reported in Table III. The ZS is identified as an EGAM, while the AM (mainly peaked around $s≃0.5$⁠) has a frequency sitting below the SAW continuum branch $(m,n)=(2,1)$ (see Fig. 6 on the left).

In the simulation where both the toroidal modes are present, we observe the EGAM dynamics (blue curve in Fig. 3 on the left) to be practically unaffected by the presence of the AM in the linear phase. On the contrary, the AM dynamics (orange curve in Fig. 3 on the left) appears to be driven by the EGAM. This results in an increase in the AM drive. Also, comparing Fig. 6(left) with Fig. 6(right) we note that in the simulation with both toroidal mode numbers the original AM (⁠$n={1}$ only) is still present, but it is subdominant respect to the new dominant AM, which is now radially localized and has a higher frequency lying now close to the Alfvén continuum (see Fig. 6 on the right)

These values of growth rate and frequency have been measured in the temporal domain $t[ωci−1]∈[28 000;37 200]$ corresponding to the exponential growth phase of the mode dynamics. They have been labeled with the superscript NL (standing for nonlinear), to underline that these values arise from the nonlinear interaction between the AM and the EGAM, as opposed to what happened in the simulations where a single toroidal mode was present. Therefore, the measured growth rates and frequencies are presented without superscript. In Fig. 5, we show the observed mode structure for each toroidal mode number in simulations where both are present $n={0,1}$⁠. The chosen times are the same at which the mode structures in Fig. 4 have been displayed. In Fig. 5 on the left, we observe the EGAM to have a similar mode structure as that shown in Fig. 4 on the left but to be slightly decreased in amplitude. The AM (see Fig. 5 on the right) appears instead to have a different mode structure respect to that shown in Fig. 4 on the right and to have larger amplitude. The modification of the AM frequency and mode structure passing from $n={1}$ to $n={0,1}$ [compare Fig. 4 with Fig. 5 and Fig. 6 (left) with Fig. 6 (right)] proves that the driven AM is not the original one (present in $n={1}$⁠).

In the regime at higher EP concentrations (Fig. 3 on the right), we observe the AM to be more unstable than the EGAM. The measured growth rates and frequencies related to this regime are reported in Table IV. We highlight that, unlike the AMs, the EGAM growth rate does not grow linearly with the EPs concentration (see Refs. 23, 28, and 36).

When the two modes are together, we observe the AM dynamics to be practically unchanged; instead, the mode $(m,n)=(0,0)$⁠, that here we identify as a ZFZF, has a growth rate

The AM pumps the ZS, with a force-driven mechanism, that was analytically derived in Ref. 17. In Ref. 17, a pumping TAE with its complex conjugate was found to be responsible for the drive of the ZFZF, growing with $γZFZFNL=2γAM$⁠.

In Sec. IV B, we extend the analytical calculation of Ref. 18 to show that a general AM can force-drive a ZS and to explain the inverse mechanism, namely, the excitation of an AM by a ZS, which is what is observed in the regime at low EPs concentration.

The theoretical explanation is given in terms of a three-wave coupling of the ZS and AM, mediated by the curvature–pressure coupling term of the EPs. We refer to Ref. 17, where the study of ZFZF driven unstable by TAE was carried out. Therefore, the most unstable mode (a TAE) was shown to trigger the ZFZF through wave–wave coupling and the ZFZF grows with twice the growth rate of the TAE: $δϕZFZF∼e2γTAEt$⁠. Following the derivation there presented, we want to generalize it to a general AM driving a ZS and to investigate the inverse problem, that is, the triggering of an AM driven unstable through wave–wave coupling by an EGAM (being in this case the most unstable mode: $γEGAM>γAM$⁠), beating with the linear AM. To do so, we follow the theoretical derivation exposed in Ref. 17. However, here we start with the vorticity equation for a general (m, n) mode (and not for the ZS), still neglecting the Reynolds and Maxwell stress, which express the nonlinearity coming from the background plasma (in the present paper, only the nonlinearity in the EPs is retained)

In Eq. (16), the first term represents the field line bending (FLB), the second term represents the Inertia Term (IT), and the third the Curvature Coupling Term (CCT). In Eq. (16), e_s is the charge of the sth species and $∂/∂l$ is derivative parallel to the background magnetic field. $δHm,n$ is the nonadiabatic response of the EPs associated with the perturbed scalar potential with polarization (m, n): $δϕm,n$⁠. The field variable $δψm,n=ωm,nδA∥,m,n/(c k∥)$ is adopted, where $δA∥,m,n$ is the perturbed parallel magnetic potential with polarization (m, n), $k∥$ is the wave-vector component parallel to the background magnetic field, and $ωm,n$ is the oscillation frequency of the (m, n) mode. In the present section, we use $⟨…⟩$ to denote the velocity–space integration and $Jk=J0(Γk)$ is the Bessel function of zero-order with argument $Γk=k⊥ρs$⁠, being $k⊥$ the module of the wave-vector component perpendicular to the background magnetic field and ρ_s the gyroradius of the sth species. In Eq. (16), the frequency drift ω_d appears. It satisfies

where θ is the gyroangle, $ωtr=v∥/(q R0)$ is the transit frequency, and $λd,m,n$ defines the coordinate transformation from the drift center to the particle gyrocenter through the following relations:

In Eq. (18), $δHd,n,m$ is the nonadiabatic response of the EPs to the (m, n) mode, in drift-orbit center. The nonadiabatic response of the EPs is calculated through the nonlinear gyrokinetic equation

In Eq. (19), e is the charge of the EPs, the subscript k is the wave-vector corresponding to the (m, n) helicity, and $δL=δϕ−v∥cδA∥$ and $Λk=∑k=k′+k″b·k′×k″$ are the coupling term. The first term on the right-hand side of Eq. (19) represents the linear response, while the second gives the nonlinear response.

The nonadiabatic response implicitly contains all the information concerning the redistribution of the EPs in phase-space. Such redistribution, responsible for the flattening of the distribution function, causes a gradual decrease in the growth rate [see Eq. (13)] with a consequent saturation of the mode. This is the basic process responsible for the mode saturation in simulations where a single toroidal mode number is retained. In the present paper, nonlinear effects are retained only in the EPs dynamics. Through a careful inspection, we have compared the EPs redistribution in phase-space between the simulation with a single toroidal mode with that with both toroidal modes. That is, at low EPs concentration we have compared the EPs redistribution observed in simulations where $n={0}$ and $n={0,1}$ were retained. At high EPs concentration, we have compared the EPs redistribution observed in simulations where $n={1}$ and $n={0,1}$ were retained. This analysis does not provide evidence of a substantial modification of the EPs redistribution, passing from single to multi-n simulations that could be retained responsible of the observed driven dynamics and, in the end, of the modification of the saturation levels. Because of this, we consider the nonlinear term in Eq. (19) (describing the coupling between two different modes) to be maily responsible for the observed driven dynamics. In summary, the nonadiabatic response δH, given in Eq. (19), describes the redistribution of EPs in phase-space. When a strong EGAM is present and is driving an AM, we do not observe a significant modification of the EPs redistribution in phase-space passing from a simulation with $n={0}$ alone, to a simulation with both $n={0,1}$⁠. The interested reader can find details of the EPs redistribution of an EGAM (⁠$n={0}$⁠) for instance in Ref. 38. We proceed, following the scheme proposed in Ref. 17, considering the (m, n) component of the CCT

This can be rewritten as

Using Eq. (17), we can express Eq. (21) as

The nonlinear nonadiabatic response of the EPs to the mode (m, n) is obtained from

On the right-hand side of Eq. (23), the term $δHk″$ is the linear response of the EPs to the mode with wave-vector $k″$⁠. Equation (23) couples two mode contributions, expressed in $δLk′$ and $δHk″$⁠, and the following relations must be satisfied:

Following the arguments presented in Ref. 17, from Eq. (24) we find that a strong AM (coupling with its complex conjugate) drives a ZFZF,

This result extends that presented in Ref. 17, generalizing to every AM (and not only to the TAE), explaining the interaction described at high EPs concentration in Sec. IV A. In the opposite regime (EGAM, stronger than the AM), the following relation has to be satisfied:

For the case investigated in the present paper, we want to study how an AM (2, 1) is driven unstable by an EGAM, characterized by a scalar potential dominated by its (0, 0) component, but having δH dominated by the mode numbers $(±1,0)$⁠. Note that the mechanism described here is valid, in general, for an AM with any helicity $(mAM,nAM)$⁠. Finally, with a similar calculation as that provided in Ref. 17, we obtain the following growth rate for the force-driven AM:

Here, γ_AM and γ_EGAM are the growth rates of the AM and of the EGAM measured in the exponential growth phase in single-toroidal mode simulations. The scalar potential on the left-hand side of Eq. (27) refers to the nonlinearly generated AM, so

This explains how a marginally unstable AM can be nonlinearly driven with a growth rate similar to that of the EGAM. The presented theoretical estimates are supported by the measurements shown in Fig. 7 in a scan against the EPs that extends the behavior found at the reference case of low EPs concentration, presented in Sec. IV A.

In the present section, we show that the results obtained Secs. IV A and IV B sections apply also to a realistic ASDEX Upgrade magnetic equilibrium, AUG shot No. 31213@0.84s. The plasma parameters and the radial temperature profiles of the background species (electrons and deuterium) together with the radial electron density profile are the same used Secs. IV A and IV B. They are described in Sec. III. Here, the q-profile has the shape shown in Fig. 1 by the blue line.

This shot belongs to a series of discharges performed in 2017 in ASDEX Upgrade.²¹ The plasma scenario has been obtained through plasma parameters previously unexplored in ASDEX Upgrade: $βEP/βth≃0.2−1$ and $EEP/Tth≃150$ (the subscripts th refers to the thermal plasma). Through these, a scenario where the mode dynamics is dominated by the EPs has been obtained and the turbulence contribution minimized. In these discharges, the interaction between EGAMs and AMs has been observed. In order to choose the proper discharge, the beam injection angle dependence of the EGAM has been investigated. The discharge $#31213$ has been selected to maximize the nonlinear EPs dynamics. This case, named the NLED-AUG, has become the base case for linear and nonlinear EP simulations within three European theory projects.²⁰ The NLED-AUG case is obtained with early off-axis NB heating, injected with an angle between the horizontal axis and the beam-line of 7.13°. As it is shown in the portion of spectrogram in Fig. 8, an EPM-TAE burst is observed to trigger the EGAM. The experimental magnetic equilibrium measured at the time $t=0.84 s$ is considered in the simulations of the present section.

EPs have an off-axis radial density profile (see green line in Fig. 2 on the right). A scan against the EPs concentration is presented as $5%≤⟨nEP⟩/⟨nD⟩≤10%$⁠. $⟨..⟩$ indicates the volume average. In this section, the electrons have a realistic mass: $me≃mi/3676$⁠.

In Fig. 9, the time evolution of the scalar potential is presented, for simulations with low and high EPs concentration. The observed mode structures are shown in Figs. 10–13. The dominant AM is an EPM, having the scalar potential dominant in its components $(m,n)=(2,1)$ and being mainly located around the radial position $s≃0.2$ (as in Ref. 24).

In the simulations where only a single toroidal mode is retained, the frequencies of the dominant modes in the exponential growth phase are

corresponding, respectively, to $νZS≃45 kHz$ and $νAM≃96 kHz$⁠. The ZS is identified with an EGAM, and it results to be in good agreement with what was detailed in Ref. 23. The frequency of the AM lies below the branch of the continuum spectrum $(m,n)=(2,1)$⁠, and its frequency is in agreement with the experimental results, lying in the EPM-TAE burst.

When both toroidal mode numbers are considered in the simulations, we observe the same tendency that was detailed in Sec. IV A. For high EPs concentration, we still observe the AM to drive the ZS, with a force-driven mechanism, before the saturation is reached. After that the first saturation of the AM is reached, the mode with the frequency of the EGAM is observed to develop. Below $⟨nEP⟩/⟨nD⟩=0.09$⁠, instead, we observe the inverse process. An EGAM drives the AM, as it was discussed in Sec. IV B.

In the regime below $⟨nEP⟩/⟨nD⟩=0.09$⁠, the difference of the saturation level of the AM with and without ZS is more pronounced (see Fig. 14) with respect to higher concentrations of EPs. We note that in the former regime, the ZSs are more unstable than the AMs, whereas in the latter the AMs are more unstable than the ZSs. The ZSs in the former regime are identified as EGAMs.

The need to understand the nonlinear dynamics of EP-driven instabilities is motivated by the aim to be predictive about future scenarios that will be met in next-generation fusion-relevant machines. In these scenarios, an EP population is present due to external heating mechanisms or as a product of fusion reactions. Most of the previous theoretical works on the nonlinear dynamics of EP-driven instabilities has been focused on the independent evolution of single modes. The main saturation mechanisms considered in most of the previous studies was the EPs redistribution in real space or in phase space.

In Ref. 18, fully nonlinear magnetohydrodynamic simulations showed that the saturation amplitude of a toroidal Alfvén eigenmode (TAE) is reduced by the nonlinear generation of zonal modes. This was analytically explained in Ref. 17, where the formation of a zero frequency zonal structure (ZFZS) was found to be force-driven by a pumping TAE, which couples with its complex conjugate. The induced EPs nonlinearities were found to be dominant over those of the background plasma, which were expressed in the Maxwell and Reynolds stress of the vorticity equation. In both these studies, the ZSs were not excited directly by the EPs.

This paper has been dedicated to the study of the interaction between two different types of modes driven unstable by the EPs: EGAMs and AMs. The results of numerical simulations obtained with the nonlinear, gyrokinetic, electromagnetic, PIC code ORB5 have been presented. In the present work only, the EPs have been allowed to follow their full trajectories, redistributing in the phase space. The nonlinearities arising from the background plasma have been neglected here. The EPs are, in this case, mainly responsible for the saturation and mode interactions. Since the thermal plasma nonlinearities are not considered, the same arguments exposed in Ref. 17 have been used to describe the interplay between the modes. The analytical derivation there exposed has been extended to the case here examined.

We observed that, in numerical simulations explained with analytical theory, through three-wave coupling mediated by the EPs via the curvature pressure term in the vorticity equation, the nonlinear growth rate of the pumped mode is given by the sum of the growth rates of the pumping modes: $γn,mNL=γn′,m′+γn″,m″$⁠. EGAMs are found to channel energy into the AMs in the regime where the EGAM growth rate is higher than the growth rate of the AM. Realistic ASDEX Upgrade plasma profiles have been considered. These results have been found both in analytical and experimental magnetic equilibrium. Despite the approximation of the EPs distribution function, these theoretical findings are qualitatively consistent with the experimental observations. This work opens the path to a new field of research, where the individual modes can saturate not only due to the EPs redistribution but also due to the mutual interaction.

In future works, simulations with more realistic distribution functions will be presented. Also, all plasma nonlinearities will be taken into account, allowing to make comparisons with experimental data of ASDEX Upgrade, where an increased AM activity in the presence of strong EGAMs has also been observed experimentally.

Simulations presented in this work were performed on the CINECA Marconi supercomputer within the ORBFAST and OrbZONE projects.

One of the authors, F. Vannini, thanks Zhyiong Qiu for useful, interesting discussions and for great help provided in understanding the interaction between A.M. and Z.S. and for the help in deriving the exposed analytical theory. Also F. Vannini wants to thank Zhixin Lu and Omar Maj for the help provided in understanding aspects of the gyrokinetic theory. The authors acknowledge stimulating discussions with F. Zonca, I. Novikau, and A. Di Siena. This work was partly performed in the frame of the “Multi-scale Energetic Particle Transport in Fusion Devices” ER project.

This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training program 2014–2018 and 2019–2020 under Grant Agreement No. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

The data that support the findings of this study are available from the corresponding author upon reasonable request.