###### Abstract

In the present paper, the evolution of the Alfvén modes is studied in a
realistic ASDEX Upgrade equilibrium by analyzing the results of simulations
with the global, electromagnetic, gyrokinetic particle-in-cell code ORB5. The
energetic particles are modelled both via the newly implemented isotropic
slowing-down and with Maxwellian distribution functions. The comparison of the
numerical results shows that modelling the energetic particles with the
equivalent Maxwellian rather than with the slowing-down, does not affect the
frequency of the driven Alfvén mode, while its growth rate appears to be
underestimated with a quantitative difference as large as almost $30\,\%$.
Additionally the choice of the isotropic slowing-down allows a better
description of the nonlinear modification of the dominant Alfvén mode
frequency, while an equivalent Maxwellian underestimates it. A good comparison
with the experimental spectrogram is found.

Gyrokinetic modelling of the Alfvén mode activity in ASDEX Upgrade with an
isotropic slowing-down fast-particle distribution
F. Vannini1, A. Biancalani4,1, A. Bottino1, T. Hayward-Schneider1, P. Lauber1,
A. Mishchenko2, E. Poli1, B. Rettino1, G. Vlad3, X. Wang1 and the ASDEX
Upgrade team5.

1Max-Planck-Institut für Plasmaphysik, 85748 Garching, Germany

2Max-Planck-Institut für Plasmaphysik, 17491 Greifswald, Germany

3ENEA, Fusion and Nuclear Safety Department, 00044 Frascati (Roma), Italy

4 Léonard de Vinci Pôle Universitaire, Research Center, Paris la Défense,
France

5 See author list of H. Meyer et al. 2019 Nucl. Fusion 59 112014
francesco.vannini@ ipp.mpg.de

##  1 Introduction

Next-generation fusion machines like ITER [1] and DEMO [2] will be
characterized by a significant population of energetic particles (EPs). With
this term we refer to fast ions in general that are going to be present in the
confining machines as fusion products ($\alpha$-particles, ${}^{4}_{2}He$) or
as products of auxiliary heating sources like neutral beam injection (NBI) or
ion cyclotron resonance heating (ICRH). The typical EP velocity $v_{EP}$ is
intermediate between the thermal bulk ions velocity $v_{th,i}$ and the thermal
electron velocity $v_{th,e}$. EPs’ gyroradius $\rho_{EP}$ is instead much
larger than the bulk plasma species gyroradii ($\rho_{i}$ and $\rho_{e}$) [3]:

$v_{th,i}\ll v_{EP}\ll
v_{th,e}\quad,\quad\rho_{EP}\gg\rho_{i}\gg\rho_{e}\quad.$ (1)

Typically in Tokamak machines the EPs’ characteristic dynamical frequencies
$\omega$ associated with their guiding-center motion (transit, bounce and
precessional) fall inside the magnetohydrodynamic (MHD) regime $\omega\sim
10^{-2}\,\omega_{ci}$, with $\omega_{ci}$ the ion cyclotron frequency. Because
of this EPs can resonate with MHD instabilities driving them unstable. Among
the MHD instabilities, the most detrimental and easily excited are the nearly
incompressible shear Alfvén wave (SAWs). These are transverse electromagnetic
waves that propagate along the magnetic field lines with group velocity equal
to the Alfvén speed $v_{A}$. Their dispersion relation is [4, 5]:

$\omega_{SAW}=k_{\parallel}v_{A}\quad,\quad
k_{\parallel}=\vec{k}\cdot\hat{\bm{b}}\quad,\quad
v_{A}=\frac{B_{0}}{\sqrt{4\pi\rho_{m,0}}}\quad,$ (2)

with $\hat{\bm{b}}=\bm{B}_{0}/B_{0}$ a unit vector pointing in the direction
of the background magnetic field $\bm{B}_{0}$ with modulus $B_{0}$,
$\rho_{m,0}$ the plasma mass density and $\bm{k}$ the wave-number of the
perturbation. In a straight cylinder [6, 7, 8]:

$k_{\parallel}=\frac{1}{R_{0}}\left(n-\frac{m}{q(r)}\right)\quad,$ (3)

with $q(r)$ the safety factor profile, $R_{0}$ the major radius of the Tokamak
and $m$ and $n$ the poloidal and toroidal mode numbers respectively. SAWs are
characterized by a frequency spectrum that varies continuously across the
radial domain because of the radial dependence contained in both the Alfvén
speed and the safety factor profile.
Different kinds of Alfvén instabilities, that we call for simplicity Alfvén
modes (AMs) [9, 10, 11], can be present in Tokamak machines. These include
both energetic-particle continuum modes (EPMs) [12] and Alfvén eigenmodes
(AEs) [10].
EPMs represent forced oscillations of the SAW continuum that arise when the EP
pressure is comparable to that of the bulk plasma. They emerge as discrete
oscillations at the frequency of the continuum where the wave-EP power
exchange is maximized, above the threshold of the continuum damping.
The AEs, on the other hand, are normal modes of the bulk plasma that can be
classified into two types. The first type includes modes that arise in
correspondence of the radial position $r_{0}$ where different branches of the
continuum, see eq. (3), cross. There two counterpropagating waves interfere
destructively, opening a gap in the continuum where modes reside. An important
example of gap modes is represented by the toroidal Alfvén eigenmodes (TAE
[13, 14, 15, 5]) that arise because of the coupling between modes with close
poloidal harmonics: $(m_{0},n)$ and $(m_{0}+1,n)$. The second type of AEs
arise in correspondence of an extremum of the continuum spectrum. There
$\partial\omega/\partial r$ vanishes and an effective potential well is
established that traps the wave. An example of these kind of modes is
represented by the reversed shear Alfvén eigenmodes (RSAE) [16, 17] that are
present in correspondence of a minimum in the safety factor profile.
The interest connected with the study of the driven AMs relies on their
potentially detrimental role. In fact the driven AMs can interact with the EPs
redistributing them in phase-space, enhancing the EP transport. In this way
EPs can be expelled before they can thermalize with the bulk plasma leading to
a less effective heating [18, 19, 20]. Additionally AMs are held responsible
for the presence of the so-called abrupt-large-events (ALE [21, 22, 23, 24]).
These are violent and rapid variations of the fluctuating magnetic field in
correspondence of which an intense migration of EPs from the core to the
periphery is observed, representing a possible threat for the safety of the
machine. Therefore, the study of the AM dynamics is of primary importance both
for the safety of the machine and because only a small fraction of EP losses
can be tolerated to achieve ignition. That is why is mandatory to gain the
necessary insight in the AM dynamics to become predictive about the scenarios
that will be met in ITER and DEMO.
In this paper we want to contribute to this task studying the AM dynamics with
the code ORB5 [25] in numerical simulations where the so-called NLED-AUG case
[26] is considered. With this term we refer to the plasma conditions of the
ASDEX Upgrade discharge $\\#31213$ at time $t=0.84\,$s. The uniqueness of this
scenario relies on the fact that it has been obtained tuning the plasma
parameters so that the EPs injected through a neutral beam (NB) have an energy
with respect to the bulk plasma temperature $\mathcal{E}_{EP}/T_{Bulk}\sim
10^{2}$ and an induced fast-ion $\beta$ comparable to that of the bulk plasma.
In such a way there is a strong Alfvén activity, the stabilizing effects of
the bulk plasma being minimized. In particular, TAE-EPM bursts are observed in
correspondence of which energetic particle driven geodesic acoustic modes
(EGAMs [27]) appear exhibiting the typical chirping.
Studies of the AM activity with the numerical tool ORB5 [25] using the NLED-
AUG scenario have already been presented in Ref. [28, 29]. There the EPs have
been modelled, respectively, with Maxwellian distribution functions and with
bi-shifted Maxwellians (also known as double-bump-on-tail). These were the two
kinds of equilibrium distribution functions available at that time in ORB5 and
the choice of one over the other was motivated by the kind of physics we were
interested to reproduce. In particular the double-bump-on-tail was chosen
because of the need to have an anisotropy in velocity space to drive unstable
an EGAM and study its interaction with the AMs. In the works already published
a good qualitative comparison with the experiments has been obtained. In this
work we want to make a further step showing that, through the slowing-down
distribution function [30] newly implemented in ORB5 we are able to go closer
to the experimental conditions obtaining an even better quantitative
comparison with the experiment.
The present paper is structured as follows. In Sec. 2 the main features and
the model of ORB5 are described. In Sec. 3 the NLED-AUG case is briefly
described. In Sec. 4 an analytical derivation of the slowing-down distribution
function together with details about its implementation in ORB5 are provided.
In Sec. 5 the results of numerical simulations obtained taking into account
the NLED-AUG case are described. Finally in Sec. 6 the conclusions of this
paper are presented.

##  2 The numerical model

ORB5 [25] is a nonlinear, global, electromagnetic, gyrokinetic, particle-in-
cell (PIC) code that can take into account collisions and sources (neglected
in this work).
The code uses a system of straight-field line coordinates:
$(r,\theta^{*},\varphi)$. As radial coordinate the code takes into account the
square root of the poloidal flux $\psi$ normalized at its value at the edge
$\psi_{0}$: $r=\sqrt{\psi/\psi_{0}}$ with $0\leq r\leq 1$ that labels the
magnetic surfaces. $\varphi$ is the toroidal angular coordinate and
$\theta^{*}$ is the poloidal magnetic angle:

$\theta^{*}=\frac{1}{2\,q(r)}\int_{0}^{2\pi}d\theta^{{}^{\prime}}\frac{\bm{B}_{0}\cdot\nabla\varphi}{\bm{B}_{0}\cdot\nabla\theta^{{}^{\prime}}}\quad,\quad
0\leq\theta^{*}\leq 2\pi$ (4)

with $\theta^{{}^{\prime}}$ the geometric poloidal angle and $q(r)$ the safety
factor profile.
In ORB5 all the physical quantities are normalized with respect to reference
parameters. The masses of the considered particle species $m_{sp}$ are given
in units of the mass of the main bulk ion species $m_{i}$, the velocities are
given in units of the sound velocity $c_{s}$, the lengths are given in units
of the sound Larmor radius $\rho_{s}$ and the time is given in units of the
inverse bulk ion cyclotron frequency:

$\omega_{ci}=\frac{q_{i}B_{0}}{m_{i}c}\quad,\quad q_{i}=e\,Z_{i}\quad$ (5)

