Gyrokinetic modelling of non-linear interaction of Alfvén waves and EGAMs in ASDEX-Upgrade

Energetic particle (EP) dynamics and excitation of EP driven instabilities is an important topic of study for the physics of fusion reactors. In this paper we consider EPs injected in the plasma by neutral beams at high energies to heat it. EP species exist far from thermal equilibrium in the form of anisotropic non-Maxwellian distribution functions. EP driven modes, such as Alfvén waves (AWs) and EP driven geodesic acoustic modes (EGAMs), can redistribute EPs in phase-space and harm confinement. We examine the effects of experimental-like anisotropic EP distribution functions on the excitation and the non-linear coupling of such instabilities with the gyrokinetic particle-in-cell code ORB5. The growth rate of EGAMs is found to be sensitively dependent on the phase-space shape of the distribution function as well as on the non-linear wave-wave coupling with AWs. Experimental findings are compared with numerical results.


Introduction
Energetic particle (EP) study and confinement is of crucial importance for future fusion reactors.In the so called burning plasma regime, energetic fusion products are supposed to be the dominant heating source for the plasma, provided they can be confined long enough to transfer energy to the bulk species.In present machines, EPs are generated with resonant frequency heating (ICRH) or neutral beam injectors Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
(NBIs).Such particles are in the weakly collisional regime and can exist far from thermal equilibrium.Therefore, EP species have anisotropic non-Maxwellian phase-space distribution functions.Our final goal is being able to predict the EP non-linear dynamics related to experimental-like EP distribution functions.To this purpose the present work is a validation exercise for our model before moving to the realistic simulation of reactor relevant burning plasmas.EP driven modes, such as Alfvén waves (AWs) [1] and EP driven geodesic acoustic modes (EGAMs) [2], are excited by EPs and can redistribute them in phase and real space harming the confinement [3,4].In particular EGAMs (toroidal mode number n = 0) can be linearly excited exclusively by positive gradients of f 0 in velocity space [5], therefore anisotropic distribution functions are needed to excite such instabilities linearly [6][7][8].AWs instead (n ̸ = 0) can be excited both by phase-space and real-space gradients [5], therefore isotropic distribution functions in velocity space as Maxwellians or isotropic slowing down with a radial profile can drive them unstable [9,10].It has been found that the radial position of such modes change with the profile of EP density [3,11].In fact, on-axis profiles trigger Alfvén modes (AMs) in the outer radial domain, while, off-axis peaked EP profiles excite AMs in the core [12,13].In both cases, the AM is located where the EP radial profile has a steep radial gradient.The eigenvalue code LIGKA [14] is useful to evaluate the nature of the AW, by calculating the local kinetic shear Alfvén continuum.LIGKA can thus determine if the mode lies in the continuum Alfvén spectrum and if it is an EP driven mode (EPM) [3,11], or if, lying in the gaps of the Alfvén continuum due to toroidicity, it is an eigenfunction of the system, namely a toroidal Alfvén eigenmode (TAE) [9].It was found, theoretically [15] and numerically [13], that AW and EGAMs can interact non-linearly via wave-wave coupling.The EP non-linear response must be retained in the simulations in order to let these modes interact non-linearly.In this paper, our goal is to enquire into the non-linear driving effect of pump AW on EGAMs [13] using experimental-like EP distribution functions [8] in the experimental relevant scenarios of NLED-AUG case [16,17].For this purpose, the gyrokinetic particle-in-cell (PIC) code ORB5 [18] was run non-linearly with experimental like EP distribution functions in a reference scenario from ASDEX-Upgrade (AUG).
This paper is structured in the following way.In section 2 the theoretical background and the model of ORB5 are presented.Section 3.2 sums up the previous results from non-linear studies with a double bump on tail [10] and the linear simulations with experimental-like distribution functions [8].From section 3.3 we present the non-linear electromagnetic results for both the analytical pitch dependent anisotropic slowing down (ASD) and in section 4 the results with the RABBIT [19] numerical f 0 .The coupling between higher mode number instability and EGAMs is shown to be fundamental for the excitation of the latter.Finally, in section 5 we compare the numerical results with the experimental observations.In section 6 we draw our conclusions and outline the path for future work.

