# Design of the IMAS interface data structure for flux-tube gyrokinetic simulations

Yann Camenen, Frédéric Imbeaux, Guillaume Fuhr, Tobias Görler, Anass Najlaoui

► To cite this version:

Yann Camenen, Frédéric Imbeaux, Guillaume Fuhr, Tobias Görler, Anass Najlaoui. Design of the IMAS interface data structure for flux-tube gyrokinetic simulations. 2026. hal-05490777

HAL Id: hal-05490777

https://hal.science/hal-05490777v1

Preprint submitted on 3 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Design of the IMAS interface data structure for
flux-tube gyrokinetic simulations

Y. Camenen¹, F. Imbeaux², G. Fuhr¹, T. Görler³, and A. Najlaoui¹

¹Aix-Marseille Univ., CNRS, PIIM UMR7345, Marseille, France

²CEA, IRFM, 13108 Saint-Paul-lez-Durance, France

³Max Planck Institute for Plasma Physics, Boltzmannstr. 2, Garching,
Germany

February 3, 2026

Abstract

A unified set of conventions and normalizations for the inputs and
outputs of delta-f flux-tube gyrokinetic simulations has been developed
within the ITER Integrated Modelling and Analysis Suite (IMAS). This
new standard aims at facilitating benchmarks and data exchange between
different code users, increasing the modularity of frameworks manipulat-
ing inputs and outputs from gyrokinetic simulations and serving as the
basis for databases storing simulations from multiple gyrokinetic codes.

This work is distributed under a Creative Commons
Attribution 4.0 International licence.

https://creativecommons.org/licenses/by/4.0/

<figure><img src="image_2.png" /></figure>
# Contents

<table>
  <tbody>
    <tr>
      <td>1 Introduction</td>
      <td>3</td>
    </tr>
    <tr>
      <td>2 The local delta-f gyrokinetic system of equations</td>
      <td>3</td>
    </tr>
    <tr>
      <td>2.1 Magnetic equilibrium</td>
      <td>3</td>
    </tr>
    <tr>
      <td>2.2 Toroidal field and plasma current direction</td>
      <td>4</td>
    </tr>
    <tr>
      <td>2.3 Flux coordinate system (r, θ, φ)</td>
      <td>4</td>
    </tr>
    <tr>
      <td>2.4 Flux-tube coordinates and Fourier representation</td>
      <td>5</td>
    </tr>
    <tr>
      <td>2.5 Chosen form of the gyrokinetic equation</td>
      <td>6</td>
    </tr>
    <tr>
      <td>2.5.1 Frame of reference</td>
      <td>6</td>
    </tr>
    <tr>
      <td>2.5.2 Velocity coordinates</td>
      <td>6</td>
    </tr>
    <tr>
      <td>2.5.3 Centrifugal effects</td>
      <td>6</td>
    </tr>
    <tr>
      <td>2.5.4 Background distribution function</td>
      <td>7</td>
    </tr>
    <tr>
      <td>2.6 The gyrokinetic equation</td>
      <td>7</td>
    </tr>
    <tr>
      <td>2.7 The field equations</td>
      <td>8</td>
    </tr>
    <tr>
      <td>2.8 Moments of the particle, guiding-center and gyrocenter distribution functions</td>
      <td>9</td>
    </tr>
    <tr>
      <td>2.9 Fluxes in the rotating frame</td>
      <td>10</td>
    </tr>
    <tr>
      <td>2.9.1 Particle flux</td>
      <td>10</td>
    </tr>
    <tr>
      <td>2.9.2 Energy flux</td>
      <td>10</td>
    </tr>
    <tr>
      <td>2.9.3 Momentum flux</td>
      <td>11</td>
    </tr>
    <tr>
      <td>2.10 Fluxes in the Laboratory frame</td>
      <td>11</td>
    </tr>
    <tr>
      <td>2.11 Transport equations</td>
      <td>12</td>
    </tr>
    <tr>
      <td>3 The IMAS gyrokinetics_local IDS</td>
      <td>12</td>
    </tr>
    <tr>
      <td>3.1 Architecture</td>
      <td>12</td>
    </tr>
    <tr>
      <td>3.2 Why a data model with normalized quantities?</td>
      <td>14</td>
    </tr>
    <tr>
      <td>3.3 normalizations</td>
      <td>14</td>
    </tr>
    <tr>
      <td>3.4 Meta-data</td>
      <td>15</td>
    </tr>
    <tr>
      <td>3.4.1 IDS properties</td>
      <td>15</td>
    </tr>
    <tr>
      <td>3.4.2 Code</td>
      <td>15</td>
    </tr>
    <tr>
      <td>3.4.3 Normalizing quantities</td>
      <td>16</td>
    </tr>
    <tr>
      <td>3.5 Model</td>
      <td>16</td>
    </tr>
    <tr>
      <td>3.6 Species parameters</td>
      <td>17</td>
    </tr>
    <tr>
      <td>3.6.1 Species all</td>
      <td>17</td>
    </tr>
    <tr>
      <td>3.6.2 Species</td>
      <td>18</td>
    </tr>
    <tr>
      <td>3.6.3 Collisions</td>
      <td>19</td>
    </tr>
    <tr>
      <td>3.7 Magnetic equilibrium</td>
      <td>19</td>
    </tr>
    <tr>
      <td>3.8 Linear runs</td>
      <td>21</td>
    </tr>
    <tr>
      <td>3.8.1 Mode amplitude</td>
      <td>22</td>
    </tr>
    <tr>
      <td>3.8.2 Eigenvalues</td>
      <td>22</td>
    </tr>
    <tr>
      <td>3.8.3 Eigenfunctions</td>
      <td>23</td>
    </tr>
    <tr>
      <td>3.8.4 Moments</td>
      <td>24</td>
    </tr>
    <tr>
      <td>3.8.5 Linear weights</td>
      <td>25</td>
    </tr>
    <tr>
      <td>3.9 Non-linear runs</td>
      <td>26</td>
    </tr>
    <tr>
      <td>3.9.1 Fields</td>
      <td>26</td>
    </tr>
    <tr>
      <td>3.9.2 Moments</td>
      <td>27</td>
    </tr>
    <tr>
      <td>3.9.3 Fluxes</td>
      <td>27</td>
    </tr>
    <tr>
      <td>4 Application examples</td>
      <td>28</td>
    </tr>
    <tr>
      <td>5 Summary and perspectives</td>
      <td>30</td>
    </tr>
    <tr>
      <td>A Catalog of Bessel function integrals</td>
      <td>31</td>
    </tr>
  </tbody>
</table>
# 1 Introduction

The gyrokinetic theory [1] addresses the description of turbulence and of the associated transport in strongly magnetised plasmas. It exploits the scale separation between the (fast) gyro-motion and (slow) turbulence evolution, to perform a reduction of the Vlasov-Maxwell system of equations and eliminate the gyro-motion from the dynamics. Gyrokinetic theory is to date the most advanced theory to model turbulent transport in tokamak plasmas. While solving the complete gyrokinetic system of equations is possible on present largest supercomputers, it still requires extreme computational resources. Reduced versions, with various trade-off between accuracy and tractability, have therefore been developed. The current reference for the construction of transport models in tokamak simulation frameworks is the delta-f flux-tube limit of the gyrokinetic system of equations [2]. In this limit, spatial and temporal scale separation is assumed between fluctuating quantities, varying on the scale of the Larmor radius and evolving on the scale of the transit/drift frequencies, and background quantities, varying on the scale of the plasma minor radius and assumed static. The gyrokinetic system of equations is solved in the vicinity of a selected flux surface in a radially and poloidally bounded domain that follows a magnetic field line (a flux-tube). Equilibrium quantities (e.g. magnetic field, temperature, density) are assumed constant across the flux-tube and turbulence is considered homogeneous in the plane perpendicular to the equilibrium magnetic field. This is the local approximation. The distribution function of gyrocenters of each species is split into an equilibrium part, usually a Maxwellian, and a fluctuating part. This is the delta-f approximation. Over the years, a number of gyrokinetic codes have been developed to solve the delta-f local gyrokinetic system of equations, e.g. CGYRO [3], GENE [4], GKV [5], GKW [6], GS2 [7] or GYRO [8]. These multiple numerical implementations are essential to check that the equations are solved correctly and code benchmarks are regularly conducted to that purpose, e.g. [9], [10], [11]. Different refinements of the physics description are also explored in parallel and, once validated, implemented in other codes. One recurrent difficulty encountered when performing benchmarks or exchanging input files between different gyrokinetic codes is that each code has its own definitions and normalizations of the physical quantities. The present work defines a common set of conventions and data normalizations for the inputs and outputs of delta-f flux-tube simulations. The goal of this common format is to ease the data exchange between groups using different gyrokinetic codes, facilitate benchmarks and allow for their automation. It is also intended to serve as the basis for the deployment of databases populated with simulations from multiple gyrokinetic codes. This new standard has been created within the ITER Integrated Modelling and Analysis Suite (IMAS) [12] which includes a standardized data model used to couple physics components and save data in a common format. In the following, we introduce the set of equations solved in delta-f flux-tube gyrokinetic codes in Sec. 2 and then detail the conventions and normalizations adopted in IMAS to describe the simulation results in Sec. 3. As an illustration, in Sec. 4, a simple benchmark is performed between GENE, GKW and the gyro-fluid code TGLF [13] using the newly developed standard.

# 2 The local delta-f gyrokinetic system of equations

## 2.1 Magnetic equilibrium

The equilibrium magnetic field $B$ is assumed to be axisymmetric and described by:

$$B = F \nabla \varphi - \frac{1}{2\pi} \nabla \varphi \times \nabla \Psi \quad (1)$$

where $F = RB_t$ is constant on a flux surface and $\Psi$ is the poloidal magnetic flux (not divided by $2\pi$). Following the ITER coordinate convention, which for IMAS Data Dictionary versions 4.0.0 and above corresponds to COCOS=17 in the COordinate COnventionS system [14], the cylindrical coordinate
<figure><img src="image_3.png" /><figcaption>Figure 1: Coordinate conventions used in the IMAS gyrokinetics_local Data Dictionary.</figcaption></figure>

system $(R, \varphi, Z)$ is chosen right handed, with $R$ the major radius, $Z$ the elevation and $\varphi$ the toroidal angle (increasing when anticlockwise from above), see Fig. 1. The poloidal flux $\Psi$ is maximum at the magnetic axis for a plasma current flowing in the direction of $\nabla\varphi$ and minimum at the magnetic axis for a plasma current flowing in the direction opposite to $\nabla\varphi$.

## 2.2 Toroidal field and plasma current direction

The toroidal field and plasma current are counted positive in the direction of $\nabla\varphi$. Their respective signs, $s_b$ and $s_j$ are defined as

$$s_b = \operatorname{sign}[\mathbf{B} \cdot \nabla\varphi] \qquad (2)$$

$$s_j = \operatorname{sign}[\mathbf{j} \cdot \nabla\varphi] \qquad (3)$$

where $\mathbf{j}$ is the plasma current density. With these conventions, a positive plasma current implies a positive poloidal magnetic field.

## 2.3 Flux coordinate system $(r, \theta, \varphi)$

The gyrokinetic system of equations is solved locally, in a flux-tube domain, wrapping around the flux surface labeled by $\Psi_0$. The poloidal contour of this flux surface is defined by $\Psi(R, Z) = \Psi_0$ and noted $\{R_{\Psi_0}, Z_{\Psi_0}\}$. A radial coordinate $r$ is first introduced:

$$r = \frac{1}{2} [\max[R_{\Psi}] - \min[R_{\Psi}]] \qquad (4)$$

It has the dimension of a length and is constant on a flux surface, by definition. The center of each flux surface is defined as:

$$R_c = \frac{1}{2} [\max[R_{\Psi}] + \min[R_{\Psi}]] \qquad (5)$$

$$Z_c = \frac{1}{2} [\max[Z_{\Psi}] + \min[Z_{\Psi}]] \qquad (6)$$

The minor radius and the center of the flux surface labeled by $\Psi_0$ are noted:

$$r_0 = r(\Psi_0) \qquad (7)$$

$$(R_0, Z_0) = (R_c(\Psi_0), Z_c(\Psi_0)) \qquad (8)$$
Note that $(R_0, Z_0)$ is not the position of the magnetic axis, although it can be in some particular cases. The poloidal angle is defined with respect to the reference point $(R_0, Z_0)$:

$$ \cos \theta = \frac{R_{\Psi_0} - R_0}{[(R_{\Psi_0} - R_0)^2 + (Z_{\Psi_0} - Z_0)^2]^{1/2}} \quad (9) $$

$$ \sin \theta = \frac{Z_{\Psi_0} - Z_0}{[(R_{\Psi_0} - R_0)^2 + (Z_{\Psi_0} - Z_0)^2]^{1/2}} \quad (10) $$

This definition is local to the flux surface $\Psi = \Psi_0$ since it depends on the flux surface center $(R_0, Z_0)$. The poloidal angle is zero at the low field side midplane, i.e. when $Z = Z_0$ and $R_{\Psi_0} > R_0$, and increases clockwise when the tokamak vertical axis is on the left of the flux surface, see Fig. 1.

## 2.4 Flux-tube coordinates and Fourier representation

We now introduce a flux-tube coordinate system $(x, y, \theta)$, where:

* $x = \frac{r-r_0}{R_0}$ is the radial coordinate. It is a flux surface label and increases from the magnetic axis to the plasma edge

* $y$ is the binormal coordinate. It is left unspecified but assumed to increase in the direction of $\varphi$, i.e. $\nabla y \cdot \nabla \varphi > 0$

* $\theta$ is the parallel coordinate (i.e. the extended ballooning angle)

By definition, $\mathbf{B} \cdot \nabla x = 0$ and $\mathbf{B} \cdot \nabla y = 0$ where $\mathbf{B}$ is the equilibrium magnetic field. Note that the sign of $\nabla y \cdot \nabla \theta$ depends on the sign of the poloidal and toroidal magnetic field. For a toroidal field and plasma current flowing in the same direction, we have $\nabla y \cdot \nabla \theta < 0$. When they flow in opposite direction, we have $\nabla y \cdot \nabla \theta > 0$.

Since for any coordinate system $(x^i, x^j, x^k)$, a vector $\mathbf{A}$ can be written as:

$$ \mathbf{A} = \frac{1}{2} \epsilon_{ijk} A^i \mathcal{J}_{x^i x^j x^k} \nabla x^j \times \nabla x^k \quad (11) $$

