Gyrokinetic investigation of the damping channels of Alfv\'en modes in ASDEX Upgrade

The linear destabilization and nonlinear saturation of energetic-particle driven Alfv\'enic instabilities in tokamaks strongly depend on the damping channels. In this work, the collisionless damping mechanisms of Alfv\'enic modes are investigated within a gyrokinetic framework, by means of global simulations with the particle-in-cell code ORB5, and compared with the eigenvalue code LIGKA and reduced models. In particular, the continuum damping and the Landau damping (of ions and electrons) are considered. The electron Landau damping is found to be dominant on the ion Landau damping for experimentally relevant cases. As an application, the linear and nonlinear dynamics of toroidicity induced Alfv\'en eigenmodes and energetic-particle driven modes in ASDEX Upgrade is investigated theoretically and compared with experimental measurements.


Introduction
In burning plasmas relevant for magnetic fusion energy (MFE) research, an important role is played by energetic particles (EPs). With the term EPs, we refer to fusion reaction products (like alpha particles) or super-thermal ions or electrons resulting from plasma heating. Such particles possess higher velocities compared to those typical of the background plasma. In typical tokamaks the frequency associated with the gyrocenter motion of the EPs falls inside the magnetohydrodynamic (MHD) domain (O(10 −2 Ω ci )), being Ω ci the ion cyclotron frequency. Because of that, through resonant wave-particle interactions, one of the three solution satisfying the MHD dispersion relation can be excited. Among them, the most detrimental are the shear Alfvén waves (SAWs), having a group velocity parallel to the equilibrium magnetic field and satisfying the dispersion relation: ω = k || v A , where the Alfvén speed has been defined v A = B/(4 π ρ m0 ) (being ρ m,0 the background plasma mass density and B the background magnetic field strength). The excitation of these modes creates a transport channel for the EPs, which can lead to loss of EPs before their thermalization, causing a less effective heating and also possibly damaging the vessel of the machine. they are also believed to be responsible of large abrupt events (ALE) observed in the Japanese tokamak (JT-60U) (see Ref. [1]). We can simply outline the whole zoology of SAWs basically in discrete Alfvén eigenmodes (AEs) and energetic particle continuum modes (EPMs), (non-normal modes of the SAW continuum spectra, merging as discrete fluctuations at the frequency that maximizes the wave-EP power exchange, above the threshold condition of the continuum damping [2]).
In this paper the attention will be mainly focused on toroidal Alfvén eigenmodes (TAE, Alfvén eignemodes lying in the frequency gap caused by the tokamak toroidicity) and EPM. The main goal of this paper will be to understand what are the damping mechanisms of the modes of interest comparing (when possible) the simulations results with the prediction of MHD theory. The proper domain to take into account all the nonlinear effects, like wave-wave and waveparticle interaction, as well as finite-Larmor-radius and finite-orbit-width effects, is represented by the gyrokinetic theory. Because of that, the simulations have been principally carried with the global, nonlinear, electromagnetic, gyrokinetic, PIC code ORB5 [3,4] whose model, if properly set, contains the MHD equations as a subset. Some comparisons with the linear gyrokinetic code LIGKA [5] have also been performed.
The paper is structured as follows. In Section 2 a description of the model implemented in the code ORB5 is given. Section 3 and Section 4 will be dedicated to the description of two different damping mechanisms affecting the SAWs: the continuum damping and the Landau damping. They will be briefly explained analytically, starting from the MHD equations. The theory in use will be compared with the results from the simulations carried with ORB5. In Section3 the numerical simulations will be performed in the cylinder limit using simplified profiles. In the simulations in Section 4, a small but finite inverse aspect ratio will be considered, using the equilibrium profiles of the International Tokamak Physics Activity (ITPA, see Ref. [6]). In Section 5 the studies on the linear and nonlinear growth rate and frequency spectra conducted considering experimental profiles from the NLED-AUG case [7] will be presented. Finally, summary and future outlook will be shown in Section 6.

