Nonlinear dynamics of energetic-particle driven geodesic acoustic modes in ASDEX Upgrade

Turbulence in tokamaks generates radially sheared zonal flows (ZFs). Their oscillatory counterparts, geodesic acoustic modes (GAMs), appear due to the action of the magnetic field curvature. The GAMs can also be driven unstable by an anisotropic energetic particle (EP) population leading to the formation of global radial structures, called EGAMs. The EGAMs might play the role of an intermediate agent between the EPs and thermal plasma, by redistributing EP energy to the bulk plasma through collisionless wave-particle interaction. In such a way, the EGAMs might contribute to the plasma heating. Thus, investigation of EGAM properties, especially in the velocity space, is necessary for precise understanding of the transport phenomena in tokamak plasmas. In this work, the nonlinear dynamics of EGAMs is investigated with the help of a Mode-Particle-Resonance (MPR) diagnostic recently implemented in the global gyrokinetic (GK) particle-in-cell code ORB5. This enables to investigate the relative importance and the evolution of the resonances responsible for the ion and electron Landau damping, and of the EP drive. An ASDEX Upgrade discharge is chosen as a reference case for this investigation due to its rich EP nonlinear dynamics.


I. INTRODUCTION
An energetic particle (EP) beam, injected in tokamak plasma, can excite a variety of modes.
One of these modes, called energetic-particle-driven geodesic acoustic mode (EGAM) 1,2 , can have a significant influence on the EP dynamics and on the plasma confinement. First of all, being excited, this mode provides an additional mechanism of the energy exchange between the energetic particles and the thermal plasma [3][4][5][6] due to the wave-particle interaction. Moreover, because of its interaction with the turbulence 7,8 , this mode might be used as an additional knob in the turbulence regulation. Significant progress has been made in the last decade in building a theoretical model which can explain the main nonlinear physics of the EGAM 9-11 . Yet, when one wants to quantitatively compare with experimental measurements, some of the approximations in the previous models can be limiting. For example, due to the importance of the wave-particle resonances with the thermal and energetic species, a kinetic treatment should be used as in Ref. 12-14. Apart from that, effect of the magnetic surface shape 15 can be considered by employing an experimental magnetic equilibrium. In this work, the importance of these effects on an experimental ASDEX-Upgrade (AUG) case by means of numerical simulations with the nonlinear gyrokinetic code ORB5 [16][17][18] is investigated. ORB5 has been previously used for nonlinear studies of EGAMs in simplified configurations 19,20 , and it is used here to compare the nonlinear EGAM dynamics with experimental measurements.
In this work, only the EGAM dynamics, leaving aside its interaction with the turbulence is going to be considered. The nonlinear behaviour of this mode, such as the mode chirping 21,22 , plasma heating by the EGAMs, and its dependence on the plasma parameters are of interest here.
It should be noted here that in contrast to finite-n modes (where n is a toroidal number), which also can transfer energy from the EPs to the thermal plasma due to the wave-particle interaction, the EGAM dynamics is practically not associated with additional particle loss, at least if one does not consider the mode propagation and topological orbital changes 23 .
The EGAM dynamics is going to be investigated in the plasma configuration of the ASDEX Upgrade (AUG) discharge #31213@0.84 s (referred to as the NLED-AUG case) 24,25 , where a rich EP nonlinear dynamics is present 26 . In this discharge, various types of EP-driven instabilities were identified, among which there are Alfvén instabilities (see, for example, Ref. 27) and EGAMs.
In particular, a clear EGAM chirping was observed. More precisely, at the experimental EGAM frequency spectrogram (Fig. 1a) one can see at least two different branches of the EGAM evolution in time. A "short" one is with a ∼ 2.5 ms period with the frequency rise from ∼ 45 kHz till ∼ 55 kHz. A "long" branch is characterised by an ascending frequency from ∼ 45 kHz till ∼ 58 kHz and frequency saturation with a period of change around 13 ms. The paper structure is as follows. The ASDEX Upgrade configuration is presented in Section II.
The Mode-Particle-Resonance (MPR) diagnostic 14 , implemented in ORB5 to study field-particle interaction in velocity space, is briefly described in Section II A.