Non-linear coupling between AW and EGAMs
The purpose of this section is offering an analytical explanation of the non-linear interactions between the AWs (n = 1) and the zonal structures (ZSs) (n = 0).TAEs, driven unstable by EPs in the gaps of the shear AW due to toroidicity where damping is minimised, can interact non-linearly with zonal flows (ZFs) via EP non-linear response, in particular with the curvature coupling terms (CCTs) [13,20,21].We briefly summarise here the steps presented in [13,15].We consider δϕ and δA ∥ as our field variables, namely the electrostatic potential and the parallel vector potential component.We take the electrostatic potential to be linear combination of the two different contributions (ZS and TAE): δϕ = δϕ Z + δϕ T .The TAE contribution can be further divided into the pump TAE and its complex conjugate: δϕ T = δϕ 0 + δϕ * 0 .In order to find the ZF behaviour we write the non-linear vorticity equation for ZF [1,3,22]: ( Equation ( 1) is composed by the following terms: the inertia term (IT), the CCT, the Maxwell (MX) and Reynolds (RE) stresses.Here we miss the field line bending term because, being the vorticity equation for ZF, it lacks electromagnetic contributions.In the equation J k = J 0 (k ⊥ ρ) is the zero-order Bessel function accounting for the finite Larmor radius effect (FLR).δH is the non-adiabatic response of the distribution function of the species s.For both δϕ and δH the subscript Z refers to the ZS contribution and k refers to the different waves contributing to the three waves interplay where k is the wave number.In particular k is referred to the ZS, while k ′ and k ′′ are referred to the pumping waves namely the TAE and its complex conjugate.
the operator ⟨. ..⟩ denotes velocity-space integration, for circular cross section large aspect ratio tokamaks , l is the coordinate along the magnetic field line so that ∂ l is the derivative along the field lines and other notation is standard.The three wave theory [23] poses the following constraint: k = k ′ + k ′′ .Since MX and RE contribution to equation ( 1) are negligible, we can solve it by δϕ Z taking the expression of the non adiabatic EP response from the nonlinear gyrokinetic equation [24]: ( In equation ( 2) Z is found to be 0 in the CCT component.Therefore we use the modified varying component of equation (2) and plug it into the expression of the CCT component of equation ( 1): with Ĥ = k θ (k r,0 + k r,0 * ), Â0 | being the amplitude of the usual 'meso'-scale structure, Φ 0 being the fine radial structure and Ĝ ⟩, here the FLR effects have been neglected.
Considering that the contribution of MX and RE stresses in equation ( 1) are negligible, plugging equation (3) into the nonlinear vorticity equation (1) and doing some substitutions we yield: where ). Trivially integrating we find that a pump TAE can drive a ZS unstable with a growth rate double with respect to the TAE's one: ( This result has been obtained in both hybrid [21] and gyrokinetic [10] simulations.Vannini et al demonstrated that the opposite mechanism is also possible [13], that is a strong ZS or EGAM can drive non-linearly an AM unstable.The theory followed by Vannini is parallel to that used by Qiu et al [15] but with the electromagnetic contributions due to the Alfvén activity.
He yielded that the non-linearly driven AM by pump-EGAM has a growth rate that follows the following relation [13]: Hence, so far theory predicted both the cases where a ZS is non-linearly driven by a pump AM in case γ ZS < γ AM [15] and the case where an AM is driven by a pump EGAM if γ ZS > γ AM [13].It is important to notice that the theory shown above was derived without making assumption on the nature of the EP distribution function, because the theory dealt only with the perturbations of the distribution function f 0 .The derivations found for the variational quantities are valid for any starting f 0 .
The main purpose of this paper is to demonstrate that the theories above are also applicable to the case of experimental like distribution functions [8], trying to push the simulations to more realistic scenarios.Finally, this theory does not include estimates of threshold values for the EP density that make the non-linear interactions between AWs and ZSs possible.Furthermore, the complex experimental equilibrium and the very steep q profile (figure 1) make agreement between theory and simulations even more feeble.Therefore, we cannot expect a perfect quantitative match of numerical results with the predicted values of growth rates or with experimental observations, we rather aim for a qualitative agreement.Before moving to the results, we give a brief description of the nonlinear gyrokinetic model in ORB5.

Non-linear gyrokinetic model in ORB5
ORB5 [18] is a gyrokinetic, global, non-linear, electromagnetic, PIC code.It solves the Vlasov-Maxwell system of equations [25][26][27] sampling distribution functions with markers.Fields are solved with a finite element method.The details about the code implementation can be found in [18,27], here only the main aspects will be outlined.As in every gyrokinetic code, the ordering of the different terms plays an important role in the implementation of the gyrokinetic theory.In particular, in ORB5 perturbations due to the stationary magnetic field geometry is considered to be of an order of magnitude smaller than the time varying perturbations of the fields.The smallness parameter associated to the geometrical variations is ε B = ρ th /L B , with ρ th thermal Larmor radius and L B = |∇B/B| −1 the characteristic scale of the spatial gradient of the magnetic field.Meanwhile, the smallness parameter associated to temporal perturbations is with E the electric field, v th thermal velocity, k ⊥ the perpendicular wave number and ϕ 1 the electrostatic potential perturbation.As in other many gyrokinetic codes, also in ORB5 ε B = ε 2 δ .Therefore, ORB5 equations will retain up to the first order terms, O(ε B ), the geometric contributions to the spatial perturbations and up to the second order, O(ε 2 δ ), for the temporal perturbations.The distribution function of each species f s is divided into a timeindependent distribution function F 0,s and a time dependent part δf s .Considering df/dt = 0 and ∂F 0,s /∂t = 0, the kinetic Vlasov equation is solved for the time depending δf s using the characteristics of the particles: where R are the gyrocenters' coordinates, ε = v 2 ∥ /2 + µB is the particle energy and µ = v 2 ⊥ /(2B) is the magnetic momentum.Equation (7) is valid in the general case for non-Maxwellian distribution functions, the isotropic Maxwellian case being F 0.s = F 0.s (R, ε) as for the bulk ions.The gyrokinetic characteristics of the particles can be written as summation of different orders contributions: Ṙ = Ṙ(0) + ε δ Ṙ(1) , v∥ = v∥ (0) + ε δ v∥ (1) , ε = ε(0) + ε δ ε(1) , μ = 0.The 0 th -order derivatives are the so called unperturbed trajectories: whereas, the 1st-order derivatives are the perturbed components due to electromagnetic time perturbations: where ) b.In the linear limit, when solving equation ( 7) the characteristic will just retain the 0 th order components neglecting the terms proportional to ε δ .In the non-linear simulations presented here, the non-linear wave-particle interactions are retained for the EP species, allowing the particles to follow the perturbed trajectories, whereas we let the thermal species follow the unperturbed trajectories.Poisson and Ampére equations are derived from a variational formulation of the problem [25] and are coupled to the systems of equation previously shown: where n s = ´F0,s dW, with dW = B * ∥ dv ∥ dµdα being the phase space volume, β s = µ 0 nsTs B 2 0 is the fluid to magnetic pressure ratio, ρ s = √ m s T s /(q s B) is the thermal gyroradius and the symplectic part of the magnetic potential A s ∥ is found with the ideal Ohm's law: The modes of our interest are mainly aligned along the magnetic surfaces where m ≃ nq(s).Therefore a Fourier filter is applied to retain only the modes within the interval n ∈ [n min , n max ] for the toroidal and m ∈ [−nq(s) ± ∆m] for the poloidal mode number, where ∆m is the filter width for the m mode number.For the simulations in this paper the nonlinearities in the solution of Vlasov equation (7) have been retained just for the EP species, bulk ions and electrons have been treated linearly instead.

NLED-AUG case
Before showing the non-linear numerical results we introduce briefly the experimental relevant configuration used to run the simulations with ORB5.The so called NLED-AUG [16] case is based on the AUG shot #31213 at 0.84 s.Here, EPs at 93 keV were produced through injection in the plasma driven by an off-axis neutral beam injection (NBI) at an angle of 7.13 • with respect to the magnetic axis.The characteristic magnetic equilibrium is shown in figure 1 (left) and the safety factor q in figure 1 (right), the off-axis EP density profile is presented in figure 2. A weakly reversed shear (figure 1 (right)) can be noticed in this equilibrium, as well as a very steep q shear in the outer domain (s ∼ 1).Additionally, the EP pressure is retained in the computation of the equilibrium.Finally, NLED-AUG case was designed for the study of EP effects, relevant in reactor conditions.To this purpose the EP pressure is particularly high, reaching the same order of magnitude of thermal ions: β EP ∼ β th,i .This characteristic lets a very rich zoo of nonlinear EP driven instabilities develop during the shots, making it exceptionally interesting for our studies.Table 1 reports some of the main parameters characterising AUG shot #31213.Further details about NLED-AUG case can be found in [17].

Linear results with experimental-like f 0 and non-linear simulations with DBoT
In this section we sum up the results of the work done so far in the direction of achieving simulations able to predict experimental-relevant scenarios for NLED-AUG case [16] with realistic anisotropic distribution functions of EP.In [8] new realistic distribution functions, namely one analytical and one numerical, were used to study anisotropy effects of the distribution function of EPs on the linear excitation of EGAMs.Firstly, an ASD distribution function was built in order to analytically represent anisotropy through two parameters: ξ 0 and σ ξ (equation ( 3) in [8]).These parameters are representative of the preferred pitch angle of the distribution function (ξ 0 ) and of the distribution of the particles around this preferred angle (σ ξ ).The same distribution function will be used in section 3.3.A scan, varying the two parameters in the ranges ξ 0 ∈ [−1, 0] and σ ξ ∈ [0.1, 0.6], was run using the configuration of the experimental set of shot #31213 in NLED-AUG case.The results show a clear dependence of the growth rate on the two parameters, proving that the phase space shape is essential to the drive of EGAMs (figure 6 in [8]).
In the last part of [8] a numerical experimental-like distribution function was obtained for shots #31213-6 in AUG by the Fokker-Planck heating solver RABBIT [28].These distribution functions were fed as input into linear electrostatic simulations in ORB5.The results show that there is no positive growth rate for any of the four different f 0 in NLED-AUG case (figure 18 in [8]).This result conflicts with the experimental measurements (figure 16) that clearly show that, even if with different growth rates, EGAMs are excited in all the four shots of NLED-AUG case.The incomplete physics of the linear electrostatic model was advanced as an explanation of this discrepancy.In fact, considering the marginal stability of the modes with f 0 with very high pitch angles, the non-linear interactions with n = 1 Alfvén Waves could drive EGAMs unstable.The very purpose of this paper is to prove this hypothesis through non-linear electromagnetic simulations retaining n = 1 electromagnetic modes with experimental-like f 0 s.
In [13], Vannini et al examined the non-linear interactions of AWs and ZFZSs in NLED-AUG case using a double bumpon-tail f 0 at different EP concentrations, namely ⟨n EP ⟩/⟨n e ⟩ = 0.053 in the low and 0.103 in the high concentration cases.In this work the theory advanced in [15] was proved for high EP concentration, namely when the growth rate of AWs is higher than the growth rate of EGAMs (figure 9 in [13], right).The opposite mechanism was also hypothesised and demonstrated with non-linear simulations, that is EGAMs pumping TAEs in  case the growth rate of the former being the highest at lower EP concentration (figure 9 in [13], left).In this paper we want to show that non-linear interactions can drive n = 0 modes unstable as predicted in [15] and shown in [13], even if they are linearly stable with ASD f 0 and numerical experimentallike EP f 0 [8].Furthermore, in [13] the reciprocal pumping mechanism between AWs and EGAMs was proved also in a simplified circular geometry.The results yielded are similar to those obtained in the case of the experimental geometry in AUG.In the next sections we will show the results yielded from non-linear electromagnetic simulations retaining either the mode n = {0} or n = {1} or both n = {0, 1}.In section 3.3 we present the results obtained using the ASD distribution function.Scans in the preferred pitch ξ 0 and the standard deviation σ ξ have been run in order to study the effects of the anisotropy in the non-linear case.In section 4 we perform the same non-linear analysis using this time the experimental-like numerical distribution function from RABBIT.
The purpose of the present work is combining the studies summarised in this section.We aim to evaluate the most experimental-like scenarios, running multi-mode non-linear electromagnetic simulations with experimental-like, anisotropic, pitch dependent distribution functions for EPs.

Analytical slowing down f 0 results
We are interested in analysing the effects of anisotropy of realistic distribution functions on the non-linear excitation of n = 0 modes via wave-wave coupling with AWs.To this purpose we ran non-linear simulations using the analytical pitchdependent ASD f 0 for the EPs, defined as in [8].The simulations were run using the equilibrium and the profiles of NLED-AUG case, shot #31213 in AUG, presented in section 3.1.We ran scans varying the pitch ξ 0 and the scattering of the particles σ ξ around the pitch angle.These parameters varied between their extremes.ξ 0 ranges between 0, in this case EPs' preferred pitch angle is perpendicular to the magnetic field, and −1, EPs' preferred pitch angle is anti-parallel to the magnetic field.The standard deviation σ ξ ranges between 0.1, which corresponds to high anisotropy due to high gradients in v ∥ , and 0.6, which moves closer to isotropy.Such distribution function and its two parameters are represented in [8], figure 1, the scan of growth rates with varying ξ 0 and σ ξ is displayed for the linear case in [8] figure 6.The simulations were run retaining both n = {0, 1} modes and letting them interact through non-linear EP response.Other species (ions and electrons) distribution functions were computed along unperturbed trajectories.In table 2 some parameters used in the simulations are reported.
3.3.1.ASD result, ξ 0 = −0.9 and σ ξ = 0.2.We picked some reference cases for which we ran the simulations only retaining either n = {0} or n = {1} as well as both n = {0, 1}.In this way we are able to compare the growth rates of the different modes with and without the non-linear interactions between them.In particular we analyse the case where the EPs f 0 was set with ξ 0 = −0.9 and σ ξ = 0.2.It is important to remind that this is one of the closest analytical cases to the 2.48 1.672×10 19  1.6018×10 19  6.98×10 17 experimental distribution function reconstructed by RABBIT (section 2.3 in [8]), therefore the results have some experimental relevance.The results are shown in figure 3, the plot shows the radial peak of the scalar potential signal in time of the dominant m mode of either the n = 0 mode (green and blue lines) or the n = 1 mode (red and orange lines) from the simulations retaining either n = {0}, n = {1} or n = {0, 1} modes.Furthermore, we also plot the (m, n) = (1, 0) sideband (SB) for the simulation retaining both n = {0, 1} modes (dashed cyan line) as a marker of the GAM, which is characterised by this m = 1 side band [29].The results (figure 3) show how the growth rate of the ZS is greatly affected by the non-linear interactions with AWs as predicted and shown in [13,15].In fact, the n = {0} mode, run without retaining the AWs (green line in figure 3), shows a marginally stable behaviour, just like what found in the linear case [8].If we keep both modes and let them interact non-linearly, in the growing stage the ZS (blue line) is now pumped through EP nonlinear response by AW (orange line), with a growth rate γ ZS ∼ 2γ AM as prescribed by equation (5).In particular, in this case γ ZS = 5.8 • 10 −2 ω A = 1.75γAM ≃ 2γ AM , considering the time interval t = [7300, 9300]ω −1 ci , with γ AM = 3.31 • 10 −2 ω A .The growth rate of the AM seems not to be much affected by the non-linear interaction with the ZS: γ AM,n={1} ≃ γ AM,n={0,1} .The frequency of damped GAM in the n = {0} simulation is ω GAM = 0.03ω A .The AW alone has instead ω AW = 0.084ω A .The frequencies of the modes in the simulation retaining both ZS and AM are similar: ω ZS,n={0,1} = 0.078 ω A ≃ ω AW = 0.083 ω A .This frequency corresponds to ∼63 kHz, not far from the GAM frequencies measured in NLED-AUG case (green lines near 50 kHz in figure 16).
As we have seen, the damped GAM frequency is much lower than the non-linearly excited ZS, and also of the GAM frequency shown in figure 16.The reason for this it two-fold.Firstly, as we know from basic GAM theory [2], the GAM frequency is proportional to the sound frequency [8,30]: which drops at high radius because of the ion temperature profile (see figure 2).Furthermore, GAM damping rate is proportional to ∝ exp(−q 2 ) [30].Therefore, the GAM will be quickly damped in the core and will survive only where the safety factor q is high: close to the edge.This is why we see the damped GAM localised at s ∼ 0.7 in figure 7.But because of the ω GAM (r) profile, for this high values of radius the frequency will be much lower than what we have with an EGAM excited near the core (as in the blue line of figure 3).Secondly, the fact that the non-linearly excited ZS has a frequency higher than the linear GAM one, highlights that the linear and nonlinear driving mechanisms are inherently different.The ZS in the blue line is not excited linearly by the EP population, it arises from non-linear interactions with the AW.The saturation levels for the n = {1} and the n = {0, 1} AWs are in ORB5 normalised units SAT AW,n={1} ≃ 6.3 and SAT AW,n={0,1} ≃ 1.15.We see that the non-linear interaction with the n = {0} mode reduces the saturation level.The saturation level for the ZS in the n = {0, 1} case is similar to the AW one in the same simulation: SAT ZS,n={0,1} ≃ 1.5.Finally, we notice that there is a second linear phase for the growing modes in figure 3. It seems this is rather a numerical instability which, anyhow, does not affect the parts we are interested in: the first growing phase and the first saturation phase.
The ZS can contain both contributions from zero frequency zonal flows (ZFZFs) and finite frequency GAMs.from the figure, we notice that in the growing phase, where the non-linear interaction between the AW and the ZS takes place, we see that a finite frequency signal is dominant.We can infer from this that the AW is actually non-linearly driving a GAM rather than a ZFZF.Instead, in the plot on the right in figure 4, the spectrogram has been computed in the deep nonlinear saturation phase (time interval t = [13 000 19 000]ω −1 ci ).In this case, we yield very different results: the dominant component has actually zero frequency, showing that eventually after saturation the ZFZF generated by the redistribution of particle in real-space by the saturated AW is dominant in the n = 0 mode.Furthermore, from [29], we know that a GAM is characterised by the presence of an (m, n) = (1, 0) SB, which is expected to have a growth rate similar to the (m, n) = (0, 0) GAM.This is what we actually observe from the trend of the dashed line in figure 3. Finally, we also notice that as we reach the saturation phase, the blue and cyan lines diverge as a symptom of a larger ZFZF developing, as mentioned above.
It is interesting to notice the radial structures of the two modes interacting together.As in [12,13], the dominant AM mode is a (m, n) = (2, 1) mode, placed in the core at s ≃ 0.25 (figure 5 left), where the spatial gradient of the EP density profile has the highest positive gradient (figure 2).The ZS has the maximum radial gradient ∂ r ϕ, corresponding to the highest radial electric field, at s ≃ 0.2 (figure 5

right).
It is possible to analyse the nature of the AW plotting the LIGKA [14] Alfvén continuum on top of the simulation spectrogram (figure 6).We can see that the mode is located at its radial peak (s ∼ 0.25) at a frequency very close to the Alfvén continuum (red lines in figure 6).This hints that the AM, during the linear phase, is actually an EPM in the BAE (betainduced) gap.Note that the mode does not peak at or close the BAE continumm minimum at s = 0.47, but further inside at s = 0.25.Thus the mode properties are mainly determined by the radial EP gradient, justifying the identification as an EPM.After the saturation is reached, the AM's frequency increases: ω AW,NL = 0.12ω A .This hints that also experimental distribution functions are able to reconstruct phenomena related to frequency chirping in NLED-AUG case [10].This phenomenon will not be analysed thoroughly in this paper and it will be enquired in future work.
If we compare the same time slice from the simulations where the n = {0} and n = {1} modes were run independently from one another, we see that the AW is nearly unaffected by the non-linear coupling with the ZS.The electrostatic potential level, in fact, is almost the same between the AW in the n = {0, 1} simulation (figure 5    (right)).In particular, we can notice two things in picture 7 (right).Firstly, as said above, the (m, n) = (0, 0) GAM survives close to the edge of the domain where the damping is minimum because of the q profile.Furthermore, we notice that a strong m = 1 mode is of a magnitude comparable to the m = 0 GAM.A thorough analysis of this mode goes beyond the aim of this work, however it may be a global Alfvén eigenmode which is marginally stable.We do not exclude the possibility that it is also a marginally stable numerical instability.In both cases, this mode does not affect the physical results highlighted above.

ξ 0 and σ ξ scan
Figure 8 shows the dependence on ξ 0 and σ ξ of the growth rates of n = {0} modes in the growing phase of non-linear simulations, where both n = {0} and n = {1} were retained.The main difference with the linear results where only the n = {0} modes were retained (figure 6 and section 4.1 in [8]) is that in the non-linear simulations all the modes have γ ZS > 0, whereas in the linear case we had some values of γ ZS < 0 [8].This is caused by the non-linear wave-wave coupling (section 2) with the AWs.As we can see, especially in the left part of the plane (high values of ξ 0 ), the trend of ZS growth rates is dominated by the trend of AM growth rates.In fact, such modes are particularly excited by deeply passing EPs.This is the case for ASD with low pitch angles with respect to the magnetic field, namely in the case where the preferred normalised parallel velocity ξ is close to −1.For these simulations, as shown in figure 3, the AM is more unstable than the n = {0} mode.This is due to the fact that EP distributions with pitch close to −1 have higher parallel velocity and therefore they can resonate better with the characteristic Alfvén velocity.For the NBIs in NLED-AUG case (93 keV, subalfvénic) at ξ = −1 we yield: Other values of ξ will lead to lower values of v EP,∥ and therefore to less resonance with AWs.Therefore, pumping AWs excite non-linearly ZS via wave-wave coupling, with γ ZS ∼ 2γ AM .At low values of ξ 0 , despite the AM is weakly excited, ZS are driven unstable by wave-wave coupling.For low values of σ ξ and ξ 0 we yield cases where the EGAM is more unstable than the AMs.These cases correspond to the most unstable conditions for the linear drive of EGAMs [8].Here, the EGAMs conserve their growth rate and it is the AM which is driven more unstable by the EGAM [13].We have to warn the reader that the results for the simulations with values of ξ 0 close to 0 are liable to numerical errors.Because of a singularity for the computation of gradients near v ∥ ∼ 0 in the case of anisotropic asymmetric EP f 0 , the results in such cases are not fully reliable.We could expect the values of γ ZS displayed in the bottom right area of the plot (ξ 0 , σ ξ ∼ 0) in figure 8 to be an upper limit.

Dependency of growth rates and saturation levels on the pitch angle
In this section we analyse and compare the growth rates and the saturation level trends from AMs and the ZSs from non-linear simulations retaining either n = {0}, n = {1} or n = {0, 1} modes from a scan using the ASD EP f 0 where the parameter σ ξ was kept constant at 0.2 and ξ 0 varied between its extremes: −1 and 0. The results display more in particular the phenomena described above.If we consider the growth rate trend for the n = {0} case (blue crosses in figure 9) we see that the results are basically the same as those obtained in section 4.3 of [8] (figures 8 and 9 in [8]).We see that the highest growth rates are found for ξ 0 − 0.4, for values of ξ 0 0; −1 the growth rates are negative.In both these cases the EP are not able to resonate with the characteristic velocity of the GAM (∼ v th,i ) because in one case they are trapped (ξ 0 = 0) or too deeply passing (ξ 0 = −1).The growth rates of AMs for the n = {1} cases (red crosses in figure 9) are always positive.Observing the radial structure we notice that for values of ξ 0 ⩾ −0.4 we obtain an n = {1} mode which is a TAE.Meanwhile values of ξ 0 ⩽ −0.4 yield AMs which are EPM.We see that in general EPMs trigger higher growth rates, where the closer ξ 0 gets to -1 the higher the growth rate.This happens because the closer the particles' pitch is to |ξ| = 1 the more deeply passing the particles will be.This enables them to resonate better with the Alfvén velocity (figure 9): as previously highlighted in equation (18).Considering the simulation with n = {0, 1} we see that the AWs (red dots in figure 9) do not change much their growth rates with respect to the n = {1} simulations (red crosses).For this reason red crosses and red dots are partially superimposed.Instead, the growth rates of ZS is greatly affected, especially for those simulations whose n = {0} growth rate was much less than the AW one.We notice that for ξ 0 ⩾ −0.6 the growth rates of the non-linearly excited ZS reaches the AW growth rate.Simultaneously, for ξ 0 ⩽ −0.6, beyond a certain threshold of the AW growth rate, the γ ZS starts to approach the trend foreseen by Qiu et al in [15].In fact, in figure 9 we see that the blue dots (γ ZS,n={0,1} ) start to approach the dashed black line that represents the Qiu's prediction for non-linearly driven ZS by EPM-dominated pump AWs [15]: It is interesting to look at the saturation levels and how they are modified with non-linear coupling between modes, displayed in figure 10.The saturation level has been found identifying the amplitude of the mode at the point in time where the mode at first abandons the linear behavior.Overall, we note that for any ξ 0 the difference in saturation levels between either n = {0} AW or n = {1} ZS with their n = {0, 1} counterpart is always within a factor of 10.Therefore, we see that non-linear interactions can play an important role in the change of the saturation level, and typically they reduce the AM SAT level.Of course, the non-linearly excited ZS have a saturation level different from the null value of the damped modes.The trend observed is rather related to the AWs, and consequently to the non-linearly excited ZS.We notice that for values of ξ 0 ⩾ −0.4 the saturation levels are generally more than one order of magnitude larger than those for ξ 0 < −0.4.As mentioned before in this range of ξ 0 values the TAEs are predominant and for these simulations the AWs manage to reach higher levels of saturation.Meanwhile, the simulations with ξ 0 < −0.4,which are dominated by EPMs in the linear phase, reach saturation at lower levels.It is interesting to observe that for these simulations, after the first EPM saturation there is a second growing phase, which has not been enquired thoroughly, whose final saturation level is comparable to that of the simulations with ξ 0 ⩽ −0.4.This phenomenon can be seen for example in figure 3.In conclusion, in this section we studied the different effects of non-linear interactions between AWs and EGAMs on growth rates, saturation, frequency, using a range of different experimental like analytical ASD EP f 0 .Most of these parameters were found to be greatly affected by this non-linear coupling.This hints that for an accurate description of the EP dynamics in scenarios like the NLED-AUG case we need to retain such effects and run such multi-mode, electromagnetic, non-linear simulations.In the next section we run similar simulations using experimental-like numerical EP f 0 from RABBIT, aiming to match the experimental observation from the NLED-AUG case [16].

Non-linear simulations of shot #31213 in AUG
Once we have analysed the non-linear effects in simulations using analytical ASD f 0 for EPs, we can move to the more realistic case using the numerical experimental-like distribution function obtained from the Fokker-Planck solver RABBIT [8,19,28] (figure 11).
We ran electromagnetic, non-linear simulations using again the equilibrium, profiles and parameters of NLED-AUG case #31213 (section 3.1).We introduced the RABBIT f 0 as input distribution function, and we ran the simulations as done previously retaining either only n = {0} or n = {1} or both.The results are plot in figure 12.In this plot we yield results similar to the case shown before for the analytical ASD with ξ 0 = −0.9 and σ ξ = 0.2.As found in the linear case [8], the case with only n = {0} yields a marginally stable GAM (green line in figure 12) with a growth rate γ GAM,n={0} = −3.8• 10 −3 ω A .The AM has instead a positive growth rate γ AM,n={1} = 1.9 • 10 −2 ω A (red line ), such growth rate is conserved almost identically in the case with both modes (orange line in figure 12) γ AM,n={1} = γ AM,n={0,1} .In the simulation with non-linear interactions, the ZS mode (blue line ) is driven unstable via non-linear wave-wave coupling by the growing AM.The growth rate of this mode is γ NL ZS = 2.9 • 10 −2 ω A = 1.52γAM ≃ 2γ AM (in particular in the time interval t ≃ [17 •   ci ) as foreseen by theory [13,15].Also in this case we added to the plot in figure 12(a) dashed cyan line highlighting the (m, n) = (1, 0) SB characteristic of the GAM.In this case this signal is affected by noise from the n = 1 modes.Nevertheless, for the (m, n) = (0, 0) ZS mode (blue line) we still can detect the strongly AW-coupled non-linear driving mechanism described above.As a final remark, we notice that as we get closer to the saturation phase, the cyan line diverges from the blue one, as observed in the case of figure 3.This hints that a ZF may be developing only close to this saturation phase.
The saturation level is almost the same for the AW in both n = {0} and n = {0, 1} simulations: SAT AW,n={1} = SAT AW,n={0,1} = 0.36 in ORB5 normalised units.While the saturation level for the non-linearly driven ZS in the n = {0, 1} simulation is almost the same as the AW one: SAT ZS,n={0,1} = 0.27.These results demonstrate that it is possible to reconstruct the non-linear EP dynamics of an experimental scenario with good qualitative agreement running non-linear electromagnetic simulations with ORB5 using experimental-like distribution functions.In fact, as it can be seen in figure 16, shot #31213 in NLED-AUG case yielded an excited n = {0} mode.Linear numerical simulations with experimental-like EP f 0 were not able to reproduce such result [8].

NLED-AUG case simulations of shots 31213-6
We ran non-linear, electromagnetic simulations retaining either n = {0}, n = {1} or both n = {0, 1} modes using the experimental density and temperature profiles as in [8] and the EP experimental-like f 0 from rabbit for all the 4 NLED-AUG cases, namely shot #31213-6 in AUG [16].These AUG shots present different NBI injection angles with respect to the magnetic axis, ranging from 6.05 • in #31214 to 7.15 • in #31213; details in section 5 and in figures 14 and 15.The results of these simulations are shown in the following plots.In figure 13 the n = {0} ZS modes from the n = {0, 1} simulations of the four shots are displayed.Figures 14 and 15 show respectively the growth rate and the saturation level trends for the four shots, or for the different NBI injection angles, for the ZSs and AWs in the simulations where either n = {0}, n = {1} or n = {0, 1} modes where retained.
In figure 13 we see the ZSs are all excited.In [8], electrostatic linear physics was unable to reconstruct this result.Therefore, it is clear that the non-linear interactions with the pump TAE drive the ZS unstable.In fact, as we can observe in figure 14, the growth rates of non-linear simulations retaining only n = {0} modes are still negative (blue crosses in figure 14).Instead, in simulations retaining both modes, in all the shots we can see that initially GAMs are damped.Only in a second phase, when the AMs overtake in magnitude the n = {0} modes, they start to grow, non-linearly excited by wave-wave coupling as described by Qiu et al in [15].We can notice that the coupling reaches full effect after some time in the coupled-physics phase (i.e. between t = 18500ω −1 ci  and t = 21000ω −1 ci in shot #31213), here the ZSs grow with a growth rate double with respect to the growth rates of the pump-TAE, as it is clearly shown in figure 14.Furthermore, both from plots 13 and 14, we can observe that the ZSs growth rate is maximum for the most on-axis case, shot #31214, where the NBI had an angle of 6.05 • with respect to the magnetic axis.The growth rate decreases as the NBI angle increases till the most off-axis case, namely shot #31213.This is due to the physics of AWs, which are driven by deeply passing EPs.Therefore, the more the NBI angle is close to the magnetic axis the more strongly the AW will be driven.In fact, particles with higher pitches are closer to the Alfvén velocity and can resonate better with AWs (equation ( 18)).Such a result is clearly shown by the trend of γ AM,n={1} (red dots) in figure 14.As we already found out above, in such multi-mode configurations the non-linear coupled dynamics makes the ZS dynamics dominated by the AWs' one.Therefore the growth rate trend of the ZS (blue dots in figure 14) follows the same pattern of the AWs (the dashed black line in the figure represents the estimated growth rate of non-linearly coupled ZS according to Qiu's formula γ NL ZS = 2 • γ AW [15]).It should also be noticed that there is almost no difference in growth rate between the AWs from the n = {1} simulations (red crosses in figure 14) and n = {0, 1} ones (red dots).In fact, in figure 14 red dots and crosses are strongly superimposed.
In figure 15, the saturation levels for all the modes in the different cases are displayed.We see also here a trend similar to that presented above.In fact, the higher saturation levels are reached for the most on-axis cases.Null saturation levels for the damped n = {0} GAM modes are driven unstable in multimode n = {0, 1} simulations up to saturation levels close to those of the AWs (blue dots in figure 15).Given the negative growth rates of GAMs in n = {0} simulations, we were not able to assign saturation levels (no blue crosses in figure 15).Saturation levels for AWs vary little between the n = {1} and the n = {0, 1} simulations, therefore red crosses and dots are partially superimposed.

Comparison with experimental measurements
The aim of this paper is to prove that numerical simulations with the gyrokinetic code ORB5 are capable of reconstructing the non-linear dynamics of EPs in realistic configurations as observed in experimental campaigns in AUG.As shown in section 3.2, linear electrostatic simulations [8] were incapable of reconstructing growth rates as observed in NLED-AUG case.Nevertheless, a dependency of the growth rate on the phase-space shape was demonstrated for realistic geometries [8].In section 3.3 we showed that, using experimental-like analytical EP f 0 , the non-linear coupling between n = {0} and n = {1} modes can excite ZSs, which would otherwise be stable.In section 4 we used experimental-like, numerical EP f 0 s from RABBIT trying to reconstruct experimental findings.In figure 16, the pick-up coil measurements are presented for the four shots of NLED-AUG case: shots #31213 − 6 in AUG [16].From the plots, it is possible to see that in all the shots the growth rates of n = {0} modes are positive.In fact, the green lines in the plots in figure 16 at ∼50 kHz are n = {0} EGAMs (as marked by the yellow arrows in the plot), that, even if with different γ, ω and saturation levels, are clearly unstable.It is interesting to notice that what is being detected here is the m = 1 side bands of the GAM.In fact, being the (n, m) = (0, 0) GAM an electrostatic perturbation, this is undetectable by magnetic coils.This is also another hint that what are we dealing with is GAMs rather than ZFZF, which have no electromagnetic side bands.We notice that in each spectrogram, along with the GAMs at ∼50 kHz, there are AWs at frequencies between 100 kHz and 150 kHz (highlighted by green arrows in figure 16).The reader may notice other signals at lower frequencies (green lines close to ∼10 kHz), these are believed to be MHD instabilities.In fact, mode number analysis reveals that these modes have finite n and are associated with rational surfaces in the plasma, i.e. they are very likely internal kink modes that are not further studied in this paper.We see a dependence of the saturation level of ZS, therefore of the excitation drive of the modes, on the injection angle of the NBI.As described in [8] and above, the four shots in NLED-AUG case differed in the NBI angle with respect to the magnetic axis: as reported in the figure below, shot #31213 has the most off-axis angle: 7.15 • (top-left in figure 16), shot #31214 the most on-axis: 6.05 • (bottom-left plot), shots #31215-6 a mid range angle: 6.65 • (respectively top and bottom-right plots).The difference between #31215/6 are short beam blips (12 ms) with a second on-axis beam in order to obtain more accurate ion temperature measurement.For the modelling time point chosen, these on-axis blips are not expected to modify the EP distribution function.The plots clearly show that the higher the NBI angle the stronger the mode.Such quantitative agreement was not found in our multi-mode, non-linear, electromagnetic simulations, using experimental temperature and density profiles and experimental-like numerical EP f 0 .Indeed, as previously shown in section 4.2, the trend found in the simulations was quite the opposite.In the multi-mode simulations the non-linear coupling with AW was the dominant driving mechanism for ZS, which, as shown above, is strongest for mainly on-axis beams.This causes a mismatch with experimental measurements.Therefore, the linear drive mechanism, which was found to be stronger for more off-axis beams [8], was less prominent in this non-linear simulations.A synthesis of these two driving mechanism agreeing with experimental observations was not found.
A number of aspects may account for this quantitative difference.First of all, RABBIT calculates f 0 as a series of Legendre polynomials which can lead to (unphysical) oscillations in v ∥ /v.This is the case here especially at the ends of the radial domain (s ∼ 0.0 and s ∼ 1.0).To improve the phase-space profiles, the RABBIT distribution function phasespace shape was approximated in the following way: f 0 (ψ < 0.4)/n ≡ f 0 (ψ = 0.4)/n and f 0 (ψ > 0.6)/n ≡ f 0 (ψ = 0.6)/n.Also, RABBIT gives only the flux-surface average distribution function (due to the approximations which make it a fast code).Monte-Carlo codes such as NUBEAM [31] would give a fully 4D (R, z, ε, v ∥ /v) distribution function, while having other drawbacks such as numeric noise and the associated challenge evaluating gradients.Secondly, there may be inaccuracies in the measurements of temperature, density and safety factor profiles.In fact, for example, the temperature profile was measured only in the last shot of NLED-AUG case, #31216, with an NBI blip.The other profiles have been inferred from that and the other plasma parameters.Also, the safety factor q profile is very steep in the outer radial domain s > 0.9, this implies that the representation of its gradients in this area may be faulty and lead to inaccuracies in the reconstruction of the modes.Finally, non-linear simulations are very sensitive to ∂f 0 /∂v ∥ gradients.Such sensitiveness coupled with the inaccuracies of RABBIT representations of experimentallike f 0 in v ∥ may easily account for the quantitative mismatch of the simulations and the experimental observations.We are confident that with a more accurate description of the experiment, with more precise profiles and better representation of the experimental f 0 we can increase our quantitative agreement with the experiments.Hopefully, such accuracy will be achieved in future studies using the numerical tools developed so far in the present work and in the context of the gyrokinetic code ORB5 [32].

Conclusion and outlook
The work in this paper represents an important step toward the capability of simulating a realistic scenario for a burning plasma fusion device.In particular in this paper it has been shown that through non-linear electromagnetic multimode simulations it is possible to reconstruct the dynamics of anisotropic EP distribution functions in a way consistent with experimental observations.In particular experimentallike EP distribution functions were adopted in ORB5 for the first time to simulate the experimental setup of NLED-AUG case.In this work non-linear electromagnetic simulations were run with the gyrokinetic code ORB5, aiming to reconstruct the experimental, coupled dynamics of GAMs and AWs, which linear studies [8] failed to reproduce.In particular, non-linear simulations retaining both n = {0, 1} modes, using the analytical ASD f 0 , showed that it is possible to excite a ZS via non-linear wave-wave coupling between AMs and ZSs also if the linear drive of the EPs on the n = {0} is not enough to drive them unstable (section 3.3).Good agreement with theory and previous simulations [13,15] was found in spite of the realistic geometry, scenario and safety factor profile.In particular, for the whole ξ 0 and σ ξ scan of the ASD f 0 it was found that the growth rates of n = {0} modes were always positive thanks to non-linear coupling with the AWs.As it was shown in section 3.5, we have found that beyond a certain threshold for γ AW,n={0,1} , i.e. for ξ 0 < −0.5 at σ ξ = 0.2, the γ ZS,n={0,1} follow the theoretical expectations from Qiu et al [15]: γ ZS,n={0,1} = 2 • γ AW,n={0,1} .Therefore, since the growth rate of AW has been found to be higher for pitch angles closer to the magnetic axis (ξ 0 −→ −1), the trend for non-linearly driven ZSs is similar.Depending on the preferred pitch angle of the particles, the AW can be a TAE (ξ 0 ⩾ −0.4 at σ ξ = 0.2), or an EPM (ξ 0 < −0.4).As seen in section 3.5, EPMs have saturation levels lower than TAEs' ones.Anyhow, such simulations were found to go through a second growing phase after the first saturation and eventually reaching levels comparable with those of TAE simulations.
We yield important results also for the case of experimental-like numerical distribution functions, namely those obtained from the heating solver code RABBIT.In section 4 we found that the non-linear drive of AWs can excite ZS unstable in those very cases that were found to be stable, considering only a linear drive from EP.In particular, in all the shots of NLED-AUG case, namely shots #31213-6 in AUG, ORB5 simulations reconstructed unstable dynamics for n = {0} ZSs as observed in the experimental findings.Growth rate measurements of the different modes from the simulations show also for these experimental-like f 0 a very good agreement with Qiu's theory (figure 14).Though, the simulations failed to quantitatively reproduce the same growth rate trend observed in pick-up coil measurements (figure 16).Indeed, the simulations show a higher growth rate for lower NBI angle, hinting that the dominant driving effect comes from the non-linear coupling with AMs.Meanwhile, the magnetic measurements show an opposite trend: the most unstable cases correspond to the most off-axis NBI angles, showing that somehow the linear EGAM drive [8] must have a quantitative effect on the growth rate too.Some explanations for such quantitative disagreement are advanced in section 5.
In conclusion, through the new numerical tools implemented in ORB5 we managed to get a very good qualitative agreement between the experimental scenario simulations and the experimental findings.In spite of some quantitative disagreements, due to numerical defects in the representations of profiles and distribution functions, we are confident that this work and those carried out recently in the ORB5 context, represents a very important step toward the prediction of non-linear EP dynamics in experimentally relevant scenarios [33].Future studies with these new numerical tools will hopefully be able to match quantitatively experimental measurements and observations, enabling us to predict EP dynamics and transport [34] also in realistic, and eventually burning plasma, scenarios.

Figure 3 .
Figure 3. Modification of ZS growth rate in presence of AMs, the plot shows the radial peak of the amplitude of the dominant m mode for each of the n = {0}, {1}, {0, 1} modes in time.For the simulation retaining both n = {0, 1} modes the (m, n) = (1, 0) side band (SB) was also plotted.These simulations used the ASD f 0 with ξ 0 = −0.9 and σ ξ = 0.2.
It is not trivial to separate the two in the signal by frequency alone.Nevertheless, from the n = 0 frequency spectra we can indirectly infer which component is dominant on the other.To do so we can look at figure 4. The spectrograms of the n = 0 ZS in the n = {0, 1} simulations have been shown for two time intervals.The one on the left has been performed in the growing stage (time interval t = [7500 11 000]ω −1 ci ).As we can see

Figure 4 .
Figure 4. Frequency spectrogram of the ZS in the linear t = [7500 11 000]ω −1 ci (left) and non-linear t = [13 000 19 000]ω −1 ci (right) phases of the simulation with both n = {0, 1} modes.We can notice the two different modes combining the ZS.In the growing phase there is a dominant GAM at finite frequency, in the non-linear phase the ZFZF becomes dominant.

Figure 6 .
Figure 6.Frequency spectrogram of the AW in the linear (left) and non-linear (right) phase of the simulation with n = {0, 1} modes, LIGKA Alfvén continuum [14] for n = {1} mode in AUG shot #31213 is superimposed with red lines to the spectrogram.

Figure 12 .
Figure 12.Modification of ZS growth rate in presence of AMs, the plot shows the radial peak of the amplitude of the dominant m mode for each of the n = {0}, {1}, {0, 1} modes in time.For the simulation retaining both n = {0, 1} modes the (m, n) = (1, 0) side band (SB) was also plotted.These simulations used the RABBIT f 0 from AUG shot #31213.
Taking the nonlinear EP response in terms of the drift orbit centre coordinates δH NL Z = e iλ dZ δH NL dZ , with λ dZ = λdZ cosθ = k Z ρd cosθ.With this notation it is possible to rewrite equation (2), we can split it into the surface averaged and the poloidally varying components.Being ω d = ω tr ∂ θ λ dZ the surface averaged component of δH NL