with $e$ the electron charge in absolute value and $Z_{i}$ the atomic number
of the main ion species.
The total distribution function of the $sp$-particle species $f_{sp}$ is
divided into a time-independent part (equilibrium or background distribution
function, $F_{0,sp}$) and a time dependent component $\delta f_{sp}$, so that
$f_{sp}=F_{0,sp}+\epsilon_{\delta}\delta f_{sp}$ with $\epsilon_{\delta}$ a
small parameter. Only the time-dependent component is discretized as will be
later discussed. The gyrokinetic Vlasov equation for the perturbed (time-
dependent) distribution function is:

$\frac{d}{dt}\delta f_{sp}=-\dot{\vec{X}}\cdot\frac{\partial
F_{0,sp}}{\partial\vec{X}}\bigg{\rvert}_{\mathcal{E},\,v_{\parallel}}-\dot{\mathcal{E}}\frac{\partial
F_{0,sp}}{\partial\mathcal{E}}\bigg{\rvert}_{\bm{X},\,v_{\parallel}}-\dot{v}_{\parallel}\frac{\partial
F_{0,sp}}{\partial
v_{\parallel}}\bigg{\rvert}_{\bm{X},\,\mathcal{E}}\quad,\quad\mathcal{E}=\frac{v_{\parallel}^{2}}{2}+\mu
B\quad,\quad\mu=\frac{v_{\perp}^{2}}{2B}$ (6)

where $\bm{X}$ is the gyrocenter position, while $v_{\parallel}$ and
$v_{\perp}$ are, respectively, the parallel and perpendicular components of
the particle velocity. The equation of motions in mixed-variable formulation
[31] of the gyrocenter characteristics of the particle species are [32]:

$\vec{\dot{X}}=v_{\parallel}\hat{\bm{b}}^{*}+\frac{1}{q_{sp}B^{*}_{\parallel}}\hat{\bm{b}}\times\mu\nabla
B+\epsilon_{\delta}\left[\frac{\hat{\vec{b}}}{B_{\parallel}^{*}}\crossproduct\nabla\langle\delta\phi-
v_{\parallel}\delta A_{\parallel}^{h}-v_{\parallel}\delta
A_{\parallel}^{s}\rangle_{\alpha}-\frac{q_{sp}}{m_{sp}}\langle\delta
A_{\parallel}^{h}\rangle_{\alpha}\hat{\vec{b}}^{*}\right]$ (7)
$\dot{v}_{\parallel}=-\frac{\mu}{m_{sp}}\hat{\bm{b}}^{*}\cdot
B-\epsilon_{\delta}\bigg{\\{}\mu\frac{\vec{\hat{b}}\crossproduct\nabla
B}{B_{\parallel}^{*}}\cdot\nabla\langle\delta
A_{\parallel}^{s}\rangle_{\alpha}+\frac{q_{sp}}{m_{sp}}\left[\hat{\vec{b}}^{*}\cdot\nabla\langle\delta\phi-
v_{\parallel}\delta A_{\parallel}^{h}\rangle_{\alpha}+\frac{\partial}{\partial
t}\langle\delta A_{\parallel}^{s}\rangle_{\alpha}\right]\bigg{\\}}$ (8)
$\dot{\mathcal{E}}=v_{\parallel}\dot{v}_{\parallel}+\mu\nabla
B\cdot\dot{\vec{X}}\quad,\quad\dot{\mu}=0\quad.$ (9)

In eqs. 7, 8 and 9, $m_{sp}$ and $q_{sp}=e\,Z_{sp}$ are respectively the mass
and charge of the $sp$-particle species with atomic number $Z_{sp}$. Note that
for the electron species $Z_{e}=-1$. Still in eqs. 7, 8 and 9 the
gyroaveraging operator appears:

$\langle
F\rangle_{\alpha}=\frac{1}{2\pi}\int_{0}^{2\pi}d\alpha\,F(\bm{X}+\bm{\rho}_{0})=\frac{1}{2\pi}\int_{0}^{2\pi}d\alpha\int
d^{3}\bm{r}\,F(\bm{r})\,\delta^{3}(\bm{X}+\bm{\rho}_{0}-\bm{r})=\mathcal{J}^{gc}_{0}(F)$
(10)

that removes the fast gyroangle $\alpha$ dependence into a general quantity
$F$, with $\bm{\rho}_{0}(\bm{X},\alpha)$ the gyroradius. $\delta\phi$ is the
perturbed scalar potential and $\delta A_{\parallel}$ is the perturbed
magnetic parallel potential decomposed into symplectic and Hamiltonian parts:
$\delta A_{\parallel}=\delta A_{\parallel}^{h}+\delta A_{\parallel}^{s}$.
$B^{*}_{\parallel}$ is the parallel component of the symplectic magnetic field
$\bm{B}^{*}$, linked to the symplectic magnetic potential $\bm{A}^{*}$ by the
following relation:

$\bm{A}^{*}=\bm{A}+\frac{m_{sp}}{q_{sp}}v_{\parallel}\hat{\bm{b}}\quad,\quad\bm{B}^{*}=\nabla\times\bm{A}^{*}\quad,\quad\hat{\bm{b}}^{*}=\frac{\bm{B}^{*}}{B^{*}_{\parallel}}=\frac{\bm{B}^{*}}{\hat{\bm{b}}\cdot\bm{B}^{*}}\quad.$
(11)

In ORB5 only the perpendicular component of the perturbed magnetic field
$\delta\bm{B}_{\perp}$ is implemented:
$\delta\bm{B}_{\perp}=\hat{\bm{b}}\times\nabla\delta A_{\parallel}$. In eq.
(11) $\bm{A}$ is the background magnetic potential:
$\bm{B}_{0}=\nabla\times\bm{A}$ and $\hat{\bm{b}}$ a unit vector pointing in
the direction of the background magnetic field. The characteristic equations
eqs. 7, 8 and 9 are coupled to the field equations. These are the gyrokinetic
quasi-neutrality equation:

$-\nabla\cdot\left[\left(\sum_{sp=i,f}\frac{q_{sp}^{2}n_{sp}}{T_{sp}}\rho_{sp}^{2}\right)\nabla_{\perp}\delta\phi\right]=\sum_{sp=i,e,f}q_{sp}\,\delta
n_{sp}\quad,\quad\delta n_{sp}=\int dW\langle\delta f_{sp}\rangle$ (12)

the parallel Ampère’s law

$\left(\sum_{sp=i,e,f}\frac{\beta_{sp}}{\rho_{sp}^{2}}-\nabla_{\perp}^{2}\right)\delta
A_{\parallel}^{h}=\mu_{0}\sum_{sp=i,e,f}\delta
j_{\parallel,sp}+\nabla^{2}_{\perp}\delta A_{\parallel}^{s}\quad,\quad\delta
j_{\parallel,sp}=q_{sp}\int dW\,v_{\parallel}\,\langle\delta f_{s}\rangle$
(13)

and the ideal Ohm’s law:

$\frac{\partial}{\partial t}\delta
A_{\parallel}^{s}+\hat{\vec{b}}\cdot\nabla\delta\phi=0$ (14)

where:

$n_{sp}=\int
dWF_{0,sp}\quad,\quad\beta_{sp}=\mu_{0}\frac{n_{sp}T_{sp}}{B_{0}^{2}}$ (15)

where $dW=B^{*}_{\parallel}\,dv_{\parallel}\,d\mu\,d\alpha$. Equations 6, 7,
8, 9, 12, 13 and 14 constitute the gyrokinetic Vlasov-Maxwell system of
equations solved by ORB5.
The discretization is achieved sampling the phase-space with a set of super-
particles called markers. The $k$-th marker of the $sp$-particle species is
associated to a weight $\omega_{sp,k}(t)$. This is a time-dependent quantity
that represents at the time $t$ the variation of the number of physical
particles contained in a small volume in phase-space $\Omega_{k}$, centered
around the $k$-th marker:

$\delta N_{sp}^{physical\,particles}=\int_{\Omega_{k}}\,d^{6}\bm{Z}\,\delta
f_{sp}(\bm{Z}(t),t)=\omega_{sp,k}(t)\quad.$ (16)

Each marker is pushed along its orbit slowing the characteristics equations
through a Runge-Kutta method at $4$th order.
The perturbed fields $\Psi=\\{\delta\phi,\delta A_{\parallel}\\}$ are
discretized using the finite-elements Galerkin approximation. The fields are
represented as linear combinations of finite dimensional function space
$\Lambda_{\mu}(\bm{x})$:

$\Psi(\bm{x},t)=\sum_{\mu}\Psi_{\mu}(t)\Lambda_{\mu}(\bm{x})\quad$ (17)

where the basis $\Lambda_{\mu}(\bm{x})$ are tensor product of $1$D polynomials
(B-splines) of degree $p=1,2,3$:

$\Lambda_{\mu}(\bm{x})=\Lambda_{j}^{p}(r)\Lambda_{k}^{p}(\theta^{*})\Lambda_{l}^{p}(\varphi)\quad\text{with}\quad\mu=(j,k,l)\quad.$
(18)

In the present work $p=3$. The number of B-splines in each direction is
connected to the number of grid points (knots)
$(N_{r},N_{\theta^{*}},N_{\varphi})$ associated to the B-splines and to the
chosen degree $p$ of the polynomials. Through the decomposition expressed in
eq. (17) the fields equations become a set of linear equations to which is
applied the double discrete Fourier transformation $\mathcal{F}$ in both the
poloidal and toroidal directions. Inverting the obtained system of equations,
it is then possible to calculate the Fourier coefficients of the
perturbations, to which a Fourier filter is applied. This is given by a field
aligned filter:

$m\in\underbrace{[n\,q(r)-\Delta m,n\,q(r)+\Delta m]}_{\text{field aligned
filter}}\quad,\quad n\in[nfilt1,nfilt2]\quad.$ (19)

