Implementation of energy transfer technique in ORB5 to study collisionless wave-particle interactions in phase-space

TL;DR

This paper introduces a new diagnostic technique for analyzing wave-particle energy transfer in phase-space within gyrokinetic simulations, enabling detailed study of species contributions and resonances in plasma dynamics.

Contribution

A novel Mode-Plasma-Resonance diagnostic implemented in ORB5 that localizes energy exchange in velocity space and applies to both linear and nonlinear plasma phenomena.

Findings

Localized energy transfer in velocity space identified

Species-specific contributions to wave dynamics quantified

Application to EGAMs in tokamak plasma demonstrated

Abstract

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.

Equations

$$\frac{\partial ϵ}{\partial t}=-\nabla \cdot \mathbf{S}-\sum _{sp}{\mathbf{J}}_{sp}\cdot \mathbf{E},$$

$$\mathcal{E}(t)=\int ϵ\mathrm{d}V,$$

$$\hat{\omega }=\omega +i\gamma ,$$

$$\mathcal{E}(t)={\mathcal{E}}_0\exp (-2i\hat{\omega }t),$$

$$-2i\hat{\omega }=\frac{1}{\mathcal{E}}\frac{\mathrm{d}\mathcal{E}}{\mathrm{d}t},$$

$$\gamma =\sum _{sp}{\gamma }_{sp}=-\frac{1}{2}Re\left[{⟨\frac{\mathcal{P}}{\mathcal{E}}⟩}_t\right],$$

$$\mathcal{P}=\sum _{sp}{\mathcal{P}}_{sp}=\sum _{sp}\int {\mathbf{J}}_{sp}\cdot \mathbf{E}\mathrm{d}V.$$

$${\mathbf{B}}_{sp}^{\ast }=\mathbf{B}+\frac{c{p}_{z,sp}}{Z_{sp}e}\nabla \times \mathbf{b},$$

$$B_{\parallel ,sp}^{\ast }=\mathbf{b}\cdot {\mathbf{B}}_{sp}^{\ast }.$$

$${\dot{\mathbf{R}}}_{\mathbf{sp}}=$$

$$\frac{c}{Z_{sp}e{B}_{\parallel ,sp}^{\ast }}\mathbf{b}\times \left[{\mu }_{sp}\nabla B+Z_{sp}e\nabla (J_{0,sp}{\Psi }_{sp})\right],$$

$${\dot{p}}_{z,sp}=$$

$${\dot{\mu }}_{sp}=$$

$${\Psi }_{sp}=\Phi -\frac{p_{z,sp}}{m_{sp}c}{A}_{\parallel },$$

$$\frac{\mathrm{d}{f}_{sp}}{\mathrm{d}t}=\frac{\partial f_{sp}}{\partial t}+{\dot{\mathbf{R}}}_{\mathbf{sp}}\cdot \nabla f_{sp}+{\dot{p}}_{z,sp}\frac{\partial f_{sp}}{\partial p_{z,sp}}=0.$$

$$\frac{\partial \delta f_{sp}}{\partial t}+{\dot{\mathbf{R}}}_{0,\mathbf{sp}}\cdot \nabla \delta f_{sp}+{\dot{p}}_{0,z,sp}\frac{\partial \delta f_{sp}}{\partial p_{z,sp}}=$$

$$-$$

$$\displaystyle\displaystyle{\frac{\mathop{}\!\mathrm{d}f_{0,sp}}{\mathop{}\!\mathrm{d}t}}\bigg{|}_{0}=\displaystyle{\frac{\partial f_{0,sp}}{\partial t}}+\mathbf{\dot{R}_{0,sp}}\cdot\boldsymbol{\nabla}f_{0,sp}+\dot{p}_{0,z,sp}\displaystyle{\frac{\partial f_{0,sp}}{\partial p_{z,sp}}}=0.$$

$$\displaystyle\displaystyle{\frac{\mathop{}\!\mathrm{d}\delta f_{sp}}{\mathop{}\!\mathrm{d}t}}\bigg{|}_{0}=-\displaystyle{\frac{\mathop{}\!\mathrm{d}f_{0,sp}}{\mathop{}\!\mathrm{d}t}}\bigg{|}_{1},$$