Model
Since the Alfvén waves have a frequency much smaller then the typical ion cyclotron frequency (Ω ci ) and their amplitude in the core is small compared to the background quantities, a good description of their propagation and interaction with the bulk plasma can be given through the gyrokinetic theory. This allows to retain a kinetic description of the events under consideration, reducing the 6D problem to a 5D one, by averaging the fast gyromotion. In this way the numerical costs are sensibly reduced.
ORB5 is a global, nonliner, gyrokinetic, electromagnetic, PIC code, which can take into account collisions and sources [3,8]. The gyrokinetic model of ORB5 [9] contains the reduced MHD equations as a subset [10]. In this section we give a brief description of the gyrokinetic model implemented in ORB5 and briefly show how the implemented equations are solved. We refer for more exhaustive explanations to Ref. [8], which also give a more complete description of the recent updates in ORB5. Concerning the magnetic equilibrium in use, ORB5 can either upload an ideal-MHD equilibrium (solution of the Grad-Shafranov equation) from the CHEASE code [11], or consider ad-hoc equilibria constituted by circular, concentric magnetic surfaces. It deals with a straight-field line set of coordinates. The magnetic surfaces are labeled by s = ψ/ψ edge , which plays the role of radial coordinate. Here ψ is the poloidal magnetic flux function. The angular dependence is given by the toroidal coordinate ϕ and by the poloidal magnetic angle: being θ the geometrical poloidal angle and q(s) the safety factor, defined as: All the quantities in the code are normalized through four reference parameters: the ion mass (m i ), the ion charge (q i = e Z i , being e the electric charge and Z i the atomic number), the value of the magnetic field strength on axis (B 0 = |B(s = 0)|) and the value of the electron temperature at a specified reference position s 0 , T e (s 0 ). All other normalized quantities are obtained through these: the time units are provided in the inverse of the ion-cyclotron frequency, Ω ci = q i B 0 /(m i c), the velocity units are normalized through the ions sound velocity (c s = q s T e (s 0 )/m i , being the temperature measured in keV ), the length units through the ion sound Larmor radius (ρ s = c s /Ω ci ) and the densities are normalized by means of their average in space. The Vlasov-Maxwell gyrokinetic equations are derived through variational principles from a discrete gyrokinetic Lagrangian. This choice allows to take into account all the simplification needed in the model under consideration, directly in the Lagrangian and then derive the gyrokinetic equation. An immediate consequence is that it is possible to consistently derive conserved quantities (like the energy, see Ref. [12]), that are also used in ORB5 to test the quality of the simulations performed. This choice also allows to derive naturally the weak gyrokinetic form of the field equations. In the Lagrangian an ordering is present, separating the effects given from the geometry of the non-uniform magnetic field, from those related to the fluctuations of the electromagnetic perturbation. This means that (as can be derived, see Ref. [9,13]) the small parameter related to the variation of the background magnetic field B = ρ th L B , (being ρ th the thermal Larmor radius and L B the typical variation length of the magnetic field) and the small parameter related to the fluctuating electromagnetic field ( δ ) are related through: In this way the action functional, written in "p z -formulation" appears to be the following: where α = 0 gives the electrostatic model, while α = 1 the electromagnetic one. In Eq.4 dΩ = dV dW , being dW = B * || dµ dp z . A sum over the species "s" also appears. The symplectic magnetic field is defined through the symplectic magnetic potential A * = A + (c/q s )p zb beingb the unitary vector parallel to the background magnetic field. The canonical gyrocenter momentum is p z = m s v + α δ (q s /c) A 1 . In the action functional, some approximations have been done. The quasi-neutrality allows to consider in Eq.4 only the contribution given from the magnetic potential, neglecting the one given from the perturbed electric field. Also the incompressibility of the parallel perturbed magnetic field is assumed B 1,|| = o(B 1,⊥ ) and only the perpendicular component of the perturbed magnetic potential is retained: In Eq.4 it must be noted that while H 0 ,H 1 multiply the total distribution functions f s , f e , H 2 is related only to the equilibrium distribution function. Thanks to this choice nonlinear second order terms do not appear in the gyrocenter dynamics and the field equations are linear. The gyrocenter hamiltonians appearing are: where the gyroaveraging operator has been introduced f = 1 2π 2π 0 dθf . The gyroaveraging is removed for the electrons that are treated as drift-kinetic: For the distribution of the species s the linear gyrokinetic Vlasov equation is: where the gyrokinetic characteristics can be derived from Eq.4 and are: The field equations, quasineutrality and Ampère, are both derived from Eq.4 via functional derivatives on the perturbed field. ORB5 splits the total distribution function in a background distribution function f 0 and in a time dependent one δf and discretize this latter through numerical particles (markers) used to sample the phase space. Through an operator splitting approach the code solves first the conlisionless dynamics (using a 4th-order Runge-Kutta method) and then treats the collisions with a Lanngevin approach. The quasineutrality and Ampère equations are solved using the Galerkin methods and the perturbed fields are discretized through cubic B-splines finite elements defined on a grid (N s , N χ , N φ ). Finally it is important to mention that from the numerical side, recently the mixed-representation ("pullback" scheme [14]) has solved the so-called "cancellation problem" for electromagnetic simulations.