In eq. (19) $nfilt1$, $nfilt2$ are the range boundaries of the rectangular
filter for the toroidal coefficients. $\Delta m$, instead, is the width of the
field aligned filter.
In Sec. 5 results of ORB5 simulations will be discussed. There we will
indicate some important parameters used in the simulations. We will specify:
the number of grid points associated to the B-splines
$(N_{r},N_{\theta^{*}},N_{\varphi})$, the number of markers used for every
particle species $nptot_{sp}$, the time step in use $\Delta t$ and the width
of the field aligned filter $\Delta\,m$. In this work only the AMs in the
NLED-AUG scenario are investigated by retaining, as in previous works [28,
33], only the modes with toroidal mode number equal to one: $nfilt1=nfilt2=1$.
Additionally only the EPs will be allowed to redistribute in phase-space,
following their full trajectories. The bulk plasma species instead will follow
their unperturbed trajectories. This means that in eqs. 7, 8 and 9 the small
parameter is taken, with abuse of notation, equal to one only for the EPs
while for the electrons and the bulk ions species it is identically zero.
In the following sections we will also specify if a drift-kinetic model is
considered and finite-Larmor-radius (FLR) are neglected, or if a gyrokinetic
model is considered and FLR effects are taken into account.

##  3 NLED-AUG scenario

As stated in the introduction, with the term NLED-AUG case [34] we refer to
the plasma conditions present at $t=0.84\,$s in the discharge number $31213$
<EMAIL_ADDRESS>performed in the Tokamak ASDEX Upgrade (AUG). EPs, in this
experimental case, are due to a NB launched with an injection angle of
$7.13^{\circ}$ with respect to the horizontal plane. This experimental case
exhibits, as anticipated in Sec. 1, plasma parameters previously unexplored in
ASDEX Upgrade. In fact, the NB-induced fast-ion $\beta$ is comparable to that
of the bulk plasma and EPs’ injection energy is approximately $100$ times
higher than the bulk species temperatures:

$\beta_{EP}/\beta_{Bulk}\sim 1\quad,\quad\mathcal{E}_{EP}/T_{Bulk}\approx
93\,keV/1\,keV\sim 10^{2}\quad.$ (20)

By doing so, the realistic ratios of plasma parameters that are going to be
met in future fusion machines have been achieved in this scenario (in
ITER/DEMO: $\mathcal{E}_{EP}/T_{Bulk}\approx 3.5MeV/30keV$). In such a way an
intense EP-driven activity is observed, with the stabilizing effects of the
bulk plasma minimized. That is why this is a very interesting scenario to
understand the EP-driven dynamics and to validate the codes against the
experiments. In particular this case presents a rich nonlinear physics of
interest. Around $t\approx 0.84\,$s a TAE-EPM burst is present after which, at
lower frequencies, an EGAM exhibit the typical chirping behaviour (cf. Fig.
1). According to the pick-up-coils measurements, in this scenario the most
unstable AM has toroidal mode number $n=1$.

Figure 1: Experimental spectrogram. At $t=0.84\,$s the plasma parameters,
profiles and the reconstructed magnetic equilibrium have been selected to
perform numerical simulations. At this time a TAE-EPM burst, covering the
frequency range $80\,$kHz$\,\leq\nu\leq 250\,$kHz, is observed. At lower
frequencies $\nu\approx 50\,$kHz, an EGAM exhibits the typical chirping
behaviour.

In Fig. 2 the temperature profiles of the bulk plasma species (left) and the
radial dependence of the electron density profile (right) are shown, while

Figure 2: Radial temperature profiles (left) of the bulk plasma species
(electrons and deuterium) and radial electron density profile (right) of the
NLED-AUG case.

in Fig. 3 left the radial dependence of the EP density profile $n_{EP}$
modelled by TRANSP [35] is shown. The EPs have an off-axis radial density
profile with a reference concentration of $\langle n_{EP}\rangle/\langle
n_{e}\rangle=0.0949$, where $\langle...\rangle$ indicates volume average. In
Fig. 3 right the radial dependence of the safety factor profile $q$ is shown.
It has a reversed shear, with a minimum around $r\approx 0.5$. Additionally
the TAE position $r\approx 0.738$ is indicated (black point). The TAE is
formed in the gap created by the interaction between the two most unstable
scalar potential Fourier components: $(m,n)=(2,1)$ and $(m,n)=(3,1)$.

Figure 3: Left: EPs off-axis radial density profile as modelled by TRANSP. It
has a reference concentration of $\langle n_{EP}\rangle/\langle
n_{e}\rangle=0.0949$. Right: Safety factor profile of the NLED-AUG case. It
exhibits a non-monotonic profile with a minimum around $r\approx 0.5$ where
$q\approx 2.28$ and a maximum at $r=1$ in the amount of $q\approx 8.5$. The
inset figure marks the TAE location at $r_{0}\approx 0.738$, corresponding to
$q_{0}=q(r_{0})=(2\,m_{0}+1)/2=2.5$, with $n=1$ and $m_{0}=2$.

In Tab. 1 the values of some important constants that will be used in the
simulations are shown. In this work frequencies and growth rates will be
provided in $\omega_{A0}$-units, that is the Alfvén frequency on-axis, at
$r=0$: $\omega_{A0}=v_{A}(r=0)/R_{0}$, with $v_{A}$ the Alfvén speed (cf. eq.
2) and $R_{0}$ the major radius of the tokamak.

$a_{0}\,[$m$]$ |  $R_{0}\,[$m$]$ |  $B_{0}\,[$T$]$ | $\beta_{e}$ | $L_{x}$ |  $\omega_{ci}\,[$rad/s$]$ |  $\omega_{A0}\,[$rad/s$]$ | $\omega_{ci}/\omega_{A0}$
---|---|---|---|---|---|---|---
$0.482$ | $1.666$ | $2.202$ | $2.7\cdot 10^{-4}$ | $551.6$ | $1.055\cdot 10^{8}$ | $4.98\cdot 10^{6}$ | $\approx 21$
Table 1: Constants in use: averaged minor radius $a_{0}$, major radius
$R_{0}$, amplitude of the background magnetic field on axis $B_{0}$.
Normalized electron pressure
$\beta_{e}=4\pi\,n_{e}(r_{0})T_{e}(r_{0})/B_{0}^{2}$ with the $r_{0}$ a
reference radial position. Normalized size of the plasma system that is
$L_{x}=2\,a_{0}/\rho_{s}$, with $\rho_{s}=c_{s}/\omega_{ci}$ and $c_{s}$ the
sound velocity and $\omega_{ci}$ the ion cyclotron frequency. $\omega_{A0}$ is
the value of the Alfvén frequency on axis.

Many works have already been dedicated to the investigation of the mode
dynamics of this interesting experimental case. The linear AM dynamics has
been investigated in Ref. [28, 33] in simulations where both on-axis and off-
axis radial density profiles have been considered for the EPs. These have been
modelled via Maxwellian distribution functions associated to temperature
profiles constant against the radius. In Ref. [36] numerical studies of the
nonlinear EGAM dynamics have been conducted. In Ref. [37] bicoherence studies
have suggested that nonlinear coupling exists between the observed TAE-EPM
burst triggering the EGAM. In Ref. [29] the interaction between AMs and EGAMs
has been studied both with numerical simulations and through an analytical
model where only the EPs have been allowed to redistribute in phase-space,
following their full trajectories. Additionally, there the EPs have been
modelled via a double-bump-on-tail distribution function, since an anisotropy
in velocity space was needed to drive an EGAM unstable [38, 39, 40].
In the present work we focus on the AM dynamics, modelling the bulk plasma
species with Maxwellian distribution functions and retaining the
nonlinearities only in the EP dynamics, as in [29, 33]. That is, only the EPs
gyrocenter characteristics eqs. 7, 8 and 9 have $\epsilon_{\delta}=1$. The
other particle species (electrons and bulk ions) follow their unperturbed
trajectories and for them $\epsilon_{\delta}=0$. The novelty of this work is
represented by the fact that we go closer to the experimental conditions
modelling the EPs via an isotropic slowing-down distribution function. We
emphasize here that EPs are present in the NLED-AUG case as injected through
an NB. This implies that their distribution function is intrinsically
anisotropic in phase-space. In the studies carried out in this paper we
neglect the effects of the anisotropy in velocity space that are, on the other
hand, important to detail the interaction between AMs and EGAMs (that here we
do not study in contrast to what was done in Ref. [29]). Nevertheless the
isotropic slowing-down represents an improvement with respect to the
Maxwellian distribution function, as will be motivated in Sec. 4. Additionally
through its choice we will observe a good quantitative comparison with the
experiment making a further step with respect to previously published works
where already a good qualitative agreement has been obtained. In the NLED-AUG
scenario the bulk ions and the EPs are constituted by deuterium ions. In all
the simulations discussed in this work the realistic electron mass will be
considered: $m_{e}\approx m_{D}/3676$.
In Sec. 4 we report a simplified analytical derivation and discuss why the
slowing-down distribution function represents a better description of the EPs
displacement in phase-space and why this is an improvement with respect to the
Maxwellian and the double-bump-on-tail distribution functions. Later we will
describe its implementation in ORB5.
In all the simulations presented in Sec. 5 the quasi-neutrality condition will
always be satisfied $n_{e}=Z_{D}\,n_{D}+Z_{EP}\,n_{EP}$ with $Z_{D}=Z_{EP}=1$,
by keeping fixed the electron radial density profile (cf. Fig. 2 right) and
varying accordingly the EP and bulk deuterium profiles.

##  4 Slowing-down distribution function

We want to investigate the AM dynamics in the NLED-AUG case with ORB5.
Approximations are present in the performed simulations. The main
approximations here concern the choice of not retaining the bulk plasma
nonlinearities and the choice for the equilibrium distribution function of the
particle species. The importance of modelling the particle species with a
proper equilibrium distribution function can be understood by looking at the
gyrokinetic Vlasov equation eq. 6. In fact $\delta f_{sp}$ can only be evolved
once the equilibrium distribution function and its derivatives are known.
Obviously changing the equilibrium distribution function implies a
modification of the time evolution of the perturbed distribution function,
thus catching different aspects of physics.
While particles of the bulk plasma are usually well described by Maxwellian
distribution functions (that from now on we denote with $F_{M}$), the question
here is what is the shape of the distribution function that better represents
the arrangements of the EPs in phase-space. For this, we choose a non-
Maxwellian equilibrium distribution, namely the slowing-down [41, 30, 42, 43,
44] that, from now on, we simply indicate with $F_{0}$.
To derive its form we summarize the main steps presented in Ref. [30] where
the analytical form of the slowing-down is derived. There the author studies
how the fast-ions with mass $m_{EP}$ and charge $q_{EP}$ injected in plasma
through a beam are displaced in phase-space because of the collisions with the
particles of the bulk plasma. The starting point is represented by the Fokker-
Planck equation [45] for the EPs:

$\partial_{t}F_{0}+\vec{v}\cdot\nabla_{\vec{x}}F_{0}+\frac{q_{EP}}{m_{EP}}\left(\vec{E}_{0}+\frac{\bm{v}}{c}\times\vec{B}_{0}\right)\cdot\nabla_{\vec{v}}F_{0}=C(F_{0})+S\quad.$
(21)

In eq. (21) $S$ is the EPs source, while $C(F_{0})$ [46] is the collision
operator, composed by the sum of two contributions

$C(F_{0})=C(F_{0},F_{M,e})+C(F_{0},F_{M,i})\quad.$ (22)

$C(F_{0},F_{M,e})$ takes into account the effects of the collisions of the EPs
with the electrons while $C(F_{0},F_{M,i})$ considers the collisions of the
EPs with the bulk ions. Considering an EP beam axisymmetrically distributed
around the background magnetic field and writing explicitly the analytical
form of the collision operator, eq. (21) becomes [30]:

$\partial_{t}F_{0}=\frac{1}{\tau_{s}v^{3}}\Bigg{\\{}v\frac{\partial}{\partial
v}\left[\left(v^{3}+v_{c}^{3}\right)F_{0}\right]+Z_{2}\frac{v_{c}^{3}}{2}\frac{\partial}{\partial\xi}\left[\left(1-\xi^{2}\right)\frac{\partial
F_{0}}{\partial\xi}\right]\Bigg{\\}}+S(v,\xi)\quad.$ (23)

In eq. (23) $v$ is the modulus of the beam particle velocity $v=|\bm{v}|$ and
$\xi=v_{\parallel}/v=\cos\chi$, where $v_{\parallel}$ is the component of the
particle velocity parallel to the background magnetic field and $\chi$ the
pitch angle. $v_{c}$ is the crossover velocity or slowing-down critical
velocity:

$v_{c}=\left(\frac{3\sqrt{\pi}}{4}\frac{m_{e}}{m_{EP}}Z_{1}\right)^{1/3}v_{th,e}\quad,\quad
Z_{1}=\sum_{i={Bulk\,ions\,species}}\frac{n_{i}}{n_{e}}\frac{m_{EP}}{m_{i}}Z_{i}^{2}\quad$
(24)

$v_{th,e}$ is the electron thermal velocity and

$Z_{2}=\sum_{i={Bulk\,ions\,species}}\frac{n_{i}\,Z_{i}^{2}}{n_{e}Z_{1}}\quad.$
(25)

The crossover velocity is a property of the bulk plasma and represents a
critical speed as for $v>v_{c}$ the electron drag dominates over the ion drag
and the EPs are mainly decelerated by friction with the electrons. For
$v<v_{c}$ the opposite trend is met and the EPs are mainly decelerated by
friction with the bulk ions. $\tau_{s}$ [47] is the slowing-down time that
describes the time after which a particle of the energetic ion beam is
decelerated by collisions with the other particle species.
We take into account a monoenergetic EP beam. In this way the EP source in eq.
(23) becomes:

$S(v,\xi)=\frac{S_{0}}{v^{2}}\delta(v-v_{EP})\delta(\xi-\xi_{EP})\quad.$ (26)

$S_{0}$ is the intensity of the EP source, $\xi_{EP}$ is the injection
parallel velocity of the beam normalized to the modulus of the total velocity
$v_{EP}=\sqrt{2\mathcal{E}_{EP}/m_{EP}}$ with $\mathcal{E}_{EP}$ the injection
energy of the beam.
In Ref. [30] the solution of eq. (23) is nicely derived. This is quite complex
and contains the dependence in $\xi$ (or $\chi$, the pitch angle). Here, and
in the results presented below, we disregard the $\xi$-dependence. It is then
straightforward to obtain the steady state solution for an isotropic,
monoenergetic, EP beam:

$F_{0}(v)=\frac{S_{0}\tau_{s}}{4\pi}\frac{\theta(v_{EP}-v)}{v_{c}^{3}+v^{3}}\quad.$
(27)

This $\xi$-independent slowing-down distribution function is appropriate to
describe fusion-born products [43] since $\alpha$ particles are not born in
any preferential direction (they are isotropic). On the other hand eq. (27)
does not appropriately describe the arrangement in phase-space of the EPs
belonging to an injected ion beam. In this case, in fact, the EPs are not
isotropic. Eq. (27) represents nevertheless an improvement with respect to the
Maxwellian distribution function and it is widely used in literature to model
EPs [43]. In fact, through the Heaviside step function $\theta(v_{EP}-v)$, it
provides a constraint in the velocity distribution of the EPs present
(born/injected) in plasma with $\mathcal{E}_{EP}$, since $0\leq v\leq v_{EP}$.
Additionally its dependence in velocity-space is governed by the crossover
velocity that, as previously stated, is a property of the bulk plasma.
By taking the zero and second order moments of eq. (27) we find [42]:

$\int
d^{3}\vec{v}\,F_{0}=S_{0}\tau_{s}I_{2}\left(\frac{v_{c}}{v_{EP}}\right)\quad,\quad\int
d^{3}\vec{v}\,v^{2}F_{0}=S_{0}\tau_{s}I_{4}\left(\frac{v_{c}}{v_{EP}}\right)v_{EP}^{2}\quad,$
(28)

where:

$I_{n}(a)=\int_{0}^{1}dx\frac{x^{n}}{x^{3}+a^{3}}=\begin{cases}n=2\quad,\quad\frac{1}{3}\ln(1+a^{-3})\\\
\\\
n=4\quad,\quad\frac{1}{2}-a^{2}\Big{\\{}\frac{1}{6}\ln\left(\frac{1-a+a^{2}}{(1+a)^{2}}\right)+\frac{1}{\sqrt{3}}\left[\tan^{-1}\left(\frac{2-a}{a\sqrt{3}}\right)+\frac{\pi}{6}\right]\Big{\\}}\quad.\end{cases}$
(29)

To ensure the zero order moment in eq. 28 to give the EPs’ radial density
profile $n_{EP}$, we impose:

$S_{0}=\frac{n_{EP}}{\tau_{s}\,I_{2}\left(\frac{v_{c}}{v_{EP}}\right)}\quad.$
(30)

With this choice it follows:

$F_{0}=\frac{n_{EP}}{\frac{4\pi}{3}\log\left[1+\left(\frac{v_{EP}}{v_{c}}\right)^{3}\right]}\frac{\theta(v_{EP}-v)}{v_{c}^{3}+v^{3}}\quad$
(31)

that represents the analytical form of the slowing-down distribution function
properly normalized. Equation 31 has been implemented in ORB5 together with
its derivatives that are needed to evolve the perturbed EP distribution
functions $\delta f_{EP}$ eq. 6:

$\begin{cases}\frac{d\,F_{0}}{d\,\psi}=\Bigg{\\{}\frac{d}{d\,\psi}\log
n_{EP}+\left(\frac{(v_{EP}/v_{c})^{3}}{\left[1+\left(\frac{v_{EP}}{v_{c}}\right)^{3}\right]\log\left[1+\left(\frac{v_{EP}}{v_{c}}\right)^{3}\right]}-\frac{v_{c}^{3}}{v_{c}^{3}+v^{3}}\right)\cdot
3\frac{d}{d\,\psi}\log v_{c}\Bigg{\\}}F_{0}\\\ \\\
\frac{d}{d(m_{EP}\mathcal{E})}F_{0}=-\frac{3v}{v_{c}^{3}+v^{3}}\frac{F_{0}}{m_{EP}}\\\
\\\ \frac{d}{dv_{\parallel}}F_{0}=0\end{cases}$ (32)

where $v=\sqrt{2\mathcal{E}}$ and:

$\frac{dv_{c}}{d\,\psi}=v_{c}\left[\frac{d}{d\,\psi}\log
v_{th,e}-\frac{1}{3}\frac{d}{d\,\psi}\log
n_{e}+\frac{1}{3}\left(\sum_{l}\frac{z_{l}^{2}}{m_{l}}\frac{dn_{l}}{d\,\psi}\right)\cdot\left(\sum_{p}\frac{z_{p}^{2}}{m_{p}n_{p}}\right)\right]\quad,$
(33)

where $\psi$ is the normalized poloidal flux. In the notation in use:

$\mathcal{E}=\frac{v^{2}}{2}=\frac{v^{2}_{\parallel}}{2}+\mu
B\quad\text{with}\quad\mu=\frac{v^{2}_{\perp}}{2B}\quad\text{and}\quad\mathcal{E}_{EP}=m_{EP}\,\mathcal{E}\quad.$
(34)

We conclude this section introducing the concept of equivalent temperature
$T_{EP}$ that is the temperature of a Maxwellian distribution function $F_{M}$
for the EPs that satisfies the following requirement:

$\int\,d^{3}\vec{v}\,v^{2}\,F_{M}=\int\,d^{3}\vec{v}\,v^{2}\,F_{0}\quad\text{with}\quad\int\,d^{3}\vec{v}\,v^{2}\,F_{M}=3\,n_{EP}(r)\frac{T_{EP}(r)}{m_{EP}}\quad$
(35)

from which it follows:

$T_{EP}(r)=\frac{2}{3}\frac{I_{4}\left(\frac{v_{c}(r)}{v_{EP}}\right)}{I_{2}\left(\frac{v_{c}(r)}{v_{EP}}\right)}\mathcal{E}_{EP}\quad.$
(36)

###  4.1 Implementation in ORB5

As it has been already anticipated in Sec. 2, in ORB5 all the quantities are
normalized with respect to reference parameters. In particular, velocities are
given in units of the sound velocity $c_{s}$:

$c_{s}=\sqrt{\frac{T_{e}(r_{0})}{m_{D}}}$ (37)

defined as the square root of the electron temperature considered at a chosen
reference position $r_{0}$, divided by the mass of the main ion species that,
for the NLED-AUG case, is constituted by deuterium ions. Here the reference
radial position is on-axis $r_{0}=0$ and $T_{e}(r_{0})=0.7088\,$keV . In Fig.
4 the velocity dependence in $(v_{\parallel},\mu\,B)$-space of the slowing-
down distribution function at $r=0.5$ is shown. This has been taken from an
ORB5-simulation. It has been obtained considering the bulk plasma parameters
of the NLED-AUG case (see Sec. 3) as well as the reference EPs parameters:
off-axis radial density profile in Fig. 3 left having a concentration of
$\langle n_{EP}\rangle/\langle n_{e}\rangle=0.0949$, injected with injection
energy $\mathcal{E}_{EP}=93\,$keV.

