Cross-code gyrokinetic verification and benchmark on the linear collisionless dynamics of the geodesic acoustic mode

TL;DR

This paper compares analytical theory with gyrokinetic simulations from three different codes to study the linear collisionless geodesic acoustic mode dynamics in tokamaks, considering various plasma parameters and regimes.

Contribution

It provides a detailed benchmark and comparison of gyrokinetic codes against analytical theory for linear GAM behavior in tokamaks, clarifying code performance and validity regimes.

Findings

Good agreement between codes and theory in certain regimes

Identification of code-specific behaviors and limitations

Insights into the linear GAM dynamics under various parameters

Abstract

The linear properties of the geodesic acoustic modes (GAM) in tokamaks are investigated by means of the comparison of analytical theory and gyrokinetic numerical simulations. The dependence on the value of the safety factor, finite-orbit-width of the ions in relation to the radial mode width, magnetic-flux-surface shaping, and electron/ion mass ratio are considered. Nonuniformities in the plasma profiles (such as density, temperature, or safety factor), electro-magnetic effects, collisions and presence of minority species are neglected. Also, only linear simulations are considered, focusing on the local dynamics. We use three different gyrokinetic codes: the lagrangian (particle-in-cell) code ORB5, the eulerian code GENE and semi-lagrangian code GYSELA. One of the main aims of this paper is to provide a detailed comparison of the numerical results and analytical theory, in the regimes where this is possible. This helps understanding better the behavior of the linear GAM dynamics in these different regimes, the behavior of the codes, which is crucial in the view of a future work where more physics is present, and the regimes of validity of each specific analytical dispersion relation.

Tables

Table 1. Table 1: 11 11 11 simulations performed for q = 2 𝑞 2 q=2 with N φ = 8 subscript 𝑁 𝜑 8 N_{\varphi}=8 by varying the number of points in r 𝑟 r , θ 𝜃 \theta , v ∥ subscript 𝑣 parallel-to v_{\parallel} and μ 𝜇 \mu directions. Results are compared on the radial wave number k r ​ ρ i subscript 𝑘 𝑟 subscript 𝜌 𝑖 k_{r}\rho_{i} , the damping rate γ 𝛾 \gamma and the frequency ω 𝜔 \omega of the zonal component of the electrostatic potential.

case	$N_r$	$N_{\theta }$	$N_{v_{\parallel }}$	$N_{\mu }$	$\Delta t{\Omega }_i$	$k_r{\rho }_i$	$\gamma {\Omega }_i^{-1}$	$\omega {\Omega }_i^{-1}$
1	$256$	$64$	$128$	$32$	$5.$	$0.05630867$	$0.02664329$	$1.81126121$
2	$256$	$64$	$128$	$32$	$10.$	$0.05630867$	$0.02663145$	$1.81126121$
3	$256$	$64$	$128$	$32$	$25.$	$0.05630868$	$0.02657794$	$1.81126121$
4	$256$	$64$	$128$	$32$	$50.$	$0.05630874$	$0.02653612$	$1.81196116$
5	$256$	$64$	$128$	$16$	$25.$	$0.05630058$	$0.02657711$	$1.81126121$
6	$256$	$64$	$128$	$8$	$25.$	$0.05601206$	$0.0265876$	$1.81126121$
7	$128$	$64$	$128$	$32$	$25.$	$0.05649958$	$0.0266591$	$1.81126121$
8	$512$	$64$	$128$	$32$	$25.$	$0.0562041$	$0.02655708$	$1.81126121$
9	$256$	$256$	$128$	$32$	$25.$	$0.05630868$	$0.02657654$	$1.81126121$
10	$256$	$64$	$64$	$32$	$25.$	$0.05630868$	$0.02655779$	$1.81196116$
11	$128$	$64$	$64$	$8$	$50.$	$0.05620426$	$0.02661145$	$1.81126121$

Equations

$$\dot{\mathbf{R}}$$

$$\dot{p_{\Vert }}$$

$$\dot{\mu }$$

$$\dot{\mathbf{R}}$$

$$\dot{p_{\Vert }}$$

$$n_{e}({\bf R},t)=n_{e0}+\frac{q_{s}n_{e0}}{T_{e0}}\big{(}\phi-\bar{\phi}\big{)}$$

$$-\nabla \cdot \frac{n_0{m}_i{c}^2}{B^2}{\nabla }_{\perp }\varphi =\sum _s\int \mathrm{d}W{q}_s\widetilde{\delta f_s}$$

$$\frac{\partial f_s}{\partial t}+\frac{d\mathbf{R}}{dt}\cdot \nabla f_s+\frac{d{v}_{\Vert }}{dt}\frac{\partial f_s}{\partial v_{\Vert }}=0$$

$$\frac{d\mathbf{R}}{dt}$$

$$\frac{d{v}_{\Vert }}{dt}$$

$$\sum _s{q}_s\delta n_s=0$$

$$\sum _s\frac{2\pi q_s}{m_s}\int B_{G\Vert }^{\ast }\widetilde{\delta f_s}\mathrm{d}{v}_{\parallel }\mathrm{d}\mu -\frac{q_s^2{n}_{0s}}{T_{0s}}\left[\varphi -\frac{B_0}{T_{0s}}\int \widetilde{\stackrel{~}{\varphi }}\exp \left(-\frac{\mu B_0}{T_{0s}}\right)\mathrm{d}\mu \right]=0,$$

$$-\frac{e^2{n}_{0e}}{T_{0e}}\left(\varphi -\bar{\varphi }\right)+\sum _{s\ne e}\frac{2\pi q_s}{m_s}\int B_{G\Vert }^{\ast }\widetilde{\delta f_s}\mathrm{d}{v}_{\parallel }\mathrm{d}\mu -\frac{q_s^2{n}_{0s}}{T_{0s}}\left[\varphi -\frac{B_0}{T_{0s}}\int \widetilde{\stackrel{~}{\varphi }}\exp \left(-\frac{\mu B_0}{T_{0s}}\right)\mathrm{d}\mu \right]=0.$$

$$-\frac{1}{n_{e_0}}\sum _s{q}_s{\nabla }_{\perp }\cdot \left(\frac{n_{s0}}{B{\Omega }_s}{\nabla }_{\perp }\varphi (\mathbf{R},t)\right)+e\left(\frac{\varphi -\bar{\varphi }}{T_e}\right)=\frac{1}{n_{e_0}}\sum _s{q}_s(n_s-n_{s0})$$

$$\Lambda^{2}(z)=z^{2}+q^{2}\omega_{ti}^{2}z\Big{[}F(z)-\frac{N^{2}(z)}{D^{2}(z)}\Big{]}=0$$

$$F(z)$$

$$N(z)$$

$$D(z)$$

$$Z(z)={\pi }^{-1/2}{\int }_{-\infty }^{+\infty }\frac{e^{-y^2}}{(y-z)}dy$$

$$\frac{\omega }{q{\omega }_{ti}}$$

$$\frac{\gamma }{q{\omega }_{ti}}$$

$$\frac{\omega }{q{\omega }_{ti}}$$

$$\frac{\gamma }{q{\omega }_{ti}}$$

$$\frac{\omega }{q{\omega }_{ti}}$$

$$\frac{\gamma }{q{\omega }_{ti}}$$

$$\frac{\omega }{q{\omega }_{ti}}$$

$$\frac{\gamma }{q{\omega }_{ti}}$$

$$Q=\frac{f_{Q2}}{f_{T}^{2}}\Big{(}\frac{e^{2}+1}{2}\Big{)}-\frac{f_{G1}}{f_{T}}e^{2}-\frac{f_{G2}}{f_{T}}+\frac{f_{G3}(e^{4}+1)+f_{G4}e^{2}}{e^{2}(e^{2}+1)}$$

$$\frac{\omega }{q{\omega }_{ti}}$$

$$\frac{\gamma }{q{\omega }_{ti}}$$

$$k_{r}\rho_{i}=\rho_{i}\sqrt{\sum_{i=0}^{N_{r}}\left(\bar{\phi}_{\rm norm}(r_{i},t)-\langle\bar{\phi}_{\rm norm}\rangle_{r}\right)^{2}}\bigg{/}\sqrt{\sum_{i=0}^{N_{r}}\left(\frac{d}{dr}\bar{\phi}_{\rm norm}(r_{i},t)\right)^{2}}$$

Full text

1 Introduction

Turbulence in tokamak plasmas is often observed accompanied by zonal, i.e. axisymmetric, radial electric fields, giving rise to zonal poloidal flows. Two kinds of zonal flows are identified: zero-frequency zonal flows (ZFZF) [1, 2, 3] and oscillating zonal flows, named geodesic acoustic modes (GAM) [4, 5] GAMs have mainly zonal polarization of the perturbed electric field, i.e. n=0, m=0, with n and m being respectively the toroidal and poloidal magnetic field, and n=0, m=1 perturbed density. Their characteristic frequency is of the order of the sound frequency $\omega\sim c_{s}/R_{0}$ (where $c_{s}=\sqrt{T_{e}/m_{i}}$ is the sound speed and $R_{0}$ is the tokamak major radius, with $m_{i}$ being the ion mass and $T_{e}$ the electron temperature). The importance of understanding the dynamics of these zonal structures in tokamaks is due to their nonlinear interaction with turbulence, being crucial for its saturation [6, 7, 8, 9, 10]. As an example of some recent experimental observations of GAMs in ASDEX Upgrade, see Ref. [11].