Continuum damping
In the present section, the continuum damping will be studied. The tokamak configuration will be selected in order to have the continuum damping as main damping mechanism affecting the Alfvén waves. In order to do so, it is first important to understand what are the equations governing the Alfvén waves. These will be obtained under the validity of the ideal magnetohydrodynamic (MHD) theory by treating MHD equations with a perturbative approach. The Alfvén wave's dynamics can be expressed starting from the quasi-neutrality condition ∇ · δJ = 0, (being δJ the perturbed current) that rewritten in terms of its components parallel and perpendicular to the background magnetic field (b = B/B) reads: Following Ref. [15], in order to obtain a simplified but relevant set of equations, modes with k ⊥ k are considered, so that the time scale between incompressible shear Alfvén waves and compressional waves can be separated. To further simplify the problem, we consider a pressureless plasma (P = 0) obtaining the following vorticity equation: A differential equation for the perturbed scalar potential δφ is thus obtained. It is linked to the perturbed magnetic potential (δA ≈ δAb) through the condition δE = 0, derived from the ideal Ohm's law. In this section a non-uniform plasma equilibrium with cylindrical limit, will be considered. a will denote the typical length scale perpendicular to the equilibrium magnetic field while R 0 will represent the typical length scale parallel to it. The equilibrium magnetic field, in a coordinate system (r, θ, z) will be assumed to be B = (0, B 0,θ (r), B 0,z (r)). The geometrical radius r can be used in the cylindrical limit instead of the flux coordinate s. By assuming a shear Alfvén oscillation of the scalar potential δφ(r, θ, φ, t) of the form: where m is the poloidal mode number, we can now write Eq.10 in cylindrical coordinates: where the local safety-factor profile has been defined: The shear Alfvén wave dispersion relation is then found to be: k the local nature of the continuum plasma oscillation can be exploited by reducing Eq.12 to: which integrated in the radial domain, becomes a differential equation for the radial electric field E r : Assuming now a dispersion relation of the form ω A (r) = ω A0 + ω A (r − r 0 ) and by Fourier transforming the radial electric field in the radial coordinate the following relation is obtained: The obtained linear dependence in time of the radial wave number is a proof of the phase mixing together with the fact that, being E r (r, t) = −i k r (t)φ(r, t), the scalar potential exhibits the characteristic decay called continuum damping: as it was proved in Ref. [15] (see also Ref. [16,17,18] for the application to Geodesic Acoustic modes, GAMs). To see evidence of this mechanism, the simulations presented in this section have been run with simplified geometries and profiles, without the presence of EPs. In order to be in the cylinder limit the inverse aspect ratio has been chosen to be = 0.01 and flat density (n e = n i = 2.22 · 10 20 m −3 ) and temperature (T e = T i = 0.01 keV ) profiles have been taken into account. This leaves all the radial dependence of the dispersion relation in the safety factor profile. Moreover, in this temperature regime the Landau damping can be neglected [19]. In the simulations under examination only one axysimmetric perturbation n = 0 , m = 1 peaked at the radial position r = 0.6, has been considered, together with a linear safety-factor profile q = q 0 + q 1 · r so that: The other important parameters in the simulations (minor and major radius, value on axis of the equilibrium magnetic field, ion cyclotron frequency, Alfvén frequency on-axis and the ratio of the last two), are reported in Tab.1.
In Fig.1 on the left, the measured values for the wave numbers k r are shown for a simulation having q 0 = 1.75 and q 1 = 0.5. They have been measured interpolating the mode structure with a sinusoidal function at times where a maximum has been reached at the radial position r = 0.6. By linearly inetrpolating the measured wave numbers it is possible to calculate the coefficient k r,1 , which is found to be in reasonable agreement with the theoretical expectations Eq.17. In Fig.1 the dynamic of the scalar potential at some radial positions is shown, together with the predicted decay, Eq.18. The theoretical prediction for this simulation (q 0 = 1.75 and q 1 = 0.5) is that k r,1 = 0.123 ω A0 , while the measured value is k r,1 = 0.107 ω A0 . Right: Perturbation amplitude dependence on time. Analytical estimation are given by the dashed lines (curves decaying in time as φ ∼ |ω A t| −1 ). The scalar potential measured at different radial positions is given by continuous lines. No EPs are present here.
In Fig.2 finally the obtained values of the coefficients k r,1 have been plotted against different values of the slope of the safety factor profiles (q 1 ) in use in the different simulations and compared with Eq.17. Given the reasonable agreement found between the results of the numerical simulations and the theory, we can say to have verified the relevance of the continuum damping as main damping mechanism for this specific case and to have observed the presence of phase mixing.