$$f_{0,sp}^{therm}=\frac{n_{sp}(\psi )}{(2\pi {)}^{3/2}{u}_{th,sp}^3(\psi )}\exp \left[-\frac{m_{sp}}{T_{sp}(\psi )}\left(\frac{1}{2}{\left(\frac{p_{z,sp}}{m_{sp}}\right)}^2+\frac{{\mu }_{sp}B}{m_{sp}}\right)\right],$$

$$u_{th,sp}(\psi )=\sqrt{\frac{T_{sp}(\psi )}{m_{sp}}},$$

$$f_{0,sp}^{fast}=$$

$$\cosh \left(\frac{p_{z,sp}}{m_{sp}}\frac{u_{H,sp}}{T_{H,sp}}\right),$$

$$A(\psi )=$$

$$⟨\delta f{⟩}_{{\Omega }_p}=\frac{1}{{\Omega }_p}{\int }_{{\Omega }_p}\delta f\mathrm{d}{\Omega }_p=\frac{1}{{\Omega }_p}{\int }_{{\Omega }_p}{w}_p\delta (\mathbf{R}-{\mathbf{R}}_p)\delta (p_z-p_{p,z})\mathrm{d}{\Omega }_p$$

$$w_p(t)=⟨\delta f{⟩}_{{\Omega }_p}{\Omega }_p,$$

$$\underset{{\Omega }_p\to 0}{\lim }⟨\delta f{⟩}_{{\Omega }_p}\to \delta f$$

$${\Omega }_p=\frac{B_{\parallel ,p}^{\ast }}{B}{v}_{\perp ,p}(\pi {\kappa }_v{u}_{th}(s){)}^2{\int }_0^1\bar{J}(s)\mathrm{d}s,$$

$$E[\zeta (\tilde{Z})]=\int \zeta (z)f(z)dz,$$

$$E[W(Z)\zeta (Z)]=\int \zeta (z)\frac{f(z)}{g(z)}g(z)dz.$$

$$E[\zeta (\tilde{Z})]=E[W(Z)\zeta (Z)]=\frac{1}{N}\sum _{i=1}^{N}w(Z_i)\zeta (Z_i),$$

$$w(Z_i)=\frac{f(Z_i)}{g(Z_i)}=f(Z_i)\Omega (Z_i),$$

$${\dot{\mathbf{R}}}_{\mathbf{sp}}={\mathbf{v}}_{\parallel ,sp}+{\mathbf{v}}_{\nabla B,sp}+{\mathbf{v}}_{curvB,sp}+{\mathbf{v}}_{\nabla p,sp}+{\mathbf{v}}_{E\times B,sp}+{\mathbf{v}}_{A_{\parallel },sp},$$

$${\mathbf{v}}_{\parallel ,sp}=\frac{p_{z,sp}}{m_{sp}}\mathbf{b},$$

$${\mathbf{v}}_{\nabla B,sp}={\mu }_{sp}B\frac{1}{Z_{sp}e{B}_{\parallel ,sp}^{\ast }}\mathbf{b}\times \frac{\nabla B}{B},$$

$${\mathbf{v}}_{curvB,sp}={\left(\frac{p_{z,sp}}{m_{sp}}\right)}^2\frac{m_{sp}}{Z_{sp}e{B}_{\parallel ,sp}^{\ast }}\mathbf{b}\times \frac{\nabla B}{B},$$

$${\mathbf{v}}_{\nabla p,sp}=-{\left(\frac{p_{z,sp}}{m_{sp}}\right)}^2\frac{m_{sp}}{Z_{sp}e{B}_{\parallel ,sp}^{\ast }}\mathbf{b}\times \left(\mathbf{b}\times \frac{\nabla \times \mathbf{B}}{B}\right),$$

$${\mathbf{v}}_{E\times B,sp}=-\frac{\nabla (J_{0,sp}\Phi )\times \mathbf{b}}{B_{\parallel ,sp}^{\ast }},$$

$${\mathbf{v}}_{A_{\parallel },sp}=\frac{p_{z,sp}}{m_{sp}}\frac{J_{0,sp}{A}_{\parallel }}{B_{\parallel ,sp}^{\ast }}\mathbf{b}\times (\mathbf{b}\times (\nabla \times \mathbf{b}))-\frac{Z_{sp}e}{m_{sp}}{J}_{0,sp}{A}_{\parallel }\mathbf{b},$$

Full text