Figure 4: Slowing-down distribution implemented in ORB5 for the EPs. The
dependence in velocity-space at a fixed radial position is shown.

In Fig. 5 slices obtained from Fig. 4 are shown. There we can observe the cuts
in velocity-space introduced because of the presence of the Heaviside function
in eq. 31:

$v_{\parallel}\approx 16.2\,c_{s}\quad,\quad\mu\,B\approx
131.1\,c_{s}^{2}\quad.$ (38)

Using the appropriate conversion factor the interested reader can easily
verify that these cuts correspond to an injection energy of
$\mathcal{E}_{EP}=93\,$keV. In Fig. 6 we show the EP radial density profile
(on the top, corresponding to Fig. 3 left) and the equivalent temperature (on
the bottom, cf. eq. 36). The analytical expressions (black dashed lines) are
compared with the corresponding moments (blue lines) calculated through
numerical integration of the distribution function implemented in ORB5. The
observed good agreement between the analytical expressions and the numerical
calculations is a further proof of the correct implementation of the new non-
Maxwellian background distribution function in the code.

Figure 5: Slices of the equilibrium distribution function obtained from Fig.
4.
Left: Dependence against $v_{\parallel}$: $F_{0}(r=0.5,v_{\parallel},\mu
B=0)$. Right: Dependence against $\mu B$: $F_{0}(r=0.5,v_{\parallel}=0,\mu
B)$.
The vertical black dashed lines correspond to the cuts in velocity space
present because of the Heaviside theta contained in eq. (31). Figure 6:  Top:
EP Radial density profile corresponding to Fig. 3 left. Bottom: EP equivalent
temperature, cf. eq. (36).
$\langle...\rangle_{\bm{v}}$ indicates integration in velocity space. The
moments of the equilibrium distribution function $F_{0}$ implemented in ORB5
and calculated through numerical integration (blue curves) are compared with
the expected values of the analytical expressions (black dashed lines).

In Fig. 7 slices of slowing-down distribution functions obtained considering
different EP concentration, electron temperature and injection energy are
shown. The corresponding moments are shown in Fig. 8. Also there a good
agreement is observed between the analytical expressions (black dashed lines)
and the numerical integrals of the implemented slowing-down.

Figure 7: Slices of equilibrium distribution functions obtained considering
different EP concentration, values of electron temperature and injection
energy. $T_{e,ref}$ indicates the reference electron temperature of the NLED-
AUG case (cf. Fig. 2 left).
Left: Dependence against $v_{\parallel}$. Right: Dependence against $\mu B$.
Figure 8:  Top: EP Radial density profiles. Bottom: EP equivalent temperature
profiles.
The black-dashed lines correspond to the analytical expressions, while the
coloured lines are calculated through numerical integration of the implemented
distribution functions. The coloured lines are obtained with the plasma
parameters considered to compute the slices of the same colour in Fig. 7

##  5 Comparison with the experiments: application to the NLED-AUG case

We discuss here the mode dynamics observed in numerical simulations where the
NLED-AUG scenario (see Sec. 3) has been taken into account. In Tab. 2 the main
parameters considered in the simulations are reported. Their choice is
motivated by convergence studies (see Appendix A). As it is indicated in Tab.
2, we will investigate here only the dynamics of modes with toroidal mode
number $n=1$. This because, as anticipated in Sec. 3, from experimental
measurements this is known to be the most unstable AM. The dynamics of zonal
structures ($n=0$) is then not treated in the present work.

$\Delta t\,[\omega_{ci}^{-1}]$ | $nptot_{D,e,EP}\cdot 10^{7}$ | $N_{r}$ | $N_{\theta^{*}}$ | $N_{\varphi}$ | $nfilt$ | $\Delta m$
---|---|---|---|---|---|---
$4$ | $3,12,3$ | $3000$ | $64$ | $32$ | $1$ | $7$
Table 2: Main simulation parameters: time step ($\Delta t$), number of markers
for particle species ($nptot$). Radial/poloidal/toroidal grid points
$(N_{r},N_{\theta^{*}},N_{\varphi})$, toroidal ($nfilt$) and poloidal width
($\Delta m$) of the field aligned Fourier filter (see Sec. 2).

We begin studying the mode dynamics observed in the exponential growth phase
of the dominant mode and compare the results obtained with a Maxwellian and a
slowing-down distribution. For the moment we run ORB5 neglecting the FLR
effects of all the particle species, i.e. considering a drift-kinetic model.
In all the simulations analyzed the dominant mode, in the exponential growth
phase, is an Alfvén mode with scalar potential dominated by its components
$(m,n)=(2,1)$ and peaked around $r=0.2$.

Figure 9: Growth rates determined in drift-kinetic simulations (no FLR
effects). Left: Scan against the EP concentration for Maxwellian ($F_{M}$) and
slowing-down ($F_{0}$) distribution functions. $T_{EP}=25\,$keV represents the
value on the center of the equivalent temperature (cf. eq. 36) corresponding
to an injection energy of $\mathcal{E}_{EP}=93\,$keV. Right: Scan against the
EP temperature in simulations where the EPs have been modelled through
Maxwellian distribution functions.

In Fig. 9 left the growth rates determined in a scan against the EP
concentration are shown. There we compare the growth rates obtained in
simulations where the EPs have been modelled with Maxwellian distribution
functions ($F_{M}$) with those obtained modelling the EPs with the slowing-
down distribution function ($F_{0}$). In Fig. 9 right the growth rate
dependence against the EPs temperature (with constant radial profile) is shown
in a scan where the EPs have been modelled via Maxwellian distribution
functions. In Fig. 10 we display the real frequencies of the dominant modes
observed in the corresponding simulations that have produced the growth rates
reported in Fig. 9. Only in Fig. 10 the error bars are present, since the
errors committed determining the growth rates in Fig. 9 are negligible.

Figure 10: Real frequencies determined in the growing exponential phase of the
dominant modes in drift-kinetic simulations (no FLR effects considered). Left:
Scan against the EP concentration for Maxwellian ($F_{M}$) and slowing-down
($F_{0}$) distribution functions.
Right: Scan against the EP temperature in simulations where the EPs have
Maxwellian distribution functions.

In the scans presented in Fig. 9 and Fig. 10, the case closer to the
experimental conditions is that where the EPs have been modelled with slowing-
down distribution function and concentration of $\langle n_{EP}\rangle/\langle
n_{e}\rangle=0.0949$. Its values of growth rate and frequency are reported in
Tab. 3. There, they are compared with those obtained in a simulation where the
EPs have been modelled with a Maxwellian distribution function with
$T_{EP}=25\,$keV. This corresponds to the value of the equivalent temperature
around $r\approx 0.5$ (cf. with Fig. 6).

Distribution function | $\gamma\,[\omega_{A0}]$ | $\omega\,[\omega_{A0}]$
---|---|---
Slowing-down, $\mathcal{E}_{EP}=93\,$keV | $0.023$ | $0.116\pm 0.008$
Maxwellian, $T_{EP}=25\,$keV | $0.016$ | $0.128\pm 0.005$
Table 3: Growth rates $\gamma$ and frequencies $\omega$ of the dominant modes,
determined in drift-kinetic simulations (no FLR effects included). The
reference EP concentration of $\langle n_{EP}\rangle/\langle
n_{e}\rangle=0.0949$ is here considered.

Referring to Tab. 3, we can then conclude that the choice of the slowing-down
distribution function over a Maxwellian distribution function with temperature
close to the equivalent temperature, does not affect significantly the
frequency of driven AM, although quantitative differences as large as almost
$40\,\%$ are found in the growth rate.

Figure 11: Temporal evolution of the dominant poloidal harmonics ($m$-numbers
indicated in the legends) in the simulations. Left: EPs modelled with slowing-
down distribution and $\mathcal{E}_{EP}=93\,$keV. Right: EPs modelled with the
equivalent Maxwellian.

We describe now in detail the mode dynamics of the AM observed in simulations
where FLR effects are included for bulk ions and EPs (the model is
gyrokinetic) and the reference EP concentration of $\langle
n_{EP}\rangle/\langle n_{e}\rangle=0.0949$ is taken into account. Here we will
not just limit ourselves to the investigation of the exponential growth phase
of the AM, but also to its nonlinear phase. There the AM mode structure
evolves radially, presenting peaks close to the edge. Therefore we have
increased the radial resolution to $N_{r}=6000$ in order to correctly resolve
those peaks. Simultaneously since in ORB5 the numerical noise is proportional
to the number of radial grid points $N_{r}$ [48], we have increased the number
of markers to $nptot_{E,e,EP}=(18,72,18)\cdot 10^{7}$.
The temporal dynamics observed in a simulation where the EPs have been
modelled with the slowing-down distribution function with
$\mathcal{E}_{EP}=93\,$keV, is compared with that observed in a simulation
where the EPs are modelled with the equivalent Maxwellian (the equivalent
temperature in Fig. 6 is taken now into account). The former is shown in Fig.
11 left, while the latter is shown in Fig. 11 right. In the two plots in Fig.
11 we have highlighted the temporal intervals corresponding to: the
exponential growth phase (yellow), the early nonlinear phase (pink) and the
deep nonlinear phase (blue). Note that both the simulations cover the same
temporal range: $t\in[0;100000]\,\omega_{ci}^{-1}$.

Distribution function | $\gamma\,[\omega_{A0}]$ | $\omega\,[\omega_{A0}]$
---|---|---
Slowing-down, $\mathcal{E}_{EP}=93\,$keV | $0.019$ | $0.133\pm 0.007$
Equivalent Maxwellian | $0.015$ | $0.138\pm 0.005$
Table 4: Growth rates $\gamma$ and frequencies $\omega$ of the dominant modes,
determined in gyrokinetic simulations (FLR effects of the EPs are retained).
The reference EP concentration of $\langle n_{EP}\rangle/\langle
n_{e}\rangle=0.0949$ is here considered.