Landau damping
In the present section bulk species temperatures high enough to make the continuum damping negligible with respect to the Landau damping will be considered, so that the latter becomes the main damping mechanisms.
The attention will be focused also on a particular Alfvén eigenmode, the toroidal Alfvén eignemode (TAE). Its characteristic frequency lies in the gap created in the continuum spectra by two close poloidal modes (m, m+1) which are coupled because of the finite tokamak toroidicity, [20]. A TAE is located at a radial position r 0 satisfying: q(r 0 ) = 2 m+1 2 n . The theoretical derivation exposed in Ref. [20,21] will be now followed. Here, a kinetic transverse part of the wave-induced current δJ k ⊥ is added to the ideal MHD current, so that Eq.9 becomes: Equation 19 is then multiplied by δφ and integrated in the overall plasma volume obtaining: where it was assumed as boundary conditions dx · δJ δφ = 0. Calling ω 0 the frequency of the wave solution of the ideal MHD vorticity equation, we can consider ω = ω 0 + δω the solution of the new vorticity equation Eq.19, being δω ω 0 . Following a perturbative approach, an expression for γ = Im{ω} is obtained from Eq.20: where all the appearing perturbed quantities have been decomposed in Fourier components in the poloidal plane. In order to obtain a simplified equation for γ, some further calculations have been done and will be now described. Eq.21 is then written in cylindrical coordinates, after writing the perturbed current in terms of the perturbed distribution function and this in terms of the unperturbed distribution function F 0 . Assuming a Maxwellian distribution function F 0 and focusing our attention on TAE (that is assuming to have a perturbation δφ strongly peaked at the radial position where we expect to have a TAE), we obtain: Being: And: In Eq.22, γ has been decomposed in the species contributions (the sum over j). It is formally identical to the one derived in Ref. [22]. The difference lies in the fact that in Ref. [22] the authors have obtained the estimation for γ starting from energy principles, while here everything has been done by adding a correction to the MHD quasi-neutrality equation and thus to the Alfvén dynamics [23]. Since we are interested in the study of the Landau damping, we will not consider the EPs contribution, which actually drives the mode unstable. Eq.22 depends on the ratio between the Alfvén speed and the thermal velocity of the considered species. This means that λ j ∼ m 1/2 j and, since the bulk ion mass is bigger than the electron mass (for Hydrogens m H ∼ 2000 m e ) one can understand that the ion contribution is negligible with respect to the electron contributions (because of the presence of the exponential terms in the polinomia). Because of that we will focus our attention on the electron Landau damping in this analytical derivation.
In this section an equilibrium with small, but finite value of inverse aspect ratio will be considered, = 0.1. The temperature profiles are flat. When considered, the fast particles have a density profile peaked on axis (see Fig.3 on the left). The magnetic equilibrium and profiles are those of the ITPA-TAE international benchmark case [6] and the safety factor profile is shown in Fig.3 on the right. The main results that will be displayed in this section, have been obtained considering heavier electrons: m e = m H /200. This has been checked to be at convergence. In Tab.2 other important details of the simulations are presented.   The chosen initial potential perturbation is peaked around r = 0.5 and is constituted by one single toroidal mode number n = 6 while the poloidal mode numbers 9 ≤ m ≤ 12 are considered (see Fig.4 on the left). A TAE is located at r = 0.5 in the gap of the continuum spectra created by the coupling of the poloidal modes m = 10 and m = 11, as it is shown in Fig.4 on the right. In Fig.5 the dependence of the damping rate against the value of the (flat) electron temperature is shown. The damping rate is found to increase with the increasing electron temperature. This is an evidence that the dominant damping is the electron Landau damping. The errorbar of the measured points correspond to 20% of their value. This because, as it is shown in Fig.6 on the left, the damping rate value has a dependence on the chosen width of the perturbation. For completeness in Fig.5 the approximated analytical electron Landau damping formula is also shown (dashed line). A reasonable qualitative agreement is found between the predicted decay and the simulation results. Finally in Fig.6 on the right, the dependence of the measured damping rate of ORB5 simulations against the electron mass has been shown. For decreasing electron masses, the absolute value of the damping rate is shown to decrease, consistently with theory of the electron Landau damping. In summary, it has been proved that the bulk electrons provide the main damping mechanism of the observed Alfvén modes in this particular regime. Several approximations have been done in the analytical theory, inter alia only passing particles are considered thus neglecting the contributions of barely trapped electrons, which are thought to be important, and which are included in our numerical simulations.