In this work, we focus on the linear properties of GAMs, and investigate the dependence on the safety factor, the effect of the finite-larmor-radius (FLR) and finite-orbit-width (FOW) of the ions in relation to the GAM radial width, the effect of the magnetic-flux-surface shaping, and of the electron/ion mass ratio. Nonuniformities in the plasma profiles (such as density, temperature, or safety factor), magnetic effects, collisions and effects of plasma minorities (such as bulk ion impurities or energetic ions) are neglected here. Our investigation is done by means of analytical theory and numerical simulations.

A great progress in the analytical investigation of the linear GAM dynamics has been achieved, starting with its first estimate in ideal MHD [4], then in kinetic theory by neglecting the effect of the FOWs of the passing ions [12, 5], then including it to the first order [13, 14, 15, 5], then including it to higher orders [16, 17], and then including the effect of the flux surface shape [18, 19]. All such analytical models neglect the effect of the finite mass of the electrons. In fact, all analytical theories derived so far treat the $m\neq 0$ component of the electrons as adiabatic (whereas the m=0 component of the electron density perturbation is imposed to zero). We will refer to this model for treating the electrons as “adiabatic”. The importance of having an analytical description is twofold. On the one hand, it allows a direct understanding of the physical mechanisms leading to the each different effect under investigation. On the other hand, it allows a detailed linear verification process of the numerical tools, which is at the basis of the development of gyrokinetic codes aimed at a rigorous turbulence investigation.

Many numerical investigations of the linear GAM dynamics and comparison with analytical theory or benchmark among codes have been carried out in the last few decades, most of which treating the electrons as adiabatic. As a non-extensive list of example, we mention here simulations performed with the gyrokinetic codes GTC [20, 21], ORB5 [22] (where the effect of the elongation was studied), TEMPEST [16] (where the effect of high-order terms of the finite ion orbit width was studied), GYRO [23] (where the effect of finite orbit width and the application to the radial velocity in the large-q limit was studied), GYSELA [24], ELMFIRE [25], and GENE with GKW [26]. In particular, a first verification of ORB5 against analytical theories, for circular geometries and low values of $k_{r}\rho_{i}$, was started in Ref. [27]. A numerical study of the effect of kinetic electrons in circular plasmas has been described in Ref. [28].

In this paper, we aim at performing a comprehensive cross-code verification and benchmark of several gyrokinetic codes, in different regimes. We perform numerical simulations with three different gyrokinetic codes, adopting equivalent physical models for the dynamics of the ions (which is the most basic species to be investigated for the physics of sound waves, and therefore of GAMs), but solving the model equations in three different ways: the Lagrangian (i.e. particle-in-cell) code ORB5 [29, 30, 31], the Eulerian code GENE [32, 33], and the semi-Lagrangian code GYSELA [34, 35]. All these codes solve the ion dynamics based on the gyrokinetic equations (see for example Ref. [36, 37, 38, 39] for some early derivations, or Ref. [40] for a recent comprehensive review). Some of the differences in the physical models of these codes are the treatment of the electrons and the possibility of investigating the dynamics in elongated magnetic equilibria.

The gyrokinetic model has been also adopted in the past for deriving analytical dispersion relations for linear collisionless GAMs in various regimes (small or large values of safety factors, negligible, small or moderate values of normalized radial wave-number, small or moderate elongation of the equilibrium magnetic flux surfaces). Testing a numerical code, by using it in a particular limit where the analytical solution is known, takes the name of code “verification”. Testing different numerical codes, which solve the basic equations in different ways, by using them in particular regimes where the hypothesis of the physical models are the same, takes the name of code “benchmark”. The verification and benchmark of gyrokinetic codes has been the main goal of an European effort developed in the last three years as an Eurofusion project. This project focuses on different physical phenomena observed in tokamaks, like for example GAMs, and ion-temperature-gradient instabilities, both crucial actors to be understood in the view of the prediction and control of the turbulent transport in existing tokamak devices and future fusion reactors. This paper in particular, summarises the results of the verification and benchmark on GAMs (see Ref. [41] for a description of the benchmark on ITGs).

The paper is organized as follows. The numerical models of the three codes adopted for this study are described in Sec. 2. The dependence of the GAM dynamics on the safety factor, the effect of ion FOW and flux-surface elongation is discussed in Sec. 3, where the electrons are treated as adiabatic. The effect of the finite mass of electrons is described in Sec. 4, where the results of numerical simulations with kinetic electrons are shown. The conclusions are summarised in Sec. 5. The investigation of the numerical convergence of the codes is presented in the appendices.

2 The models

The choice of the model for the investigation of the dynamics of GAMs is dictated by their specific spatial and temporal characteristic scales. In particular, GAMs are zonal oscillations (i.e. with toroidal and poloidal mode numbers equal to zero) with radial wavelength bigger than (or in same case of the order of) the ion Larmor radius, and time scales of the order of the sound time $\sim 2\pi R/c_{s}$. The basic physics is the one of the sound waves, therefore the MHD description is sufficient for estimating the order of magnitude of the frequency of the mode [4]. Nevertheless, such spatio-temporal scales make the need of a kinetic treatment clear. This is due to the importance of resonances with passing ions, which can determine frequency and damping rate of the modes. Considering the requirements for the spatial scales, and the frequencies being lower than the ion cyclotron frequency, we can easily see that the gyrokinetic model [36, 37, 38, 39, 42, 43, 40] is the most appropriate tool.

2.1 The numerical model of ORB5

ORB5 is a nonlinear gyrokinetic code based on a particle-in-cell (PIC) algorithm. The basic discretization scheme of a PIC code (also known as “Lagrangian” code), for the Vlasov-Maxwell problem, is presented in Ref. [44]. A PIC code discretizes the distribution function with macro-particles, also known as markers, associated with weights. In a gyrokinetic PIC code, the markers are pushed along the trajectories derived from the gyrokinetic model while the fields are known on a spatial grid and evolved by solving the gyrokinetic field equations either with finite differences or with finite element methods. The sources (charge density and current density) needed for solving the field equations are calculated by projecting the marker weights on a spatial grid. In ORB5, the distribution function is decoupled in a background distribution analytically known, while only the perturbation is solved using markers, with a control-variate Monte Carlo method, hystorically known as $\delta$f PIC [45] (see Ref. [31] for a recent overview).

ORB5 was originally developed for electrostatic turbulence studies [29]. In the last few years, it has been extended to the electromagnetic, multi-species version within the NEMORB project [30, 31]. In this paper, only the linearized electrostatic model of ORB5 is used. Only collisionless simulations are considered. Although only the local GAM dynamics is of interest in this paper, no flux-tube version of ORB5 exists, therefore only global simulations are considered, and the global effects are neglected (see Ref. [46] for an investigation of the global effects with ORB5). The model equations of ORB5 are derived in a Lagrangian formulation [31], based on the gyrokinetic Vlasov-Maxwell equations of Sugama, Brizard and Hahm [43, 40].

The gyrocenter trajectories describe the motion of the markers of the kinetic species in phase-space coordinates written in $p_{\|}$-formalism, ${\bf Z}_{p_{\|}}=({\bf R},p_{\|},\mu)$, i.e. respectively the gyrocenter position, canonical parallel momentum $p_{\|}=m_{s}v_{\|}+(q_{s}/c)\tilde{A}_{\|}$ and magnetic momentum $\mu=m_{s}v_{\perp}^{2}/(2B)$ (with $m_{s}$ and $q_{s}$ being the mass and charge of the species). $v_{\|}$ and $v_{\perp}$ are respectively the parallel and perpendicular component of the particle velocity. The gyroaverage operator is labeled here by the tilde symbol $\tilde{}$ . The gyroaverage operator reduces to the Bessel function $J_{0}$ if we transform into Fourier space. In all simulations with ORB5 shown in this paper, the gyroaverage is always calculated by considering non-vanishing Larmor radius for the ions, whereas it is calculated with zero argument for the electrons. In other words, finite-Larmor-radius (FLR) effects are retained for ions and neglected for electrons. The code ORB5 is based on straight-field-line tokamak coordinates. Dirichlet boundary conditions are imposed in the radial direction, while periodicity is assumed in the two angles. The nonlinear electromagnetic version of the trajectories is [31]:

[TABLE]

Here, the time-dependent fields are the scalar potential $\phi$ and the parallel component of the vector potential $A_{\|}$, and ${\bf B}^{*}={\bf B}+(c/q_{s}){\bf\nabla}\times({\bf b}\,p_{\|})$, where ${\bf B}$ and ${\bf b}$ are the equilibrium magnetic field and magnetic unitary vector. The linearization of the Vlasov equation is performed by pushing the markers along unperturbed trajectories:

[TABLE]