the magnetic field can be expressed as

$$ \mathbf{B} = B^{\theta} \mathcal{J}_{xy\theta} \nabla x \times \nabla y \quad (12) $$

where the Jacobian is given by

$$ \mathcal{J}_{xy\theta} = \frac{1}{\nabla x \cdot \nabla y \times \nabla \theta} \quad (13) $$

Since, in the local approximation, the turbulence is homogeneous in the plane perpendicular to the magnetic field, the perturbed distribution function and the fields can be expressed as a sum of Fourier modes. For any function $h$, we note $\hat{h}$ its Fourier components:

$$ h(x, y, \theta) = \sum_{k_x, k_y} \hat{h}(k_x, k_y, \theta) \exp[ik_x x + ik_y y] \quad (14) $$

Note that for another coordinate system $(x_1, y_1)$ where the binormal coordinate is of opposite sign, for instance $x_1 = x$ and $y_1 = -y$, then:

$$ \hat{h}(k_x, k_y, \theta) = \hat{h}(k_{x_1}, -k_{y_1}, \theta) = \hat{h}(-k_{x_1}, k_{y_1}, \theta)^{\dagger} \quad (15) $$

where $\dagger$ represents the complex conjugate and hermitian symmetry has been used for the last equality. This, together with the convention $\nabla y \cdot \nabla \varphi > 0$ for the binormal coordinate needs to be kept in mind when coupling a gyrokinetic code to IMAS.

The binormal wavevector is related to the toroidal mode number $n$ by

$$ k_y = \frac{2\pi n}{\Delta y} \quad (16) $$
where $\Delta y$ is the increment in the $y$ coordinate over one toroidal turn. In a flux-tube, the perpendicular wavevector $k_{\perp}$ is a function of the parallel coordinate:

$$k_{\perp}^2(\theta) = g^{xx}(\theta)k_x^2 + 2g^{xy}(\theta)k_x k_y + g^{yy}(\theta)k_y^2 \quad (17)$$

where the metric tensors are defined as $g^{xx} = \nabla x \cdot \nabla x$, $g^{xy} = \nabla x \cdot \nabla y$ and $g^{yy} = \nabla y \cdot \nabla y$.

## 2.5 Chosen form of the gyrokinetic equation

In the literature, the gyrokinetic equation is expressed in various forms, depending, for instance, on the choice of the velocity space coordinates, of the frame of reference, or whether the equation is solved for the gyro-center or guiding-center distribution function. Some forms are more suited than others to specific studies or numerical implementations, but the differences in formulation do not affect the physics. We choose here to formulate the gyrokinetic system of equations that describes the evolution of the gyrocenter distribution function in the toroidally rotating frame and use the parallel velocity and magnetic moment as velocity coordinates. The gyrokinetic equation is introduced to help clarifying the IMAS definitions of Sec. 3 but these definitions are general and independent of the chosen form of the gyrokinetic equation.

### 2.5.1 Frame of reference

The frame rotation, $\Omega$, is chosen to be the zeroth-order flow at the flux surface of interest, $r = r_0$, so that the background $E \times B$ drift cancels in the gyrokinetic equation. The zeroth-order plasma flow is purely toroidal and species independent [15]:

$$v_0 = R^2 \omega_{\Phi} \nabla \varphi \quad (18)$$

with the toroidal angular frequency $\omega_{\Phi}$ proportional to the background radial electric field:

$$\omega_{\Phi} = 2\pi \frac{\partial \Phi_0}{\partial \Psi} \quad (19)$$

and the background potential $\Phi_0$ constant on a flux surface. The frame toroidal angular frequency is therefore given by $\Omega = \omega_{\Phi}(r_0)$.

### 2.5.2 Velocity coordinates

The parallel velocity $v_{\parallel}$ and the magnetic moment $\mu_s = \frac{m_s v_{\perp}^2}{2B}$ are used as velocity coordinates in the rotating frame. Here, $m_s$ is the species mass, $B$ the norm of the magnetic field and $v_{\perp}$ the perpendicular velocity coordinate. The volume element for integrals over velocity space is given by:

$$d\mathbf{v} = 2\pi v_{\perp} dv_{\perp} dv_{\parallel} = \frac{2\pi B}{m_s} d\mu_s dv_{\parallel} \quad (20)$$

### 2.5.3 Centrifugal effects

In terms of physics, including centrifugal effects induced by the background plasma flow adds several terms to the gyrokinetic equation. These terms, related to centrifugal drifts and electrostatic trapping have important implications on stability and transport when the background toroidal rotation is comparable to the thermal velocity. Centrifugal effects have been introduced into the gyrokinetic code GKW [16, 17] based on a formulation of the equations in the frame moving toroidally with the background flow [18]. They have also been implemented in the gyrokinetic code CGYRO [19] based on a formulation of the equations in the Laboratory frame [20]. Both implementations have been shown to be equivalent [19]. All terms related to centrifugal effects are included in the present description.
### 2.5.4 Background distribution function

The background distribution function of gyrocenters $F_{Ms}$ is assumed to be a Maxwellian:

$$F_{Ms} = \frac{n_s}{\pi^{3/2} v_{ths}^3} \exp \left[ -\frac{m_s v_{\parallel}^2}{2T_s} - \frac{\mu_s B}{T_s} \right] \quad (21)$$

where $n_s$, $T_s$ and $v_{ths} = \sqrt{2T_s/m_s}$ are the species density, temperature and thermal velocity, respectively. In the presence of sonic toroidal flows, the plasma density is no longer constant on a flux surface and is given by:

$$n_s(\theta) = n_{s0} \exp \left[ -\frac{\mathcal{E}_s(\theta)}{T_s} \right] \quad (22)$$

where $n_{s0} = n_s(\theta = 0)$. The term $\mathcal{E}_s$ represents the centrifugal energy and the potential energy associated to the poloidally varying component of the electric potential, see e.g. [16, 19]. It is quadratic in the toroidal flow and vanishes at $\theta = 0$ by choice of $n_{s0}$.

## 2.6 The gyrokinetic equation

In the toroidally rotating frame with angular frequency $\Omega$, the temporal evolution of the perturbed distribution function of gyrocenters, $f_s$, is given by:

$$\begin{aligned} \frac{\partial f_s}{\partial t} = & -[v_{\parallel} \mathbf{b} + v_{\Phi_0} + v_{\chi} + v_D] \cdot \nabla f_s + \frac{\mathbf{b}}{m_s} \cdot [\nabla \mathcal{E}_s - \mu_s \nabla B] \frac{\partial f_s}{\partial v_{\parallel}} \\ & -[v_{\chi} + v_D] \cdot \nabla_S F_{Ms} - \frac{Z_s e v_{\parallel}}{T_s} \frac{\partial \bar{A}_{\parallel}}{\partial t} F_{Ms} \\ & -[v_{\parallel} \mathbf{b} + v_D + v_{\bar{A}_{\parallel}}] \cdot [Z_s e \nabla \bar{\phi} + \mu_s \nabla \bar{B}_{\parallel}] \frac{F_{Ms}}{T_s} \\ & +C(f_s) \end{aligned} \quad (23)$$

In this equation, $\mathbf{b} = \mathbf{B}/|\mathbf{B}|$ is the unit vector in the direction of the background magnetic field and $e$ is the elementary charge. The background electric potential $\Phi_0$, the perturbed electric potential $\phi$, the parallel component of the perturbed vector potential $A_{\parallel}$ (magnetic flutter) and the parallel component of the perturbed magnetic field $B_{\parallel}$ (magnetic compression) lead to a perpendicular motion given by

$$v_{\Phi_0} + v_{\chi} = \frac{\mathbf{b} \times \nabla \Phi_0}{B} + \frac{\mathbf{b} \times \nabla \chi}{B} \quad (24)$$

where

$$\chi = \bar{\phi} - v_{\parallel} \bar{A}_{\parallel} + \frac{\mu_s}{Z_s e} \bar{B}_{\parallel} \quad (25)$$

$Z_s$ is the species charge number and the gyro-averaged fields are given by:

$$\bar{\phi} = \sum_{k_x, k_y} J_0(k_{\perp} \rho_s) \hat{\phi} \exp[ik_x x + ik_y y] \quad (26)$$

$$\bar{A}_{\parallel s} = \sum_{k_x, k_y} J_0(k_{\perp} \rho_s) \hat{A}_{\parallel} \exp[ik_x x + ik_y y] \quad (27)$$

$$\bar{B}_{\parallel s} = \sum_{k_x, k_y} \hat{J}_1(k_{\perp} \rho_s) \hat{B}_{\parallel} \exp[ik_x x + ik_y y] \quad (28)$$

where $k_{\perp}$ is the perpendicular wavevector defined in Eq. (17), $\rho_s = m_s v_{\perp} / (Z_s e B)$ the species Larmor radius and $J_n$ the n-th order Bessel function of the first kind with $\hat{J}_1$ defined as:

$$\hat{J}_1(k_{\perp} \rho_s) = \frac{2 J_1(k_{\perp} \rho_s)}{k_{\perp} \rho_s} \quad (29)$$
The component of $v_\chi$ due to magnetic flutter is noted $v_{\bar{A}_\parallel}$. The curvature drift $v_{cu}$, the $\nabla B$ drift $v_{\nabla B}$, the Coriolis drift $v_{co}$ and the centrifugal drift $v_{cf}$ are grouped into:

$$v_D = v_{cu} + v_{\nabla B} + v_{co} + v_{cf} \qquad (30)$$

where

$$v_{cu} = \frac{m_s v_{\parallel}^2}{Z_s e B} \mathbf{b} \times (\mathbf{b} \cdot \nabla)\mathbf{b} \qquad (31)$$

$$v_{\nabla B} = \frac{\mu_s \mathbf{b} \times \nabla B}{Z_s e B} \qquad (32)$$

$$v_{co} = \frac{2m_s v_{\parallel}}{Z_s e B} \Omega_{\perp} \qquad (33)$$

$$v_{cf} = \frac{1}{Z_s e B} \mathbf{b} \times \nabla \mathcal{E}_s \qquad (34)$$

with $\Omega_\perp = \Omega - (\Omega \cdot \mathbf{b})\mathbf{b}$. In most, if not all, gyrokinetic codes, the radial force balance is invoked to split the curvature drift into a term proportional to $\mathbf{b} \times \nabla B$ and another term proportional to the radial pressure gradient:

$$v_{cu} = \frac{m_s v_{\parallel}^2}{Z_s e B} \left[ \frac{\mathbf{b} \times \nabla B}{B} + \frac{\mu_0}{B^2} \frac{\partial p}{\partial r} \mathbf{b} \times \nabla r \right] \qquad (35)$$

where $p$ is the total plasma pressure entering the Grad-Shafranov equation. The source term due to the radial gradient of the maxwellian background in Eq. (23) is given by:

$$\nabla_S F_{Ms} = \left[ \frac{\nabla n_{s0}}{n_{s0}} + \frac{1}{T_s} \frac{\partial \mathcal{E}_s}{\partial r} + \left( \frac{m_s v_{\parallel}^2}{2T_s} + \frac{\mu_s B}{T_s} + \frac{\varepsilon_s}{T_s} - \frac{3}{2} \right) \frac{\nabla T_s}{T_s} + \frac{m_s}{T_s} \left( \frac{RB_t}{B} v_{\parallel} + \Omega(R^2 - R_{\theta=0}^2) \right) \nabla \omega_{\Phi} \right] F_{Ms}$$

Finally, $\mathcal{C}(f_s)$ is the collision operator.

## 2.7 The field equations

The Fourier components of the perturbed electromagnetic fields, $\hat{\phi}$, $\hat{A}_\parallel$ and $\hat{B}_\parallel$ entering the gyrokinetic equation are obtained from the Poisson and Ampère equations expressed in terms of the gyrocenter distribution function. The gyrokinetic Poisson equation is given by:

$$\left[ k_{\perp}^2 \lambda_D^2 + \sum_s Z_s^2 \frac{T_e}{T_s} \frac{n_s}{n_e} [1 - \Gamma_0(b_s)] \frac{e \hat{\phi}}{T_e} = \sum_s \frac{Z_s}{n_e} \int J_0(k_{\perp} \rho_s) \hat{f}_s d\mathbf{v} + \sum_s Z_s \frac{n_s}{n_e} [\Gamma_0(b_s) - \Gamma_1(b_s)] \frac{\hat{B}_\parallel}{B} \right] \qquad (36)$$

where the species density $n_s$ varies poloidally in the presence of centrifugal effects and is given by Eq. (22). The Debye length is $\lambda_D = \sqrt{\frac{\varepsilon_0}{T_e} n_e e^2}$ with $T_e$ and $n_e$ the electron temperature and density, respectively. The notation $\Gamma_n(b_s) = I_n(b_s) e^{-b_s}$ is used where $I_n$ is the n-th order modified Bessel function of the first kind and the parameter $b_s$ is defined as

$$b_s = \frac{1}{2} (k_{\perp} \rho_{ths})^2 \qquad (37)$$

where $\rho_{ths} = m_s v_{ths} / (Z_s e B)$ is the Larmor radius of species 's' at the thermal velocity. The parallel component of the gyrokinetic Ampère equation is given by:

$$k_{\perp}^2 \hat{A}_{\parallel} = \mu_0 \int Z_s e v_{\parallel} J_0(k_{\perp} \rho_s) \hat{f}_s d\mathbf{v} \qquad (38)$$
Finally, the perpendicular component of the gyrokinetic Ampère equation is given by:

$$ \left[ 1 + \sum_s \frac{2\mu_0 n_s T_s}{B^2} [\Gamma_0(b_s) - \Gamma_1(b_s)] \right] \hat{B}_\parallel = -\mu_0 \sum_s \int \mu_s \hat{J}_1(k_\perp \rho_s) \hat{f}_s d\mathbf{v} - \sum_s \frac{2\mu_0 n_s T_s}{B^2} [\Gamma_0(b_s) - \Gamma_1(b_s)] \frac{Z_s e \hat{\phi}}{2T_s} \quad (39) $$

## 2.8 Moments of the particle, guiding-center and gyrocenter distribution functions

Equation (23) evolves the perturbed distribution of gyrocenters, $f_s$, whereas the calculation of the perturbed electromagnetic fields and that of the turbulent fluxes involve the perturbed distribution of particles, $f_s^p$. The transformations required to obtain velocity space moments of the particle distribution function from that of the gyrocenter distribution function are detailed below. Generic velocity moments of the particle distribution function are first introduced as:

$$ \mathcal{M}_s^p(\mathbf{x}, \alpha_s) = \int \alpha_s(v_\parallel, v_\perp) f_s^p(\mathbf{x}, v_\parallel, v_\perp, \xi) \frac{d\mathbf{v}}{2\pi} d\xi \quad (40) $$

where $\mathbf{x}$ is the position vector, $\alpha_s$ is an arbitrary function of $v_\parallel$ and $v_\perp$, and $\xi$ is the gyro-angle. Moments of the particle distribution function are related to that of the guiding-center distribution function $f_s^g$ by:

$$ \mathcal{M}_s^p(\mathbf{x}, \alpha_s) = \int \alpha_s f_s^{gc}(\mathbf{X}, v_\parallel, v_\perp) \delta(\mathbf{X} + \rho_s - \mathbf{x}) \frac{d\mathbf{v}}{2\pi} d\xi d\mathbf{X} \quad (41) $$

where $\rho_s$ is the Larmor radius vector pointing from the guiding-center position to the corresponding particle position and $\delta$ is the Dirac function. The guiding-center distribution function $f_s^g$ is itself related to the gyrocenter distribution function $f_s$ by:

$$ f_s^{gc}(\mathbf{X}) = f_s(\mathbf{X}) - \frac{F_{Ms}(\mathbf{X})}{T_s} Z_s e [\phi(\mathbf{X} + \rho_s) - \bar{\phi}(\mathbf{X})] + \frac{F_{Ms}(\mathbf{X})}{T_s} \mu_s \bar{B}_\parallel(\mathbf{X}) \quad (42) $$

The terms involving $\phi$ and $B_\parallel$ account for the effect of polarisation and magnetisation of the background plasma, respectively. Inserting Eq. (42) into Eq. (41) and focusing on a single Fourier mode, we get:

$$ \hat{\mathcal{M}}_s^p(\alpha_s) = \hat{\mathcal{M}}_s^{J_0}(\alpha_s) - \frac{Z_s e \hat{\phi}}{T_s} \hat{\mathcal{M}}_s^{pol}(\alpha_s) + \frac{\hat{B}_\parallel}{T_s} \hat{\mathcal{M}}_s^{mag}(\alpha_s) \quad (43) $$

On the right hand side, the first term is the contribution from the perturbed distribution function of gyrocenters, the second stems from polarisation effects and the third from magnetisation effects:

$$ \hat{\mathcal{M}}_s^{J_0}(\alpha_s) = \int \alpha_s J_0(k_\perp \rho_s) \hat{f}_s d\mathbf{v} \quad (44) $$

$$ \hat{\mathcal{M}}_s^{pol}(\alpha_s) = \int \alpha_s [1 - J_0^2(k_\perp \rho_s)] F_{Ms} d\mathbf{v} \quad (45) $$

$$ \hat{\mathcal{M}}_s^{mag}(\alpha_s) = \int \alpha_s \mu_s J_0(k_\perp \rho_s) \hat{J}_1(k_\perp \rho_s) F_{Ms} d\mathbf{v} \quad (46) $$

(47)

We also introduce the $\hat{J}_1$ moments which will be useful to compute the fluxes driven by magnetic compression in the next section:

$$ \hat{\mathcal{M}}_s^{J_1}(\alpha_s) = \int \alpha_s \hat{J}_1(k_\perp \rho_s) \hat{f}_s d\mathbf{v} \quad (48) $$

and the moments of the gyrocenter distribution function which is part of the IMAS data model:

$$ \hat{\mathcal{M}}_s(\alpha_s) = \int \alpha_s \hat{f}_s d\mathbf{v} \quad (49) $$
## 2.9 Fluxes in the rotating frame

For each species, the flux-surface averaged contravariant radial component of the particle flux, energy flux and toroidal angular momentum flux are expressed in the rotating frame as a function of moments of the gyrocenter distribution function. The fluxes are counted positive when radially outward (for a transported quantity that is positive). The definitions given in [19, 20] are used, keeping only the non-vanishing contributions. The resulting fluxes are listed below and split as a function of the transporting field.

### 2.9.1 Particle flux

As detailed in Appendix B, the particle flux is given by:

$$ \Gamma_s^r = \sum_{k_x, k_y > 0} \left\langle 2k_y \mathcal{E}^{ry} \mathrm{Im} \left[ \int \hat{f}_s \hat{\chi}^\dagger \mathrm{d}v \right] \right\rangle \qquad (50) $$

where

$$ \hat{\chi} = J_0(k_{\perp} \rho_s)(\hat{\phi} - v_{\parallel} \hat{A}_{\parallel}) + \frac{\mu_s}{Z_s e} \hat{J}_1(k_{\perp} \rho_s) \hat{B}_{\parallel} \qquad (51) $$

is the generalised field potential, $\langle \cdot \rangle$ is the flux surface average, $\dagger$ denotes the complex conjugate and

$$ \mathcal{E}^{ry} = \frac{\mathbf{b} \times \nabla y}{B} \cdot \nabla r \qquad (52) $$

Separating the contributions to the particle flux as a function of the Fourier components and transporting fields, we get:

$$ \Gamma_s^r = \sum_{k_x, k_y > 0} \Gamma_{s, \phi}^r(k_x, k_y) + \Gamma_{s, A_{\parallel}}^r(k_x, k_y) + \Gamma_{s, B_{\parallel}}^r(k_x, k_y) $$

where

$$ \Gamma_{s, \phi}^r = \left\langle 2k_y \mathcal{E}^{ry} \mathrm{Im} \left[ \hat{M}_s^{J_0}(1) \hat{\phi}^\dagger \right] \right\rangle $$

$$ \Gamma_{s, A_{\parallel}}^r = - \left\langle 2k_y \mathcal{E}^{ry} \mathrm{Im} \left[ \hat{M}_s^{J_0}(v_{\parallel}) \hat{A}_{\parallel}^{\dagger} \right] \right\rangle $$

$$ \Gamma_{s, B_{\parallel}}^r = \left\langle 2k_y \mathcal{E}^{ry} \mathrm{Im} \left[ \hat{M}_s^{J_1} \left( \frac{\mu_s}{Z_s e} \right) \hat{B}_{\parallel}^{\dagger} \right] \right\rangle $$

### 2.9.2 Energy flux

Similarly, the energy flux is given by

$$ Q_s^r = \sum_{k_x, k_y > 0} \left\langle 2k_y \mathcal{E}^{ry} \mathrm{Im} \left[ \int \epsilon \hat{f}_s \hat{\chi}^\dagger \mathrm{d}v \right] \right\rangle \qquad (53) $$

with $\epsilon = \frac{1}{2} m_s v^2 + \mathcal{E}_s$, and its components are

$$ Q_{s, \phi}^r = \left\langle 2k_y \mathcal{E}^{ry} \mathrm{Im} \left[ \hat{M}_s^{J_0}(\epsilon) \hat{\phi}^\dagger \right] \right\rangle $$

$$ Q_{s, A_{\parallel}}^r = - \left\langle 2k_y \mathcal{E}^{ry} \mathrm{Im} \left[ \hat{M}_s^{J_0}(v_{\parallel} \epsilon) \hat{A}_{\parallel}^{\dagger} \right] \right\rangle $$

$$ Q_{s, B_{\parallel}}^r = \left\langle 2k_y \mathcal{E}^{ry} \mathrm{Im} \left[ \hat{M}_s^{J_1} \left( \frac{\mu_s}{Z_s e} \right) \hat{B}_{\parallel}^{\dagger} \right] \right\rangle $$
### 2.9.3 Momentum flux

For the momentum flux, the transported quantity is the toroidal angular momentum, defined positive for a flow in the direction of $\nabla\varphi$. The radial flux of toroidal angular momentum is split into two parts: one grouping the contributions arising from the parallel dynamics $\Pi_s^{\parallel,r}$ and the other one from the perpendicular dynamics $\Pi_s^{\perp,r}$.

For the parallel dynamics, we have:

$$ \Pi_{s,\phi}^{\parallel,r} = \left\langle 2k_y \mathcal{E}^{ry} s_b m_s R^{\frac{|B_t|}{B}} \mathrm{Im} \left[ \hat{\mathcal{M}}_s^{J_0}(v_{\parallel}) \hat{\phi}^{\dagger} \right] \right\rangle $$

$$ \Pi_{s,A_{\parallel}}^{\parallel,r} = - \left\langle 2k_y \mathcal{E}^{ry} s_b m_s R^{\frac{|B_t|}{B}} \mathrm{Im} \left[ \hat{\mathcal{M}}_s^{J_0}(v_{\parallel}^2) \hat{A}_{\parallel}^{\dagger} \right. \right. \\ \left. \left. - \frac{n_s}{2} v_{ths}^2 \left[ \frac{Z_s e \hat{\phi}}{T_s} \Gamma_0 + \frac{\hat{B}_{\parallel}}{B} [\Gamma_0 - \Gamma_1] \right] \hat{A}_{\parallel}^{\dagger} \right] \right\rangle $$

$$ \Pi_{s,B_{\parallel}}^{\parallel,r} = \left\langle 2k_y \mathcal{E}^{ry} s_b m_s R^{\frac{|B_t|}{B}} \mathrm{Im} \left[ \hat{\mathcal{M}}_s^{J_1} \left( \frac{\mu_s}{Z_s e} v_{\parallel} \right) \hat{B}_{\parallel}^{\dagger} \right] \right\rangle $$

where the arguments of $\Gamma_0(b_s)$ and $\Gamma_1(b_s)$ have been omitted for brevity. In contrast to the particle and energy fluxes, polarisation and magnetisation effects contribute to the parallel momentum flux (caused by magnetic flutter).

The contributions arising from the perpendicular dynamics are given by:

$$ \Pi_{s,\phi}^{\perp,r} = \left\langle 2nk_{\perp}^{r} \frac{m_s}{B} \left\{ \mathrm{Re} \left[ \hat{\mathcal{M}}_s^{J_1} \left( \frac{\mu_s}{Z_s e} \right) \hat{\phi}^{\dagger} \right] \right. \right. \\ \left. \left. + \frac{n_s}{B} (\Gamma_0 - \Gamma_1) |\hat{\phi}|^2 \right. \right. \\ \left. \left. + \frac{n_s T_s}{2Z_s e B^2} (\Gamma_0 - \Gamma_1) \mathrm{Re} \left[ \hat{B}_{\parallel}^{\dagger} \hat{\phi}^{\dagger} \right] \right\} \right\rangle $$

$$ \Pi_{s,A_{\parallel}}^{\perp,r} = - \left\langle 2nk_{\perp}^{r} \frac{m_s}{B} \mathrm{Re} \left[ \hat{\mathcal{M}}_s^{J_1} \left( \frac{\mu_s}{Z_s e} v_{\parallel} \right) \hat{A}_{\parallel}^{\dagger} \right] \right\rangle $$

$$ \Pi_{s,B_{\parallel}}^{\perp,r} = \left\langle 2 \frac{nk_{\perp}^{r}}{k_{\perp}^2} \left\{ \mathrm{Re} \left[ \left( \hat{\mathcal{M}}_s^{J_1}(\mu_s) - 2\hat{\mathcal{M}}_s^{J_0}(\mu_s) \right) B_{\parallel}^{\dagger} \right] \right. \right. \\ \left. \left. + \frac{n_s Z_s e}{B} \left[ (2b_s - 1)\Gamma_0 - 3b_s\Gamma_1 \right] \mathrm{Re} \left[ \hat{\phi} \hat{B}_{\parallel}^{\dagger} \right] \right. \right. \\ \left. \left. + \frac{n_s T_s}{B^2} |\hat{B}_{\parallel}|^2 \left[ \frac{1}{2}(\Gamma_0 - \Gamma_1) \right. \right. \\ \left. \left. - 2[2(1 - 4b_s)\Gamma_0 + (2b_s - 1)\Gamma_1] \right] \right\} \right\rangle \quad (54) $$

where $n$ is the toroidal mode number, related to the binormal wavevector by $2\pi n = \Delta y k_y$, and $k_{\perp}^r = k_{\perp} \cdot \nabla r = R_0(k_x g^{xx} + k_y g^{xy})$.

## 2.10 Fluxes in the Laboratory frame

In section 2.9, the fluxes were expressed in the frame rotating toroidally as a rigid body with angular frequency $\Omega$. They are related to the fluxes expressed in the Laboratory frame as follows:

$$ \Gamma_s^{r,L} = \Gamma_s^r \qquad (55) $$

$$ \Pi_s^{\parallel,r,L} = \Pi_s^{\parallel,r} + \sum_{k_x, k_y > 0} \left\langle 2k_y \mathcal{E}^{ry} m_s \Omega \left[ \frac{RB_t}{B} \right]^2 \mathrm{Im} \left[ \int \hat{f}_s \hat{\chi}^{\dagger} d\mathbf{v} \right] \right\rangle \quad (56) $$

$$ \Pi_s^{\perp,r,L} = \Pi_s^{\perp,r} + \sum_{k_x, k_y > 0} \left\langle 2k_y \mathcal{E}^{ry} m_s \Omega \left[ \frac{RB_p}{B} \right]^2 \mathrm{Im} \left[ \int \hat{f}_s \hat{\chi}^{\dagger} d\mathbf{v} \right] \right\rangle \quad (57) $$

$$ Q_s^{r,L} = Q_s^r + \Omega \left[ \Pi_s^{\parallel,r} + \Pi_s^{\perp,r} \right] + \sum_{k_x, k_y > 0} \left\langle 2k_y \mathcal{E}^{ry} \frac{1}{2} m_s R^2 \Omega^2 \mathrm{Im} \left[ \int \hat{f}_s \hat{\chi}^{\dagger} d\mathbf{v} \right] \right\rangle \quad (58) $$

where the $L$ superscript is used to label the Laboratory frame fluxes and $B_p$ is the poloidal magnetic field.
## 2.11 Transport equations

It is important to realize that the definition of the fluxes is linked to that of the transport equations in which they enter. The fluxes defined in Sec. 2.10 follow the definitions of $\Gamma_a^A$, $\Pi_a^A$ and $Q_a^A$ given by Eq. (53) in [20]. More precisely:

$$ \Gamma_a^A = \frac{1}{2\pi} \left| \frac{\partial \Psi}{\partial r} \right| \Gamma_s^{r,L} \qquad (59) $$

$$ \Pi_a^A = \frac{1}{2\pi} \left| \frac{\partial \Psi}{\partial r} \right| \left[ \Pi_s^{\parallel,r,L} + \Pi_s^{\perp,r,L} \right] \qquad (60) $$

$$ Q_a^A = \frac{1}{2\pi} \left| \frac{\partial \Psi}{\partial r} \right| Q_s^{r,L} \qquad (61) $$

The corresponding transport equations are Eqs. (55-57) in [20], which must be combined with Eqs. (60-63) of the same paper to reveal $\Gamma_a^A$, $\Pi_a^A$ and $Q_a^A$. As an example, the resulting transport equation for momentum transport, retaining only turbulent contributions and without sources, reads:

$$ \frac{\partial}{\partial t} \left\langle \sum_s m_s n_s R^2 \mathbf{v}_0 \cdot \nabla \varphi \right\rangle + \frac{1}{V'} \frac{\partial}{\partial r} \left[ V' \sum_s \Pi_s^{\parallel,r,L} + \Pi_s^{\perp,r,L} \right] + \frac{1}{V'} \frac{\partial}{\partial r} \left[ V' \Pi_M^r \right] = 0 \qquad (62) $$

where $V' = \partial V / \partial r$, with $V$ the volume enclosed by a flux surface, and $\Pi_M^r$ the part of the Maxwell stress not included in $\Pi_s^{r,L}$:

$$ \Pi_M^r = \sum_{k_x, k_y > 0} \frac{1}{\mu_0} \left\langle 2nk_{\perp}^r \left[ |\hat{A}_{\parallel}|^2 + \frac{1}{k_{\perp}^2} |\hat{B}_{\parallel}|^2 \right] \right\rangle \qquad (63) $$

# 3 The IMAS gyrokinetics_local IDS

The IMAS gyrokinetics_local Interface Data Structure (IDS) has been created specifically to provide a unified standard for flux-tube gyrokinetic simulations. This IDS is a subset of the IMAS Data Dictionary, which describes magnetic fusion devices physics and subsystems in a machine generic way. The gyrokinetics_local IDS (GK IDS) describes the inputs (density, temperature, etc.) and outputs (growth rate, fluxes, etc.) of numerical simulations solving the system of delta-f flux-tube gyrokinetic equations presented in section 2. It also contains all the metadata required to repeat a simulation. The architecture of the GK IDS, the naming and the definition of each quantity, including sign conventions and normalizations, are introduced in this section.

## 3.1 Architecture

An overview of the gyrokinetics_local IDS is given in Table 1. It shows the overall structure, together with the name, type and dimension of the data. It corresponds to tag 4.2.0 of the open source IMAS Data Dictionary [21]. Some generic elements which are mandatory in any IMAS IDS but not relevant in the present context are omitted. Their description can be found in the online documentation of the IMAS Data Dictionary [22]. The content of a few substructures omitted, for brevity, in the overview table is detailed in Table 2. Based on the data model, an instance of the GK IDS can be generated for various programming languages and filled with results from gyrokinetic simulations. A single GK IDS instance can be used to store a set of linear simulations and one non-linear simulation provided they are performed with the same code, including the same physics and sharing the same physical inputs. Simulations from different codes or with different inputs are stored in different IDS instances. In an IDS instance, not all the elements of the data model needs to be filled: when data is unavailable, they are left at their default "empty" value (in case of scalars) or are not allocated (in case of arrays).
<table><caption>**Table 2:** Content of the moment, field and fluxes/linear weights structures. The size of the corresponding arrays is specified in the overview table.</caption><thead><tr><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td><strong>moments_norm_*</strong></td><td>struct</td><td><strong>linear_weights_* fluxes_*</strong></td><td>struct</td></tr><tr><td>density</td><td>cpx</td><td>particlesier.potential</td><td>flt</td></tr><tr><td>j_parallel</td><td>cpx</td><td>particlesierfield_parallel</td><td>flt</td></tr><tr><td>pressure_parallel</td><td>cpx</td><td>particlesierfield_parallel</td><td>flt</td></tr><tr><td>pressure_perpendicular</td><td>cpx</td><td>energyier.potential</td><td>flt</td></tr><tr><td>pressure</td><td>cpx</td><td>energyierfield_parallel</td><td>flt</td></tr><tr><td>heat flux_parallel</td><td>cpx</td><td>energyierfield_parallel</td><td>flt</td></tr><tr><td>v_parallel_energy_perpendicular</td><td>cpx</td><td>momentumier_ parallel_ phi potential</td><td>flt</td></tr><tr><td>v_perpendicular_square_energy</td><td>cpx</td><td>momentumier_ parallel_ a field_parallel</td><td>flt</td></tr><tr><td><strong>fields_*</strong></td><td>struct</td><td>momentumier_ parallel_ b field_parallel</td><td>flt</td></tr><tr><td>phi potential_ perturbed norm</td><td>flt/cpx</td><td>momentumier_ perpendicular_ phi potential</td><td>flt</td></tr><tr><td>a field_parallel_ perturbed norm</td><td>flt/cpx</td><td>momentumier_ perpendicular_ a field_parallel</td><td>flt</td></tr><tr><td>b field_parallel_ perturbed norm</td><td>flt/cpx</td><td>momentumier_ perpendicular_ b field_parallel</td><td>flt</td></tr></tbody></table>

## 3.2 Why a data model with normalized quantities?

Most of the IMAS Data Dictionary describes dimensional quantities expressed in SI units. For the gyrokinetics_local IDS, however, it was chosen to use normalized quantities. This is because the local version of the gyrokinetic system of equations is based on an expansion in the ratio of the ion Larmor radius to the plasma size, $\rho_*$, which is assumed to be small and does not appear anymore in the normalized equations. The solution is then valid for any value of $\rho_*$ (provided it is small) and this value only needs to be specified to convert the solution to physical units. In practice, the normalized version therefore provides the solution for an infinite number of $\rho_*$ values. Said differently, there is an infinite set of physical inputs that yield the same normalized inputs and hence the same normalized solution. This is why it is more efficient to use the normalized version in the IMAS gyrokinetics_local IDS.

## 3.3 Normalizations

All quantities in the data model are normalized using six reference quantities. The definition of these six fundamental reference quantities are:

* $q_{ref} = e = 1.602176634 \times 10^{-19}$ C, the positive elementary charge

* $m_{ref} = m_D = 3.3435837724 \times 10^{-27}$ kg, the Deuterium mass

* $T_{ref} = T_e(r_0, \theta = 0) = T_{e0}$, the electron temperature at $\theta = 0$

* $n_{ref} = n_e(r_0, \theta = 0) = n_{e0}$, the electron density at $\theta = 0$. When several electron species are included, e.g. thermal and fast electrons, $n_{e0}$ is the total electron density.

* $L_{ref} = R_0$, the major radius of the flux surface centre

* $B_{ref} = |B_t(R_0)| > 0$, the toroidal magnetic field at $R = R_0$, for the flux surface under consideration

Following IMAS conventions, the values of the fundamental physical constants are the CODATA recommended values by the National Institute of Standards and Technology [23]. The reference temperature, density, major radius and magnetic field are local to the flux surface of interest. For convenience, the following derived reference quantities are also defined:

* $v_{thref} = \sqrt{\frac{2T_{ref}}{m_{ref}}}$, the reference thermal velocity

* $\rho_{ref} = \frac{m_{ref}v_{thref}}{q_{ref}B_{ref}}$, the reference Larmor radius

If needed to convert the inputs of a gyrokinetic code to IMAS, the electron mass is taken to be $m_e = 9.1093837015 \times 10^{-31}$ kg.
## 3.4 Meta-data

The three structures `ids_properties`, `code` and `normalizing quantities` give information about the origin of the data.

### 3.4.1 IDS properties

In IMAS, each IDS contains a generic `ids_properties` structure which provides information about how the instance of the IDS was generated. A complete description of the `ids_properties` structure can be found in the IMAS documentation [22]. We only describe here the four elements relevant in the present context:

* **provider**

Name of the person who has produced the data.

* **creation_date**

Date and time at which the IDS instance has been produced with format YYYY-MM-DDTHH:MM:SSZ.
Example : 2020-07-24T14:19:00Z

* **provenance**

This generic IMAS structure can be used to keep track of the original files used to generate this IDS instance, i.e. the inputs and outputs of a gyrokinetic code, including the date and time at which they were generated.

* **comment**

Any comment describing the content of this IDS instance.

* **homogeneous_time**

Set to 2, since the `gyrokinetics_local` IDS is static (it does not describe the dynamics of a plasma pulse).

### 3.4.2 Code

A `code` structure appears in three places: at the root, under the `linear` structure and under the `non_linearity` structure. The one located at the root contains information that pertains to the whole GK IDS. It is also a generic IMAS structure and we document here a few of its elements, see [22] for a comprehensive description.

* **name**

Name of the gyrokinetic code used.

* **repository**

URL of the gyrokinetic code repository.

* **commit**

Unique commit reference for the source of the gyrokinetic code.

* **version**

Unique version (tag) of the gyrokinetic code.

* **library**

List of external libraries used by the gyrokinetic code. If the conversion to IMAS is done in post-processing, the name, description, repository and commit version of the post-processing scripts can be stored here.

The two additional `code` structure are partial code structures. They only contains the input files used of the simulations (stored in the form of an XML string) and a flag indicating whether the simulation was successful or not.
* **parameters**
  List of the code inputs in XML format.

* **output_flag**
  Set to 0 if the run was successful. Other values mean some difficulty has been encountered, with the exact meaning code specific. Negative values mean that the result shall not be used.

There is one partial code structure for each simulation stored in the IDS instance.

### 3.4.3 Normalizing quantities

In principle, there is no need to store the value of the normalizing temperature, density, magnetic field and length, only their definition is required. When modelling an experimental case, it can nevertheless be convenient to store the value of the normalizing quantities used to prepare the inputs for the simulations. This can be done in the normalizing-quantities structure:

* **t_e**
  Normalizing temperature, $T_{e0}$ in [eV].

* **n_e**
  Normalizing density, $n_{e0}$ in [m$^{-3}$].

* **b_field**
  Normalizing magnetic field, $B_t(R_0)$ in [T].

* **r**
  Normalizing length, $R_0$ in [m].

## 3.5 Model

The model structure gives an overview of the physical effects included in the simulation. It contains the following switches:

* **adiabatic_electrons**
  Set to 1 if electrons are assumed to be adiabatic, 0 otherwise.

* **include_a_field_parallel**
  Set to 1 if the fluctuations of the parallel vector potential $A_\parallel$ are included, 0 otherwise.

* **include_b_field_parallel**
  Set to 1 if the fluctuations of the parallel magnetic field $B_\parallel$ are included, 0 otherwise.

* **use_mhd_rule**
  Set to 1 if the geometrical dependence of the $\nabla B$ drift is approximated to be that of the curvature drift, 0 otherwise.

$$v_{\nabla B} = \frac{\mu_s \mathbf{b} \times \nabla B}{Z_s e B} \sim \frac{\mu_s}{Z_s e} \mathbf{b} \times (\mathbf{b} \cdot \nabla)\mathbf{b}, \quad (64)$$

Using `use_mhd_rule = 1` only makes sense in combination with the neglect of $B_\parallel$ fluctuations. It can be used to approximate the effect of retaining $B_\parallel$ fluctuations thanks to a cancellation at the MHD level [24].

* **include_coriolis_drift**
  Set to 1 if the Coriolis drift, see Eq. (33), is included, 0 otherwise.

* **include_centrifugal_effects**
  Set to 1 if centrifugal effects are included, 0 otherwise.
* `collisions pitch_only`

Set to 1 if pitch-angle scattering only is retained in the collisions, 0 otherwise.

* `collisions momentum_conservation`

Set to 1 if momentum conservation is enforced in the collision operator, 0 otherwise.

* `collisions energy_conservation`

Set to 1 if energy conservation is enforced in the collision operator, 0 otherwise.

* `collisions finite larmor radius`

Set to 1 if finite Larmor radius effects are retained in the collision operator, 0 otherwise.

## 3.6 Species parameters

The parameters relative to the plasma species are grouped into three structures: `species_all`, `species` and `collisions`. The content of these structures is detailed below.

### 3.6.1 Species all

All species are assumed to have a purely toroidal flow with a common toroidal angular frequency given by Eq. (19). This species independent toroidal flow and the associated $E \times B$ shearing rate are specified in the `species_all` structure, together with the plasma beta (ratio of plasma pressure to magnetic energy) and the Debye length.

* `velocity phi norm`

$$u_N = \frac{R_0}{v_{thref}} \omega_{\Phi}$$

Positive for a flow in the direction of $\nabla\varphi$.

* `shearing rate norm`

$$\gamma_{E_N} = -\frac{R_0}{v_{thref}} \frac{1}{B_{ref}} \left(\frac{\partial\Psi}{\partial r}\right)^2 \frac{\partial^2\Phi_0}{\partial\Psi^2}$$

* `beta reference`

$$\beta_{eN} = 2\mu_0 \frac{n_{ref} T_{ref}}{B_{ref}^2}$$

The plasma beta entering the Ampère's equations (38) and (39) is obtained from the reference value above combined with the species parameters (normalized density and temperature). When $\beta_{eN} = 0$, the run is electrostatic.

* `debye length norm`

$$\lambda_{DN} = \frac{\lambda_D}{\rho_{ref}} = \frac{1}{\rho_{ref}} \sqrt{\frac{\epsilon_0 T_{ref}}{n_{ref} q_{ref}^2}}$$

The Debye length enters the gyrokinetic Poisson equation (36). When neglected, this term is set to zero. In practice, it only matters for small scale ETG modes.

* `angle pol equilibrium`

This is a poloidal angle grid, $\theta_{cf}$ from $-\pi$ to $\pi$, common to all species. It is used to specify the poloidal variation of the species effective potential energy and of its gradient in the presence of centrifugal effects.
### 3.6.2 Species

For each species s, the background distribution function is assumed to be a local Maxwellian, see Eq. (21). Poloidal asymmetry in the background density induced by centrifugal effects is allowed following Eq. (22). Poloidally asymmetric quantities are given in input at $\theta = 0$. An arbitrary number of kinetic species is allowed, including several instances of the same species. The species parameters are stored under **species** which is an array of $n_s$ structures where $n_s$ is the number of species.