Gyrokinetics; PIC; Wave-particle interaction; Zonal flows; GAMs; 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 Brizard07 . 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 Hatzky07 ; Mishchenko17 , 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) Winsor68 ; Diamond05 ; Qiu18 or energetic-particle driven GAMs, called EGAMs Boswell06 ; Fu08 ; Horvath16 ; Ido18 . There are different kind of techniques to investigate dynamics of modes in the phase-space. Correlation techniques Howes17 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 Keiling09 . 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 Jolliet07 ; Bottino15 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 Bottino04 ; Tronko16 . 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 Ido15 ; Osakabe14 . 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 Zhang10 ; Biancalani17 ; Novikau17 ; Ehrlacher18 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 turbulence-transport processes. The GAMs are an oscillating branch of zonal flows Rosenbluth98 ; Diamond05 (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 Manz12 ; Schmitz12 ; Chen00 ; Medvedeva17 . Contrarily, the role of the GAMs Zonca08 ; Conway11 ; Kobayashi18 ; Liang18 in the turbulence suppression is still unclear and even contradictory Scott03 ; Silva17 . 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 Zonca08 and arising from the ZFs due to the magnetic curvature, or they can return the energy back to the instabilities Scott03 . Due to this complex dynamics, the role of the EGAMs in the turbulence suppression is still a subject of study Zarzoso13 ; Zarzoso17 ; Sasaki17 ; MSasaki17 ; Biancalani18 . 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 Osakabe14 . 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.

II Formulation and implementation of the Mode-Particle-Resonance diagnostic

II.1 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 Brambilla98 of electromagnetics:

[TABLE]

where $\mathcal{E}$ is the mode energy, $\mathbf{S}$ is the Poynting’s vector, $\mathbf{J}_{sp}$ is the current, produced by a species $sp$, and $\mathbf{E}$ is the electric field. By integrating over the whole real space, we get $\int\boldsymbol{\nabla}\cdot\mathbf{S}\mathop{}\!\mathrm{d}V=0$ using the Gauss’s theorem. Introducing the mode complex frequency:

[TABLE]

the mode energy evolves in time as

[TABLE]

from where we can get an expression for the mode complex frequency:

[TABLE]

Finally, using Eq. 1, we get an expression for the mode damping rate:

[TABLE]

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 $\mathcal{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 $\gamma<0$ corresponds to a positive signal $\mathcal{P}$, indicating the energy transfer from a wave to plasma particles. On the other hand, a positive rate $\gamma>0$ corresponds to the growth of the wave.

II.2 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 particle-in-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 Bottino15 . Every marker has a particular magnetic moment $\mu_{sp}=m_{sp}v_{\perp,sp}^{2}/(2B)$, a position in real space $\mathbf{R}_{sp}$, a parallel canonical momentum $p_{z,sp}$ and it is moving in a background magnetic field $\mathbf{B}=\mathbf{b}B$ with

[TABLE]

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=(\mathbf{R}_{sp},p_{z,sp},\mu_{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:

[TABLE]

which are obtained by varying a GK Lagrangian with respect to the phase-space coordinates $Z=(\mathbf{R}_{sp},p_{z,sp},\mu_{sp})$ Tronko16 ; Tronko18 . The orbits are perturbed by the field perturbation

[TABLE]

with $\Phi$ and $A_{\parallel}$ being electric and parallel magnetic potential perturbations respectively, where only $\Phi$ 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:

[TABLE]

Considering perturbations of the distribution function and of the particle orbits till the first order, one can linearize the Vlasov equation:

[TABLE]

Assuming that $f_{0,sp}$ is an equilibrium distribution function, it should be conserved along unperturbed particle trajectories $(\mathbf{\dot{R}_{0,sp}},\dot{p}_{0,z,sp})$:

[TABLE]

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:

[TABLE]

where \bigg{|}_{1} indicates that it is necessary to take derivatives along the perturbed parts of species orbits $(\mathbf{\dot{R}_{1,sp}},\dot{p}_{1,z,sp})$. Thermal species have an equilibrium distribution function in a form of the Maxwellian one:

[TABLE]

where $n_{sp}(\psi)$, $T_{sp}(\psi)$ are species density and temperature profiles along the radial coordinate $\psi$, which is the poloidal flux. A symmetric two-bumps-on-tail distribution function has been used in this work for the fast species Biancalani14 ; Zarzoso14 . This distribution assumes a flat temperature profile of the fast species:

[TABLE]

where $u_{H,sp}$, $T_{H,sp}$ are constant input parameters, which specify a shift and width of the bumps respectively.

The perturbation $\delta f$ is discretized in the $Z=(\mathbf{R},p_{z},\mu)$ 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. Jolliet07 ; Jolliet10 ; Bottino15 . A marker weight can be associated to a phase space volume $\Omega_{p}$ and correspondent averaged perturbation distribution function $\langle\delta f\rangle_{\Omega_{p}}$:

[TABLE]

Considering uniform spreading of the markers in real space and Maxwellian distribution in the velocity space, it can be shown Jolliet10 that the phase space volume $\Omega_{p}$, associated to a marker $p$, is

[TABLE]

where $\bar{J}(s)$ is the flux-surface-averaged Jacobian, $\kappa_{v}$ defines maximum value of the species parallel and perpendicular velocities, normalized to a species thermal speed $u_{th}(s)=\sqrt{T/m}$, at every radial point $s=\sqrt{\psi/\psi_{edge}}$.

The meaning of the variable $\Omega_{p}$ can be explained proceeding directly from the Monte Carlo integration Bottino15 . The expectation value of an arbitrary function $\zeta(\tilde{Z})$ is

[TABLE]

where $\tilde{Z}$ is a random variable, distributed according to the function $f$. To minimize the variance of the function $\zeta$, one can chose another distribution function $g(\tilde{Z})$, which does not vanish in the support of the distribution function $f$ (so-called importance sampling):

[TABLE]

In this case, speaking in terms of marker weights and using random variable $Z$, distributed with density $g$, the expectation value of the function $\zeta(\tilde{Z})$ is calculated as

[TABLE]

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. Bottino15 .

To clarify different terms in Eq. 10, the characteristic $\mathbf{\dot{R}_{sp}}$ can been split on several terms:

[TABLE]

where $\mathbf{b}\times(\boldsymbol{\nabla}\times\mathbf{B})/B=\boldsymbol{\nabla}p/B^{2}$ in Eq. 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. Tronko16 :

[TABLE]

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 $\Delta W_{sp}$, related to a particular velocity bin, we project the energy transfer signal to the velocity space of a particular species:

[TABLE]

where the sum $\sum_{i\in V,\Delta W_{sp}}$ is taken on all markers in the phase volume $V\Delta W_{sp}$. The gyro-averaged electric field $-\boldsymbol{\nabla}(J_{0,sp}\Phi)|_{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 $\mathbf{\dot{R}_{sp}}$ is taken into account (Eq. 31 - 34). Since the GK model, that is used in ORB5, is based on the Hamiltonian formulation Tronko16 , $p_{z,sp}$ is used as one of the velocity variables:

[TABLE]

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:

[TABLE]

In the ES case, the variables $u_{sp}$ and $v_{\parallel,sp}$ are identical $u_{sp}=v_{\parallel,sp}$, and in EM simulations with low $\beta$ they are close $u_{sp}\approx v_{\parallel,sp}$. With the rise of $\beta$, the difference between these two variables increases because of the contribution of the parallel magnetic potential $A_{\parallel}$. A proper transition from the variable $p_{z,sp}$ to the variable $v_{\parallel,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. The temperature is defined by the value of $\rho^{*}=1/205$, where $\rho^{*}=\rho_{s}/a$, with $\rho_{s}=c_{s}/\omega_{ci}$ and $c_{s}=\sqrt{T_{e}/m_{i}}$ being the sound speed and $\omega_{ci}=Z_{i}eB_{0}/m_{i}$ being the ion cyclotron frequency. The simulation has been performed with the electrostatic version of ORB5 with adiabatic electrons. Since here we are interested only in the GAM dynamics, the simulation has been done without energetic species. Non-zonal modes, i.e. modes with toroidal numbers $n\neq 0$, have been filtered out to keep only the physics of the zonal modes. Background magnetic field at the magnetic axis is $B_{0}=2.0$ T, the minor and major radii are $a_{0}=0.5$ m, $R_{0}=1.65$ m respectively. For simplicity, a radial domain $s=[0.5,1.0]$ has been simulated. The radial coordinate is $s=\sqrt{\psi/\psi_{edge}}$, where $\psi$ is the poloidal flux coordinate. The real space has been discretized with $n_{s}=300$ grid points along the radial direction, with $n_{\chi}=64$ along poloidal direction and $n_{\phi}=4$ along toroidal direction. A time step $dt[\omega^{-1}_{ci}]=10$ has been chosen, where the time is normalised to the inverse deuterium cyclotron frequency $\omega_{ci}$. The number of the ion markers is $N_{i}=10^{8}$.

To simulate the GAM dynamics, so-called Rosenbluth-Hinton test Rosenbluth98 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 $\mathcal{P}(v_{\parallel},\mu,t)$ (Eq. 38) as a function of the velocity variables $(v_{\parallel},\mu)$ 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:

[TABLE]

where $\omega_{GAM}$ is the GAM frequency, that can be found directly from the radial zonal electric field $\overline{E}_{r}$. Since the perturbation of the plasma distribution function related to the GAM dynamics can have higher poloidal modes $m\geq 1$, the GAM-particle interaction can be observed at smaller parallel velocities as well

[TABLE]

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_{\parallel,res}$. By integrating the signal in the whole velocity domain, one gets the time evolution of $\mathcal{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). The reason is that the GAM energy periodically oscillates between field and plasma components. Because of that, the total GAM energy can be estimated as an interpolation of the maxima of the ES field energy.

Eq. 6 involves an integration in time. Varying and choosing different time intervals, one can estimate an errorbar of the GAM damping rate by building a distribution (or histogram) of the damping rate values. Every chosen time interval has to contain a whole number of GAM periods. The result histogram can be fitted with the normal distribution function, that gives a mean value of the damping rate and a $95\%$ confidence interval as $3\sigma$, where $\sigma$ is the standard deviation, found from the distribution function. One can use a different distribution function to take into account a non-zero skewness of a histogram. For example, the Generalized-Extreme-Value (Ref.  Kotz00 ) distribution function has been used as well. Both distribution functions give quite close results for the mean value and the standard deviation of the GAM damping rate. A result value of the GAM damping rate, found from the MPR diagnostic is the following:

[TABLE]

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 $\overline{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 $\overline{E}_{r}$, that can be seen from the left plot of Fig. 1. First of all, a zero-frequency component of $\overline{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 $\overline{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

[TABLE]

is used as a test one, which is fitted to the time evolution of the $\overline{E}_{r}(s=0.74)$. This method has been used previously in Ref. Novikau17 to study the influence of the drift-kinetic electrons on the GAM dynamics in linear global GK simulations. But here, as in the MPR diagnostic, an opportunity to estimate errorbars of the frequency and especially of the damping (or growth) rate by varying time intervals has been added as well.

Finally, the GAM frequency and damping rate, found using the non-linear fitting of $\overline{E}_{r}(s=0.74)$ to the test function Eq. 44, are the following:

[TABLE]

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. Zonca96 ; Zonca08 ) is

[TABLE]

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:

[TABLE]

with a radial wavenumber $k$ and a complex frequency $\hat{\omega}=\omega^{D}+i\gamma^{D}$, that verifies the dispersion relation Eq. 47. The corresponding perturbation of the deuterium distribution function is:

[TABLE]

where $c$ is the light speed, $\omega_{tr}=v_{\parallel}/(qR_{0})$ is the passing frequency, $\theta_{p}$ is the poloidal angle in a simplified circular geometry, $F_{0}$ being the deuterium equilibrium distribution function:

[TABLE]

and $v_{d}$ being the amplitude of the radial drift composed by the curvature drift and grad-B drift:

[TABLE]

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:

[TABLE]

where $\mathbf{B}^{*}_{sp}$ and $B^{*}_{\parallel,sp}$ are defined in Eq. 8 and 9 respectively. Considering a low-pressure plasma ($\mathbf{J}_{0}\times\mathbf{B}\ll 1$ with a plasma current $\mathbf{J}_{0}$) in a circular plasma cross-section with a curvature $\boldsymbol{\kappa}$, one gets the following simplifications:

[TABLE]

Applying the introduced approximations, we get that

[TABLE]

using which together with Eq. 55 and by putting everything to Eq. 37, an expression for the energy transfer signal can be derived:

[TABLE]

Evaluating the velocity integrals $I_{1}$ and $I_{2}$, one gets the expression

[TABLE]

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)$:

[TABLE]

As a result, we have the following complex expression for the plasma-field energy exchange signal:

[TABLE]

In the GK model of ORB5, the ES field energy is given (Ref. Bottino15 ) by:

[TABLE]

with $\Phi$ 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

[TABLE]

Finally, taking into account only the perturbation of the zonal radial electric field $|\nabla_{\perp}\Phi|^{2}\approx\overline{E}_{r}^{2}$ and since $\int\mathop{}\!\mathrm{d}WF_{0}=1$, we get

[TABLE]

Since $\int\mathop{}\!\mathrm{d}V\overline{E}_{r}^{2}\sin\theta_{p}/\int\mathop{}\!\mathrm{d}V\overline{E}_{r}^{2}=1/2$ and applying Eq. 68, 71 to Eq. 6, the consistency between the GAM dispersion relation and the MPR method can be proved:

[TABLE]

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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

In this case, the analytical energy transfer signal takes the following form:

[TABLE]

If we take the GAM frequency (Eq. 45) and damping rate (Eq. 43 or 46) from the simulation, described in Sec. III, we can see that Eq. 73 and especially 76 give a quite similar behaviour in time (Fig. 4) with the same frequency of $\mathcal{P}$ as the one from the numerical simulation. And the frequency of $\mathcal{P}$ signal is double of that of $\overline{E}_{r}$. The amplitude modulation, that is observed in the numerical energy transfer signal, can be explained by the constant component of the zonal radial electric field, that is emphasized in Eq. 76, and as one can see from the comparison of the blue and green lines in Fig. 4.

V Application to EGAMs in AUG shot #31213

V.1 Equilibrium and definition of the numerical simulation

The AUG shot #31213 at time 0.84 s has been selected within the Non-Linear Energetic-particle Dynamics (NLED) Eurofusion enabling research project LauberSitee ; Horvath16 . It has been chosen to study the effect of the energetic particles (EPs) on the dynamics of EGAMs. That is why, in these simulations we have three species: gyro-kinetic thermal deuterium, gyro-kinetic energetic (fast) deuterium, and thermal electrons, either adiabatic (AE) or drift-kinetic (KE). The linear dynamics of EGAMs in this NLED-AUG case has been recently investigated with the gyrokinetic codes GENE and ORB5 by considering adiabatic electrons DiSiena_NF2018 . Here, we extend the previous study by investigating the effect of kinetic electrons and describing the contribution of the resonances of all species in phase space. The simulation with the AE is performed in the electrostatic limit, while the simulation with the KE has been done including dynamics of the magnetic potential perturbation as well. In this latter case the pullback method Mishchenko17 has been used for the mitigation of the cancellation problem in EM simulations Chen01 ; Hatzky07 . 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 $R_{0}=1.62$ m, $a=0.482$ m respectively.

The real space of the system has been discretized using the following parameters: $n_{s}=256$, $n_{\chi}=256$, $n_{\phi}=32$. In the ES simulation the time grid has a step $dt[\omega_{ci}^{-1}]=20$ with $N_{i}=5\cdot 10^{8}$ being a number of markers for the thermal ions, and $N_{f}=5\cdot 10^{8}$ for the fast ions. In the EM case, the time step and number of markers have been changed: $dt[\omega_{ci}^{-1}]=5$, $N_{i}=N_{f}=10^{8}$, $N_{e}=4\cdot 10^{8}$. Such a high number of markers is needed to provide at least several thousands of numerical markers in every velocity bin, 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 $\beta_{e}=\langle n_{e}\rangle T_{e}/(B_{0}^{2}/(2\mu_{0}))=2.7\cdot 10^{-4}$, where $\langle n_{e}\rangle$ is the electron density, averaged in a tokamak volume, $\mu_{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 by the expression Eq. 20 with $u_{H,f}=8$ and $T_{H,f}=1$. The ORB5 simulation with such parameters of the fast species results in one of the biggest EGAM growth rate for the given plasma configuration.

V.2 Numerical investigation of the wave-particle resonances in the EGAM dynamics

First of all, one can notice from Fig. 7 that the radial structure of the EGAMs slightly changes when the dynamics of the drift-kinetic electrons is switched on. The position of the crest in the EGAM radial structure shifts inwards from around $s=0.48$ to $s=0.40$.

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 $\overline{E}_{r}$. They appear to be consistent within the error bars of the measurements:

[TABLE]

From here on, only the radial point $s=0.40$ is considered in the following calculations. Consistency between the EGAM growth rate, calculated directly from $\overline{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:

[TABLE]

The consistency between both methods is observed in the EM case with KE as well:

[TABLE]

From Eq. 78 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:

[TABLE]

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. Zarzoso14 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:

[TABLE]

From the above equations one can see that the drive on the fast particles does not change (Eq. 86 and 89). On the other hand, there is a significant increase of the EGAM damping on the thermal deuterium plasma (Eq. 85 and 87). Since this increase is comparable with the electron contribution (Eq. 88), one can not claim from these results that the decrease of the EGAM growth rate occurs only directly due to the additional damping on electrons. But apart from the direct damping, inclusion of drift-kinetic electrons changes the position of the EGAM crest (Fig. 7) and slightly changes the EGAM frequency. These changes can lead to the increase of the EGAM damping on thermal deuterium. On the other hand, the corresponding errorbars of the ion contributions become higher in the simulation with KE in comparison with AE. Nevertheless, it is clearly shown here that in the experimentally relevant plasma conditions the inclusion of the drift-kinetic electrons significantly decreases the EGAM growth rate of about a factor 2.

We now want to investigate the role of the different resonances in phase space. In Fig. 8(b) 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:

[TABLE]

where $\epsilon$ 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 Zhang10 . We can separate two velocity domains $e11$ and $e21$, shown in Fig. 8(b). 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. 8(a)), 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:

[TABLE]

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. 8(c)):

[TABLE]

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. 8(c)), 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. 8(d)):

[TABLE]

One can see that there is an EGAM damping even on the energetic particles (Eq. 103 and 106). 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 verify some of the results, obtained in Sec. V, we have performed a comparison with the gyrokinetic GENE code, that has a similar diagnostic. The Gyrokinetic Electromagnetic Numerical Experiment (GENE) GENE is an Eulerian code, which solves the Vlasov-Maxwell system of coupled equations on the phase-space grid $(\mathbf{R},v_{\parallel},\mu)$ at each time step. Here, $\mathbf{R}$ denotes the gyrocenter position, $v_{\parallel}$ the velocity component parallel to the magnetic field and $\mu$ the magnetic moment. The gyrokinetic description employs an approach based on the study of a distribution function $f_{s}(\mathbf{R},v_{\parallel},\mu)$ for each plasma species $(s)$, which contrarily as it is done in a particle-in-cell code as ORB5, is not discretized with markers. The distribution function is split, accordingly to the so-called $\delta$-f approach, into a background component $f_{0,s}$ and in a small fluctuating part $f_{1,s}$, i.e. $f_{s}=f_{0,s}+f_{1,s}$. The equilibrium distribution function $f_{0,s}$ is usually modelled with a Maxwellian distribution. However, recently, this assumption has been relaxed and more flexible equilibrium distributions can be considered Varenna ; PoP ; EGAM . In particular, different analytic choices, e.g. slowing down, bi-Maxwellian and bump-on-tail, as well as numerical distributions as obtained from numerical models are supported. While the equilibrium distributions are considered time independent on the turbulent time scales, their perturbed components evolve in time accordingly to the Vlasov equation, which in the linear and electrostatic limit employed throughout this paper reads as (for more details one can see Ref. Tobias ; Dannert_1 ; Dannert_2 )

[TABLE]

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\bar{\phi}_{1}/(B_{0}T_{0,s})\partial f_{0,s}/\partial\mu$. 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, $\bar{\phi}_{1}$ denotes the gyro-averaged electrostatic potential, $\mathcal{J}$ the phase-space jacobian, $\mathcal{K}_{x}\sim-\partial_{y}B_{0}-\partial zB_{0}$ and $\mathcal{K}_{y}\sim\partial_{x}B_{0}-\partial_{z}B_{0}$, respectively, the radial and bi-normal curvature terms and $\mathcal{C}^{2}=\textbf{B}_{0}\cdot\textbf{B}_{0}$. Finally, the Poisson brackets are defined as

[TABLE]

Eq. 107 needs to be solved self-consistently with the Poisson field equation. The full plasma dynamic can be investigated in GENE either in a flux-tube (local assumption) GENE or in a full-global radial domain Tobias . The local approximation allows the radial direction to be Fourier transformed by assuming periodic boundary conditions. GENE is able to study the contribution of each plasma species to the overall more unstable mode-dynamic through the study of the time evolution of the potential energy of the system $E_{w}$ Banon_Navarro_PoP2011 ; Navarro_PRL2011 . It is defined only in Fourier space (only in the local flux-tube limit) for each wave vector $\textbf{k}=(k_{x},k_{y})$ as follows

[TABLE]

Here, the bracket represents the field-aligned z-average, namely

[TABLE]

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

[TABLE]

From the energy relation of Eq. 111 it is to compute the more unstable linear growth rate $\gamma$ through the time variation of the potential energy, as shown in details in Ref DiSiena_NF2018 ; Hatzky_PoP2002 ; Banon_Navarro_PoP2011 ; Manas_PoP2015 , by the relation

[TABLE]

Eq. 112 allows us to distinguish between the contribution of each species to the total growth rate, by removing the sum over all the species and studying each term separately. Positive (negative) values of $\partial E_{k,s}/\partial t$ indicate that the plasma species considered is giving (taking) energy to (from) the electrostatic field component with a consequent growth (damping) of the mode. Moreover, by studying $\gamma_{s}$ in phase-space, i.e. $(v_{\parallel},\mu)$ for each plasma species, velocity resonances, which are the main drive term of the EGAMs studied in this paper, can be investigated in details.

The same AUG shot has been simulated in GENE in case with adiabatic electrons (one can see also Ref. DiSiena_NF2018 ), using the flux-tube version of the code at $s=0.5$. In Fig. 9 one can see that both ORB5 and GENE give the same positions of the resonances of the EGAM - fast deuterium plasma interaction. According to chosen parameters of the fast deuterium distribution function, peaks of the energetic bumps are located at $|v_{\parallel}|=8$. The opposite signs of the same resonances in two codes is explained by the fact that different signs are used in the MPR diagnostic in ORB5 (Eq. 6) and in the corresponding diagnostic in GENE (Eq. 112). On the other hand, one of the possible explanation of the opposite nature of the resonance asymmetry for positive and negative parallel velocities is the opposite direction of the background magnetic field, used by ORB5 and GENE.

As a benchmark, comparison of the EGAM frequency and total growth rate has been done as well. One can see that both codes give the same values of the mode frequency:

[TABLE]

On the other hand, there is $18\%$ consistency between the codes for the EGAM total growth rate:

[TABLE]

and for the contributions of different plasma species to the mode dynamics:

[TABLE]

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 Zonca96 ; Zonca08 , 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_{\parallel,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. 118) and energetic deuterium (Eq. 120) 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 Berk97 ; Berk99 ; Berk06 ; Berk10 ; HWang13 ; BiancalaniJPP17 , 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 $\beta$, by performing a proper transition from the variable $p_{z,sp}$ to the velocity variable $v_{\parallel,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, resulting in collisionless damping of the turbulent fluctuations. More precisely, a particular challenge in tokamak plasma and plasma physics in general is to identify the physical mechanisms by which the EM field and plasma flow fluctuations are damped and how their energy is converted to plasma heat, or some other energization of particles. For example, in astrophysical plasmas dissipation of the turbulence energy through the Landau damping of the Alfvén waves can take place Li16 . It would be interesting to investigate the influence of the plasma $\beta$ on the energy channeling Parashar18 , especially on the contributions of different species in the plasma heating by EGAMs and Alfvén waves. There are also physical phenomena, which are specific to the space plasma, such as particle acceleration by the magnetic energy released during collisions of the magnetic islands in solar and heliospheric environments Du18 . It might be interesting to investigate role of such processes in tokamak plasmas as well.

Acknowledgement

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

Simulations, presented in this work, have been performed on the CINECA Marconi supercomputer within the framework of the OrbZONE and ORBFAST projects.

Stimulating discussions with X. Garbet and V. Grandgirard on the role of kinetic electrons in the GAM dynamics are kindly acknowledged. One of the authors, I. Novikau, would like to thank F. Vannini for useful conversations. One of the authors, A. Biancalani, also wishes to acknowledge stimulating discussions with F. Zonca and Z. Qiu on the GAM/EGAM analytical theory.