In this paper, the trajectories given by Eq. 4, 5 are always calculated for the ions (which are always treated kinetically), whereas the electrons can be either treated kinetically (by considering $J_{0}=1$ and neglecting the electron polarization) or with an “adiabatic” model, where the electron gyrocenter density is calculated directly from the value of the scalar potential as [31]:

[TABLE]

where $\bar{\phi}$ is the flux-surface averaged potential. The quantities with subscript “0” refer to the equilibrium, and therefore are functions of the radial coordinate $\rho$ only.

The equation for solving the scalar potential is the gyrokinetic Poisson equation, also known as polarization equation. This is derived from the gyrokinetic Lagrangian of ORB5, using the variational derivation, and imposing that the ExB drift energy of the particles is larger than the field energy (quasi-neutrality condition) [31]. The gyrokinetic Poisson equation is [31]:

[TABLE]

with $n_{0}m_{i}$ being here the total plasma mass density (approximated as the ion mass density), and the summation over the species is performed when the electrons are treated as kinetic, otherwise the electron contribute is given by $-n_{e}({\bf R},t)$. Here $\delta f=f-f_{0}$ is the gyrocenter perturbed distribution function, with $f$ and $f_{0}$ being the total and equilibrium (i.e. independent of time, assumed here to be a Maxwellian) gyrocenter distribution functions. The integrals are over the phase space volume, with $\mathrm{d}W=(2\pi/m_{i}^{2})B_{\|}^{*}\mathrm{d}p_{\|}\mathrm{d}\mu$ being the velocity-space infinitesimal. The gyrokinetic Poisson equation is solved with a finite element method, by using B-splines in all the spatial directions.

Eqs 4, 5, 6, 7 are the constitutive equations of the model of ORB5 used in this paper for studying the collisionless electrostatic linearized dynamics of GAMs. In the electromagnetic version, the Ampère equation is also solved for calculating the time evolution of the parallel component of the vector potential $A_{\|}$, which is neglected in this paper.

2.2 The numerical model of GENE

The Gyrokinetic Electromagnetic Numerical Experiment (GENE) code, is also a nonlinear gyrokinetic code originally developed for electromagnetic turbulence studies in the flux-tube (i.e. local) limit [32], recently extended to its global representation [33]. The model of GENE is also based on the gyrokinetic Vlasov-Maxwell equations of Brizard and Hahm [40]. Intra- and inter-species collisions (both pitch angle and energy scattering) are implemented. In this paper, only the linearized electrostatic collisionless version of GENE is used.

GENE is a Eulerian code. In a Eulerian description, the distribution function is not discretized with markers, but it is discretized on a 5D fixed grid in phase-space. The gyrokinetic equation is then solved on this grid for each species $s$. The coordinate system of GENE in the 5D phase space is written in $v_{\|}$-formalism, ${\bf Z}_{v_{\|}}=({\bf R},v_{\|},\mu)$, i.e. respectively the gyrocenter position, parallel velocity and magnetic momentum. GENE adopts a field-aligned coordinate system to represent the fluctuation fields in the configuration space of ${\bf R}$. This coordinate system becomes singular at the magnetic axis which therefore cannot be simulated. In the local version of the code, the radial direction is Fourier transformed and periodic boundary conditions are applied. In the global version the radial direction is instead treated in real space and Dirichlet boundary conditions are applied. The binormal direction (i.e. perpendicular to the radial direction and to the equilibrium magnetic field) is always Fourier transformed as axisymmetry corresponds to invariance in this direction, and each linear mode corresponds to a toroidal mode number $n$.

The distribution function $f_{s}$ of each species is evolved accordingly to the gyrokinetic equation in the form [33]:

[TABLE]

where the equations of motion of the gyrocenters are given by [33]:

[TABLE]

Here ${\bf B}^{*}_{G}={\bf B}+(m_{s}c/q_{s}){\bf\nabla}\times({\bf b}\,v_{\|})$, the generalized ExB drift is ${\bf v}_{E}=(c/B^{2}){\bf B}\times{\bf\nabla}(\tilde{\phi}-(v_{\|}/c)\tilde{A}_{\|})$, the grad-B drift is ${\bf v}_{\nabla B}=(\mu c/q_{s}B^{2}){\bf B}\times{\bf\nabla}B$, and the curvature drift is ${\bf v}_{c}=(v^{2}_{\|}/\Omega_{s})(\nabla\times{\bf b})_{\perp}$. In the electrostatic version of the code, used in this paper, the ExB drift is ${\bf v}_{E}=(c/B^{2}){\bf B}\times{\bf\nabla}\tilde{\phi}$, and the second term on the right hand side of Eq. 10 is dropped, so that the equation of the time derivative of the parallel component of the velocity takes the form:

[TABLE]

The linearization in GENE is done by plugging in the equation of motion in the Vlasov equation and neglecting all the nonlinear terms. Only linear simulations are considered in this paper.

Equation (8) is then solved self-consistently with the gyrokinetic Maxwell equations, which in the cases considered here reduce to the gyrokinetic Poisson equation, Eq. 7, which is solved for obtaining the scalar potential (whereas the Ampère equation can also solved in case of electromagnetic simulations). As in the ORB5 code, different models are available for describing each species dynamics. In this paper ions are always assumed to be fully gyrokinetic whereas electrons, depending on the particular case being simulated, are treated either as a second kinetic species or assumed to respond adiabatically. For typical tokamak parameters the Debye length is much smaller than the characteristic wave-length of microinstabilities. The gyrokinetic Poisson equation can thus be reduced to the quasi-neutrality condition, which, having assumed a quasi-neutral background, reads

[TABLE]

where $\delta n_{s}$ indicates the perturbed gyrocenter density of the $s$-th species, obtained from the gyrokinetic model. When all species are treated kinetically, the equation for quasi-neutrality, Eq. (12) can be rewritten as

[TABLE]

while in case of adiabatic electrons it reduces to

[TABLE]

As mentioned above, the linear physical models of ORB5 and GENE are equivalent (see Ref. [50] for a detailed discussion on the comparison of the two models), and no difference in the results is expected for the linear collisionless GAM dynamics, depending on that. Nevertheless, the numerical schemes are different. Moreover, the existence of the two representations of GENE, namely the global and the local (i.e. flux-tube) representations, offers the possibility to solve the model equations in two independent ways. As shown in the following sections, no difference is found in the results, for the chosen tests. This means that, for these particular cases, the local dynamics is dominant and independent on the adopted numerical scheme.

2.3 The numerical model of GYSELA

Like ORB5 and GENE, the GYrokinetic SEmiLAgrangian code (GYSELA) is also a nonlinear 5D gyrokinetic code [35]. No linear version exists, therefore nonlinear simulations are considered in this paper, but with sufficiently small initial perturbation, in order to focus on the linear dynamics. The GYSELA code is dedicated to electrostatic Ion Temperature Gradient (ITG) turbulence with possibility to address transport of impurities. Electrons are presently assumed adiabatic but a kinetic version is under development. GYSELA is a global full-$f$ flux-driven code which addresses turbulent and neoclassical transports on an equal footing.

GYSELA is a global code with a toroidal geometry with a simplified concentric circular magnetic configuration. Its coordinate system in the 5D space is written as GENE in $v_{\|}$-formalism, ${\bf Z}_{v_{\|}}=({\bf R},v_{\|},\mu)$ but where $R=(r,\theta,\varphi)$ with $r$ the radial direction and $\theta$ and $\varphi$ the poloidal and toroidal geometric angles. Boundary conditions are periodic in $\theta$ and $\varphi$ directions. Non-axisymmetric fluctuations of the electric potential and of the distribution function - i.e $(m,n)\neq(0,0)$ modes - are forced to zero at both radial boundaries of the simulated domain. As far as the axisymmetric component is concerned, the value of the potential is prescribed at the outer boundary, while the radial electric field is set to zero at the inner boundary. No flux-tube version of GYSELA exists, but since in this paper only local physics is concerned, density, temperature and safety factor profiles will be considered constants to minimize the global effects. GYSELA is a full-$f$ code, namely the back reaction of turbulent transport is accounted for in the time evolution of the equilibrium. In such a framework, the turbulence regime is evanescent if no free energy is injected in the system. A flux-driven version of the code is available since 2009 [52], where the system can be driven by a prescribed volumetric source, versatile enough to allow for separate injection of heat, parallel momentum and vorticity. However in this paper, the temperature and density profiles are constant and therefore we only use the forcing governed by the two equal thermal baths at the two radial boundaries. A linearized multi-species collision operator is implemented in the code [51] but here only collisionless simulations are considered. No filters in the toroidal mode number are imposed in these simulations with GYSELA, therefore all components are allowed to develop. As shown in the following sections, the results of GYSELA are found to be in good agreement of those obtained with codes which use a linearized version of the model equations and filter out the non-zonal component. This means that, for the tests chosen in this paper, the nonlinear excitation of non-zonal components is negligible and does not sensibly modify the evolution of the zonal component.