II. NUMERICAL AUG CONFIGURATION
The AUG plasma is simulated in the gyrokinetic (GK) particle-in-cell (PIC) global code ORB5 [16][17][18] , where a realistic magnetic configuration up to the last closed flux surface, constructed by the Grad-Shafranov solver code CHEASE 29 from experimental data, is used (Fig. 1b). The magnetic field at the magnetic axis is B 0 = 2.2 T. The major radius at the axis is R 0 = 1.67 m.
The geometrical major and minor radii are R 0 = 1.62 m and a 0 = 0.482 m respectively. Realistic temperature and density profiles of the thermal species (deuterium and electrons) are used as well (Figs. 2a-2b). Under "energetic particles" we understand here a deuterium beam, which actually drives EGAMs in a plasma system. This species is described by a two-bumps-on-tail distribution function in velocity space (Fig. 2c) which leads to the EGAM formation through the wave-particle interaction (due to the mechanism of the inverse Landau damping). Here, p z is the canonical parallel momentum, µ is the magnetic moment, B is the background magnetic field, while v ,EP , T EP being constant input parameters, which specify a shift and width of the bumps respectively. It means, that the EPs are described by space is explained by the necessity to avoid the input of an additional momentum into the plasma system, which might have some effect on the EGAM dynamics.
In the nonlinear (NL) electrostatic simulation, performed in Section III to study the EGAM up-chirping, the EPs have a realistic density profile, shown in Fig. 2a. In other sections, the EPs are described by an axisymmetric Gaussian distribution function: and are localised at a radial point s EP = 0.50, with s = ψ/ψ edge being a radial coordinate, where ψ is the poloidal flux coordinate. Such a simplification has been done to have more freedom in the variation of the EP parameters. The radial width of the EP Gaussian σ EP , used in this work, is equal to 0.10 to have a localised EP beam. An EP concentration is defined as n EP /n e , where n EP , n e are the EP and electron densities, averaged in volume, respectively.
In the simulations, presented in this work, three species are taken into account: gyro-kinetic thermal deuterium, gyro-kinetic energetic deuterium (EP), and electrons, either adiabatic (AE) or drift-kinetic (KE). Simulations with AE have been performed electrostatically, while the ones with KE have been made electromagnetically. The experimental pressure ratio at the reference magnetic surface has been adopted: Electromagnetic simulations have the advantage of the absence of the unstable ω H -mode 30

A. Mode-Particle-Resonance diagnostic
To localise velocity domains of the most intensive EGAM-particle interaction, the power balance diagnostic, implemented in ORB5 33 , has been recently extended 14 by keeping information about the wave-plasma interaction in velocity space. Using this diagnostic, a damping/growth rate γ of an ES wave can be calculated as 14,33 where P is the energy transfer signal (heat rate), E f is the field energy, f sp (r, p z , µ,t) = F 0,sp (r, p z , µ)+ δ f sp (r, p z , µ,t) is the species distribution function, Φ, J 0,sp Φ are the ES potential perturbation and the gyroaveraged one,Ṙ 0,sp is the species equations of motion, unperturbed by the field perturbations, Z sp e is the species charge. The time integration in Eq. 5 is performed on several wave periods n w T w . Finally, the integration in Eq. 6 and Eq. 7 is performed over the real V and species velocity W sp domains. For the equilibrium distribution functions F 0,sp , considered in this work, the part of P related to F 0,sp is negligible. The energy transfer is normalized in the following way:

III. GK SIMULATION OF EGAM CHIRPING IN THE AUG DISCHARGE
To model the nonlinear evolution of the EGAM frequency, the EP parameters are taken consistent with previous works (Refs. 14 and 15), where the order of the mode frequency has been reproduced in linear GK simulations. The following EP velocity v ,EP = 8.0, temperature T EP = 1.0 and concentration n EP /n e = 0.095 are applied. In Fig. 3a, one can see a zoom of the 'long' branch, Since the EGAM is a strongly driven mode, it significantly depends on EP parameters. Moreover, the considered simulation is performed electrostatically without taking into account the dynamics of the drift-kinetic electrons. As it is shown in Section VI, electrons can significantly modify the EGAM behaviour and, as a result, might have some effect on the mode chirping as well.  (Fig. 4a) and linear growth rate (Fig. 4b) on the temperature and velocity shift of the EP beam.

IV. INFLUENCE OF EP PARAMETERS ON THE EGAM DYNAMICS
To investigate the plasma heating by the EGAMs, we start from the study of the mode behaviour dependence on EP parameters. First of all, the linear growth rate is strongly modified by EP frequency close to the original GAM frequency. The low velocity energetic beam (Fig. 6a), which interacts with the geodesic mode through a high order m = 2 resonance, drives an EGAM with a frequency close to half the GAM frequency.
As it was mentioned in Refs. 6 and 23, interactions between the EPs and thermal species significantly benefit from the existence of high order resonances, which is confirmed here as well. The EGAM-thermal ion energy exchange occurs at the higher order (m = 2) resonance in both cases, independently of the EP parallel velocity (Fig. 5b). In other words, even when the EGAM is driven by an EP beam with a high parallel speed, the EGAM is still damped by the thermal deuterium through resonances at higher order. This can be simply explained by the fact that the position of these resonances are closer to the bulk of the thermal ion distribution function. Since the thermal ion energy transfer signal (Fig. 5b) is positive, and the EP energy transfer signal (Fig. 5a) is mostly