* charge_norm

$$Z_s = \frac{Z_s e}{q_{ref}}$$

* mass_norm

$$m_{sN} = \frac{m_s}{m_{ref}}$$

* density_norm

$$n_{sN} = \frac{n_{s0}}{n_{ref}}$$

with $n_{s0} = n_s(\theta = 0)$.

* density_loggradient_norm

$$R/L_{n_{s0}} = -\frac{R_0}{n_{s0}} \frac{\partial n_{s0}}{\partial r}$$

Positive for peaked profiles.

* temperature_norm

$$T_{sN} = \frac{T_s}{T_{ref}}$$

* temperature_loggradient_norm

$$R/L_{T_s} = -\frac{R_0}{T_s} \frac{\partial T_s}{\partial r}$$

Positive for peaked profiles.

* potential_energy_norm

$$\mathcal{E}_{sN}(\theta_{cf}) = \frac{\mathcal{E}_s(\theta_{cf})}{T_{ref}}$$

The effective potential energy responsible for the poloidal density asymmetry in Eq. (22), with $\mathcal{E}_s(\theta_{cf} = 0) = 0$ by definition.

* potential_energygradient_norm

$$\mathcal{E}'_{sN}(\theta_{cf}) = \frac{R_0}{T_{ref}} \frac{\partial \mathcal{E}_s}{\partial r}(\theta_{cf})$$

* velocity_p__ gradient_norm

$$u'_{sN} = -\frac{R_0^2}{v_{thref}} \frac{\partial \omega_{\varphi,s}}{\partial r}$$

Positive for peaked profiles of a flow in the direction of $\nabla\varphi$. This input is used to specify the parallel flow shear entering the radial gradient of the background Maxwellian distribution function, Eq. (36). It is assumed to arise entirely from the radial shear in a toroidal flow of constant angular frequency on a flux surface. The possibility to have a species dependent toroidal velocity gradient with $\omega_{\varphi,s} \neq \omega_{\Phi}$ is retained. Strictly speaking, this is inconsistent with the local limit assumption $\rho_* \to 0$ but can be interesting for the interpretation of tokamak experiments. Note that in the strict application of the local limit, i.e. for $\omega_{\varphi,s}(r) = \omega_{\Phi}(r)$, one has:

$$\gamma_{E_N} = \frac{1}{2\pi R_0 B_{ref}} \frac{\partial \Psi}{\partial r} u'_{sN} = -\frac{RB_p}{R_0 B_{ref}} u'_{sN}$$
### 3.6.3 Collisions

To help capturing the variety of collisions implementation in gyrokinetic codes, the normalized collisionality is given for each species couple a/b in an array of size $n_s \times n_s$ (lines for species a, columns for species b):

* `collisionality_norm`

$$ \nu_N^{a/b} = \frac{R_0}{v_{thref}} \nu^{a/b} $$

where $\nu^{a/b}$ is the collisionality frequency of species a against species b. For collision frequencies consistent with the species temperature and density, we have:

$$ \nu_N^{a/b} = \frac{R_0}{v_{thref}} \frac{n_{b0} Z_a^2 Z_b^2 e^4 \ln \Lambda^{a/b}}{4\pi \epsilon_0^2 m_a^2 v_{tha}^3} $$

where $n_{b0}$ is the density of species b at $\theta = 0$ and $\ln \Lambda^{a/b}$ the Coulomb logarithm. For electron-ion collisions, an extra $Z_{eff}$ factor can be added to the formula above to mimic the impact of heavier impurities, as is frequently done in gyrokinetics codes.

In addition, specific collisionality switches are described in the `model` structure, see Sec. 3.5.

## 3.7 Magnetic equilibrium

The magnetic equilibrium in the vicinity of the flux surface of interest is specified by the parameters stored in the `fluxsurface` structure, given at $r = r_0$.

* `r_minor_norm`

$$ r_N = \frac{r_0}{R_0} $$

* `b_fieldphi_sign`

$$ s_b = \text{sign}[\mathbf{B} \cdot \nabla \varphi] $$

* `ip_sign`

$$ s_j = \text{sign}[\mathbf{j} \cdot \nabla \varphi] $$

with j the plasma current density.

* `q`

$$ q = \frac{1}{2\pi} \int_{0}^{2\pi} \frac{\mathbf{B} \cdot \nabla \varphi}{\mathbf{B} \cdot \nabla \theta} d\theta $$

The safety factor q can be positive or negative, its sign depends on the sign of the toroidal plasma current and magnetic field:

$$ q = s_b s_j |q| $$

* `magnetic_shear_r_minor`

$$ \hat{s} = \frac{r}{q} \frac{\partial q}{\partial r} $$

* `pressureGradient_norm`

$$ p'_{N} = - \frac{2 \mu_{0} R_{0} \partial p}{B_{ref}^{2} \partial r} $$

with p the total plasma pressure (i.e. the pressure term entering the Grad-Shafranov equation). $p'_{N}$ (often called $\beta'$) is a quantity related to the magnetic equilibrium, which for instance affects the curvature drift, see Eq (35). It does not need to be consistent with the gradients ($R/L_{ns}$, $R/L_{Ts}$) specified for the kinetic species.
In addition to the parameters above, the shape of the flux surface is specified using an extension of the Miller parametrisation of local magnetic equilibria proposed in [25] and particularly convenient for local gyrokinetic codes. The plasma shape is described by a parametrisation of the constant poloidal flux contours:

$$R_{\Psi}(r, \theta) = R_{c}(r) + r \cos \theta_{R} \qquad (65)$$

$$Z_{\Psi}(r, \theta) = Z_{c}(r) - r\kappa(r) \sin \theta \qquad (66)$$

with the angle $\theta_R$ defined as:

$$\theta_R(r, \theta) = \theta + \sum_{n=0}^{n_{sh}-1} c_n(r) \cos n\theta + s_n(r) \sin n\theta \qquad (67)$$

The flux surface center ($R_c, Z_c$) and radial coordinate $r$ are defined in Sec. 2.3. Note that the flux surface center ($R_c, Z_c$) is defined here for each flux surface and therefore depends on the radial coordinate $r$. Since $Z_0 = Z_c(r_0)$ does not enter the gyrokinetic equations, it is assumed to be zero. The elongation $\kappa$ is defined as:

$$\kappa(r) = \frac{\max[Z_{\Psi}] - \min[Z_{\Psi}]}{\max[R_{\Psi}] - \min[R_{\Psi}]} \qquad (68)$$

The specification of the local magnetic equilibrium then requires $4n_{sh}+3$ terms ($s_0$ and its derivative are never needed), given at $r = r_0$, where $n_{sh}$ is the degree of the expansion:

* dgeometric_axis_r_dr_minor
  $$\frac{\partial R_c}{\partial r}(r_0)$$

* dgeometric_axis_z_dr_minor
  $$\frac{\partial Z_c}{\partial r}(r_0)$$

* elongation
  $$\kappa(r_0)$$

* delongation_dr_minor_norm
  $$R_0 \frac{\partial \kappa}{\partial r}(r_0)$$

* shape Coefficients_c
  $$c_n(r_0)$$

* shape Coefficients_s
  $$s_n(r_0)$$

* dc_dr_minor_norm
  $$R_0 \frac{\partial c_n}{\partial r}(r_0)$$

* ds_dr_minor_norm
  $$R_0 \frac{\partial s_n}{\partial r}(r_0)$$

Typically $n_{sh} \le 4$ is sufficient to capture most of the plasma shapes but higher values can be used when needed (for instance close to the edge).
## 3.8 Linear runs

The outputs of linear simulations are stored in the **linear** structure, together with a few specific inputs. The array of structures **wavevector** is placed at the root. It contains the radial and poloidal wavevectors of the considered Fourier component, see Eq. (14), normalized and expressed in a coordinate independent form:

* radialWavevector_norm

$$k_r \rho_{ref} = k_x \sqrt{g^{xx}(\theta=0)} \rho_{ref} \qquad (69)$$

* binormalWavevector_norm

$$k_{\theta} \rho_{ref} = k_y \sqrt{g^{yy}(\theta=0)} \rho_{ref} \qquad (70)$$

Remembering Eq. (17), $k_r$ is the dimensional perpendicular wave vector at $\theta = 0$ obtained when $k_y = 0$ and similarly, $k_{\theta}$ is the dimensional perpendicular wave vector at $\theta = 0$ obtained when $k_x = 0$. By convention, only positive binormal wavevector are stored in the GK IDS since this is sufficient to describe any real field owing to hermitian symmetry. The ($k_r\rho_{ref}, k_{\theta}\rho_{ref}$) pairs stored in the **wavevector** array of structure are expected to be unique: no duplicated wavevectors are allowed. For each wavevector, multiple eigenmodes are solution of the linear problem. These eigenmodes are stored under the **eigenmode** array of structures. To capture the case where each eigenmode is the result of a dedicated simulation, the corresponding code specific inputs are stored in the **code** structure under **eigenmode**, (see Sec. 3.4.2) together with the following quantities:

* angle_pol

This is the extended ballooning angle, $\theta$, used as the parallel coordinate in the flux-tube coordinate system $(x, y, \theta)$ introduced in Sec. 2.4.

* poloidal Turns

The number $N_p$ of poloidal turns (or coupled radial modes). In the spectral representation of a flux-tube domain, the periodicity of the torus in the poloidal direction leads to the coupling of radial modes in the parallel direction:

$$\hat{f}_s(k_x, k_y, \theta - \pi) = \hat{f}_s(k_x + \Delta k_x, k_y, \theta + \pi)$$

with $\Delta k_x$ depending on the magnetic shear. A flux-tube extending over an odd number $N_p$ of poloidal turns with 1 radial mode $k_x$ is therefore equivalent to a flux-tube extending over 1 poloidal turn with $N_p$ coupled radial modes $k_x + p\Delta k_x$ with $p$ an integer varying from $-(N_p - 1)/2$ to $(N_p + 1)/2$. In practice, running with a sufficiently high number of poloidal turns (or coupled radial modes) is necessary to get a result independent of the boundary condition imposed at the end of the field line (when the number of poloidal turns increases, $k_{\perp}$ at the end of the field line increases and the perturbations are damped by finite Larmor radius effects).

* time_norm

$$t_N = t \frac{v_{thref}}{R_0}$$

This is the normalized time at which the outputs are given. It can contain a single value when the outputs are given at the end of the linear simulation or an array of values in case the temporal evolution is stored. The second option can be useful to store the results of linear zonal-flow damping simulations (Rosenbluth-Hinton test).

* initial_value_run

Set to 1 for an initial value run and to 0 for an eigenvalue run.
### 3.8.1 Mode amplitude

In linear runs, unstable modes are exponentially growing in time and so are the associated fluxes. The exponentially growing outputs are therefore normalized to a dimensionless mode amplitude defined as:

$$ \hat{A}_f = \sqrt{\int [|\hat{\phi}_N(\theta)|^2 + |\hat{A}_{\parallel N}(\theta)|^2 + |\hat{B}_{\parallel N}(\theta)|^2] \frac{d\theta}{2\pi}} \quad (71) $$

where the Fourier components of the electrostatic potential and vector potential are normalized as follows:

$$ \hat{\phi}_N = \frac{1}{\rho_*} \frac{e}{T_{e0}} \hat{\phi} \quad (72) $$

$$ \hat{A}_{\parallel N} = \frac{1}{\rho_*^2 R_0 B_{ref}} \hat{A}_{\parallel} \quad (73) $$

$$ \hat{B}_{\parallel N} = \frac{1}{\rho_* B_{ref}} \hat{B}_{\parallel} \quad (74) $$

The integral over $\theta$ is performed over the full length of the flux-tube (it is equivalent to an integral over one poloidal turn of the sum of all connected $k_x$ modes).

### 3.8.2 Eigenvalues

The normalized mode growth rate $\gamma$ and frequency $\omega_r$ are defined as:

*   **growth_rate_norm**

$$ \gamma_N(k_x, k_y) = \gamma(k_x, k_y) \frac{R_0}{v_{thref}} $$

Positive for unstable modes.

*   **frequency_norm**

$$ \omega_{rN}(k_x, k_y) = \omega_r(k_x, k_y) \frac{R_0}{v_{thref}} $$

Positive for modes propagating in the ion diamagnetic drift direction (and therefore negative for modes propagating in the electron diamagnetic drift direction).

For initial value runs, a tolerance on the run convergence, calculable in post-processing, is defined as:

*   **growth_rate{tolerance**

$$ \gamma_{tol} = \frac{1}{|\gamma_f|} \sqrt{\frac{1}{\Delta t} \int_{t_f - \Delta t}^{t_f} [\gamma(t) - \gamma_f]^2 dt} $$

where $t_f$ and $\gamma_f$ are the time and growth rate at the end of the simulation, respectively. The time dependent growth rate is given by

$$ \gamma(t) = \frac{1}{\delta t} [\ln \hat{A}_f(t) - \ln \hat{A}_f(t - \delta t)] $$

with $\delta t$ depending on the code setup but in any case several times smaller than the interval $\Delta t$ over which the standard deviation is computed. This time interval is taken to be

$$ \Delta t = \frac{1}{2} \max(\tau_{\parallel s}, \tau_{\perp s}) $$

where $\tau_{\parallel s}$ is the time for a thermal particle of species s to stream along the field line over one poloidal turn and $\tau_{\perp s}$ is the time for a thermal particle of species s to drift over one fluctuation wavelength:

$$ \tau_{\parallel s} = 2\pi q \sqrt{\frac{m_{sN}}{T_{sN}}} \frac{R_0}{v_{thref}} \quad (75) $$

$$ \tau_{\perp s} = \frac{2\pi}{k_{\theta*} \rho_{ref}} \frac{|Z_{sN}|}{T_{sN}} \frac{R_0}{v_{thref}} \quad (76) $$

Only non-trace species, with $n_{sN} > 0$ are retained for the computation of the time interval.
### 3.8.3 Eigenfunctions

To facilitate the comparison of the mode structure between different runs, the eigenfunctions are rotated in the complex plane by an arbitrary phase factor $\alpha$ chosen so that the real part of the normalized electric potential is positive and its imaginary part is zero at $\theta = \theta^*$, where $\theta^*$ is the $\theta$ value at which $|\hat{\phi}_N|$ is maximum. When there are several maxima, the lowest $\theta^*$ value is chosen. For a given eigenmode, all the parallel outputs (eigenfunctions and moments) are rotated in the complex plane using the same value of $\alpha$ computed at the last time slice of `time_norm` and defined as:

$$e^{i\alpha} = \frac{|\hat{\phi}_N(k_x, k_y, \theta^*)|}{\hat{\phi}_N(k_x, k_y, \theta^*)}$$

The normalized eigenfunctions and associated derived quantities are defined as:

* `phi potenti al perturbed norm`

$$\hat{\phi}_{Nf}(\theta) = \frac{\hat{\phi}_N(k_x, k_y, \theta)}{\hat{A}_f(k_x, k_y)} e^{i\alpha}$$

* `a field parallel perturbed norm`

$$\hat{A}_{\parallel Nf}(\theta) = \frac{\hat{A}_{\parallel N}(k_x, k_y, \theta)}{\hat{A}_f(k_x, k_y)} e^{i\alpha}$$

* `b field parallel perturbed norm`

$$\hat{B}_{\parallel Nf}(\theta) = \frac{\hat{B}_{\parallel N}(k_x, k_y, \theta)}{\hat{A}_f(k_x, k_y)} e^{i\alpha}$$

* `phi potential perturbed weight`

$$\hat{A}_{\phi} = \sqrt{\int |\hat{\phi}_{Nf}(\theta)|^2 \frac{d\theta}{2\pi}}$$

* `a field parallel perturbed weight`

$$\hat{A}_{A_{\parallel}} = \sqrt{\int |\hat{A}_{\parallel Nf}(\theta)|^2 \frac{d\theta}{2\pi}}$$

* `b field parallel perturbed weight`

$$\hat{A}_{B_{\parallel}} = \sqrt{\int |\hat{B}_{\parallel Nf}(\theta)|^2 \frac{d\theta}{2\pi}}$$

* `phi potential perturbed parity`

$$\pi_{\phi} = \frac{|\int \phi_{Nf}(\theta) d\theta|}{\int |\phi_{Nf}(\theta)| d\theta}$$

* `a field parallel perturbed parity`

$$\pi_{A_{\parallel}} = \frac{|\int A_{\parallel Nf}(\theta) d\theta|}{\int |A_{\parallel Nf}(\theta)| d\theta}$$

* `b field parallel perturbed parity`

$$\pi_{B_{\parallel}} = \frac{|\int B_{\parallel Nf}(\theta) d\theta|}{\int |B_{\parallel Nf}(\theta)| d\theta}$$
### 3.8.4 Moments

Moments of the particle and gyrocenter distribution functions expressed in the rotating frame are given as a function of the poloidal angle $\theta$ and time $t_N$ grids for each species and each eigenmode. Four structures are defined:

* moments_normosed.particle

* moments_normose.gyrocenter

* moments_normose.gyrocenter_bessel_0

* moments_normose.gyrocenter_bessel_1

which contains, respectively, the $\hat{M}_s^p(\alpha_s)$, $\hat{M}_s(\alpha_s)$, $\hat{M}_s^{J_0}(\alpha_s)$ and $\hat{M}_s^{J_1}(\alpha_s)$ moments introduced in Sec. 2.8. For each of the four structures, the normalized moments are computed using the values of $\alpha_s$ listed below:

* density

$$ \alpha_s = \frac{1}{n_{s0}\rho_*} \frac{e^{i\alpha}}{\hat{A}_f} $$

* j_parallel

$$ \alpha_s = Z_s v_{\parallel} \frac{1}{n_{s0} v_{thref} \rho_*} \frac{e^{i\alpha}}{\hat{A}_f} $$

* pressure_parallel

$$ \alpha_s = \frac{1}{3} m_s v_{\parallel}^2 \frac{1}{n_{s0} T_{e0} \rho_*} \frac{e^{i\alpha}}{\hat{A}_f} $$

* pressure_perpendicular

$$ \alpha_s = \frac{1}{3} m_s v_{\perp}^2 \frac{1}{n_{s0} T_{e0} \rho_*} \frac{e^{i\alpha}}{\hat{A}_f} $$

* pressure

$$ \alpha_s = \frac{1}{3} m_s (v_{\parallel}^2 + v_{\perp}^2) \frac{1}{n_{s0} T_{e0} \rho_*} \frac{e^{i\alpha}}{\hat{A}_f} $$

* heat flux_parallel

$$ \alpha_s = \frac{1}{2} m_s v_{\parallel} v_{\perp}^2 \frac{1}{n_{s0} v_{thref} T_{e0} \rho_*} \frac{e^{i\alpha}}{\hat{A}_f} $$

* v_parallel_energy_perpendicular

$$ \alpha_s = \frac{1}{2} m_s v_{\parallel} v_{\perp}^2 \frac{1}{n_{s0} v_{thref} T_{e0} \rho_*} \frac{e^{i\alpha}}{\hat{A}_f} $$

* v_perpendicular_square_energy

$$ \alpha_s = \frac{1}{2} m_s v_{\perp}^2 v_{\parallel}^2 \frac{1}{n_{s0} v_{thref}^2 T_{e0} \rho_*} \frac{e^{i\alpha}}{\hat{A}_f} $$

With these definitions, the perturbed pressure is simply the sum of the perturbed parallel and perpendicular pressures. The species density $n_{s0} = n_s(\theta = 0)$ rather than $n_{ref}$ is used for the normalization to allow for trace species with zero density.
### 3.8.5 Linear weights

For linear runs, the fluxes defined in Sec. 2.9 and 2.10 for each wavevector, species and transporting field, are normalized to the mode amplitude $\hat{A}_f(k_x, k_y)$ to yield the so called weights used in quasi-linear transport models. The structure `linear_weights_rotate frame` contains the linear weights expressed in the rotating frame using the definitions of Sec. 2.9 and made dimensionless as follows:

* `particles phi potential`

$$ \Gamma_{sN,\phi} = \frac{\Gamma_{s,\phi}^r}{n_{s0}v_{thref}\rho_*^2 \hat{A}_f^2} \frac{1}{\hat{A}_f^2} $$

* `particles a field parallel`

$$ \Gamma_{sN,A_\parallel} = \frac{\Gamma_{s,A_\parallel}^r}{n_{s0}v_{thref}\rho_*^2 \hat{A}_f^2} \frac{1}{\hat{A}_f^2} $$

* `particles b field parallel`

$$ \Gamma_{sN,B_\parallel} = \frac{\Gamma_{s,B_\parallel}^r}{n_{s0}v_{thref}\rho_*^2 \hat{A}_f^2} \frac{1}{\hat{A}_f^2} $$

* `energy phi potential`

$$ Q_{sN,\phi} = \frac{Q_{s,\phi}^r}{n_{s0}T_{e0}v_{thref}\rho_*^2 \hat{A}_f^2} \frac{1}{\hat{A}_f^2} $$

* `energy a field parallel`

$$ Q_{sN,A_\parallel} = \frac{Q_{s,A_\parallel}^r}{n_{s0}T_{e0}v_{thref}\rho_*^2 \hat{A}_f^2} \frac{1}{\hat{A}_f^2} $$

* `energy b field parallel`

$$ Q_{sN,B_\parallel} = \frac{Q_{s,B_\parallel}^r}{n_{s0}T_{e0}v_{thref}\rho_*^2 \hat{A}_f^2} \frac{1}{\hat{A}_f^2} $$

* `momentum phi parallel phi potential`

$$ \Pi_{sN,\phi}^{\parallel,r} = \frac{\Pi_{s,\phi}^{\parallel,r}}{n_{s0}m_{ref}R_0 v_{thref}^2 \rho_*^2 \hat{A}_f^2} \frac{1}{\hat{A}_f^2} $$

* `momentum phi parallel a field parallel`

$$ \Pi_{sN,A_\parallel}^{\parallel,r} = \frac{\Pi_{s,A_\parallel}^{\parallel,r}}{n_{s0}m_{ref}R_0 v_{thref}^2 \rho_*^2 \hat{A}_f^2} \frac{1}{\hat{A}_f^2} $$

* `momentum phi parallel b field parallel`

$$ \Pi_{sN,B_\parallel}^{\parallel,r} = \frac{\Pi_{s,B_\parallel}^{\parallel,r}}{n_{s0}m_{ref}R_0 v_{thref}^2 \rho_*^2 \hat{A}_f^2} \frac{1}{\hat{A}_f^2} $$

* `momentum phi perpendicular phi potential`

$$ \Pi_{sN,\phi}^{\perp,r} = \frac{\Pi_{s,\phi}^{\perp,r}}{n_{s0}m_{ref}R_0 v_{thref}^2 \rho_*^2 \hat{A}_f^2} \frac{1}{\hat{A}_f^2} $$
* momentum phi perpendicular a field parallel

$$ \Pi_{sN,A_{\parallel}}^{\perp} = \frac{\Pi_{s,A_{\parallel}}^{\perp,r}}{n_{s0}m_{ref}R_0v_{thref}^2\rho_{*}^2\hat{A}_f^2} \frac{1}{\hat{A}_f^2} $$

* momentum phi perpendicular b field parallel

$$ \Pi_{sN,B_{\parallel}}^{\perp} = \frac{\Pi_{s,B_{\parallel}}^{\perp,r}}{n_{s0}m_{ref}R_0v_{thref}^2\rho_{*}^2\hat{A}_f^2} \frac{1}{\hat{A}_f^2} $$

The structure linear weights contains the linear weights expressed in the Laboratory frame made dimensionless in the same way as above but using the unnormalized fluxes of Sec. 2.10. All linear weights are given at the last time step of time norm.

## 3.9 Non-linear runs

The outputs of non-linear runs, together with a few specific inputs, are stored in the non-linear structure. Unlike the linear structure, where the wavevectors are stored in an array of structures, a wavevector grid is used in the non-linear structure. Storing each wavevector separately is indeed desirable for linear simulations where all wavevectors are uncoupled but makes less sense and proves impractical for non-linear runs based on a grid with a large number of coupled wavevectors. The radial and binormal components of the wavevectors are therefore defined as 1D arrays given at $\theta = 0$ with the same definition as for linear runs, see Eqs. (69) and (70) in Sec. 3.8, and stored in:

* radial wavevector norm

* binormal wavevector norm

The combination of the 1D radial and binormal wavevector grids forms the 2D rectangular grid used in the non-linear simulations. The poloidal angle (parallel coordinate) and time grids

* angle pol

* time norm

have the same definition as for linear runs. Some of the outputs stored in the non-linear structure are time averaged. The time interval used for the average is given by:

* time interval norm

$$ [t_N^{start}, t_N^{stop}] $$

Time averages are computed on the interval starting at $t_N^{start}$ and ending at $t_N^{stop}$, with $t_N$ the normalized time.

Finally, a switch allowing to store the results of quasi-linear calculations is added for versatility:

* quasi-linear

Set to 1 if the fluxes and fields spectra are computed from a quasi-linear model. 0 if non-linear.

### 3.9.1 Fields

The fluctuating fields $\phi$, $A_{\parallel}$ and $B_{\parallel}$ depend on the binormal wavevector $k_y$, the radial wavevector $k_x$, the poloidal angle $\theta$ and time $t_N$. The 4D fields arrays (complex values) contain maximal information but can represent a very large amount of data. For flexibility, the IMAS gyrokinetics data model defines a hierarchy of structures to store the fluctuating fields (or their intensity) in arrays of reduced dimensionality. These structures are listed below with $\zeta$ used as a generic field used to explicit their content:
* fields_4D

$\hat{\zeta}(k_y, k_x, \theta, t_N)$

The fields themselves (complex values). Function of $k_y, k_x, \theta, t_N$.

* fields_intensity_3D

$$ \int_{t_{\text{start}}^N}^{t_{\text{stop}}^N} |\hat{\zeta}(k_y, k_x, \theta, t_N)|^2 \frac{dt_N}{\Delta t_{\text{av}}} $$

The time averaged fields intensity, with $\Delta t_{\text{av}} = t_{\text{stop}}^N - t_{\text{start}}^N$. Function of $k_y, k_x, \theta$.

* fields_intensity_2D\_surface_average

$$ \left\langle \int_{t_{\text{start}}^N}^{t_{\text{stop}}^N} |\hat{\zeta}(k_y, k_x, \theta, t_N)|^2 \frac{dt_N}{\Delta t_{\text{av}}} \right\rangle $$

The time averaged and flux surface averaged fields intensity. Function of $k_y, k_x$.

* fields_zonal_2D

$$ \langle \hat{\zeta}(k_y = 0, k_x, \theta, t_N) \rangle $$

The flux surface averaged zonal fields. Function of $k_x, t_N$.

* fields_intensity_1D

$$ \left\langle \sum_{k_x} \int_{t_{\text{start}}^N}^{t_{\text{stop}}^N} |\hat{\zeta}(k_y, k_x, \theta, t_N)|^2 \frac{dt_N}{\Delta t_{\text{av}}} \right\rangle $$

The time averaged and flux surface averaged fields intensity, summed over $k_x$. Function of $k_y$.

In all these structures, the field represented by $\hat{\zeta}$ is specified as follows:

* phi potenti al perturbed norm: $\hat{\zeta} = \hat{\phi}_N$

* a field parallel perturbed norm: $\hat{\zeta} = \hat{A}_{\parallel N}$

* b field parallel perturbed norm: $\hat{\zeta} = \hat{B}_{\parallel N}$

Note that for non-linear runs, the fields are not rotated in the complex plane and not normalized to the mode amplitude.

### 3.9.2 Moments

With validation against the experiments in mind, moments of the particle distribution functions are stored in the moments\_norm\_particle structure. These moments are listed in Table 2 and follow the definitions given in Sec. 2.8, except that they are not normalized to the mode amplitude, nor rotated in the complex plane. They are stored in a 5 dimensional array with coordinates $s, k_y, k_x, \theta, t_N$, with $s$ the species label. For now, moments of the gyrocenter distribution function are not included in the GK IDS for non-linear runs.

### 3.9.3 Fluxes

The fluxes are stored in arrays of various dimensionality grouped in a hierarchy of structures. A generic flux $\mathcal{F}$ and its primitive $\mathcal{F}_{\theta}$, with $\mathcal{F} = \langle \mathcal{F}_{\theta} \rangle$, are introduced to describe the content of the flux structures:

* fluxes_5D

$$ \mathcal{F}_{\theta}(s, k_y, k_x, \theta, t_N) $$

The fluxes primitives (non flux-surface averaged). Function of $s, k_y, k_x, \theta, t_N$, with $s$ the species label.
* fluxes_4D

$$ \int_{t_N^{start}}^{t_N^{stop}} \mathcal{F}_{\theta}(s, k_y, k_x, \theta, t_N) \frac{dt_N}{\Delta t_{av}} $$

The time averaged fluxes primitives. Function of s, k_y, k_x, θ.

* fluxes_3D

$$ \int_{t_N^{start}}^{t_N^{stop}} \mathcal{F}(s, k_y, k_x, t_N) \frac{dt_N}{\Delta t_{av}} $$

The time averaged fluxes. Function of s, k_y, k_x.

* fluxes_3D_k_x_sum

$$ \sum_{k_x} \mathcal{F}(s, k_y, k_x, t_N) $$

The time dependent fluxes summed over k_x. Function of s, k_y, t_N.

* fluxes_2D_k_x_sum

$$ \sum_{k_x} \int_{t_N^{start}}^{t_N^{stop}} \mathcal{F}(s, k_y, k_x, t_N) \frac{dt_N}{\Delta t_{av}} $$

The time averaged fluxes, summed over k_x. Function of s, k_y.

* fluxes_2D_k_x_k_y_sum

$$ \sum_{k_x, k_y > 0} \mathcal{F}(s, k_y, k_x, t_N) $$

The time dependent fluxes, summed over k_x and k_y. Function of s, t_N.

* fluxes_1D

$$ \sum_{k_x, k_y > 0} \int_{t_N^{start}}^{t_N^{stop}} \mathcal{F}(s, k_y, k_x, \theta, t_N) \frac{dt_N}{\Delta t_{av}} $$

The time averaged fluxes, summed over k_x and k_y. Function of s only.

These six structures are used to store the fluxes expressed in the Laboratory frame. An additional set of similar structures with the suffix _rotating_frame are used to store the corresponding fluxes expressed in the rotating frame. In all these structures, the fluxes represented by $\mathcal{F}$ are the ones listed in Sec. 3.8.5 with the same naming and normalizing conventions, excepted the absence of normalization with the mode amplitude ($A_f = 1$).

# 4 Application examples

The gyrofluid code TGLF and the gyrokinetic codes GENE and GKW have been interfaced with IMAS to read and write gyrokinetics_local IDSs. This is done in pre- and post-processing using Python scripts available in the following git repository [26]. Recently, the IMAS interface has been exploited to benchmark TGLF against GKW in high beta regimes prone to KBM turbulence [27]. As an additional illustration, we use the IMAS standard to compare the linear response of GENE, GKW and TGLF for a simple reference case defined in Table 3. Electric potential and perpendicular magnetic field fluctuations are included and magnetic field compression is neglected. The flux surfaces have a circular cross-section. Two kinetic species are included, electrons and deuterium ions, with a zero electron temperature gradient, no collisions and no rotation. For all simulations, the radial wavevector is set to zero. Grid convergence studies show that rather high resolution is required in GENE and GKW to obtain values of the mode growth rate converged within less than 1%. In velocity space, based on the convergence study, 24 points were chosen for the $\mu$ grid and 160 points for the $v_\|\|$ grid. The maximum velocity was set to 3 times the thermal species velocity for each species. The
<table><caption>**Table 3:** Inputs of the reference case used to compare GENE, GKW and TGLF using the IMAS standard.</caption><tbody><tr><td colspan="2"><strong>model</strong></td></tr><tr><td>adiabatic-electrons</td><td>0</td></tr><tr><td>include_a_field_parallel</td><td>1</td></tr><tr><td>include_b_field_parallel</td><td>0</td></tr><tr><td>use_mhd-approximation</td><td>0</td></tr><tr><td>include_coriolis_drift</td><td>0</td></tr><tr><td>include_centrifugal-effects</td><td>0</td></tr><tr><td>collisions-pitch_only</td><td>0</td></tr><tr><td>collisions_momentum_conservation</td><td>0</td></tr><tr><td>collisions_energy_conservation</td><td>0</td></tr><tr><td>collisionsfinite_larmor radius</td><td>0</td></tr><tr><td colspan="2"><strong>flux surface</strong></td></tr><tr><td>ip_sign</td><td>-1.0</td></tr><tr><td>b_field phi_sign</td><td>-1.0</td></tr><tr><td>r_minor_norm</td><td>1/6</td></tr><tr><td>q</td><td>2.0</td></tr><tr><td>magnetic_shear r_minor</td><td>1.0</td></tr><tr><td>pressuregradient_norm</td><td>0.0</td></tr><tr><td>dgeometric_axis r_dr_minor</td><td>0.0</td></tr><tr><td>dgeometric_axis z_dr_minor</td><td>0.0</td></tr><tr><td>elongation</td><td>1.0</td></tr><tr><td>delongation_dr_minor_norm</td><td>0.0</td></tr><tr><td>shapeCoefficients_c</td><td>[0.0,0.0]</td></tr><tr><td>dc_dr_minor_norm</td><td>[0.0,0.0]</td></tr><tr><td>shapeCoefficients_s</td><td>[0.0,0.0]</td></tr><tr><td>ds_dr_minor_norm</td><td>[0.0,0.0]</td></tr></tbody></table><table><tbody><tr><td colspan="2"><strong>species all</strong></td></tr><tr><td>velocity phi_norm</td><td>0.0</td></tr><tr><td>shearing rate_norm</td><td>0.0</td></tr><tr><td>beta_reference</td><td>5e-4</td></tr><tr><td>debye length_norm</td><td>0.0</td></tr><tr><td colspan="2"><strong>species</strong></td></tr><tr><td colspan="2">D</td></tr><tr><td>charge_norm</td><td>1.0</td></tr><tr><td>mass_norm</td><td>1.0</td></tr><tr><td>density_norm</td><td>1.0</td></tr><tr><td>density_loggradient_norm</td><td>3.0</td></tr><tr><td>temperature_norm</td><td>1.0</td></tr><tr><td>temperature_loggradient_norm</td><td>6.0</td></tr><tr><td>velocity phigradient_norm</td><td>0.0</td></tr><tr><td colspan="2"><strong>species</strong></td></tr><tr><td colspan="2">e<sup>-</sup></td></tr><tr><td>charge_norm</td><td>-1.0</td></tr><tr><td>mass_norm</td><td>2.724437e-4</td></tr><tr><td>density_norm</td><td>1.0</td></tr><tr><td>density_logGradient_norm</td><td>3.0</td></tr><tr><td>temperature_norm</td><td>1.0</td></tr><tr><td>temperature_logGradient_norm</td><td>0.0</td></tr><tr><td>velocity phiGradient_norm</td><td>0.0</td></tr><tr><td colspan="2"><strong>collisions</strong></td></tr><tr><td>collisionality_norm</td><td>[[0.0,0.0],[0.0,0.0]]</td></tr></tbody></table>

dynamics along the magnetic field line was discretised with 64 points per poloidal turns and a total of 31 poloidal turns. In GENE and GKW, numerical dissipation is applied in the $v_\parallel$ and $\theta$ directions, as is customary in finite-difference codes to prevent numerical instabilities. Default settings were used for TGLF. The linear growth rate and the mode frequency are shown as a function of the binormal wavevector in Fig. 2. The values obtained for GENE and GKW are nearly identical at all wavevectors for this Ion Temperature Gradient mode. The predictions of TGLF are within 10 to 20% from the gyrokinetic predictions, except at the highest wavevector where the real frequency computed by TGLF is significantly lower. In the left plot of Fig. 3, the real and imaginary part of the normalized electric potential fluctuations, $\dot{\phi}_{Nf}$, are shown as a function of the parallel coordinate $\theta$. Here as well, a perfect match is obtained between GENE and GKW. TGLF eigenfunction is very close to the gyrokinetic result. In the right plot of Fig. 3, the absolute value of the perturbed parallel current

<figure><img src="image_4.png" /><figcaption>**Figure 2:** Linear growth rate (left) and mode frequency (right) as a function of the binormal wave vector $k_\theta \rho_{ref}$ for GKW (blue crosses), GENE (orange circles) and TGLF (green stars).</figcaption></figure>

density carried by the electrons is shown as a function of the parallel coordinate. It reveals a very extended structure, with values still higher than 10% of the maximum after 30 poloidal turns. The overall agreement between GENE and GKW for the electron current density is excellent, with slight
<figure><img src="image_5.png" /><figcaption>**Figure 3:** Left: parallel structure of the real (full line) and imaginary (dashed line) part of the normalized electric potential $\hat{\phi}_{Nf}$ at $k_{\theta}\rho_{ref} = 0.2$ for GKW (blue), GENE (orange) and TGLF (green). Right: parallel structure of the particle current density moment for the electron species, $|\dot{M}_{eN}^{p}(-ev_{\|})|$ at $k_{\theta}\rho_{ref} = 0.7$ for GKW (blue) and GENE (orange).</figcaption></figure>

differences appearing close to the boundary of the flux tube, presumably due to a different numerical dissipation in the parallel direction. The parallel current density is not shown for TGLF since it is not a standard output for this code.

# 5 Summary and perspectives

A new Interface Data Structure, the gyrokinetics_local IDS, has been introduced in the IMAS Data Dictionary [21] to store inputs and outputs of local delta-f gyrokinetic simulations. This new standard has been documented in details, including all sign conventions, definitions and normalizations. The present version of the GK IDS, which corresponds to tag 4.2.0 of the IMAS Data Dictionary, is meant to cover basic needs and ensure reproducibility of the simulations. Many more fields could be added. For instance, one may wish to cover non-Maxwellian background distribution functions, poloidal asymmetries due to temperature anisotropy, some finite $\rho_*$ effects (e.g. profile shearing, parallel derivatives), storage of the perturbed distribution function or of additional moments, sonic flow modification of the magnetic equilibrium, specifics of the $E \times B$ shearing implementation, etc. However, trying to be exhaustive brings the risk of slowing down the deployment process and making the GK IDS heavier without any need. The choice has therefore been made to first deliver basic functionalities and keep additional development for future revisions. Importantly, the Data Dictionary is not meant to be frozen and extensions can be triggered by users to cover new needs.

As a first application, a public database of simulations relying on the gyrokinetics_local standard is presently being deployed [28]. It provides a practical tool to store and share gyrokinetic simulations, to perform gyrokinetic code benchmarks and to develop turbulent transport models based on the training of machine learning surrogates.

# Acknowledgments

We are grateful to Jonathan Citrin, Pierre Manas, Karel van de Plassche, Haomin Sun and Daan van Vugt for many fruitful discussions. This research was funded in part by l'Agence Nationale de la Recherche (ANR), project ANR-22-CE30-0023.

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.

## A Catalog of Bessel function integrals

The integrals of Bessel functions over velocity space used to compute fluxes and moments are listed in this section. We note $J_n$ the n-th order Bessel function of the first kind, $I_n$ the n-th order modified Bessel function of the first kind and $\Gamma_n(b_s) = I_n(b_s)e^{-b_s}$ with $b_s = \frac{1}{2}(k_\perp \rho_{ths})^2$. We also introduce $\hat{J}_1(k_\perp \rho_s) = 2J_1(k_\perp \rho_s)/(k_\perp \rho_s)$. In the following, the argument of the Bessel functions is omitted for brevity, i.e. $J_n$ is used for $J_n(k_\perp \rho_s)$ and $\Gamma_n$ for $\Gamma_n(b_s)$.

$$
\begin{align*}
\int v_{\parallel}^2 J_0^2 F_{Ms} \mathrm{d}\mathbf{v} &= n_s v_{ths}^2 \frac{1}{2} \Gamma_0 \\
\int \mu_s J_0^2 F_{Ms} \mathrm{d}\mathbf{v} &= n_s \frac{T_s}{B} [(1 - b_s)\Gamma_0 + b_s\Gamma_1] \\
\int \mu_s J_0 \hat{J}_1 F_{Ms} \mathrm{d}\mathbf{v} &= n_s \frac{T_s}{B} (\Gamma_0 - \Gamma_1) \\
\int \mu_s^2 J_0 \hat{J}_1 F_{Ms} \mathrm{d}\mathbf{v} &= n_s \frac{T_s^2}{B^2} [2(1 - b_s)\Gamma_0 + (2b_s - 1)\Gamma_1] \\
\int v_{\parallel}^2 \mu_s J_0 \hat{J}_1 F_{Ms} \mathrm{d}\mathbf{v} &= n_s v_{ths}^2 \frac{T_s}{B} \frac{1}{2} (\Gamma_0 - \Gamma_1) \\
\int \mu_s^2 \hat{J}_1^2 F_{Ms} \mathrm{d}\mathbf{v} &= n_s \frac{T_s^2}{B^2} \frac{1}{2} (\Gamma_0 - \Gamma_1)
\end{align*}
$$

## B Fluxes calculation

Following [19, 20], we first introduce the non-adiabatic distribution function of gyrocenters,

$$ \hat{h}_s = \hat{f}_s + \frac{Z_s e}{T_s} F_{Ms} J_0(k_{\perp} \rho_s) \hat{\phi} + \frac{\mu_s}{T_s} F_{Ms} \hat{J}_1(k_{\perp} \rho_s) \hat{B}_{\parallel} \quad (77) $$

and the generalised field potential,

$$ \hat{\chi} = J_0(k_{\perp} \rho_s)(\hat{\phi} - v_{\parallel} \hat{A}_{\parallel}) + \frac{\mu_s}{Z_s e} \hat{J}_1(k_{\perp} \rho_s) \hat{B}_{\parallel} \quad (78) $$

The contravariant radial component of the flux surface averaged particle flux is given by [19]:

$$ \Gamma_s^r = \sum_{k_x, k_y} \left\langle ik_y \mathcal{E}^{ry} \int \hat{h}_s^\dagger \hat{\chi} \mathrm{d}\mathbf{v} \right\rangle \quad (79) $$

where $\langle \cdot \rangle$ is the flux surface average and the factor

$$ \mathcal{E}^{ry} = \frac{\mathbf{b} \times \nabla y}{B} \cdot \nabla r \quad (80) $$

arises from

$$ \mathbf{v}_{\chi} \cdot \nabla r = \frac{\mathbf{b} \times \nabla \chi}{B} \cdot \nabla r = \sum_{k_x, k_y} ik_y \mathcal{E}^{ry} \hat{\chi}(k_x, k_y, \theta) \quad (81) $$

The Fourier components of any real field, for instance $\phi$, verify $\hat{\phi}_{k_x, -k_y}^\dagger \hat{\phi}_{k_x, -k_y} = \hat{\phi}_{-k_x, k_y} \hat{\phi}_{-k_x, k_y}^\dagger$. In addition, $k_\perp(k_x, -k_y) = k_\perp(-k_x, k_y)$, leading to:

$$
\begin{align}
\sum_{k_x, k_y} k_y J_0(k_\perp \rho_s) \hat{\phi}_{k_x, k_y}^\dagger \hat{\phi}_{k_x, k_y} &= \sum_{k_x, k_y>0} k_y J_0(k_\perp \rho_s) \hat{\phi}_{k_x, k_y}^\dagger \hat{\phi}_{k_x, k_y} \nonumber \\
&\quad - \sum_{k_x, k_y>0} k_y J_0(k_\perp \rho_s) \hat{\phi}_{k_x, -k_y}^\dagger \hat{\phi}_{k_x, -k_y} \nonumber \\
&= 0 \tag{82}
\end{align}
$$
As a consequence, all the terms of this form will vanish in Eq. (79). Noting that the integral over velocity space of $v_\parallel F_{Ms}$ is zero, the expression of the particle flux therefore reduces to

$$ \Gamma_s^r = \sum_{k_x, k_y > 0} \left\langle 2k_y \mathcal{E}^{ry} \mathrm{Im} \left[ \int \hat{f}_s \hat{\chi}^\dagger \mathrm{d}\mathbf{v} \right] \right\rangle \quad (83) $$

This means that the radial flux of particles given by Eq. (79) is equal to the radial flux of gyrocenters, Eq. (83). The contributions from polarisation and magnetisation cancel out. Similarly, the energy flux given in [19]

$$ Q_s^r = \sum_{k_x, k_y} \left\langle ik_y \mathcal{E}^{ry} \int \left( \frac{1}{2} m_s v^2 + \mathcal{E}_s \right) \hat{f}_s \hat{\chi}^\dagger \mathrm{d}\mathbf{v} \right\rangle \quad (84) $$

where $\mathcal{E}_s$ is the centrifugal energy entering Eq. (22), can be written as:

$$ Q_s^r = \sum_{k_x, k_y > 0} \left\langle 2k_y \mathcal{E}^{ry} \mathrm{Im} \left[ \int \left( \frac{1}{2} m_s v^2 + \mathcal{E}_s \right) \hat{f}_s \hat{\chi}^\dagger \mathrm{d}\mathbf{v} \right] \right\rangle \quad (85) $$

As for the particle flux, there is no contribution from polarisation or magnetisation effects. For the toroidal angular momentum flux, the situation is a bit more involved. We first consider the contribution arising from the parallel dynamics, expressed in the rotating frame:

$$ \Pi_s^{\parallel,r} = \sum_{k_x, k_y} \left\langle ik_y \mathcal{E}^{ry} \int s_b m_s R \frac{|B_t|}{B} v_{\parallel} \hat{h}_s^{\dagger} \hat{\chi} \mathrm{d}\mathbf{v} \right\rangle \quad (86) $$

Keeping only non-canceling components, we get:

$$ \Pi_s^{\parallel,r} = \sum_{k_x, k_y > 0} \left\langle 2k_y \mathcal{E}^{ry} s_b m_s R \frac{|B_t|}{B} \mathrm{Im} \left[ \int v_{\parallel} \hat{f}_s \hat{\chi}^{\dagger} \mathrm{d}\mathbf{v} \right. \right. \\ \left. - \frac{Z_s e}{T_s} \hat{\phi} \hat{A}_{\parallel}^{\dagger} \int v_{\parallel}^2 J_0^2 F_{Ms} \mathrm{d}\mathbf{v} \right. \\ \left. - \frac{1}{T_s} \hat{B}_{\parallel} \hat{A}_{\parallel}^{\dagger} \int \mu_s v_{\parallel}^2 J_0 \hat{J}_1 F_{Ms} \mathrm{d}\mathbf{v} \right] \right\rangle \quad (87) $$

Using the velocity space integrals of Bessel functions listed in Appendix A, this yields

$$ \Pi_s^{\parallel,r} = \sum_{k_x, k_y > 0} \left\langle 2k_y \mathcal{E}^{ry} s_b m_s R \frac{|B_t|}{B} \mathrm{Im} \left[ \int v_{\parallel} \hat{f}_s \hat{\chi}^{\dagger} \mathrm{d}\mathbf{v} - \frac{n_s}{2} v_{ths}^2 \left[ \frac{Z_s e \hat{\phi}}{T_s} \Gamma_0 + \frac{\hat{B}_{\parallel}}{B} [\Gamma_0 - \Gamma_1] \right] \hat{A}_{\parallel}^{\dagger} \right] \right\rangle \quad (88) $$

In addition to the radial flux of parallel momentum carried by the gyrocenters (first term under the integral in Eq. (88)), a momentum flux is driven by the interplay between polarisation and magnetisation of the background plasma and the magnetic flutter of the field lines (second and third terms under the integral in Eq. (88)).

The contribution arising from the perpendicular dynamics is given by

$$ \Pi_s^{\perp,r} = \sum_{k_x, k_y} \left\langle n \mathbf{k}_{\perp} \cdot \nabla r \frac{Z_s e}{k_{\perp}^2} \int \hat{h}_s^{\dagger} \hat{\chi}_a \mathrm{d}\mathbf{v} \right\rangle \quad (89) $$

where $n$ is the toroidal mode number, related to the binormal wavevector by $2\pi n = \Delta y k_y$ and an additional potential $\hat{\chi}_a$ has been introduced for convenience:

$$ \hat{\chi}_a = k_{\perp} \rho_s J_1(k_{\perp} \rho_s) [\hat{\phi} - v_{\parallel} \hat{A}_{\parallel}] + \frac{\mu_s}{Z_s e} [\hat{J}_1(k_{\perp} \rho_s) - 2J_0(k_{\perp} \rho_s)] \hat{B}_{\parallel} \quad (90) $$
Invoking again Hermitian symmetry, see Eq. (15), we notice that:

$$ \sum_{k_x, k_y} n \mathbf{k}_{\perp} \cdot \nabla r \hat{h}_s^{\dagger} \hat{\chi}_a = \sum_{k_x, k_y > 0} 2n \mathbf{k}_{\perp} \cdot \nabla r \operatorname{Re} [\hat{h}_s \hat{\chi}_a^{\dagger}] \quad (91) $$

We then have:

$$ \begin{aligned} \Pi_s^{\perp, r} = & \sum_{k_x, k_y > 0} \left\langle 2n\mathbf{k}_{\perp} \cdot \nabla r \left\{ \begin{aligned} & \operatorname{Re} \left[ \frac{m_s}{B} \int \frac{\mu_s}{Z_s e} \hat{J}_1 \hat{f}_s \hat{\phi}^{\dagger} d\mathbf{v} \right] \\ & + \frac{n_s m_s}{B^2} (\Gamma_0 - \Gamma_1) |\hat{\phi}|^2 \\ & + \frac{n_s m_s T_s}{2 Z_s e B^3} (\Gamma_0 - \Gamma_1) \operatorname{Re} [\hat{B}_{\parallel} \hat{\phi}^{\dagger}] \\ & - \operatorname{Re} \left[ \frac{m_s}{B} \int v_{\parallel} \frac{\mu_s}{Z_s e} \hat{J}_1 \hat{f}_s \hat{A}_{\parallel}^{\dagger} d\mathbf{v} \right] \\ & + \operatorname{Re} \left[ \frac{1}{k_{\perp}^2} \int \mu_s (\hat{J}_1 - 2J_0) \hat{f}_s \hat{B}_{\parallel}^{\dagger} d\mathbf{v} \right] \\ & + \frac{n_s Z_s e}{k_{\perp}^2 B} [(2b_s - 1)\Gamma_0 - 3b_s\Gamma_1] \operatorname{Re} [\hat{\phi} \hat{B}_{\parallel}^{\dagger}] \\ & + \frac{n_s T_s}{k_{\perp}^2 B^2} \frac{1}{2} (\Gamma_0 - \Gamma_1) |\hat{B}_{\parallel}|^2 \end{aligned} \right\rangle \\ & - \frac{2n_s T_s}{k_{\perp}^2 B^2} [2(1 - b_s)\Gamma_0 + (2b_s - 1)\Gamma_1] |\hat{B}_{\parallel}|^2 \end{aligned} \quad (92) $$

## References

[1] A. Brizard and T.S. Hahm. Foundations of nonlinear gyrokinetic theory. *Rev. Mod. Phys.*, 79:421, 2007.

[2] M. Yoshida (Chair Transport and Confinement), R.M. McDermott (co-Chair Transport and Confinement), C. Angioni, Y. Camenen, J. Citrin, M. Jakubowski, J.W. Hughes, Y. Idomura, P. Mantica, A. Mariani, S. Mordijck, E.J. Paul, T. Tala, G. Verdoolaege, A. Zocco, F.J. Casson, G. Dif-Pradalier, B. Duval, B.A. Grierson, S.M. Kaye, P. Manas, M. Maslov, T. Odstrcil, J.E. Rice, L. Schmitz, F. Sciortino, E.R. Solano, G. Staebler, M. Valovič, E. Wolfrum, J.A. Snipes (co-Chair Transport and Confinement) and the Transport and Confinement Group. Transport and confinement physics: Chapter 2 of the special issue: on the path to tokamak burning plasma operation. *Nucl. Fusion*, 65:033001, 2025.

[3] J. Candy, E.A. Belli, and R.V. Bravenec. A high-accuracy Eulerian gyrokinetic solver for collisional plasmas. *J. Comput. Phys.*, 324:73–93, 2016.

[4] F. Jenko, W. Dorland, M. Kotschenreuther, and B.N. Rogers. Electron temperature gradient driven turbulence. *Phys. Plasmas*, 7:1904, 2000.

[5] T-H Watanabe and H Sugama. Velocity-space structures of distribution function in toroidal ion temperature gradient turbulence. *Nucl. Fusion*, 46(1):24, 2005.

[6] A.G. Peeters, Y. Camenen, F.J. Casson, W.A. Hornsby, A.P. Snodin, D. Strintzi, and G. Szepesi. The non-linear gyrokinetic flux tube code GKW. *Comput. Phys. Commun.*, 180:2650, 2009.

[7] M. Kotschenreuther, G. Rewoldt, and W.M. Tang. Comparison of initial value and eigenvalue codes for kinetic toroidal plasma instabilities. *Comput. Phys. Commun.*, 88:128, 1995.

[8] J. Candy and R.E. Waltz. An Eulerian gyrokinetic-Maxwell solver. *J. Comput. Phys.*, 186:545, 2003.

[9] R.V. Bravenec, J. Candy, M. Barnes, and C. Holland. Linear and nonlinear verification of gyrokinetic microstability codes. *Phys. Plasmas*, 18:122505, 2011.
[10] G. Merlo, O. Sauter, S. Brunner, A. Burckel, Y. Camenen, F. J. Casson, W. Dorland, E. Fable, T. Goerler, F. Jenko, A. G. Peeters, D. Told, and L. Villard. Linear multispecies gyrokinetic flux tube benchmarks in shaped tokamak plasmas. *Phys. Plasmas*, 23:032104, 2016.

[11] E. Hassan, D.R. Hatch, W. Guttenfelder, Y. Chen, and S. Parker. Gyrokinetic benchmark of the electron temperature-gradient instability in the pedestal region. *Phys. Plasmas*, 28:062505, 2021.

[12] F. Imbeaux, S.D. Pinches, J.B. Lister et al. Design and first applications of the iter integrated modelling and analysis suite. *Nucl. Fusion*, 55:123006, 2015.

[13] G.M. Staebler, J.E. Kinsey, and R. E. Waltz. A theory-based transport model with comprehensive physics. *Phys. Plasmas*, 14:055909, 1995.

[14] O. Sauter and S.Yu. Medvedev. Tokamak coordinate conventions: Cocos. *Comput. Phys. Commun.*, 184:293, 2013.

[15] F.L. Hinton and S.K. Wong. Neoclassical ion transport in rotating axisymmetric plasmas. *Phys. Fluids*, 28:3082, 1985.

[16] A.G. Peeters, D. Strintzi, Y. Camenen, C. Angioni, F. J. Casson, W.A. Hornsby, and A.P. Snodin. Influence of the centrifugal force and parallel dynamics on the toroidal momentum transport due to small scale turbulence in a tokamak. *Phys. Plasmas*, 16:042310, 2009.

[17] A.G. Peeters, D. Strintzi, Y. Camenen, C. Angioni, F. J. Casson, W.A. Hornsby, and A.P. Snodin. Erratum: "influence of the centrifugal force and parallel dynamics on the toroidal momentum transport due to small scale turbulence in a tokamak" [phys. plasmas16, 042310 (2009)]. *Phys. Plasmas*, 19:099901, 2012.

[18] A.J. Brizard. Nonlinear gyrokinetic vlasov equation for toroidally rotating axisymmetric tokamaks. *Phys. Plasmas*, 2:459, 1995.

[19] E.A. Belli and J. Candy. Impact of centrifugal drifts on ion turbulent transport. *Phys. Plasmas*, 25:032301, 2018.

[20] H. Sugama and W. Horton. Nonlinear electromagnetic gyrokinetic equation for plasmas with large mean flows. *Phys. Plasmas*, 5:2560, 1998.

[21] https://github.com/iterorganization/IMAS-Data-Dictionary.

[22] https://imas-data-dictionary.readthedocs.io/en/latest/.

[23] https://physics.nist.gov/cuu/Constants/.

[24] R.E. Waltz and R.L. Miller. Ion temperature gradient turbulence simulations and plasma flux surface shape. *Phys. Plasmas*, 6:4265, 1999.

[25] R. Arbon, J. Candy, and E.A. Belli. Rapidly-convergent flux-surface shape parameterization. *Plasma Phys. Control. Fusion*, 63:012001, 2021.

[26] https://gitlab.com/gkdb/imas-gk.

[27] A. Najlaoui, Y. Camenen, C. Bourdelle, and A. Loarte. Verifying turbulence model reduction in high beta tokamak plasmas. *Plasma Phys. Control. Fusion*, 67:045016, 2025.

[28] https://gitlab.com/gkdb/gkdb.