The numerical scheme of GYSELA is based on a semi-Lagrangian method [48] (more specifically on a “backward semi-Lagrangian approach”), which is a mix between PIC and Eulerian methods exhibiting good properties of conservation [34]. In this method, the phase-space mesh grid is kept fixed in time (like in Eulerian codes) and the Vlasov equation is integrated along the trajectories (like in PIC codes) using the invariance of the distribution function along the trajectories (Liouville theorem). In GYSELA, the interpolation step is presently performed with cubic splines.

Like for ORB5 and GENE, the model equations of GYSELA are based on the gyrokinetic equations of Brizard and Hahm [40]. Then, the time evolution of the full guiding-center distribution function $F_{s}$ is governed for each species $s$, by the same form of equation as the one of GENE i.e. Eq. 8 where the characteristics, i.e. the trajectories of the gyrocenters, are given by Eqs. 9 and 11. These 5D gyrokinetic Vlasov equations are self-consistently coupled to a 3D quasi-neutrality equation defined as:

[TABLE]

where the gyrocenter density $n_{s}$ of species $s$ reads $n_{s}(\mathbf{R},t)=\int{\cal J}_{\rm v}\>{\rm d}\mu\>{\rm d}v_{\parallel}\;\tilde{F}_{s}(\mathbf{R},\mathbf{v},t)$ with ${\cal J}_{\rm v}=2\pi B_{G\|}^{*}/m_{s}$ the jacobian in velocity space. The equilibrium gyrocenter density $n_{s0}$ corresponds to the same expression as $n_{s}$ where $F_{s}$ is replaced by the equilibrium Maxwellian $F_{s0}$. The gyroaverage operator was historically approximated by a Padé expansion but in this paper the new version based on direct integration on the gyro-circles with Hermite interpolation is used. In GYSELA, the quasi-neutrality equation Eq. 15 is solved with finite differences in radial direction and Fourier projection in $\theta$ direction ($\varphi$ plays the role of a parameter). See Appendix A in Ref. [35] to see how the presence of $\bar{\phi}$ is overcome.

3 Numerical simulations with adiabatic electrons

In this section, the results of numerical simulations of GAMs with adiabatic electrons are discussed, and compared with analytical theory. The main aim here is to perform a detailed verification and benchmark of the different gyrokinetic codes. This has the triple role of: a) understanding better the behavior of the linear GAM dynamics in different regimes; b) understanding better the behavior of the codes, which is crucial in the view of a future work where more physics is present; c) understanding better the regimes of validity of each specific analytical dispersion relation. Two main regimes are considere: one where the GAM radial size is large with respect to the ion larmor radius, and therefore the FOW effects are smaller, and one where the GAM radial structure is finer, and therefore the FOW effects are larger.

3.1 GAMs with broad radial structure

3.1.1 Analytical predictions

In the case of GAMs with broad radial structure ($k_{r}\rho_{i}\ll 1$), an analytical theory neglecting the FLR and FOW corrections can be considered as a first approximation. Although an MHD theory would be sufficient for estimating the order of magnitude of the GAM frequency, nevertheless, due to the resonances with ions, a gyrokinetic treatment of the ions is necessary for a proper estimation of the GAM frequency and damping rate. Such a dispersion relation in the case of circular flux surfaces has been provided by F. Zonca in 1996 [12] in the general electro-magnetic case, for low-frequency Alfvén modes, and can be adopted for GAMs as well, when neglecting diamagnetic effects, as discussed in details in Ref. [5]. No resonances of the electrons are retained. It reads:

[TABLE]

where $z=(\omega+i\gamma)/\omega_{ti}$, $\omega_{ti}=\sqrt{2}\,v_{ti}/(qR_{0})$ is the transit ion frequency, $v_{ti}=\sqrt{T_{i}/m_{i}}$, and the functions $F$, $N$ and $D$ are defined by:

[TABLE]

where $\tau_{e}=T_{e}/T_{i}$ is the ratio of electron over ion temperatures, and $Z(z)$ is the plasma dispersion function:

[TABLE]

Eq. 16 is the desired dispersion relation. It is in implicit form, i.e. the zeroes of the function $\Lambda(z)$ must be found in the complex plane.

For shorter wavelengths and/or larger $q$, FLR/FOW effects become more important [17], and higher order transit resonances play a more important role in the Landau damping of GAMs, in addition to the modification of their real frequency. An extension of equation (16) to the case where FLR and FOW effects are also considered to the first order (still in circular geometry, and with adiabatic electrons), was made for general low-frequency Alfvén modes by F. Zonca in 1998 [13]. An approximated explicit formula for the frequency and damping rates of GAMs with $\omega=l\omega_{ti}$, $l=\pm 1,\pm 2$ transit resonances accounted for was provided by H. Sugama in 2006 [14] and 2008 [15], in the regime of moderate values of $q$:

[TABLE]

with $x=\omega/\omega_{ti}=\Re(z)$, $\rho_{i}=\sqrt{2T_{i}/m_{i}}/\Omega_{i}$, $f_{T}=7/4+\tau_{e}$ and $f_{S1}=23/8+2\tau_{e}+\tau_{e}^{2}/2$, $f_{S2}=(1+\tau_{e})/16$, and $f_{S3}=3/8+7\tau_{e}/16+5\tau_{e}^{2}/32$ (with $\Omega_{i}$ being the ion cyclotron frequency). Note that, in the limit of large values of $q$ (i.e. $q>4$) and $\tau_{e}=1$, the normalized frequency, Eq. 21, tends to $\omega/q\omega_{ti}\simeq f_{T}^{1/2}=1.66$, and $x\simeq f_{T}^{1/2}q=1.66\,q$. For cases with large enough values of $q$ to satisfy $\omega\simeq 2\omega_{ti}$, the second term in equation (22) (namely the one proportional to $k_{r}^{2}\rho_{i}^{2}$) becomes dominant since the ions with lower energy (and thus, which are present in larger number) resonate with the GAM frequency. Note that FLR corrections are not included in Eq. 22.

The effect of elongation $e$, in a gyrokinetic treatment, has been included by Z. Gao in 2009 [18], in the large aspect ratio limit, and neglecting the FLR/FOW effects. The resulting GAM frequency and damping rate, where we neglect here the effect of the radial derivative of the elongation because not taken into account in our paper, are111a typo was present in the original paper, due to a missing proper normalization of $\omega$ in the formula for the damping rate.:

[TABLE]

Note that, for $e=1$, i.e. for circular flux surfaces, and for large values of $q$, the frequency given by Gao-2009, Eq. 23, reduces to the one of Sugama-2008, Eq. 21. Also note that, for deriving the damping rate of Gao-2009, Eq. 24, one has to neglect the second of the two terms of the damping rate given by Sugama-2008, Eq. 22 (which means assuming that the values of $q$ are below a certain threshold) and at the same time assuming the limit of large values of $q$ (i.e. large values of $x$). Due to these strong approximations, we expect the formula for the damping rate of Gao-2009 to give a good qualitative comparison with the results of numerical simulations, but some divergence in the absolute values are not to be surprising.

The previous dispersion relations, namely Eq. 16 given by F. Zonca in 1996, Eqs. 21 and 22 given by H. Sugama in 2006 and 2008, and Eqs. 23 and 24 given by Z. Gao in 2007, are considered as a reference for the comparison with all the results of numerical simulations on GAMs with broad radial structure, shown in Sec. 3.1. Separate dispersion relations, where higher-order FLR/FOW effects are taken into account, are discussed in Sec. 3.2 and used for comparison with the results of numerical simulations shown in the same section.

3.1.2 Equilibrium and definition of the simulation

For our numerical test, we choose a tokamak equilibrium with high aspect ratio ($\varepsilon=a/R=0.1$), with $R_{0}=1.3$ m, $a=0.13$ m. The equilibrium magnetic field is given by ${\bf B}=(B_{0}R_{0}/R)({\bf e}_{\varphi}+(r/qR_{0}){\bf e}_{\theta})$. The value of the magnetic field on axis is $B_{0}=1.9$ T. Each simulation has a different $q$ profile, flat, and each one with different value of $q$. Flat temperature and density profiles are also always considered. The value of $\rho^{*}=\rho_{s}/a$ is chosen as $\rho^{*}=1/160$ for all simulations shown in Sec. 3.1 (with $\rho_{s}=c_{s}/\Omega_{i}$ being the sound Larmor radius). The value of the density is irrelevant for electrostatic simulations.

We initialize a charge density perturbation with only zonal component (i.e. independent of the poloidal and toroidal angle), and generating a scalar potential with a sine dependence on the radius, of the form $\phi(\rho,t=0)=\sin(k_{r}a\rho)$, with $k_{r}=2\pi/a$ (where $\rho=r/a$ is the normalized minor radius, with values in [0,1]). In GYSELA, the initial perturbation of the distribution function is chosen such that $\phi(\rho,t=0)=1-\cos(k_{r}a\rho)$. This choice has been preferred because the radial profile of the zonal component stays more stable in time than for the case $\phi(\rho,t=0)=\sin(k_{r}a\rho)$, leading to a mean $k_{r}$ value in time closer to the initial one. This is particularly true for large values of $k_{r}$ as those explored in section 3.2. One explanation could be that in the case of $1-\cos$ profile the gradients are flatter at radial boundaries so that boundary conditions seem to have less impact. Anyway, this raise the delicate point of confronting global nonlinear code results with linear theory results which are based on local approximation. With this choice of $\rho^{*}$ and $k_{r}$, we obtain a relatively low value of $k_{r}\rho_{i}$, namely $k_{r}\rho_{i}=0.055$ which corresponds to a regime where ion FOW effects are relatively small, for moderate values of $q$. The perturbation is let evolve in a linear electrostatic simulation with adiabatic electrons. GAMs oscillations are observed, and we measure the scalar potential, and calculate frequency and damping rate.