NLED-AUG case
In the present section the results of numerical simulations involving a realistic scenario will be presented.
The shot number #31213 of ASDEX-Upgrade (AUG) has been selected within the Non-Linear Energetic-particle Dynamics (NLED) Eurofusion enabling research project [7]. Here an early offaxis NBI (with T EP ∼ 93 keV ) occurs with an injection angle (angle between the horizontal axis and the beam-line) of 7.13 • . The magnetic equilibrium measured at the time t = 0.84 s is considered in the present simulations (see Fig.7 on the bottom left). This case is referred to as "NLED-AUG case". Further description of this case can be found in Ref. [7]. The NLED-AUG case is found to be of great interest because of its rich linear and nonlinear dynamics arising from the interaction of the modes with the EPs.
Tab.3 contains the details of the main parameters considered in the simulations. Tab.4 shows the values of the bulk species profiles on axis in the absence of EPs. The bulk ions, as well as the EPs (when considered), are constituted by deuterium. The EP temperature will be always considered to be radially flat and equal to T EP = 113 keV . For the EPs density profiles, an offaxis density profile fitting the experimental profiles is considered with Maxwellian distribution function. For comparison we also run simulations with an on-axis EPs density profiles. Note that when EPs are included the electron density profiles is changed in order to match quasi-neutrality n e = Z i · n i + Z EP · n EP .
In Fig.7 the safety factor profile is shown, together with the temperature profiles of the bulk species. The safety factor profile has a reversed shear, with q min (r = 0.5) 2.28. The density profiles in use will be shown in the following subsection. For numerical reasons, the electron mass is chosen to be m e = m D /500, being m D the deuterium mass.
This section is divided in two subsections. In the first, the results of numerical simulations only involving the linear dynamics will be discussed. In particular, the dependence of γ against the electron temperature will be shown, with and without EPs contribution. Also the results of a benchmark with LIGKA are presented. In the second subsection instead, the results of simulations involving also the nonlinear dynamics will be presented. Finally, it is important to remind that, unless specifically written, the bulk and energetic ions will be treated as drift-kinetic. The initial perturbation considered will take into account just one toroidal mode number (n = 1) and the poloidal mode number 0 ≤ m ≤ 7.