V. PLASMA HEATING BY EGAMS
As it has been already discussed previously in Refs. 3-6, the EGAMs might provide an additional route for plasma heating by EPs. The thermal species can obtain energy from the EGAM due to the Landau damping of the mode, which in turn is driven by the energetic particles through the inverse Landau damping. Note that this is also true for Alfvén modes, with the substantial difference that the EP radial redistribution due to EGAMs is negligible with respect to Alfvén modes.
Therefore, EGAMs represent a privileged mode for this plasma heating mechanism. This process is investigated here by varying the EP parameters in the AUG plasma configuration by performing nonlinear GK simulations firstly in the ES limit with adiabatic electrons.
To calculate the amount of energy transferred from the mode to the thermal ions, the EGAMspecies energy exchange signal P sp is integrated in time. As it is shown in Fig. 7, at the beginning,  (Fig. 8a). Evolution of the thermal ion temperature in the same cases (Fig. 8b).
a rise of the heat rate due to the growth of the mode is observed. After a while, the mode saturates, the heat rate level decreases due to the relaxation of the EP distribution function, and in the deep saturated domain the energy transfer from the mode to the thermal ions is practically suppressed. It means that the total energy transferred to the bulk plasma (red points in Fig. 7) remains practically unchanged after the mode saturation. The amount of this energy depends on the EP parameters ( Fig. 8a), which is directly reflected in the bulk ion temperature evolution (Fig. 8b).
By varying the EP parameters in NL electrostatic simulations, one can observe a general correlation between the EGAM saturation levels and the amount of energy transferred from the mode to the bulk deuterium plasma (Fig. 9). However, for example, the case with v ,EP = 6.0, T EP = 0.4, n EP /n e = 0.01 (the most right red star) has practically the same saturation level as the case with v ,EP = 3.5, T EP = 0.15, n EP /n e = 0.09 (the most right blue point). Having said that, we should notice that the corresponding EGAM-bulk ion energy exchange for these cases is significantly different. This means, that having the same mode level, one can achieve higher plasma heating by varying EP parameters. A possible explanation is that the EP beam with a lower velocity v ,EP = 3.5 drives the EGAM through the resonances of the same order as those responsible for the mode interaction with the thermal ions (Fig. 5). In other words, the plasma heating by the

VI. INFLUENCE OF ELECTRON DYNAMICS
In this section, the effect of the drift-kinetic electrons on the EGAM dynamics is examined, by considering electrons with a realistic mass m d /m e = 3672. In Fig. 10a, one can observe the EGAM growth rate significantly decreases when drift-kinetic electrons are included. The interaction between the EGAMs and the electrons is investigated by the wave-particle energy transfer signal in the velocity space. In Fig. 10b, where ε is an inverse aspect ratio, the parallel velocity v is normalized to the sound speed c s , and the magnetic moment µ is normalised to m d c 2 s /(2B 0 ). The energy exchange between the EGAM and trapped electrons occurs inside of this cone. The fact that the energy transfer resonances are close to the passing-trapped boundary, means that the mode interacts mainly with the barely trapped electrons with a bounce frequency close to that of the mode.
The inclusion of the electron dynamics influences the nonlinear EGAM behaviour as well.
Consistently with the decrease of the EGAM linear growth rate, the mode saturation levels are also reduced in NL simulations (Fig. 11a) that leads to the lowering of the bulk plasma heating by the mode. This reduction is clearly observed in the time evolution of the thermal ion temperature (Fig. 12b). However, even in the case with the drift-kinetic electrons, the EGAM transfers most of its energy to the thermal ions, and not to the electrons, that can be seen in Fig. 12a. In other words, it is the ion plasma heating, which mainly benefits due to the presence of the EGAM dynamics.  (Fig. 11a). Evolution of the thermal deuterium temperature (Fig. 11b). By varying the EP parameters, general correlation between the mode level and the EP-thermal plasma energy flow had been revealed. Apart from that, it has been emphasized that the plasma heating by EGAMs benefits from the presence of high order mode-particle resonances which are often responsible for the EGAM-thermal species interaction. More precisely, it has been shown that having the same mode levels, one can have higher energy transfer from the EPs to the bulk plasma if the EGAM drive and damping occur through the resonances of the same order.
Finally, it was indicated that although the EGAM transfers most of its energy to the thermal ions, and not to the electrons, the electron dynamics might significantly reduce the plasma heating by EGAMs by lowering the mode amplitude.
Since the nonlinear GK simulations with drift-kinetic electrons have shown a significant influence of the electron dynamics on the EGAM behaviour, it would be reasonable to study effect of the kinetic electrons on the mode chirping. The results obtained here will be useful in a future work to perform investigation of the EGAM influence on the turbulent dynamics in the ASDEX Upgrade, that will be a next step in the code validation and an interesting example of a complex plasma system.

ACKNOWLEDGMENTS
The authors would like to acknowledge stimulating discussions with the Multi-scale Energetic