{'fraction_non_alphanumeric': 0.034482758620689655, 'fraction_numerical': 0.014560689069229192, 'mean_word_length': 4.424984106802289, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Implementation of energy transfer technique in ORB5 to study collisionless wave-particle interactions in phase-space

A new diagnostic has been developed to investigate the wave-particle interaction in the phase-space in gyrokinetic particle-in-cell codes. Based on the projection of energy transfer terms onto the velocity space, the technique has been implemented and tested in the global code ORB5 and it gives an opportunity to localise velocity domains of maximum wave-plasma energy exchange for separate species. Moreover, contribution of different species and resonances can be estimated as well, by integrating the energy transfer terms in corresponding velocity domains. This Mode-Plasma-Resonance (MPR) diagnostic has been applied to study the dynamics of the Energetic-particle-induced Geodesic Acoustic Modes (EGAMs) in an ASDEX Upgrade shot, by analysing the influence of different species on the mode time evolution. Since the equations on which the diagnostic is based, are valid in both linear and nonlinear cases, this approach can be applied to study nonlinear plasma effects. As a possible future application, the technique can be used, for instance, to investigate the nonlinear EGAM frequency chirping, or the plasma heating due to the damping of the EGAMs.


I. INTRODUCTION
Gyrokinetic (GK) codes have recently become standard tools for the investigation of waves and instabilities in tokamak plasmas, with frequency below the ion cyclotron frequency 1 . Although they have been traditionally considered numerically heavy, in comparison to lighter hybrid models, in the last years GK codes have become capable of providing global electromagnetic predictions of the nonlinear plasma dynamics, thanks to smart schemes improving the numerical performance 2,3 , and to the access to high-performance computers. One advantage of using GK codes is that their model includes kinetic effects such as wave-particle resonances, which are neglected in fluid descriptions.
Wave-particle interaction, such as Landau damping, can be best detected by phase space resolving diagnostics. In particular, investigating collisionless energy transfer signals as a function of particle velocity, necessary details can be provided to identify dominant collisionless processes governing the damping or growth of electrostatic (ES) zonal modes, such as geodesic acoustic modes (GAMs) 4-6 or energetic-particle driven GAMs, called EGAMs 7-10 .
There are different kind of techniques to investigate dynamics of modes in the phase-space.
Correlation techniques 11 can be used to clarify the origin of the energy-transfer process and the nature of mechanisms that lie beyond observed mode dynamics by calculating correlations of the energy transfer terms with different fields signals. Conjunction diagnostics based on the measurements at different positions along the same magnetic flux tube can be used to study the integrated effect of wave-particle interactions between the two space points 12 . The conjunction studies are particularly well suited to study the waves, that are propagating along the magnetic field lines, such as shear-Alfvén waves.
In this work we develop a Mode-Particle-Resonance (MPR) diagnostic in the code ORB5 13,14 to investigate energy transfer signals in velocity space in global gyrokinetic (GK) simulations. The previous version of this diagnostic gave only time evolution of the energy transfer terms, averaged over the whole phase space 15,16 . We extend it, by taking projection of these terms onto the velocity space, that gives an opportunity to investigate the contribution of different resonances in different velocity domains to the mode dynamics.
This technique is applied in global GK simulations of an experimental shot on ASDEX-Upgrade machine to study EGAMs. These modes are characterised by the oscillations of mainly toroidally symmetric global radial ES field with frequency comparable to that of the GAMs. The energetic particles (EPs) excite the mode through the inverse Landau damping, and EPs are displaced from higher to lower energy range 17,18 . On the other hand, the GAMs and EGAMs are mainly damped by Landau damping. In addition to ion Landau damping, GAMs have been found to be subject to the electron Landau damping [19][20][21][22] as well. Here, we show that EGAMs are also subject to electron Landau damping, which can be as important as ion Landau damping in experimentally relevant conditions. Moreover, in these simulations the MPR diagnostic provides additional details to clarify the role of different species in the EGAM-plasma interaction.
The GAMs and EGAMs can play a significant role in the regulation of the turbulencetransport processes. The GAMs are an oscillating branch of zonal flows 5,23 (ZFs). The ZFs can reduce the radial transport in tokamak plasma acting as a sink for the turbulence energy through the inverse energy cascading or/and by shearing plasma eddies [24][25][26][27] . Contrarily, the role of the GAMs [28][29][30][31] in the turbulence suppression is still unclear and even contradictory 32,33 . It could be explained by the fact that the GAMs can transfer the energy in both directions. They can either take the energy from the turbulence, being directly excited by instabilities 28 and arising from the ZFs due to the magnetic curvature, or they can return the energy back to the instabilities 32 . Due to this complex dynamics, the role of the EGAMs in the turbulence suppression is still a subject of study [34][35][36][37][38] . At the same time, EGAMs might play a role of an intermediate agent between the fluctuating fields and thermal plasma, by spreading fluctuating field energy to the bulk plasma through the collisionless wave-particle interaction 18 . In such a way, the EGAMs might be a crucial component in tokamak plasma stabilisation and be significantly helpful in the plasma heating.
Thus, investigation of EGAMs characteristics, especially in the velocity space, is necessary for precise understanding of the transport phenomena in fusion reactors, where the EPs are produced as the result of the nuclear fusion reaction or by external sources such as neutral beam injection (NBI) or ion cyclotron resonance heating (ICRH).
The remainder of the paper is structured in the following way. In section II the theoretical background and the implementation of the MPR diagnostic in ORB5 is presented. After that, the processing of the output signals from the diagnostic is demonstrated, and an example of an ES simulation of the GAMs in a circular magnetic configuration is given (Sec. III). In Sec. IV we show that the MPR diagnostic verifies the GAM dispersion relation.
Having discussed the technique, the experimental AUG shot #31213 is investigated in section V using the developed diagnostic in linear simulations. Results, calculated in ORB5, are compared with simulations of the GENE code for an ES case with adiabatic electrons in Sec. VI.

A. Theoretical background
The MPR diagnostic is based on the projection of energy transfer terms onto the velocity space. It gives an opportunity to localize velocity domains of maximum mode-plasma interactions at particular time moments or averaged in a specific time interval, that is more relevant for the modes that oscillate in time. By integrating the energy transfer terms in a chosen velocity domain and normalizing them to the mode energy, one gets a damping or growth rate of the mode. More precisely, considering the case of GAMs/EGAMs, the theoretical background of the diagnostic can be explained using the Poynting's theorem 39 of electromagnetics: where E is the mode energy, S is the Poynting's vector, J sp is the current, produced by a species sp, and E is the electric field. By integrating over the whole real space, we get ∇ · S dV = 0 using the Gauss's theorem. Introducing the mode complex frequency: the mode energy evolves in time as from where we can get an expression for the mode complex frequency: Finally, using Eq. 1, we get an expression for the mode damping rate: Since the GAMs/EGAMs oscillate in time, it is necessary to take time averaging on several GAM/EGAM periods to get a pure damping or growth rate of these modes. Moreover, the GAM/EGAM electric field energy is transformed periodically to the plasma pressure perturbation. To take into account this change we consider only the envelop of the ES field energy. Speaking about P, we are dealing with the work done by the ES field on the plasma. Since the mode energy is taken to be only positive, with a chosen convention for signs in Eq. 4-3, a negative rate γ < 0 corresponds to a positive signal P, indicating the energy transfer from a wave to plasma particles. On the other hand, a positive rate γ > 0 corresponds to the growth of the wave.

B. Discretization
To describe the implementation of the diagnostic in the code ORB5, we should start from the discretization of the plasma distribution function in the code. ORB5 is a particlein-cell (PIC) code, where the Vlasov equation is solved using a Monte Carlo algorithm, and the Maxwell equations are solved using a finite-element method. At the beginning of a simulation a finite collection of initial positions in phase space is sampled by a set of numerical markers 14 . Every marker has a particular magnetic moment µ sp = m sp v 2 ⊥,sp /(2B), a position in real space R sp , a parallel canonical momentum p z,sp and it is moving in a background magnetic field B = bB with Here, c is the speed of light, m sp and Z sp e are the species mass and charge, where for electrons Z e e = −e and e is the absolute value of the electron charge. Taking a phase-space position Z = (R sp , p z,sp , µ sp ) of a species marker as a random variable, the code distributes the markers in the phase space according to the initial particle distribution function f 0,sp .
It means, that each marker is a realisation of the random variable Z. For simplicity, a marker will be considered as a particle that is moving along a particular orbit defined by the following equations of motion: which are obtained by varying a GK Lagrangian with respect to the phase-space coordinates Z = (R sp , p z,sp , µ sp ) 16,40 . The orbits are perturbed by the field perturbation with Φ and A being electric and parallel magnetic potential perturbations respectively, where only Φ remains in ES simulations. In the gyro-kinetic approximation the code deals with the dynamics of the gyrocentres, whose orbits are perturbed by the potentials, averaged in a space domain, defined by the species Larmor radius, around a marker position.
This averaging is represented by the operator J 0,sp . In the drift-kinetic approximation, the potential perturbation is considered at a space point, where a marker is localised, without performing the gyro-averaging. In ORB5 the thermal and fast ions can be treated either gyro-kinetically or drift-kinetically, while the electrons are calculated drift-kinetically.
The time evolution of the plasma distribution function f sp is described by the Vlasov equation: Considering perturbations of the distribution function and of the particle orbits till the first order, one can linearize the Vlasov equation: ∂δf sp ∂t +Ṙ 0,sp · ∇δf sp +ṗ 0,z,sp ∂δf sp ∂p z,sp = − ∂f 0,sp ∂t +Ṙ 0,sp · ∇f 0,sp +ṗ 0,z,sp ∂f 0,sp ∂p z,sp − Ṙ 1,sp · ∇f 0,sp +ṗ 1,z,sp ∂f 0,sp ∂p z,sp Assuming that f 0,sp is an equilibrium distribution function, it should be conserved along unperturbed particle trajectories (Ṙ 0,sp ,ṗ 0,z,sp ): In other words, the first bracket on the right hand side of Eq. 15 is equal to zero. Finally, the time evolution of the perturbation of the species distribution function in linear simulations is described in the following way: where 1 indicates that it is necessary to take derivatives along the perturbed parts of species orbits (Ṙ 1,sp ,ṗ 1,z,sp ). Thermal species have an equilibrium distribution function in a form of the Maxwellian one: where n sp (ψ), T sp (ψ) are species density and temperature profiles along the radial coordinate ψ, which is the poloidal flux. A symmetric two-bumps-on-tail distribution function has been used in this work for the fast species 41,42 . This distribution assumes a flat temperature profile of the fast species: where u H,sp , T H,sp are constant input parameters, which specify a shift and width of the bumps respectively.
The perturbation δf is discretized in the Z = (R, p z , µ) phase space by N sp markers.
Apart of its location Z, every marker has a particular weight w p (t), which should evolve consistently with the GK Vlasov equation Eq. 14. Here, we omit the index sp to simplify equations and use the index p, indicating that a variable is related to a particular marker.
Detailed derivation of the weight time evolution can be found in Ref. 13,14,43 . A marker weight can be associated to a phase space volume Ω p and correspondent averaged perturbation distribution function δf Ωp : Considering uniform spreading of the markers in real space and Maxwellian distribution in the velocity space, it can be shown 43 that the phase space volume Ω p , associated to a marker p, is whereJ(s) is the flux-surface-averaged Jacobian, κ v defines maximum value of the species parallel and perpendicular velocities, normalized to a species thermal speed u th (s) = T /m, at every radial point s = ψ/ψ edge .
The meaning of the variable Ω p can be explained proceeding directly from the Monte Carlo integration 14 . The expectation value of an arbitrary function ζ(Z) is whereZ is a random variable, distributed according to the function f . To minimize the variance of the function ζ, one can chose another distribution function g(Z), which does not vanish in the support of the distribution function f (so-called importance sampling): In this case, speaking in terms of marker weights and using random variable Z, distributed with density g, the expectation value of the function ζ(Z) is calculated as that is consistent with Eq. 23. In other words, if we have a small amount of markers in a finite phase space volume, their weights will be increased in comparison to a domain where there are higher number of markers at the same phase space volume. More details can be found in Ref. 14 .
To clarify different terms in Eq. 10, the characteristicṘ sp can been split on several terms: 34 indicates the dependence on the gradient of the kinetic plasma pressure p. A precise form of the GK energy transfer signal, valide in both linear and nonlinear cases, can be derived from the GK Hamiltonian using the Noether theorem as it is shown in Ref. 16 : with V and W sp being real and velocity spaces. By integrating the signal over the whole real space V and in a small velocity domain ∆W sp , related to a particular velocity bin, we project the energy transfer signal to the velocity space of a particular species: where the sum i∈V,∆Wsp is taken on all markers in the phase volume V ∆W sp . The gyroaveraged electric field −∇(J 0,sp Φ)| i is taken at a position of a marker i. The sum is normalized to a total number of species markers N sp in the whole phase-space domain. In the current version of the diagnostic, only the electrostatic part ofṘ sp is taken into account (Eq. 31 -34). Since the GK model, that is used in ORB5, is based on the Hamiltonian formulation 16 , p z,sp is used as one of the velocity variables: This is a common choice in most of the modern GK PIC codes. In the MPR diagnostic a variable u sp is used for a parallel velocity: In the ES case, the variables u sp and v ,sp are identical u sp = v ,sp , and in EM simulations with low β they are close u sp ≈ v ,sp . With the rise of β, the difference between these two variables increases because of the contribution of the parallel magnetic potential A . A proper transition from the variable p z,sp to the variable v ,sp (instead of u sp ) is necessary for the investigation of the dynamics of EM modes and for proper analysis of EM simulations.
It is a matter of future publications.

III. POST-PROCESSING
Here, a GAM in a circular magnetic configuration is considered to show how the diagnostic is organised, and how the MPR data are treated. A circular deuterium plasma with flat safety factor q = 1.5, and flat density and temperature radial profiles is considered. Hinton test 23 has been sent by introducing an axisymmetric density perturbation designed to produce an initial electric potential field.
First of all, the MPR diagnostic provides the energy transfer signal P(v , µ, t) (Eq. 38) as a function of the velocity variables (v , µ) and time. By averaging this signal on several GAM periods, resonances of the mode-particle interaction can be localised in the velocity space.
Their location can be compared with the analytically given parallel resonance velocity: where ω GAM is the GAM frequency, that can be found directly from the radial zonal electric field E r . Since the perturbation of the plasma distribution function related to the GAM dynamics can have higher poloidal modes m ≥ 1, the GAM-particle interaction can be observed at smaller parallel velocities as well v (m) By integrating in corresponding velocity domains, one can estimate contribution of these resonances to the mode dynamics. In this particular case, it can be seen from Fig. 1, that the energy transfer occurs mainly at the first resonance v ,res . By integrating the signal in the whole velocity domain, one gets the time evolution of P. By normalizing it to the mode energy, the GAM damping rate can be estimated using Eq. 6. The mode energy has been taken as an envelop (grey dotted line at the left plot in Fig. 2) of the field energy (green line and the distribution of the GAM damping rate values is shown in Fig. 2. The result from the MPR diagnostic can be compared with the direct calculation of the GAM damping rate, by fitting the zonal radial electric field E r at a particular radial point.
Here, the point s = 0.74 has been taken, since it is very close to a crest in E r , that can be seen from the left plot of Fig. 1. First of all, a zero-frequency component of E r is filtered out for more precise calculation of the GAM characteristics. After that, a GAM frequency is estimated, for example, by the Fast Fourier Transform. On the other hand, the damping rate is estimated by the linear least-square root method from the peaks in the time evolution of E r . This preliminary processing gives the first assumption of the GAM frequency and damping rate, that are used as initial guesses in the non-linear fitting procedure. A function is used as a test one, which is fitted to the time evolution of the E r (s = 0.74). This method has been used previously in Ref. 21  Finally, the GAM frequency and damping rate, found using the non-linear fitting of E r (s = 0.74) to the test function Eq. 44, are the following: and have been calculated from the corresponding distribution functions, shown in Fig. 3.
As it can be seen here, the calculation of the GAM frequency is quite precise with an errorbar being around 1%, while the errorbar of the damping rate prediction is around 20%.
Comparing both methods (Eq. 43 and 46), one can see that the MPR diagnostic is not as precise as the non-linear fitting, at least, in case of the calculation of the GAM damping rate.
On the other hand, it provides additional information such as a position of the GAM-plasma resonances in the velocity space (Fig. 1).

IV. ANALYTICAL VERIFICATION
Here, we are going to show the consistency of the MPR diagnostic by comparing the measurements on a GAM investigated with ORB5, with the analytical dispersion relation derived in the GK framework, by neglecting the effects of the finite Larmor radius and finite orbit width, and considering adiabatic electrons. The corresponding GAM dispersion relation (Ref. 28,45 ) is We omit species indices, since all relevant plasma variables are related to the deuterium. A GAM is described by the evolution of the zonal electric field: with a radial wavenumber k and a complex frequencyω = ω D + iγ D , that verifies the dispersion relation Eq. 47. The corresponding perturbation of the deuterium distribution function is: where c is the light speed, ω tr = v /(qR 0 ) is the passing frequency, θ p is the poloidal angle in a simplified circular geometry, F 0 being the deuterium equilibrium distribution function: and v d being the amplitude of the radial drift composed by the curvature drift and grad-B To derive an expression for the energy transfer term, we need the equation of motion Eq.
10, which in a linear ES system can be rewritten as: where B * sp and B * ,sp are defined in Eq. 8 and 9 respectively. Considering a low-pressure plasma (J 0 × B ≪ 1 with a plasma current J 0 ) in a circular plasma cross-section with a curvature κ, one gets the following simplifications: Applying the introduced approximations, we get thaṫ using which together with Eq. 55 and by putting everything to Eq. 37, an expression for the energy transfer signal can be derived: Evaluating the velocity integrals I 1 and I 2 , one gets the expression which can be significantly simplified using the GAM dispersion relation (Eq. 47) to get rid of the functions N(z), D(z), and F (z): As a result, we have the following complex expression for the plasma-field energy exchange signal: In the GK model of ORB5, the ES field energy is given (Ref. 14 ) by: with Φ being an ES potential perturbation. Since we consider the drift-kinetic approximation here, the gyro-averaging operator J 0 is equal to 1. Using the GK Poisson equation, the expression can be reduced to Finally, taking into account only the perturbation of the zonal radial electric field |∇ ⊥ Φ| 2 ≈ E 2 r and since dW F 0 = 1, we get Since dV E 2 r sin θ p / dV E 2 r = 1/2 and applying Eq. 68, 71 to Eq. 6, the consistency between the GAM dispersion relation and the MPR method can be proved: The expression Eq. 72 means that using the field and plasma perturbations (Eq. 54 and Eq. 55), which verify the GAM dispersion relation Eq. 47, in the MPR diagnostic (Eq. 6), one gets a GAM damping rate that verifies the starting GAM dispersion relation.
To check the time behaviour of the analytical energy transfer signal, one should take the real part of Eq. 68: Taking into account, that all signals in the code ORB5 are real, we should also consider the time behaviour of the energy transfer, got from only real parts of the plasma perturbation and zonal radial electric field: Since the GAM damping rate is a small value in comparison with the GAM frequency, the Eq. 73 and Eq. 74 give the same time evolution. On the other hand, we should take into account the zero-frequency component of the zonal radial electric field: In this case, the analytical energy transfer signal takes the following form: If we take the GAM frequency (Eq. 45) and damping rate (Eq. 43 or 46) from the sim- Eq. 76 (green line). Here, the ratio E r,0 /E r,1 ≈ 0.25, which is used in Eq. 76, is estimated from the E r (s = 0.74), given by the numerical simulation in ORB5.

V. APPLICATION TO EGAMS IN AUG SHOT #31213
A. Equilibrium and definition of the numerical simulation The AUG shot #31213 at time 0.84 s has been selected within the Non-Linear Energetic- done including dynamics of the magnetic potential perturbation as well. In this latter case the pullback method 3 has been used for the mitigation of the cancellation problem in EM simulations 2,48 . Corresponding profiles of the safety factor, species density and temperature are shown in Fig. 5. The magnetic field is reconstructed with experimental data, including all geometrical effects (Fig. 5). 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 where the mode-plasma resonances are observed (Fig. 6). In the EM case the radial domain has been reduced to s = [0.0, 0.9] to avoid numerical instabilities due to the abrupt increase of the safety factor at the edge. The density profile, that is depicted in Fig. 5, corresponds to the case with β e = n e T e /(B 2 0 /(2µ 0 )) = 2.7 · 10 −4 , where n e is the electron density, averaged in a tokamak volume, µ 0 is the magnetic constant, and T e is measured at the radial position s = 0.0. In both cases, the velocity distribution of the fast particles is described Considering firstly the ES case with AE, we compare the EGAM frequency and growth rate, calculated at radial positions s = 0.40 and s = 0.48, using the non-linear fitting of E r . They appear to be consistent within the error bars of the measurements: Consistency between the EGAM growth rate, calculated directly from E r and by the MPR diagnostic, using Eq. 6, significantly improves in time due to the growth of the EGAM signal in comparison with the zero-frequency zonal flow. Skipping initial transient time period, the MPR diagnostic can be applied to measure the EGAM growth rate, that appears to be consistent with Eq. 78: The consistency between both methods is observed in the EM case with KE as well: From Eq. 77 and Eq. 82 one can see that the change in the EGAM frequency is small in comparison with the change in the growth rate, when dynamics of the drift-kinetic electrons is included. In particular, the EGAM growth rate decreases from Eq. 81 to Eq. 84. We now want to investigate the role of the drift-kinetic electrons in the EGAM dynamics to understand which wave-particle interactions lead to the decrease of the EGAM total growth rate, by estimation of the contribution of different species. In the simulation with adiabatic electrons: These equations show that the total EGAM growth rate is a balance between the drive on the fast species and damping on the thermal one (one can see also Ref. 42 for a similar analysis in the case of EGAMs in simplified configurations, with adiabatic electrons). Moreover, the absolute values of the species contributions are much higher than the absolute value of the EGAM total growth rate.
In case with drift-kinetic electrons, the species contributions are the following: We now want to investigate the role of the different resonances in phase space. In Fig. 8b one can see the energy transfer signal for the EGAM-electron interaction in the velocity space, averaged on several EGAM periods. The white cone there indicates an analytical estimation of the boundary between the passing-trapped electrons: where ǫ is an inverse aspect ratio. According to that figure, the EGAMs are damped by the electrons which are localised mainly near this boundary, similar to what happens for GAMs 19 .
We can separate two velocity domains e11 and e21, shown in Fig. 8b. The resonances in the domain e11 correspond to the EGAM interaction with the barely trapped electrons, while the domain e21 corresponds to the EGAM damping on the barely passing electrons. By averaging in the chosen velocity domains, one gets the time evolution of the energy transfer signal (Fig. 8a), that should be filtered for its proper use in Eq. 6. According to the MPR diagnostic, the contribution of the barely trapped electrons is much more significant than that of the barely passing electrons as one can see from Eq. 91 and 95: On the other hand, if we consider wider velocity domains, we can notice that there is still a significant contribution of electrons with higher parallel velocity to the EGAM damping.
One can see it, for instance, comparing the e13 and e14 velocity domains (Eq. 93 and 94).
The reason might be in the choice of the velocity space variables in ORB5, which has been explained in Sec. II in Eq. 39. We can estimate as well contribution of different resonances in the thermal deuterium velocity space (Fig. 8c): i31 : γ[ √ 2v th,i /R 0 ] = −6.9 · 10 −2 ± 7.8 · 10 −3 , In Eq. 101 the contributions of all considered resonances are summed up. The result value is close enough to the total contribution of the thermal ions to the EGAM damping (Eq. 87).
From one point of view, it is an additional option to verify the implemented diagnostic.
Moreover, it shows the dominant role of the m = 1 resonances (Eq. 42), which are localised in the velocity domains i12 and i21 (Fig. 8c), in the EGAM-thermal deuterium interaction.
On the other hand, by comparing Eq. 96 and 97, one can see that the parallel dynamics has a predominant contribution to the energy exchange between the EGAMs and the thermal deuterium plasma. Finally, we can consider different velocity domains in the EGAM -energetic deuterium interaction (Fig. 8d): One can see that there is an EGAM damping even on the energetic particles (Eq. 103 and 105). But it is significantly smaller than the dominant drive (Eq. 102 and 104). Sum on the resonances (Eq. 106) indicates that the EGAMs are driven by the fast species and its absolute value is close enough to the total drive, found in Eq. 89.

VI. COMPARISON WITH GENE
To Here, the function h 1,s represents the non-adiabatic part of the perturbed distribution function f 1,s . It is defined as h 1,s = f 1,s − qφ 1 /(B 0 T 0,s )∂f 0,s /∂µ. Eq. 107 is written in the field aligned coordinate system (x, y, z) with x the radial, y the bi-normal and z the field aligned directions. Moreover,φ 1 denotes the gyro-averaged electrostatic potential, J the phase- Fourier space (only in the local flux-tube limit) for each wave vector k = (k x , k y ) as follows Here, the bracket represents the field-aligned z-average, namely The time derivative of Eq. 109 determines the energy flow during the whole simulation time domain. In particular it represents the energy effectively transferred from the particles to the field. It reads as The same AUG shot has been simulated in GENE in case with adiabatic electrons (one can see also Ref. 47 ), using the flux-tube version of the code at s = 0.5. In Fig. 9 ORB5 : ω/2π = 42.7 ± 0.1 (kHz) On the other hand, there is 18% consistency between the codes for the EGAM total growth rate: and for the contributions of different plasma species to the mode dynamics: ORB5 : γ = (−87.6 ± 0.6) · 10 3 (s −1 ), f ast deuterium : GENE : γ = 115 · 10 3 (s −1 ), The difference in the values can be explained mainly by the fact that the simulation in GENE has been performed using the local flux-tube version, while the simulation in ORB5 is a global one.

VII. CONCLUSIONS
In this paper a Mode-Particle-Resonance (MPR) diagnostic has been implemented in the gyrokinetic code ORB5 to investigate mode-plasma interaction processes. The technique is based on the projection of energy transfer terms on the velocity space (Eq. 38) and gives an opportunity to localise velocity domains of maximum energy exchange between an electrostatic mode and different species. Moreover, integrating in a chosen velocity domain, a rate of the mode damping or growth can be calculated using Eq. 6 and contribution of different species to the mode dynamics can be estimated as well. Using a GAM dispersion relation, which neglects finite-Larmor-radius and finite-orbit-width effects and treats the electrons as adiabatic 28,45 , the theoretical principle, which lies behind the MPR diagnostic, has been analytically verified for an ES case. It has been shown that the GAM damping rate, derived from the energy exchange principle (Eq. 6), is identical to the GAM damping rate, given by the GAM dispersion relation Eq. 47. Analytical time evolution of the energy transfer signal, given in Eq. 68, has been found to have the same frequency as the numerical one (Fig. 4).
In Sec. V, the MPR diagnostic has been applied to the case of AUG shot #31213 (NLED AUG base case) to investigate contributions of different resonances to the EGAM dynamics.
It has been shown that inclusion of the drift-kinetic electrons significantly decreases the EGAM growth rate of about a factor 2 for the selected case (Eq. 81 and Eq. 84). It has been shown that the EGAM damping occurs at the first resonance v ,res in case of the interaction with deuterium plasma. On the other hand, in case of the electrons the EGAMs are damped mainly by the barely trapped electrons (Fig. 8). The total EGAM growth rate (Eq. 115) and contribution of the thermal (Eq. 117) and energetic deuterium (Eq. 119) to the mode dynamics has been calculated in the codes ORB5 and GENE in case with adiabatic electrons. The benchmark has shown 18%-consistency for the total growth rate and species contributions.
From the point of view of further possible application, the MPR diagnostic can be used, for instance, to study the energy exchange between energetic and thermal species indirectly through the zonal waves, such as EGAMs, that play a role of a mediator in this case. Other interesting effects, associated with the EGAM nonlinear evolution, are the EGAM frequency chirping, which consists in a fast modification of the mode frequency, and the saturation mechanisms. Since the frequency shift during the chirping is considered to occur as a result of the wave-particle interaction 60-65 , the MPR diagnostic can be used to investigate this phenomenon as well. The saturation mechanisms (wave-particle or wave-wave interactions) are important to investigate in order to build a theoretical model capable of predicting the saturation levels in experimentally relevant conditions, and as a consequence, the EP redistribution in phase space.
The current version of the diagnostic can be applied only to the case of mainly electrostatic modes, such as GAMs and EGAMs. As it has been discussed in Sec. II, the reason is in the choice of the velocity space variables in ORB5. The MPR diagnostic can be extended to work with EM simulations with arbitrary β, by performing a proper transition from the variable p z,sp to the velocity variable v ,sp . There are different possible areas of application of the EM-MPR diagnostic. A wider range of the modes whose dynamics is mainly controlled by wave-particle resonances, like energetic-particle driven MHD instabilities, can be investigated. For a turbulent plasma, the collisionless interactions between the EM fields and the plasma particles may lead to a secular transfer of energy from fields to particles,