Linear simulations
In the present subsection the results of numerical simulations involving only the wave linear dynamics will be presented. Fig.8 and Fig.10 show the frequency spectra, mode structure and poloidal view of the scalar potential φ, obtained considering respectively off-axis and on-axis density profile for the EPs and a concentration equal to 3%. The frequency spectra have been analyzed in the same temporal domain, when a clearly growing mode is observed. In both Fig.8 and Fig.10 the continuum spectra obtained with the linear gyrokinetic code LIGKA [5] is shown (red crosses), together with the analytical curve for the continuum spectra calculated in cylindrical coordinates and including the toroidicity effects, [20] (green dotted line).
When the EPs possess an off-axis profile, a mode sitting at the radial position r 0.22 is observed. The dominant poloidal component of the scalar potential appears to be that having m = 2. Due to the measured frequency lying on the branch of the continuum, this mode is identified as an EPM. The frequency is measured as: When finite Larmor radius are take into account, a slightly change in the frequency (from 129 to 131 kHz) is observed. The frequency measured in numerical simulations can be compared with the experimental measurements. In Fig.9, the spectrogram obtained with Mirnov coils is shown. A big variety of EPs driven modes can be found. At t = 0.84 s, the modes with frequencies around 50 kHz have been identified as EGAMs (see Ref. [24,25,26]). We focus here on the Alfvén modes with frequency lying in the domain between 100 and 150 kHz. The numerical result shows that these modes are indeed EPMs (see the white cross in Fig.9). Despite the approximation of the EPs distribution function we notice that the results of the numerical simulations appear to be in good agreement with the modes observed in the spectrogram.   An on-axis density profile is also considered for the EPs. The radial dependence of the EPs density profiles is expressed by the formula: n EP (1 − r α ) β . The coefficients α, β have been chosen in order to have the second derivative of n EP equal to zero at the position where an Alfvén mode is expected. The numerical analysis shows that a mode lying in the gap of the continuum spectra, created by the poloidal modes m = 2 and m = 3 is observed. It appears to be peaked at the radial position r 0.738. Due to the radial localization and frequency this is identified as a TAE.
In Fig.11 the dependence of γ against the value of the bulk species temperature is shown. In the plot on the top left, the dependence of the growth rate against the electron temperature, keeping the bulk ions temperature constant (T i (r = 0) = 3.5 keV ), is shown. In the plot on the top right, the dependence of the growth rate against the bulk ion temperature, keeping the electron temperature constant (T e (r = 0) = 0.707 keV ), is shown. The EPs temperature is flat T EP = 113 KeV and the EPs have a concentration of 3%. In the plot on the bottom instead, the dependence of the damping rate (simulations without fast particles) is shown, against the value of the electron temperature. This study of the dependence of the growth or damping rate against the electron temperature shows that the electrons are the main responsible of the damping of the Alfvén modes even in this realistic scenario, which is identified here as electron Landau damping. Since a realistic scenario is considered here, the approximate theoretical predictions for the Landau damping described in the previous sections is outside its validity regime and therefore is not shown. Figure 11: Top: respectively on the left and on the right, two scans in the electron and bulk ions temperature for a growing mode are shown. EPs concentration 3%,T EP = 113 KeV . Bottom: Scan in the electron temperature for a damped mode.
In Fig.12 a theoretical estimation of the regions in the radial domain where the Landau damping is supposed to dominate over the continuum damping is presented for simulations without EPs. In order to perform this calculation, the characteristic radial structure has been measured in a simulation with an EPM, giving the value of k r,0 . This allows us to calculate the analytical prediction for the half decay time of the continuum damping. The half decay time due to Landau damping, on the other hands, can be measured in a simulation without EPs: The radial regions where the Landau damping dominates over the continuum damping are those where the half decay time of the Landau damping is smaller than the half decay time of the continuum damping: t 1/2,continuum damping > t 1/2,Landau damping . Equivalently it can be said that the Landau damping dominates over the continumm damping in those regions where: In Fig.12 on the left the continuum spectra calculated with the code LIGKA are shown. They are used to calculate the derivative of the frequency to obtain the estimation of the regions where the Landau damping dominates over the continuum damping. To do so, the frequency has been divided in two branches, denoted as upper and lower branch. In Fig.12 on the right the regions where t 1/2,continuum damping > t 1/2,Landau damping are shown. Therefore the green regions correspond to the radial domain where the Landau damping is dominant over the continuum damping if the frequency of the mode is sitting on the lower branch. The cyan regions instead represent the domains where the Landau damping is dominating on the continuum damping if the frequency of the mode is sitting on the upper branch. Figure 12: Left: continuum spectra calculated with the code LIGKA. It has been divided into two branches (denoted as upper and down). Right: regions where t 1/2,continuum damping > t 1/2,Landau damping . Fig.13, the results of a first benchmark between ORB5 and the code LIGKA are shown. Here a scan in the EPs temperature is depicted. The EPs have an on-axis profiles and their concentration is kept constant and equal to 3%. A reasonable agreement has been found between the two codes for the measured growth rate and frequency. Figure 13: Scan in T EP . TAE growth rate and frequency, calculated with LIGKA and ORB5 for an EPs concentration equal to 3 % (same density, temperature profiles in use).