The growth rate and frequencies determined in the exponential growth phase of
the two simulations in Fig. 11 are reported in Tab. 4, while the corresponding
mode structures and frequency spectra are shown in Fig. 12. The dominant AM
here, as in the drift-kinetic simulations, is an EPM whose scalar potential is
dominated by the Fourier component $(m,n)=(2,1)$ peaked around the radial
position $r\approx 0.2$. Its frequency lies almost on the continuum spectrum
calculated with the linear gyrokinetic code LIGKA [49] (see red dashed lines
in Fig. 12 on the bottom). We can conclude also here that, as for the drift-
kinetic simulations, the choice of the equivalent Maxwellian over the slowing-
down distribution function, does not modify significantly the kind of driven
AM investigated (mode structure and frequency). However, referring to Tab. 4,
we still observe a growth rate quantitative difference that in this case is as
large as almost $30\,\%$ of the value of the growth rate obtained in the
simulation with the EPs modelled with the equivalent Maxwellian.

Figure 12: Mode structures and frequency spectra observed in the exponential
growth phases (yellow regions) of the simulations in Fig. 11, where the EPs
have been modelled with slowing-down (left) and equivalent Maxwellian (right).

Figure 13: Frequency spectra calculated in the early nonlinear phases (pink
regions) of the simulations in Fig. 11. Left: EPs modelled with slowing-down
distribution function. Right: EPs modelled via equivalent Maxwellian.

In Fig. 13 we compare the frequency spectra calculated in the early nonlinear
phases (pink regions) of the two different simulations shown in Fig. 11. In
this temporal phase the scalar potential Fourier component $(m,n)=(2,1)$ is
still the dominant one, since it is almost an order of magnitude higher than
the other Fourier components. Comparing Fig. 13 with Fig. 12, we observe that
both simulations catch correctly the nonlinear modification of the AM
frequency passing from the linear to the nonlinear phase. In this time
interval, the dominant AM frequency exhibits tiny differences modelling the
EPs with slowing-down or equivalent Maxwellian. Additionally in both
simulations, at the end of this phase, the Fourier component $(m,n)=(3,1)$
begins to grow, even though with different growth rate, as indicated in Tab.
5.

Distribution function | $\gamma\,[\omega_{A0}]$
---|---
Slowing-down, $\mathcal{E}_{EP}=93\,$keV | $0.002$
Equivalent Maxwellian | $0.00053$
Table 5: Growth rates $\gamma$ of $(m,n)=(3,1)$ in the deep nonlinear phase
(blue regions in Fig. 11). Gyrokinetic simulations (FLR effects of the EPs are
retained). The reference EP concentration of $\langle n_{EP}\rangle/\langle
n_{e}\rangle=0.0949$ is here considered.

Figure 14: Frequency spectra calculated in the deep nonlinear phases (blue
regions) of the simulations in Fig. 11. Left: EPs modelled with slowing-down
distribution function. Right: EPs modelled via equivalent Maxwellian.

In Fig. 14 the frequency spectra calculated in the deep nonlinear phases (blue
regions in Fig. 11) are shown. In these temporal phases the scalar potential
Fourier component $(m,n)=(3,1)$ reaches an amplitude comparable to that of the
dominant $(m,n)=(2,1)$. Passing from Fig. 13 to Fig. 14 we observe that in the
simulation where the EPs have been modelled with the slowing-down (left), a
higher frequency modification of the AM is detected with respect to the right
plots where the EPs had been modelled with the equivalent Maxwellian. In fact
in Fig. 14 left the AM oscillates mainly at the TAE frequency. We remind here
that this is an AE whose frequency lies inside the gap created by the branches
of the continuum $(m_{0},n)$ and $(m_{0}+1,n)$ at the radial position $r_{0}$
that labels the rational surface satisfying $q(r_{0})=(2\,m_{0}+1)/(2\,n)$
[13] where two close poloidal harmonics interact. For the case under
investigation: $m_{0}=2$, $n=1$ and $r_{0}=0.738$ (cf. Fig. 3 right). We
emphasize that in Fig. 14 left the AM oscillates almost at every radial
position at the TAE frequency while, in Fig. 14 right, such a higher frequency
modification is not observed. Finally we note in Fig. 14 the presence, in this
deep nonlinear phase, of higher frequencies in the EAE gap: $\omega\sim
0.6,-0.5\,\omega_{A0}$. We do not further investigate the emerging of these
higher frequencies as the focus of the present paper is on the study of the
TAE-EPM frequency modification observed in the experiment (cf. Fig. 1), that
correspond to the dominant frequencies also in the presented simulations.
We summarize the results of the AM frequency modification in Fig. 15. There
the spectrogram calculated in a simulation with EPs modelled with the slowing
down (left) is compared with that obtained in the simulation where the EPs
have been modelled with the equivalent Maxwellian (right). The spectrograms
presented have been obtained through the fast Fourier transform of the scalar
potential at $r=0.2$. This corresponds to the radial positions where the mode
structure of the AM has a maximum. The time interval where the Fourier
transform has been calculated $[t_{0};t_{0}+\Delta t]$ has been obtained
varying continuously $t_{0}$ (as indicated in the $x$-axis in Fig. 15) and
choosing $\Delta t=10^{4}\,\omega_{ci}^{-1}$. In Fig. 15 the time units are
provided also in milliseconds, while the frequency units are provided also in
kHz, to allow a better comparison with the experimental spectrogram in Fig. 1.
We notice, at first, that the total time width of the simulations is of
$\approx\,1\,$ms which lies well inside the width of the TAE-EPM burst at
$t=0.84\,$s (cf. Fig. 1). Finally, the choice of the slowing-down over the
equivalent Maxwellian allows a better comparison with the experiment. Through
this, in fact, a broader modification of the AM frequency in the nonlinear
phase is observed $100\,$kHz$\leq\nu\leq 180\,$kHz, resulting in a better
agreement with the experiment (cf. Fig. 1).

Figure 15: Spectrograms calculated from ORB5 simulations with the EPs modelled
with different distribution functions (as indicated in the plot titles). The
Fourier transform of the signal has been calculated in the temporal range
$t\in[t_{0},t_{0}+\Delta t]$ with $\Delta t=10^{4}\,\omega_{ci}^{-1}$. The low
boundary $t_{0}$ has been varied continuously as indicated in the $x$-axis of
the plot.

##  6 Conclusion

In this work the study of the AM instabilities driven by EPs has been
conducted for the NLED-AUG scenario. This is an interesting case obtained in
the Tokamak ASDEX Upgrade where the EPs have been injected through a NB with
energy $\mathcal{E}_{EP}\approx 93\,$keV. The peculiarity of this scenario
relies on the fact that it has been obtained through plasma ratios
($\beta_{EP}/\beta_{Bulk}$ and $\mathcal{E}_{EP}/T_{Bulk}$) close to those
that are going to be met in future fusion machines.
Studies of the AM dynamics in the NLED-AUG case have already been conducted
through numerical simulations [28, 29, 33]. The main novelty in the present
work is represented by the fact that the EPs have been modelled thorough an
equilibrium isotropic slowing-down distribution function. Through it, we have
been able to go closer to the experimental conditions. In particular we have
obtained a good quantitative agreement with the experimental spectrogram doing
a further step with respect to previous published works [28, 29, 33]. There,
in fact, a good qualitative agreement with the experiment has been obtained
and the EPs have been modelled with Maxwellian or double-bump-on-tails, whose
choice was motivated by the kind of physics we were interested to investigate.
In the present work we have proved that the isotropic slowing-down
distribution function has been correctly implemented in ORB5. We have also
shown that its choice allows a more accurate description of the frequency
modification of the AM in the nonlinear phase, while the choice of an
equivalent Maxwellian underestimate it. The key conclusion is that with the
inclusion of the FLR effects we have been able to achieve a good quantitative
agreement with the experimental spectrogram. These results are already
encouraging and prove that we are able to reproduce the relevant physics
behind the experiments in an attempt to become predictive about future fusion
scenarios and pave the path to new studies that have to be carried out to
catch more of the nonlinear aspects in the NLED-AUG case.
In future works we will retain the nonlinearities in all the particle species
and we will consider the interaction between the AM and the EGAM, to describe
the triggering of the EGAM by the TAE-EPM burst. An attempt on this has
already been done in Ref. [29] where the EPs have been modelled with a double-
bump-on-tail. Moreover the AM and EGAM interaction will be studied considering
a $\xi$-dependent slowing-down distribution function [50] in order to retain
the anisotropy in velocity space needed to drive the EGAM.

##  7 Acknowledgments

Simulations presented in this work were performed on the CINECA Marconi
supercomputer within the ORBFAST, OrbZONE and TSVV10 projects. The authors
acknowledge stimulating discussions with Z. Lu, F. Zonca, S. Briguglio, I.
Novikau and A. Di Siena. This work was partly performed in the frame of the
“Multi-scale Energetic Particle Transport in Fusion Devices” ER project. This
work has been carried out within the framework of the EUROfusion Consortium,
funded by the European Union via the Euratom Research and Training Programme
(Grant Agreement No 101052200 — EUROfusion). Views and opinions expressed are
however those of the author(s) only and do not necessarily reflect those of
the European Union or the European Commission. Neither the European Union nor
the European Commission can be held responsible for them.

##  Appendix A Slowing-down validation and convergence studies

We present here a series of convergence studies used to find the reference
simulation parameters in Tab. 2. In Fig. 16, Fig. 17, Fig. 18 and Fig. 19 we
show the discrepancies $\sigma$ of the determined quantities $X$ from the most
accurate value, that is:

$\sigma\,[\%]=\frac{X-X_{most\,accurate\,vale}}{X_{most\,accurate\,vale}}\cdot
100\quad.$ (39)

The most accurate value $X_{most\,accurate\,value}$ corresponds to the value
(almost) at convergence. Its discrepancy, because of the definition in eq.
(39), is $\sigma=0\,\%$. In Fig. 16, Fig. 17, Fig. 18 and Fig. 19 the orange
points have been obtained in simulations where the reference parameters in
Tab. 6 have been considered.

$\Delta t\,[\omega_{ci}^{-1}]$ | $nptot_{D,e,EP}\cdot 10^{7}$ | $N_{r}$
---|---|---
$4$ | $3,12,3$ | $1000$
Table 6: Reference simulation parameters: time step ($\Delta t$), number of
markers for particle species ($nptot$) and radial grid points $(N_{r})$, cf.
Sec. 2.

In all the simulations here presented FLR effects are retained and the
reference EP concentration of $\langle n_{EP}\rangle/\langle
n_{e}\rangle=0.0949$ is considered. The quantities $X$ here studied (growth
rates $\gamma$ and frequencies $\omega$) have been determined in the
exponential growth phase of the simulations. The scans here presented are
against: the time width $\Delta t$, the number of markers $nptot$ and the
number of radial grid points $N_{r}$ (cf. Sec. 2). The growth rate and
frequency obtained using the reference parameters in Tab. 6 correspond to:

$\gamma_{ref}=0.0197\,\omega_{A0}\quad,\quad\omega_{ref}=0.133\,\omega_{A0}\quad.$
(40)

Figure 16: Simulations with $\\#$ of markers${}_{D,e,EP}=(3,12,3)\cdot 10^{7}$
and $N_{r}=1000$. Scan against the time step $\Delta t$ in use in the
simulations. In the plots, the discrepancies $\sigma$ of growth rate (left)
and frequencies (right) with respect to the best resolved case (at $\Delta
t=0.25\,\omega_{ci}^{-1}$) are shown. The orange points correspond to the
results obtained with the parameters shown in Tab. 6.

Figure 17: Simulations with $\Delta t=4\,\omega_{ci}^{-1}$ and $N_{r}=1000$.
Scan against the $\\#$ of markers of Deuterium and EPs. The number of markers
for the electrons is four times higher than that of the EPs. In the plots the
discrepancies $\sigma$ of growth rate (left) and frequencies (right) with
respect to the best resolved case (corresponding to the simulation with number
of markers $=30\cdot 10^{7}$) are shown. The orange points correspond to the
results obtained with the parameters shown in Tab. 6.

Figure 18: Simulations with $\Delta t=4\,\omega_{ci}^{-1}$ and $N_{r}=1000$.
Scan against the $\\#$ of markers of the electrons. The number of markers for
EPs and Deuterium particles is fixed at $3\cdot 10^{7}$. In the plots the
discrepancies $\sigma$ of growth rate (left) and frequencies (right) with
respect to the best resolved case (here corresponding to the simulation with
number of markers $=50\cdot 10^{7}$) are shown. The orange points correspond
to the results obtained with the parameters shown in Tab. 6.

Figure 19: Simulations with $\Delta t=4\,\omega_{ci}^{-1}$ and $\\#$ of
markers${}_{D,e,EP}=(3,12,3)\cdot 10^{7}$. Scan against the $\\#$ of radial
grid points $N_{r}$. In the plots the discrepancies $\sigma$ of growth rate
(left) and frequencies (right) with respect to the best resolved case (here
corresponding to the simulation with $N_{r}=6000$) are shown. The orange
points correspond to the results obtained with the parameters shown in Tab. 6.