3.1.3 Dependence on the safety factor

For the simulations shown in this section, we have considered an analytical equilibrium with concentric circular flux surfaces. The radial electric field is measured at the radial position of its peak, namely at mid-radius, $\rho=0.5$. It is observed to oscillate in time, and be damped due to Landau damping. The frequency is measured for different simulations with different value of q, obtained with ORB5, GENE and GYSELA, and it is found to scale correctly with the theoretical dispersion relation Zonca-1996 (Eq. 16), and the explicit formula Sugama-2006 (Eq. 21), as shown in Fig. 1. In particular, note that the value of the frequency tends to $\omega_{q\rightarrow\infty}/q\omega_{ti}=1.66$ for large values of $q$, as discussed in Sec. 3.1.1 after Eq. 21. Some minor differences are found at low values of $q$ in the results of the dispersion relations (due to the hypothesis of large $q$ considered by Sugama for the calculation of the explicit formula).

The dependence of the damping rate on $q$ has also been studied, for the same simulations (see Fig. 1). All codes match well with the analytical predictions of Zonca-1996 (Eq. 16) at low values of $q$ ($q\leq 1.5$), where the FOW effects are negligible. At larger values of $q$ ($q>1.5$), the FOW effects included at the first order in the explicit formula Sugama-2008 (Eq. 22) are shown to be dominant. All codes fit well with Sugama-2008 for values of $q$ smaller than 3.5. At even larger values of $q$ ($q\geq 3.5$) the higher-order FOW effects become dominant, and deviations from the formula Sugama-2008 are observed. This regime is studied more in details in Sec. 3.2. The difference at large q between the flux-tube version of GENE (which agrees perfectly with ORB5) and the global version of GENE is due to the fact that the $k_{r}$ used for global GENE runs was slightly larger. For this value of $\rho^{*}$, the choice of $k_{r}\rho_{i}$ = 0.055 requires to simulate the entire domain in minor radius, while simulations of global GENE accounted only for 98% of it. This affects only the high q, i.e. when the damping is very small and the relative effect of $k_{r}$ is large. Values of damping rates larger than ORB5 at large values of q are also observed with GYSELA, probably because the value of $k_{r}$ has been observed to evolve in time towards values which are a bit larger than at the initial time of the simulations with GYSELA, and this increases the averaged damping rate.

3.1.4 Dependence on the elongation

The dependence on the elongation has been studied by loading magnetic equilibria with different elongation (and no triangularity) calculated with the CHEASE code [49]. These simulations have been performed with ORB5 and GENE. The safety factor has been chosen with $q=1.4$, and we have varied the elongation from $e=1$ (circular flux surfaces) to $e=1.75$ (elongated plasmas). The results are shown in Fig. 2.

The frequency measured with the two codes has been found to fit very well, falling within the error bars for the whole scan. The fit with the analytical prediction of Gao-2009 is also very good (with a maximum of 3% of difference not depending on the elongation), showing the correct decrease of the frequency with the increasing elongation.

Regarding the damping rate, a very good matching of the two codes has also been found, showing an increase of the damping rate with the elongation, as predicted by the analytical theory. A quantitative fit of the damping rate with the analytical theory has been found worse than for the frequency. This is probably due to the fact that the dependence of the dispersion relation on the safety factor $q$, the FOW effects, proportional to $k_{r}^{2}\rho_{i}^{2}$, and the elongation $e$ at the same time, forces some strong approximations to be taken when deriving an explicit analytical formula, as discussed in Sec. 3.1.1. Therefore, the analytical derivation is based on some assumptions, like the assumption of negligible FOW effects, and of large values of $q$ at the same time, which is most likely at the origin for the divergences with the results of the numerical simulations. Note that this difference, up to 40%, of the result of the numerical simulations with respect to the analytical theory is of the same order of magnitude of what found also in the previous scan, shown in Sec. 3.1.3, for the value of $q=1.4$. This confirms that the quantitative analytical prediction of the damping rate is very challenging, due to the many approximations needed in deriving explicit formulas.

3.2 GAMs with arbitrary radial structure

In this section, we want to investigate the linear collisionless dynamics of GAMs in a regime where the FOW effects are more important, therefore we push towards higher values of $k_{r}\rho_{i}$, corresponding to GAMs with finer radial structure with respect to the ones considered in Sec. 3.1. We neglect here the effect of the elongation, and we still consider only the results obtained by the analytical theory and numerical simulations with adiabatic electrons.

3.2.1 Analytical predictions

As $k_{r}\rho_{i}q^{2}$ further increases, higher and higher order transit resonances must be taken into account to properly get more accurate GAM damping rates, as it was first shown in Ref. [16], and discussed in details by Z. Qiu in 2009 [17]. The collisionless damping of GAMs for $k_{r}\rho_{i}q^{2}\rightarrow\infty$ was derived in Ref. [5] with all the transit resonances FOW and FLR properly accounted for, and the dispersion was later extended to relatively smaller $q$ parameter region in Ref. [17] to compare with numerical simulations [16]. The dispersion relation of Qiu-2009 was calculated in the limit of large values of $q$ and moderate values of $k_{r}\rho_{i}$, i.e. $1/q^{2}\ll k_{r}\rho_{i}\ll 1$. It reads:

[TABLE]