Nonlinear simulations
In this subsection results involving the nonlinear dynamics of the Alfvén waves are presented, when both on-axis (Fig.14) and off-axis (Fig.15) density profiles for the EPs are considered. With an on-axis density profiles of the EPs, a mode sitting in the frequency gap is observed (TAE), Fig.14. Its mode structure and frequency spectra are not observed to change passing from the linear to the nonlinear phase, confirming its nature of an eignemode of this system, which is only weakly perturbed by EPs. Figure 14: Mode structure and frequency spectra in the linear phase (left) and in the saturation phase (right). EPs have ah on-axis profile.
When an off-axis density profile for the EPs is considered, a mode with dominant poloidal mode number m = 2 and peaked around r 0.22 is observed (see Fig.15). This is consistent with what was observed in the previous sections, when just the linear effects in the simulations were involved. Passing to the nonlinear phase a secondary mode with m = 2 and m = 3 is observed to grow around the radial position r 0.738. This second mode is identified as the previously described TAE. This happens, because in the first linear phase the EPs drive the EPM unstable, which appears in fact to be dominant. In the nonlinear phase, the coexistence of the EPM and TAE is observed, due to an earlier saturation of the EPM (see Fig.16) .

Conclusion
The presence of Alfvén modes in burning plasma can affect negatively the energy confinement and can also cause a damage in the confining machine. Because of their importance, the present paper has dealt with the main damping mechanisms affecting the Alfvén modes, trying to outline them and to understand in which regime they are acting. Among the great zoology of existing Alfvén modes, the attention has been focused on toroidal Alfvén eigenmodes and energetic particle modes. These studies have been mainly carried by means of numerical simulations conducted with the code ORB5. The obtained results have been compared, when possible, with the presented analytic theory developed in a simplified geometry, and with the results of the linear code Ligka.
In Section 3 simulations with very small inverse aspect ratio ( = 0.01) have been considered, in order to lead the analysis in the cylinder limit. Simplified profiles have been taken into account and very low electron temperature has been considered in order to have the continuum damping dominant over the Landau damping. The developed theory has been used to analyze the results of simulations without energetic particles. The dependence of the radial wave number of the mode against the time have been observed (phase mixing). Also the scalar potential has been found to decay as δφ ∝ |ω A t| −1 (continuum damping).
In Section 4 higher bulk ion and electron temperatures have been considered, in order to observe the Landau damping to be dominant over the continuum damping. The numerical simulations have been conducted using plasma equilibrium and profiles from the ITPA-TAE international benchmark case. In order to separate the ions and electrons contributions to the damping, a kinetic term to the ideal MHD vorticity equation has been added. Following a perturbative approach, a simplified estimation for the Landau damping has been obtained using cylindrical coordinates. The developed analytical theory has been compared with the dependence found in numerical simulations of the damping rate against the bulk electron temperature. A reasonable agreement has been found and and this has also proved that the electron are the main responsible for the damping.
In Section 5 a realistic plasma equilibrium taken from a shot in ASDEX Upgrade has been considered. The results of the linear numerical simulations have shown the dependence of the damping rate against the bulk electron temperature describing, also in this case, the action of the Landau damping. A theoretical estimation of the radial regions where the Landau damping is expected to be dominant over the continuum damping has been presented, for simulations near marginal stability. A benchmark with the code LIGKA has shown good agreement for the frequency and growth rate dependence on the EPs temperature. Finally, the nonlinear simulations have shown the interaction of an EPM and a TAE in the scenario with off-axis EPs density profile.
The future works will extend the developed theory in order to find a better agreement between the predicted estimation of the decaying mode and the numerical simulations. In Ref. [2] it was suggested that all the Alfvén fluctuations can be explained within the framework of a single general fishbone-like dispersion relation (GFLDR). This could represent the starting point to improve the anaytical prediction of the damping rate and it would be a very interesting analytical and numerical task. Future and deeper benchmark with the code Ligka and the Hybrid Magnetohydrodynamics Gyrokinetic code HYMAGYC [27] would be of great interest in order to better understand the linear and nonlinear dynamics contained in the NLED-AUG case.

Acknowledgments
Simulations presented in this work were performed on the CINECA Marconi supercomputer within the ORBFAST and OrbZONE projects.
One of the authors, F. Vannini, would like to thank Xin Wang, Zhixin Lu and Gregorio Vlad for useful, interesting discussions and for great help provided in understanding Alfvén dynamics, Gyrokinetic and MHD theory. The authors wish to acknowledge stimulating discussions with F. Zonca, G. Fogaccia, A. Könies, J. Gonzalez-Martin 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 number 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.