## References

  * [1] K. Tomabechi, J. R. Gilleland, Yu. A. Sokolov, and R. Toschi.  ITER conceptual design.  Nuclear Fusion, 31(6):1135–1224, jun 1991.  doi:10.1088/0029-5515/31/6/011.
  * [2] H. Zohm, C. Angioni, E. Fable, G. Federici, G. Gantenbein, T. Hartmann, K. Lackner, E. Poli, L. Porte, O. Sauter, G. Tardini, D. Ward, and M. Wischmeier.  On the physics guidelines for a tokamak DEMO.  Nuclear Fusion, 53(7):073019, jun 2013.  doi:10.1088/0029-5515/53/7/073019.
  * [3] W. W. Heidbrink and G. J. Sadler.  The behaviour of fast ions in tokamak experiments.  Nuclear Fusion, 34(4):535–615, apr 1994.  doi:10.1088/0029-5515/34/4/i07.
  * [4] S. M. Mahajan.  Spectrum of Alfvén waves, a brief review.  Physica Scripta, T60:160–170, jan 1995.  doi:10.1088/0031-8949/1995/t60/020.
  * [5] G. Vlad, F. Zonca, and S. Briguglio.  Dynamics of Alfvén waves in tokamaks.  La Rivista del Nuovo Cimento (1978-1999), 22:1–97, 1999.
  * [6] E. M. Barston.  Electrostatic oscillations in inhomogeneous cold plasmas.  Annals of Physics, 29(2):282–303, 1964.  URL: https://www.sciencedirect.com/science/article/pii/0003491664902222, doi:https://doi.org/10.1016/0003-4916(64)90222-2.
  * [7] Z. Sedláček.  Electrostatic oscillations in cold inhomogeneous plasma I. Differential equation approach.  Journal of Plasma Physics, 5(2):239–263, 1971.  doi:10.1017/S0022377800005754.
  * [8] H. Grad.  Plasmas.  Physics Today, 22(12):34, 1969.
  * [9] P. Lauber.  Super-thermal particles in hot plasmas—Kinetic models, numerical solution strategies, and comparison to tokamak experiments.  Physics Reports, 533(2):33–68, December 2013.  doi:10.1016/j.physrep.2013.07.001.
  * [10] W. W. Heidbrink.  Basic physics of Alfvén instabilities driven by energetic particles in toroidally confined plasmas.  Physics of Plasmas, 15(5):055501, 2008.  arXiv:https://doi.org/10.1063/1.2838239, doi:10.1063/1.2838239.
  * [11] L. Chen and F. Zonca.  Physics of Alfvén waves and energetic particles in burning plasmas.  Rev. Mod. Phys., 88:015008, Mar 2016.  URL: https://link.aps.org/doi/10.1103/RevModPhys.88.015008, doi:10.1103/RevModPhys.88.015008.
  * [12] L. Chen.  Theory of magnetohydrodynamic instabilities excited by energetic particles in tokamaks.  Physics of Plasmas, 1(5):1519–1522, 1994.  arXiv:https://doi.org/10.1063/1.870702, doi:10.1063/1.870702.
  * [13] C. Z. Cheng, L. Chen, and M. S. Chance.  High-n ideal and resistive shear Alfvén waves in tokamaks.  Annals of Physics, 161(1):21–47, 1985.  URL: https://www.sciencedirect.com/science/article/pii/0003491685903355, doi:https://doi.org/10.1016/0003-4916(85)90335-5.
  * [14] C. Z. Cheng and M. S. Chance.  Low‐n shear Alfvén spectra in axisymmetric toroidal plasmas.  The Physics of Fluids, 29(11):3695–3701, 1986.  URL: https://aip.scitation.org/doi/abs/10.1063/1.865801, arXiv:https://aip.scitation.org/doi/pdf/10.1063/1.865801, doi:10.1063/1.865801.
  * [15] W. W. Heidbrink, E. J. Strait, E. Doyle, G. Sager, and R. T. Snider.  An investigation of beam driven Alfvén instabilities in the DIII-D tokamak.  Nuclear Fusion, 31(9):1635–1648, sep 1991.  doi:10.1088/0029-5515/31/9/002.
  * [16] S. E. Sharapov, B. Alper, H. L. Berk, D. N. Borba, B. N. Breizman, C. D. Challis, A. Fasoli, N. C. Hawkes, T. C. Hender, J. Mailloux, S. D. Pinches, and D. Testa.  Alfvén wave cascades in a tokamak.  Physics of Plasmas, 9(5):2027–2036, 2002.  cited By 131.  URL: https://www.scopus.com/inward/record.uri?eid=2-s2.0-0012632611&doi=10.1063%2f1.1448346&partnerID=40&md5=2d120ca22ae0f5f423bd6ece2e0e0fb6, doi:10.1063/1.1448346.
  * [17] H. L. Berk, D. N. Borba, B. N. Breizman, S. D. Pinches, and S. E. Sharapov.  Theoretical interpretation of Alfvén cascades in tokamaks with nonmonotonic q profiles.  Physical Review Letters, 87(18):1850021–1850024, 2001.  cited By 154.  URL: https://www.scopus.com/inward/record.uri?eid=2-s2.0-0035968822&partnerID=40&md5=063c3fae02c2bf6d6bc5d4e03ec128a0.
  * [18] K. McGuire, R. Goldston, M. Bell, M. Bitter, K. Bol, K. Brau, D. Buchenauer, T. Crowley, S. Davis, F. Dylla, H. Eubank, H. Fishman, R. Fonck, B. Grek, R. Grimm, R. Hawryluk, H. Hsuan, R. Hulse, R. Izzo, R. Kaita, S. Kaye, H. Kugel, D. Johnson, J. Manickam, D. Manos, D. Mansfield, E. Mazzucato, R. McCann, D. McCune, D. Monticello, R. Motley, D. Mueller, K. Oasa, M. Okabayashi, K. Owens, W. Park, M. Reusch, N. Sauthoff, G. Schmidt, S. Sesnic, J. Strachan, C. Surko, R. Slusher, H. Takahashi, F. Tenney, P. Thomas, H. Towner, J. Valley, and R. White.  Study of high-beta magnetohydrodynamic modes and fast-ion losses in PDX.  Phys. Rev. Lett., 50:891–895, Mar 1983.  URL: https://link.aps.org/doi/10.1103/PhysRevLett.50.891, doi:10.1103/PhysRevLett.50.891.
  * [19] R. B. White, R. J. Goldston, K. McGuire, Allen H. Boozer, D. A. Monticello, and W. Park.  Theory of mode‐induced beam particle loss in tokamaks.  The Physics of Fluids, 26(10):2958–2965, 1983.  URL: https://aip.scitation.org/doi/abs/10.1063/1.864060, arXiv:https://aip.scitation.org/doi/pdf/10.1063/1.864060, doi:10.1063/1.864060.
  * [20] D. J. Sigmar, C. T. Hsu, R. White, and C. Z. Cheng.  Alpha‐particle losses from toroidicity‐induced Alfvén eigenmodes. Part II: Monte Carlo simulations and anomalous alpha‐loss processes.  Physics of Fluids B: Plasma Physics, 4(6):1506–1516, 1992.  arXiv:https://doi.org/10.1063/1.860061, doi:10.1063/1.860061.
  * [21] K. Shinohara, Y. Kusama, M. Takechi, A. Morioka, M. Ishikawa, N. Oyama, K. Tobita, T. Ozeki, S. Takeji, S. Moriyama, T. Fujita, T. Oikawa, T. Suzuki, T. Nishitani, T. Kondoh, S. Lee, M. Kuriyama, JT-60 Team, G. J. Kramer, N. N. Gorelenkov, R. Nazikian, C. Z. Cheng, G. Y. Fu, and A .Fukuyama.  Alfvén eigenmodes driven by Alfvénic beam ions in JT-60U.  Nuclear Fusion, 41(5):603–612, may 2001.  doi:10.1088/0029-5515/41/5/313.
  * [22] K. Shinohara, M. Takechi, M. Ishikawa, Y. Kusama, K. Tsuzuki, K. Urata, H. Kawashima, K. Tobita, A. Fukuyama, C. Z. Cheng, D. S. Darrow, G. J. Kramer, N. N. Gorelenkov, R. Nazikian, Y. Todo, Y. Miura, and T. Ozeki.  Energetic particle physics in JT-60U and JFT-2M.  Plasma Physics and Controlled Fusion, 46(7):S31–S45, may 2004.  doi:10.1088/0741-3335/46/7/s03.
  * [23] M. Ishikawa, M. Takechi, K. Shinohara, Y. Kusama, C. Z Cheng, G. Matsunaga, Y. Todo, N. N. Gorelenkov, G. J. Kramer, R. Nazikian, A. Fukuyama, V. A. Krasilnikov, Yu Kashuck, T. Nishitani, A. Morioka, M. Sasao, M. Isobe, and the JT-60 team.  Energetic ion transport by abrupt large-amplitude event induced by negative-ion-based neutral beam injection in the JT-60U.  Nuclear Fusion, 45(12):1474–1480, nov 2005.  doi:10.1088/0029-5515/45/12/002.
  * [24] A. Bierwage, K. Shinohara, Y. Todo, N. Aiba, M. Ishikawa, G. Matsunaga, M. Takechi, and M. Yagi.  Simulation tackle abrupt massive migrations of energetic beam ions in tokamak plasma.  Nature Communications, 9, 2018.  doi:https://doi.org/10.1038/s41467-018-05779-0.
  * [25] E. Lanti, N. Ohana, N. Tronko, T. Hayward-Schneider, A. Bottino, B. F. McMillan, A. Mishchenko, A. Scheinberg, A. Biancalani, P. Angelino, S. Brunner, J. Dominski, P. Donnel, C. Gheller, R. Hatzky, A. Jocksch, S. Jolliet, Z. X. Lu, J. P. Martin Collar, I. Novikau, E. Sonnendrücker, T. Vernay, and L. Villard.  Orb5: A global electromagnetic gyrokinetic code using the PIC approach in toroidal geometry.  Computer Physics Communications, 251:107072, 2020.  URL: https://www.sciencedirect.com/science/article/pii/S0010465519303911, doi:https://doi.org/10.1016/j.cpc.2019.107072.
  * [26] P. Lauber.  Studies and details on the NLED-AUG case.  http://www2.ipp.mpg.de/~pwl/NLED_AUG/data.html.
  * [27] Z. Qiu and L. Chen.  Kinetic theories of geodesic acoustic modes: Radial structure, linear excitation by energetic particles and nonlinear saturation.  Plasma Science and Technology, 13(3):257–266, jun 2011.  doi:10.1088/1009-0630/13/3/01.
  * [28] F. Vannini, A. Biancalani, A. Bottino, T. Hayward-Schneider, P. Lauber, A. Mishchenko, I. Novikau, and E. Poli.  Gyrokinetic investigation of the damping channels of Alfvén modes in ASDEX Upgrade.  Physics of Plasmas, 27(4):042501, 2020.  arXiv:https://doi.org/10.1063/1.5134979, doi:10.1063/1.5134979.
  * [29] F. Vannini, A. Biancalani, A. Bottino, T. Hayward-Schneider, P. Lauber, A. Mishchenko, E. Poli, and G. Vlad.  Gyrokinetic investigation of the nonlinear interaction of Alfvén instabilities and energetic particle-driven geodesic acoustic modes.  Physics of Plasmas, 28(7):072504, 2021.  arXiv:https://doi.org/10.1063/5.0049588, doi:10.1063/5.0049588.
  * [30] J. D. Gaffey.  Energetic ion distribution resulting from neutral beam injection in tokamaks.  Journal of Plasma Physics, 16(2):149–169, 1976.  doi:10.1017/S0022377800020134.
  * [31] A. Mishchenko, A. Könies, R. Kleiber, and M. Cole.  Pullback transformation in gyrokinetic electromagnetic simulations.  Physics of Plasmas, 21(9):092110, 2014.  arXiv:https://doi.org/10.1063/1.4895501, doi:10.1063/1.4895501.
  * [32] A. Mishchenko, A. Bottino, A. Biancalani, R. Hatzky, T. Hayward-Schneider, N. Ohana, E. Lanti, S. Brunner, L. Villard, M. Borchardt, R. Kleiber, and A. Könies.  Pullback scheme implementation in ORB5.  Computer Physics Communications, 238:194 – 202, 2019.  URL: http://www.sciencedirect.com/science/article/pii/S0010465518304181, doi:https://doi.org/10.1016/j.cpc.2018.12.002.
  * [33] G. Vlad, X. Wang, F. Vannini, S. Briguglio, N. Carlevaro, M. V. Falessi, G. Fogaccia, V. Fusco, F. Zonca, A. Biancalani, A. Bottino, T. Hayward-Schneider, and P. Lauber.  A linear benchmark between HYMAGYC, MEGA and ORB5 codes using the NLED-AUG test case to study Alfvénic modes driven by energetic particles.  Nuclear Fusion, 61(11):116026, oct 2021.  doi:10.1088/1741-4326/ac2522.
  * [34] P. Lauber, B. Geiger, G. Papp, G. Por, L. Guimarais, P. Z. Poloskei, V. Igochine, M. Maraschek, G. Pokol, T. Hayward-Schneider, Z. Lu, X. Wang, A. Bottino, F. Palermo, I. Novikau, A. Biancalani, G. Conway, and The ASDEX Upgrade Team.  ”strongly non-linear energetic particle dynamics in ASDEX Upgrade scenarios with core impurity accumulation”, in proceedings of the 27th IAEA Fusion energy, 2018.
  * [35] R. J. Hawryluk.  An empirical approach to tokamak transport.  In Physics of plasmas close to thermonuclear conditions, pages 19–46. Elsevier, 1981.
  * [36] I. Novikau, A. Biancalani, A. Bottino, Ph. Lauber, E. Poli, P. Manz, G. D. Conway, A. Di Siena, N. Ohana, E. Lanti, and L. Villard.  Nonlinear dynamics of energetic-particle driven geodesic acoustic modes in ASDEX Upgrade.  Physics of Plasmas, 27(4):042512, 2020.  arXiv:https://doi.org/10.1063/1.5142802, doi:10.1063/1.5142802.
  * [37] P. Z. Poloskei, G. Papp, G. Pokol, P. Lauber, and L. Horvath.  ”strongly non-linear energetic particle dynamics in ASDEX Upgrade scenarios with core impurity accumulation”, in proceedings of the 44th EPS conference on plasma physics, 2017.
  * [38] G. Y. Fu.  Energetic-particle-induced geodesic acoustic mode.  Phys. Rev. Lett., 101:185002, Oct 2008.  URL: https://link.aps.org/doi/10.1103/PhysRevLett.101.185002, doi:10.1103/PhysRevLett.101.185002.
  * [39] F. Zonca and L. Chen.  Radial structures and nonlinear excitation of geodesic acoustic modes.  EPL (Europhysics Letters), 83(3):35001, jul 2008.  doi:10.1209/0295-5075/83/35001.
  * [40] D. Zarzoso, A. Biancalani, A. Bottino, P. Lauber, E. Poli, J. B. Girardo, X. Garbet, and R. J. Dumont.  Analytic dispersion relation of energetic particle driven geodesic acoustic modes and simulations with NEMORB.  Nuclear Fusion, 54(10):103006, sep 2014.  doi:10.1088/0029-5515/54/10/103006.
  * [41] J. G. Cordey and W. G. F. Core.  Energetic particle distribution in a toroidal plasma with neutral injection heating.  The Physics of Fluids, 17(8):1626–1630, 1974.  URL: https://aip.scitation.org/doi/abs/10.1063/1.1694943, arXiv:https://aip.scitation.org/doi/pdf/10.1063/1.1694943, doi:10.1063/1.1694943.
  * [42] C. Estrada-Mila, J. Candy, and R. E. Waltz.  Turbulent transport of alpha particles in reactor plasmas.  Physics of Plasmas, 13(11):112303, 2006.  arXiv:https://doi.org/10.1063/1.2364149, doi:10.1063/1.2364149.
  * [43] C. Angioni and A. G. Peeters.  Gyrokinetic calculations of diffusive and convective transport of $\alpha$ particles with a slowing-down distribution function.  Physics of Plasmas, 15(5):052307, 2008.  arXiv:https://doi.org/10.1063/1.2913610, doi:10.1063/1.2913610.
  * [44] A. Di Siena, T. Görler, H. Doerk, R. Bilato, J. Citrin, T. Johnson, M. Schneider, and E. Poli.  Non-maxwellian fast particle effects in gyrokinetic gene simulations.  Physics of Plasmas, 25(4):042304, 2018.  arXiv:https://doi.org/10.1063/1.5020122, doi:10.1063/1.5020122.
  * [45] M. N. Rosenbluth, W. M. MacDonald, and D. L. Judd.  Fokker-Planck Equation for an Inverse-Square Force.  Phys. Rev., 107:1–6, Jul 1957.  URL: https://link.aps.org/doi/10.1103/PhysRev.107.1, doi:10.1103/PhysRev.107.1.
  * [46] W. B. Thompson.  An introduction to plasma physics.  American Journal of Physics, 32(1):68–68, 1964.  arXiv:https://doi.org/10.1119/1.1970111, doi:10.1119/1.1970111.
  * [47] P. Helander and D. J. Sigmar, editors.  Collisional Transport in Magnetized Plasmas.  Cambridge University Press, 2002.  URL: https://www.amazon.com/Collisional-Transport-Magnetized-Cambridge-Monographs/dp/0521807980.
  * [48] A. Bottino and E. Sonnendrücker.  Monte carlo particle-in-cell methods for the simulation of the vlasov–maxwell gyrokinetic equations.  Journal of Plasma Physics, 81(5):435810501, 2015.  doi:10.1017/S0022377815000574.
  * [49] P. Lauber, S. Günter, A. Könies, and S. D. Pinches.  LIGKA: A linear gyrokinetic code for the description of background kinetic and fast particle effects on the mhd stability in tokamaks.  Journal of Computational Physics, 226(1):447–465, 2007.
  * [50] B. Rettino, T. Hayward-Schneider, A. Biancalani, A. Bottino, P. Lauber, I. Chavdarovski, F. Vannini, and F. Jenko.  Gyrokinetic modelling of anisotropic energetic particle driven instabilities in tokamak plasmas.  submitted to Nuclear Fusion, 2021.