where $\hat{\omega}=\omega/q\omega_{ti}=x/q$, $\hat{k}=k_{r}\rho_{i}$, $f_{Q1}=31/16+9\tau_{e}/4+\tau_{e}^{2}$, $f_{Q2}=747/32+481\tau_{e}/32+35\tau_{e}^{2}/8+\tau_{e}^{3}/2$, and $f_{Q3}=\tau_{e}^{2}+5\tau_{e}/4+1$. Note that, differently from the analytical predictions described in Sec. 3.1, the GAM frequency has a dependence on $k_{r}\rho_{i}$. The $k_{r}\rho_{i}$ in Eq. 26 comes from both FLR ($J^{2}_{0}(k_{r}\rho_{i})$) and also FOW ($J_{p}(k_{r}\rho_{d})$, with p being integers, and $k_{r}\rho_{d}\simeq k_{r}\rho_{i}q$ for circulating particles. In the limit of large values of $q$ for a fixed $\hat{k}$, then $\hat{\omega}$ tends to a constant with value $\hat{\omega}_{q\rightarrow\infty}\simeq f_{T}^{1/2}$. In this limit, the last term in the first squared bracket of the formula for the damping rate, Eq. 26, and the last two terms in the second squared bracket of the same formula, can be neglected, and the GAM damping rate tends to a constant value. Note that Eqs. 25, 26 can be used for both short/long wavelength regimes. E.g., in long wavelength limit, with $k_{r}\rho_{i}\ll 1$, FOW effects are still important if $q$ is large. So, in general, when we say $q$ is large or small, it is not compared to 1, but to $\sqrt{1/k_{r}\rho_{i}}$.

Finally, an analytic dispersion relation where the effects of the non-circular geometry are also included has been derived by Z. Gao in 2010 [19]. We report here the formulas for $\omega$ and $\gamma$ where no radial derivative of the elongation is considered, for concentric flux surfaces (i.e. with no Shafranov shift), and neglecting the effects of finite inverse aspect ratio. It reads:

[TABLE]

where

[TABLE]

and $f_{G1}=(13-2\tau_{e}-4\tau_{e}^{2})/16$, $f_{G2}=(39+50\tau_{e}+20\tau_{e}^{2})/16$, $f_{G3}=(9+4\tau_{e})/16$ and $f_{G4}=(6-8\tau_{e})/16$. In the case of no elongation, $e=1$, the dispersion relation of Gao-2010 reads:

[TABLE]

Note that the frequency reduces exactly to the one of Qiu-2009. The damping rate can be derived as an approximation of the one of Qiu-2009, in the large-$q$ regime, and when considering only the largest terms in the squared parenthesis of Eq. 26, i.e. respectively $\hat{\omega}^{4}$ and $\hat{\omega}^{2}$. Some differences are therefore expected for the damping rates of Qiu-2009 and Gao-2010.

3.2.2 Equilibrium and definition of the simulation

We choose a tokamak equilibrium with circular flux surfaces and high aspect ratio ($\varepsilon=a/R=0.1$), with $R_{0}=1.3m$, $a=0.13m$. Each simulation has a different q profile, flat, and each one with different value of q. Flat temperature and density profiles are considered. Different values of $\rho^{*}$ are considered (and $T_{i}=T_{e}$). The value of the density is irrelevant for electrostatic simulations. The initialization is done in a similar way as described in Sec. 3.1.2, but we initialize here different simulations with different value of $k_{r}$.

3.2.3 Dependence on the safety factor

A scan with $q$ has been repeated here with ORB5, GYSELA and GENE, similarly to the one reported in Sec. 3.1.3. The value of $\rho^{*}$ here has been chosen as in Sec. 3.1, i.e. $\rho^{*}=1/160$.

Frequency and damping rates of GAMs depend on the safety factor q. To the lowest order in $k_{r}\rho_{i}$, the frequency is well described by the limit of $k_{r}\rho_{i}\rightarrow 0$, and the FOW effects provide corrections which do not modify the order of magnitude of the frequency. All codes seem to follow the analytical prediction obtained without FOW effects at low values of q, whereas there is a change in trend occurring around q=2, where all codes start following the analytical predictions where the FOW effects are included in the frequency.

For the damping rates, the trend of the dependence on q is well described by the limit of small $k_{r}\rho_{i}$, where first order corrections (i.e. accounting for the 2nd harmonic resonance $v_{\parallel}=qR\omega_{GAM}/2$ of the passing ions), only at small q ($q<3$). At larger values of q ($q>3-4$), higher order corrections (i.e. accounting for the 4nd harmonic resonance and higher) are necessary for estimating analytically the GAM damping rate. Note that, in the limit of large values of $q$, the damping rate tends to a constant, as predicted by the analytical theory, Eq. 26.

As a result of this verification test, a good agreement in the scalings measured with ORB5, GYSELA and GENE (both local and global) and with the theoretical prediction of the analytical theory is found, both for the frequency and the damping rate.

3.2.4 Dependence on the radial wave number

The dependence of the frequency and the damping rate on the radial wave number is discussed here. As shown in the previous section, the frequency is well described by the limit of zero FOW to the lowest order, and the corrections of the FOW effects to the value of the frequency provide some modifications, up to 10%. The damping rate dependence on $k_{r}\rho_{i}$, on the other hand, must be considered to orders higher than the first, when $k_{r}\rho_{i}>0.1$, if realistic values of $q$ are considered as measured in tokamaks ($q>4$). Good agreement of ORB5, GYSELA and GENE (both local and global) and the analytical theory is observed for the frequency at low values of $k_{r}\rho_{i}$, while at higher values of $k_{r}\rho_{i}$, the numerical codes predict slightly lower frequencies with respect to the analytical theory. The origin of this discrepancy is thought to be the breaking of the regime of validity of the analytical predictions, derived with the hypothesis of moderate values of $k_{r}\rho_{i}$.

A very good agreement of all codes is observed for the damping rate, except at low values of $k_{r}\rho_{i}$. This is the regime where the damping rates are very small and therefore very difficult to measure, in some cases hidden below the noise (especially for PIC codes). In particular, with PIC codes the cases at very low damping rate require a very high resolution (i.e. a large number of markers) in order for the signal to overcome the statistical error. Therefore, for very low values of the damping rate, the measured numerical value is less trustable, and the error bar becomes bigger. The comparison of the gyrokinetic simulations with the three different analytical formulas of Sugama-2008, Qiu-2009 and Gao-2010 shows that the damping rate is better approximated by Sugama-2008 at very low values of $k_{r}\rho_{i}$ (although this formula still underestimates the damping rate at this large values of $q$), by Qiu-2009 at intermediate values, and by Gao-2010 at large values (see Fig. 4).

4 Numerical simulations with kinetic electrons

4.1 Effect of the finite electron mass, for radially broad modes

In this section, the same equilibrium as in Sec. 3.1.2 is adopted. The flux surfaces are circular, and the safety factor profile is flat, with $q=3.5$. We initialize a scalar potential perturbation with only zonal component, and with a sine dependence on the radius, of the form $\phi(\rho,t=0)=\sin(k_{r}a\rho)$, with $k_{r}=2\pi/a$ (corresponding to a relatively low value of $k_{r}\rho_{i}$). The perturbation is let evolve in a linear electrostatic simulation with kinetic electrons. Our simulations have a spatial grid of (s,$\theta$,$\phi$) = 64x64x4 and a time step of 2 $\Omega_{i}^{-1}$, with $10^{8}$ markers. The length of the simulations is $4\cdot 10^{4}\,\Omega_{i}^{-1}$, corresponding to 20000 time steps.

The dependence of the frequency and damping rate on the ion/electron mass ratio is depicted in Fig. 5, for simulations performed with ORB5 and the global version of GENE. We can see that for the frequency, a convergence towards the values of the adiabatic electrons is observed very soon for increasing $m_{i}/m_{e}$, whereas for the damping rate, the convergence is not observed. For realistic values of $m_{i}/m_{e}$ in deuterium plasmas, the measured damping rate is more than 10 times larger than the value given by the adiabatic electrons, for the chosen value of the safety factor ($q=$ 3.5). A good agreement of the two codes is found for both frequencies (giving results within 2% of difference) and damping rates (within 25% of difference at large mass ratios). Such a difference in the damping rate measured in simulations with kinetic electrons and with adiabatic electrons is due to the effect of the resonance with the bounce motion of barely trapped electrons [28].

5 Summary and conclusions

Zonal (i.e. axisymmetric) poloidal flows, corresponding to zonal radial electric fields, are known to develop in tokamak plasmas, as the result of nonlinear interaction with turbulence. They appear in the form of zero-frequency zonal flows (ZFZF) [1, 2, 3] and oscillating zonal flows, named geodesic acoustic modes (GAM) [4, 5, 10]. Their different behavior in time results in a different efficiency in the turbulence regulation. Both ZFZFs and GAMs are crucial to be understood (linearly and then nonlinearly) for a theoretical characterization of a turbulent plasma. In this paper, we have focused on the linear collisionless dynamics of GAMs.

The linear collisionless theory of GAMs has been developed in different regimes and several numerical investigations have been performed and compared with the theory in the past. Many gyrokinetic codes have also been developed for the study of the nonlinear interaction of turbulence and zonal structures. Nevertheless, no comprehensive linear verification and benchmark effort has been done, to test multiple gyrokinetic codes comparing them with each other and the different analytical theories derived in different limits.

In this paper, we have selected a list of tests which serve for investigating the behaviour of some of the most known gyrokinetic codes in the magnetic-confinement-fusion turbulence community, especially in comparison with each other or with analytical theory. The choice of the codes has been made in order to give an approximative representation of the big variety of turbulence codes existing in our community. The chosen codes have been ORB5 [29, 30, 31], GENE [32, 33], and GYSELA [34, 35]. These codes are based on the same basic gyrokinetic formalism for the treatment of the ion dynamics, which makes them equivalent in the linear electrostatic collisionless regime, which is the one considered here. Additional features can be optionally switched on in some codes, like for example a non-circular geometry of the magnetic flux surfaces, or non-adiabatic models for the electrons. The main basic difference of the three codes, even when circular flux surfaces are considered and the electrons are treated as adiabatic, resides in the numerical algorithm which is used to solve the model equations. In fact, the Lagrangian algorithm is used for ORB5, the Eulerian algorithm for GENE, and the Semi-Lagrangian algorithm for GYSELA. This difference of the numerical schemes, makes the detailed cross-code comparison and verification against analytical theory even more meaningful - the numerical result is controlled not to depend on the numerical approximation of the basic model, but only on the considered physics. The tests have been divided into two main classes, depending on the model used for the treatment of the electrons. In the first class, where the electrons are treated adiabatically, analytical dispersion relations exist in literature, and this makes not only a cross-code benchmark, but also a detailed verification of the codes, possible. On the other hand, when the electrons are treated kinetically, no analytical theory presently exists, and therefore a cross-code benchmark only has been performed.

The first test with adiabatic electrons has been chosen in a regime where all three codes can be compared, namely with a magnetic equilibrium with circular flux surfaces. The frequency and damping rates of GAMs have been observed to fit well among codes, in the limit of moderate-low values of the safety factor, and for small values of the wave-number normalized to the ion Larmor radius (see Sec. 3.1.3). In the same regime, a comparison with the analytical dispersion relation of Zonca-1996 [12], where no FOW effects are retained, and with the explicit formulas for the frequency and damping rate respectively of Sugama-2006 [14] and Sugama-2008 [15], where FOW effects are retained to the first-order, has also been successfully done. When introducing a non-circular geometry of the flux surfaces, the codes ORB5 and GENE have also been been benchmarked and verified against the analytical dispersion relation of Gao-2009 [18], for a scan on the flux surface elongation. The result has been a good agreement of the codes for both frequency and damping rates, a quantitative agreement with the analytical theory for the frequency, and qualitative for the damping rates (see Sec. 3.1.4).

When a regime with higher radial wave-numbers is considered, the ion FOW effects play a more important role. A comparison of ORB5, GENE and GYSELA with adiabatic electrons, with the analytical theories of Sugama-2008, Qiu-2009 and Gao-2010 has been made, scanning in the range $0<k_{r}\rho_{i}\leq 0.45$ (see Sec. 3.2). All codes have shown a very good comparison of the frequency with each other for all values of $k_{r}\rho_{i}$, and a good comparison with the analytical theory for low values of $k_{r}\rho_{i}$. The difference with the analytical theory which is found for higher values of $k_{r}\rho_{i}$, is thought to be due to the breaking of the regime of validity of the analytical theories, derived as expansions for small values of wave-numbers. The damping rate has given a very good matching of all codes, especially at moderate and large values of $k_{r}\rho_{i}$, where the theories of Qiu-2009 and Gao-2010 have been recovered. At low values of $k_{r}\rho_{i}$, corresponding to low values of the damping rate, a difference among codes has been found, due to the general difficulty to measure low damping rates. For example, for a PIC code like ORB5, a high resolution in number of markers is necessary to kill the statistical noise and properly measure a very low value of damping rate, but typically some uncertainty still remains, unless a very big number of markers is used.

Benchmark tests with kinetic electrons have also been performed, with ORB5 and GENE, in a low-$k_{r}\rho_{i}$ regime, with circular flux surfaces, and moderate value of the safety factor. The results of the two codes have been found to fit very well. No analytical theory presently exists providing the modification of the frequency and damping rate due to the effect of the kinetic electrons, therefore no verification has been possible in this regime. The scan of the frequency and damping rate in the ion to electron mass ratio has shown that a convergence of the frequencies with the analytical prediction obtained with adiabatic electrons is found when electrons are sufficiently light (when approaching realistic values of $m_{i}/m_{e}$ for hydrogen and deuterium plasmas) whereas no convergence is found for the damping rate, which stays one order of magnitude higher than the analytical prediction obtained with adiabatic electrons (in agreement with Ref. [28]).

Detailed convergence tests have been performed with all three codes in order to assess the numerical stability for the considered GAM dynamics. Convergence scans have been done for ORB5 with respect to the number of markers, which characterizes the type of discretization of a PIC code (see Appendix A). Analogously, the numerical description of the simulations performed with GENE and the convergence scans in $v_{\|}$ are reported in Appendix B. Finally the description of the numerical parameters used for simulations with GYSELA, and convergence scans in the spatial and velocity space are presented in Appendix C.

In conclusion, we have made a choice of three gyrokinetic codes and tested them for the physics of linear electrostatic collisionless GAMs in different regimes, by means of verification against analytical theory and cross-code benchmarks. These tests have shed light on the regimes of validity of the different analytical theories derived in the different limits. In particular, we have shown that there is not one approximate analytical formula, which can be applied for all the considered regimes. In fact, each considered formula has been found to match the results of the numerical simulations in a different regime of application, but to fail in other regimes. These regimes have been properly identified here, making their usage more sensible for the future. These tests have also improved the trustability of the codes. In particular, we have shown that the results of the three selected codes match very well for all simulations performed in regimes where the damping rate is moderate or high, whereas some differences have been found for very small values of the damping rates, where the numerical error can strongly affect the measurement. These tests performed on zonal structures like GAMs, and complementary tests performed on the linear dynamics of microturbulence modes (see for example Ref. [41]), serve to prepare a solid basis for a more comprehensive theoretical understading of the turbulent transport in tokamak plasmas, based on the numerical simulations with the set of available gyrokinetic codes, analytical theory, and intermediate reduced models, which is one of the major goals of our community.

Acknowledgements

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 under Grant Agreement No. 633053, for the WP15-ER-01/IPP01 project on “Verification and development of new algorithms for gyrokinetic codes”, and WP15-ER-01/IPP02 project on “Micro-turbulence properties in the core of tokamak plasmas: close comparison between experimental observations and theoretical predictions”. The views and opinions expressed herein do not necessarily reflect those of the European Commission. This work was also partly supported by a grant from the Swiss National Supercomputing Centre (CSCS) under project IDs s617 and s704. Simulations were performed on the International Fusion Energy Research Center (IFERC) CSC Helios supercomputer within the framework of the VERIGYRO project, and on the CINECA Marconi supercomputer within the framework of the OrbZONE project. The authors acknowledge discussions with F. Zonca, P. Lauber, E. Poli, D. Del Sarto, A. Ghizzo, P. Niskala, Ö. Gürcan, P. Morel. Part of this work was done while two of the authors, A. Biancalani and I. Novikau, were visiting LPP-Palaiseau (France), whose team is acknowledged for the hospitality.

Appendix A Numerical convergence tests with ORB5

The numerical stability of the codes is crucial to be investigated, in order to assess the efficiency and the regime of validity of the algorithm. A test consists in measuring the convergence of GAM frequencies $\omega_{GAM}$ and damping rates $\gamma_{GAM}$ with different number of markers. This kind of tests has been done for simulations with adiabatic electrons and repeated for simulations with kinetic electrons (see Sec. 2.1 for a description of these two models for the electrons).

For these scans a plasma configuration with a major radius $R_{0}=1.3$ m and a minor radius $a=0.13$ m was chosen with the toroidal magnetic field on axis be $B=1.9$ T and flat profiles for the safety factor, with $q=3.5$, temperature (defined from the value of $\rho^{*}=1/160=0.00625$) and density profiles (with value irrelevant for the present electrostatic simulations). Electrostatic linear simulations are evaluated with an initial electric potential perturbation with $k_{r}=2\pi/a$. This configuration corresponds to the case depicted in Fig. 1 (point with $q=3.5$), and Fig. 5 (points with $m_{i}/m_{e}=2000$ and $m_{i}/m_{e}=4000$).

For the simulations with adiabatic electrons, the typical spatial grid is $(\rho,\alpha_{1},\alpha_{2})=64\times 64\times 4$ and the time step is 40 $\Omega_{i}^{-1}$ (where the ion cyclotron frequency is evaluated here with the magnetic field on axis, i.e. at $\rho=0$, and with $\alpha_{1}$ and $\alpha_{2}$ being the two periodic coordinates, i.e. the poloidal and toroidal angles). Simulations with $1000$ time-steps are considered, with a total time length of $4\cdot 10^{4}$ $\Omega_{i}^{-1}$, where we observe about 10 GAM oscillations. A scan in the number of ion markers is performed, from $10^{5}$ to $10^{8}$. For kinetic electrons, the typical spatial grid is $(\rho,\alpha_{1},\alpha_{2})=64\times 64\times 4$, and the time-step is $2$ $\Omega_{i}^{-1}$. For the case with kinetic electrons, simulations with $2\cdot 10^{4}$ time-steps are considered, a number of ion markers of $10^{8}$, and the number of electron markers is scanned from $10^{7}$ to $5\cdot 10^{8}$.

The frequency has been calculated directly by measuring the averaged period of oscillation at one radial position. To apply other techniques, like for example the Fourier decomposition, it’s necessary to have more oscillations that increases significantly the calculation time of simulations. The damping rate and its standard deviation have been found by using the method of least squares. Results of the convergence tests are given in Fig. 6, where it can be seen that for the case of adiabatic electrons, the frequency and the damping rate converge well to the analytical value calculated by using the explicit expressions of Sugama-2006 and Sugama-2008 [14, 15], for increasing number of ion markers. On the other hand, the absolute values of the damping rates for simulations with kinetic electrons are found to stay considerably higher (as described in Sec. 4), and no convergence with the results of simulations with adiabatic electrons is observed in the range of number of electron markers considered. The GAM frequency does not change much with the number of markers, except for the cases with very small number of markers. Error bars are also reported in the values of the damping rates, to emphasize that at very low number of markers, the Monte-Carlo error becomes comparable with the physical signal damping.

In the simulations with kinetic electrons performed with ORB5 and discussed in this paper, the dynamics of passing electrons is treated kinetically, and consequently high frequency oscillations are observed on top of the lower frequency GAM oscillation (see also Ref. [27]). These high-frequency oscillations correspond to the limit of kinetic Alfvén waves for $\beta$ going to zero (electrostatic model) at fixed temperature, also known as the $\omega_{H}$-mode [47]. These high-frequency oscillations have been observed to create numerical instabilities for low number of markers (below $10^{7}$). For this reason, the results of simulations with kinetic electrons and electron markers below $10^{7}$ have not been reported in Fig. 6.

Regarding the numerical parameters of the simulations of GAMs with broad radial structure described in Sec. 3.1.2, we have used a spatial grid of ($N_{\rho}$,$N_{\theta}$,$N_{\phi}$) = 256x64x4 and a time step of $dt=$100 $\Omega_{i}^{-1}$, with $N_{i}=10^{8}$ markers. The length of the simulations is $4\cdot 10^{5}\,\Omega_{i}^{-1}$, corresponding to 400 time steps. Regarding the simulations of GAMs with fine radial structure described in Sec. 3.2.2, a typical simulation has a spatial grid of ($N_{\rho}$,$N_{\theta}$,$N_{\phi}$) = 256x64x4 and a time step of 25, 50 and 100 $\Omega_{i}^{-1}$, with $10^{7}$ and $10^{8}$ markers.

Appendix B Numerical convergence tests with GENE

GENE simulations are carried out considering an initial density perturbation with the same sinusoidal functional form as described in Sec. 3.1.2. In order to match the radial wave-number of the initial perturbation, the radial domain $L_{x}$ is adapted for each value of $k_{r}=2\pi/L_{x}$. The mid-radius location, $r/a=0.5$ is the reference position used to measure all normalization quantities and define the dimensionless machine size parameter $\rho^{*}$. The typical resolution used in the radial direction is one point per ion larmor radius, with the number of points adapted such as to have always one grid-point located at $r/a=0.5$. A high spatial resolution is used in the parallel direction, up to 96 points, which turn out to be necessary in order to correctly converge the GAM damping for the large $q$ - small $k_{r}$ cases. In velocity space we consider the domain $L_{v_{\|}}\times L_{\mu}$=$6\times 12$, a choice that will be justified in the following. A typical grid is $n_{v_{\|}}\times n_{\mu}=256\times 32$ points, where the high resolution in the parallel velocity is motivated by the need of avoiding the recurrence problem. A detailed discussion of this issue is outside the scope of this paper, and the interested reader is referred to e.g. [53], where the recurrence problem is discussed in details. With the aforementioned resolution, the recurrence time is longer than the final simulated time in all cases considered here, thus the mode frequency and damping can be easily extracted. Alternatively one could have used a small hyperdiffusion in the $v_{\|}$ direction obtaining the same result. However, in general, we prefer avoiding introducing any numerical dissipation as this might impact the residual level of zonal flows (not considered in this paper).

The properties of the GAM are evaluated by analyzing the time traces of the flux-surface-averaged electrostatic potential $\bar{\phi}$, measured at mid-radius. When comparing flux-tube and global simulations, the same time interval is used. Simulations with adiabatic electrons are run typically up to 150 $R_{0}/\sqrt{2}\,v_{ti}$, in order to collect sufficiently long statistics (for the large damping cases it suffices to run the simulation for much shorter times). The damping rate $\gamma_{\rm GAM}$ of the GAM is then evaluated by separately fitting maxima and minima of the curve $\bar{\phi}$ in time. The frequency is computed using an Hilbert transform.

In Figure 7 we plot GAM frequency and damping for two different values of the safety factor (2 and 4 respectively) varying the ion velocity space domain and resolution. These simulations have been performed with the flux-tube version of GENE with adiabatic electrons. They have been repeated for global simulations (results not shown here) obtaining, as expected, the same behavior and an almost perfect agreement with local results. We observe how the GAM frequency rapidly converges, whereas the damping is much more sensitive to resolution, and a sufficiently large velocity space must be considered in order to correctly converge the simulation results.

We remark that kinetic electron runs are instead carried out for a significantly shorter time, $\sim 30R_{0}/\sqrt{2}\,v_{ti}$, as the damping is found to be much stronger and it is therefore not necessary to simulate longer times.

Appendix C Numerical convergence tests with GYSELA

The convergence scan proposed for GYSELA has been performed with the same parameters as described in Sec. 3.1.2. Density and temperature profiles are considered flat. The flat safety factor is taken equal to $2$ for the following tests. Electrons are considered adiabatic. In GYSELA, due to its full-$f$ character, the initial condition is performed on the distribution function $F_{s}$ and consists of an equilibrium distribution function $F_{s0}$ added to a perturbation $\delta F_{s}$, namely $F_{s}=F_{s0}+\delta F_{s}$. Then, the electrostatic potential $\phi(r,\theta,\varphi)$ is computed at time $t=0$ by solving quasi-neutrality equation (15). In the present test, the perturbation part $\delta F_{s}$ reads $\delta F_{s}=F_{s0}\>g(r)$ with $g(r)=-\frac{1}{r}\left(k_{r}\cos(k_{r}r)-k_{r}^{2}r\sin(k_{r}r)\right)$ where $k_{r}=(k+1)\pi/L_{r}$ with $k\in\mathbb{N}$, and $L_{r}=160\rho_{s}$. The corresponding radial profile of the zonal component $\phi_{00}(r)\sim\sin(k_{r}r)$ is plotted in Fig. 8 (black line) for $k=1$ which is the value used for the following simulations.

In GYSELA, the $5$D space $(r,\theta,\varphi,v_{\parallel},\mu)$ is uniformly discretized with $N_{r}\times N_{\theta}\times N_{\varphi}$ points in the 3D real space and $N_{v_{\parallel}}\times N_{\mu}$ points in the $2$D velocity space. This mesh grid is fixed in time with $r\in[0,L_{r}]$, $\theta\in[0,2\pi]$, $\varphi\in[0,2\pi]$, $v_{\parallel}\in[-\alpha\,v_{T_{i}},\alpha\,v_{T_{i}}]$ and $\mu\in[0,L_{\mu}]$. Due to the toroidal axisymmetry of the test the number of toroidal points $N_{\varphi}$ is fixed to $N_{\varphi}=8$. A comparison (not presented here) with $N_{\varphi}=16$ has shown really good agreement with $N_{\varphi}=8$. Simulations with $N_{\varphi}=4$ would be probably close to those with $N_{\varphi}=8$ but are not possible in the code due to parallelization constraint. This technical constraint could be removed. However simulations with so little number of points in toroidal direction are not standard simulations, so choice has been made to run with $N_{\varphi}=8$ and to postpone the required modification of the code for now. The maximum of thermal velocities in parallel velocity space is fixed at $\alpha=7$. A simulation with $\alpha=5$ has been performed (not presented here) showing very small departure ($<2\%$) compared to the case $\alpha=7$. However, as this value could have more impact for larger $q$ values due to resonance position the value $\alpha=7$ has been preferred for the following tests. $L_{\mu}$ is fixed to $L_{\mu}=12\,T_{i}/B_{0}$ (with $B_{0}=1$). All simulations have been performed for a flat safety factor profile equal to $2$ and until $t=50000\,\Omega_{i}^{-1}$. Flat density and temperature profiles are also considered with $\tau_{e}=T_{e}/T_{1}=1$.

Parameters and results are summarized in Table 1. Comparisons are performed on the three quantities: (i) the radial wave number $k_{r}\rho_{i}$, (ii) the damping rate $\gamma$ and (iii) frequency $\omega$ of the zonal component of the electrostatic potential $\phi_{00}$. The radial wave number is computed with the following formula:

[TABLE]

with $\bar{\phi}_{\rm norm}(r_{i},t)=\bar{\phi}(r_{i},t)/\max_{r_{i}}\bar{\phi}(r_{i},t)$. The values reported in Table 1 correspond to the mean values of $k_{r}\rho_{i}$ computed at times where $\log(\bar{\phi}(r_{p},t))$ is maximum with $r_{p}$ the radial position of the maximum value of $\bar{\phi}(r)$ at initial time. The damping rate is estimated by using the method of least squares also on the maximum values of $\log(\bar{\phi}(r_{p},t))$. $\gamma$ values reported in Table 1 are computed with 6 maximums (see red circles in Figure 9). Four first simulations (cases 1 to 4 in Table 1) have been performed for the same 5D mesh of $\sim 536.8$ millions of points $(N_{r},N_{\theta},N_{\varphi},N_{v_{\parallel}},N_{\mu})=(256,64,8,128,32)$ but varying the time step $\Delta t$ from $\Delta t=5\,\Omega_{i}^{-1}$ to $\Delta t=50\,\Omega_{i}^{-1}$. All the other simulations except the last one (cases 5 to 10) have been performed with $\Delta t\Omega_{i}=25$ varying: (i) the number of points in $\mu$ direction (case 5: $N_{\mu}=16$, case 6:$N_{\mu}=8$); (ii) the number of points in radial direction (case 7: $N_{r}=128$, case 8: $N_{r}=512$); (iii) the number of points in poloidal direction (case 9: $N_{r}\times N_{\theta}=256^{2}$) and (iv) finally the number of points in parallel velocity space (case 10: $N_{v_{\parallel}}=64$). The last case (case 11) corresponds to a simulation where all varying parameters have been taken to their smaller tested value, namely $\Delta t=50\Omega_{i}^{-1}$, $N_{r}=128$, $N_{\theta}=64$, $N_{v_{\parallel}}=64$ and $N_{\mu}=8$.

Considering case 1 as the reference case, the maximum relative error is less than $1\%$ for $k_{r}\rho_{i}$ and $\omega$ estimations and less than $2\%$ for $\gamma$ (see Table 1). As conclusion all these simulations, even the coarse grained one (case 11), are fully accurate. However considering that these tests have been performed for a small safety factor value $q=2$ and a small radial wave number $k_{r}\rho_{i}\sim 0.056$ we could suggest to avoid parameters where we observe a small departure from the reference case, namely $\Delta t=50\Omega_{i}^{-1}$ and $N_{\mu}=8$. Then, more secure parameters for larger $q$ values or larger $k_{r}\rho_{i}$ values could correspond to those of case 5, namely a mesh $(N_{r},N_{\theta},N_{\varphi},N_{v_{\parallel}},N_{\mu})=(256,64,8,128,16)$ of $268.4$ millions of points with a time step of $\Delta t=25\,\Omega_{i}^{-1}$. Such a simulation requires $2$ hours on $256$ cores for $2000$ time iterations compared to the coarse grained simulation which takes around $1$ hour on $64$ cores ($1000$ iterations).