# Electromagnetic version of the flux-driven full-F global gyrokinetic code GYSELA

Zhisong Qu, Youngwoo Cho, Chang Kai Chai, Michael Go, Roméo Bigué, Peter Donnel, Kévin Obrejan, Yi Wen Ban, Ludovica de Gianni, G Dif-Pradalier, et al.

► To cite this version:

Zhisong Qu, Youngwoo Cho, Chang Kai Chai, Michael Go, Roméo Bigué, et al. Electromagnetic version of the flux-driven full-F global gyrokinetic code GYSELA. 2025. hal-05342386

HAL Id: hal-05342386

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

Preprint submitted on 2 Nov 2025

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 Attribution 4.0 International License
Electromagnetic version of the flux-driven full-$F$ global gyrokinetic code GYSELA

Z. S. Quª, Y. W. Choª, C. K. Chaiª, M. Goª, R. Biguéb,c, P. Donnelb, K. Obrejanb, Y. W. Banª, L. De Giannib, G. Dif-Pradalierb, X. Garbeta,b, V. Grandgirardb, K. Limª, J. Y. Ngª, S. Rajª, Y. Sarazinb, E. Sonnendrückerd,e, R. Varennesa, C. G. Wanª

ªSchool of Physical and Mathematical Sciences, Nanyang Technological University, Singapore, 637371, Singapore

bCEA, IRFM, Saint-Paul-lez-Durance, F-13108, France

cAix-Marseille University, CNRS, PIIM, Marseille, France

dMax Planck Institute for Plasma Physics, Garching, D-85748, Germany

eSchool of Computation Information and Technology, Technical University of Munich, Garching, D-85748, Germany

**Abstract**

We present the electromagnetic (EM) extension of the flux-driven full-$F$ global gyrokinetic code GYSELA (for GYrokinetic SEMi-LAgrangian), enabling self-consistent simulations of turbulence and transport in magnetised plasmas beyond the electrostatic limit. A mixed-variable pullback scheme has been implemented and enhanced with a predictor for the non-ideal component of the parallel vector potential, reducing numerical stiffness and preserving conservation properties. The resulting formulation allows stable timestepping with realistic parameters. Verification is performed against theoretical shear Alfvén wave dispersion and linear benchmarks in the Cyclone Base Case, showing excellent agreement. Furthermore, nonlinear flux-driven electromagnetic ion-temperature-gradient turbulence simulations demonstrate

Email address: zhisong.qu@ntu.edu.sg (Z. S. Qu)

Preprint submitted to Elsevier

November 2, 2025
the robustness of the scheme, with conservation of particles, energy, and
momentum below 1% error. These developments establish the EM version
of GYSELA as a reliable tool for investigating multi-scale electromagnetic
plasma dynamics and transport in fusion-relevant regimes.

Keywords: Tokamak, Gyrokinetics, Turbulence, Semi-Lagrangian

# 1. Introduction

Electromagnetic effects are becoming increasingly significant in turbu-
lent transport within fusion plasmas, particularly as current and forthcom-
ing experiments operate at high plasma beta (β). In addition, energetic
particle-driven instabilities such as Alfvén eigenmodes interact strongly with
turbulence, adding considerable complexity to the system dynamics [1, 2]. A
comprehensive understanding therefore requires a self-consistent treatment
of both the particle distribution function and the electromagnetic fields. Gy-
rokinetic theory [3, 4] serves as the principal framework for investigating
electromagnetic turbulence and its interactions with energetic particles.

A number of gyrokinetic codes, especially those capable of simulating the
global plasma domain, have already incorporated electromagnetic effects, in-
cluding EUTERPE [5], GENE [6], GENE-X [7], GENE-3D [8], GKNET [9],
gKPSP [10], GKW [11], GTC [12], NLT [13], ORB5 [14], TRIMEG-GKX [15],
and XGC [16], among others. The GYSELA (GYrokinetic SEmi-LAgrangian)
code [17] is a global full-F nonlinear gyrokinetic code that has historically
been devoted to simulating electrostatic plasma turbulence and has achieved
numerous successes, particularly in the study of flux-driven physics and non-
local transport [18, 19]. A preliminary implementation of an electromagnetic
scheme was previously proposed by Gillot [20]. In the present work, we rigorously derive the governing equations, implement a numerically robust scheme with second-order accuracy in time, and carry out detailed benchmarking tests, thereby establishing GYSELA as a practical tool for studying electromagnetic turbulence.

At the core of a successful electromagnetic gyrokinetic code lies the numerical scheme used to evolve the perturbed parallel vector potential $A_\|$ and the treatment of kinetic electrons. The former is closely associated with the so-called cancellation problem in the Hamiltonian formalism of gyrokinetics [21], which can lead to unphysical results at high plasma beta, while the latter imposes a stringent constraint on the time step size due to the rapid streaming of electrons along the magnetic field line [22]. In recent years, the mixed-variable pull-back scheme [23, 24, 25, 26] has emerged as a powerful approach, enabling electromagnetic simulations that are free from the cancellation problem. The parallel vector potential is rewritten into the sum of a “symplectic part” $A_\|S$, and a “Hamiltonian part” $A_\|H$, while keeping the latter small. A further advantage of this scheme is the removal of the parallel gradient of the electrostatic potential from the advection equation of the parallel momentum, which significantly reduces numerical stiffness and allows for the use of much larger time steps without compromising accuracy.

It turns out that although GYSELA does not suffer from the cancellation problem as severely as particle-in-cell (PIC) codes, owing to its grid-based formulation, the presence of a cross term between the parallel momentum and the vector potential compromises the conservation properties of the Strang splitting scheme employed in the code. It is therefore essential to minimise
$A_{\parallel H}$ as much as possible, making the mixed-variable pull-back scheme a natural and effective choice. In this work, we introduce a predictor for $A_{\parallel S}$ based on its time history, which further reduces the amplitude of $A_{\parallel H}$. In addition, we have carefully designed a numerical scheme with enhanced stability, enabling simulations with realistic time steps. In GYSELA, the rapid electron streaming is treated using a semi-Lagrangian scheme, which is known to possess good stability properties [27]. With these developments, GYSELA now demonstrates not only consistent results in numerical verifications and physical benchmarks, but also the capability to perform flux-driven nonlinear simulations with robust conservation properties.

This paper is organised as follows. The gyrokinetic Vlasov equation and Poisson-Ampère equations are introduced in Section 2.1 and 2.2, respectively, with their derivations detailed in Appendix B. Section 2.3 briefly gives the conservation laws to be tested numerically. Section 3 describes the numerical method. Various physical and numerical benchmarks are performed in Section 4. Finally, Section 5 draws the conclusion.

## 2. Equations

### 2.1. The gyrokinetic Vlasov equation

The gyrokinetic Vlasov equation is given by

$$ \partial_t(B_{\parallel s}^* \bar{F}_s) + \nabla \cdot (B_{\parallel s}^* \dot{\mathbf{X}} \bar{F}_s) + \partial_{p_{\parallel}} (B_{\parallel s}^* \dot{p}_{\parallel} \bar{F}_s) = B_{\parallel s}^* \mathcal{R}(\bar{F}_s), \quad (1) $$

where $\bar{F}_s(\mathbf{X}, p_{\parallel}, \mu, t)$ is the gyrocentre distribution function, $B_{\parallel s}^*$ the velocity space Jacobian, $\mathbf{X}$ the gyrocentre position, $p_{\parallel}$ the canonical velocity, $\mu$ the magnetic moment, $t$ the time variable, and $\mathcal{R}$ the right-hand side operator
containing the sources, Krook operator, collisions, etc. The forms of these right-hand side operators can be found in [17]. The overhead dot represents the total time derivatives along the gyrocentre trajectory, i.e. $d/dt$. Strictly speaking, an overhead bar should be added to all the gyrocentre variables above to distinguish them from guiding centre ones (e.g. $\mu \to \bar{\mu}$), but we have ignored the bar for simplicity. The species s is labelled by the subscript “s”. The definition of $p_\parallel$ is given by

$$p_\parallel = v_\parallel + \frac{q_s}{m_s} (A_{\parallel H} + A_{\parallel S} - \mathcal{J}[A_{\parallel S}]), \quad (2)$$

where $v_\parallel$ is the parallel velocity, $A_{\parallel S}$ and $A_{\parallel H}$ the symplectic and Hamiltonian parts of the parallel vector potential, respectively, $m_s$ and $q_s$ the mass and charge, respectively, and $\mathcal{J}[\cdots]$ the gyro-average operator. For large-scale modes, $p_\parallel = v_\parallel + q_s A_{\parallel H}/m_s$, and therefore $A_{\parallel H}$ is the part of the canonical momentum and the Hamiltonian as well, while $A_{\parallel S}$ appears only in the symplectic part of the Lagrangian $A_s^* \cdot \dot{\mathbf{X}}$. This definition follows Brizard and Hahm [3] Eq. (169) and differs from the definition in some other gyrokinetic codes (e.g. ORB5) in which $p_\parallel = v_\parallel + \frac{q_s}{m_s} \mathcal{J}[A_{\parallel H}]$. The equations of motion are the same when the first-order Hamiltonian is used. However, the second-order Hamiltonian differs, leading to a slightly different form of the Ampère equation in Section 2.2.

We define

$$A_s^* = A + \left( \frac{m_s}{q_s} p_\parallel + \mathcal{J}[A_{\parallel S}] \right) \mathbf{b}, \quad (3)$$

where $A$ is the equilibrium magnetic vector potential, $\mathbf{B} = \nabla \times A$ the equilibrium magnetic field, and $\mathbf{b} = \mathbf{B}/|\mathbf{B}|$ the unit vector along $\mathbf{B}$. Taking the
curl of (3), one gets

$$ \mathbf{B}_s^* = \nabla \times \mathbf{A}_s^* = B_{\parallel s}^* \mathbf{b}_s^* - \left( \frac{m_s}{q_s} p_{\parallel} + \mathcal{J}[A_{\parallel s}] \right) \frac{\nabla B \times \mathbf{B}}{B^2} + \nabla \mathcal{J}[A_{\parallel s}] \times \mathbf{b}, \quad (4) $$

where $B = |\mathbf{B}|$ and

$$ B_{\parallel s}^* = \mathbf{B}_s^* \cdot \mathbf{b} = B + \left( \frac{m_s}{q_s} p_{\parallel} + \mathcal{J}[A_{\parallel s}] \right) \frac{\mathbf{b} \cdot \nabla \times \mathbf{B}}{B}, \quad (5) $$

while

$$ \mathbf{b}_s^* = \frac{\mathbf{B}}{B_{\parallel s}^*} + \left( \frac{m_s}{q_s} p_{\parallel} + \mathcal{J}[A_{\parallel s}] \right) \frac{\nabla \times \mathbf{B}}{B_{\parallel s}^* B}. \quad (6) $$

Note that $\mathbf{b}_s^* \neq \mathbf{B}_s^* / B_{\parallel s}$.

The gyrocentre Hamiltonian [3] up to the first order of gyroradius is written as

$$ \bar{H}(\mathbf{X}, p_{\parallel}, \mu) = \bar{H}_0 + \bar{H}_1, \quad (7) $$

where the zeroth- and first-order components are given by

$$ \bar{H}_0 = \frac{m_s}{2} p_{\parallel}^2 + \mu B, \quad (8) $$

$$ \bar{H}_1 = q_s (\mathcal{J}[\Phi] - p_{\parallel} \mathcal{J}[A_{\parallel H}]), \quad (9) $$

respectively, where $\Phi$ is the electrostatic potential. The equations of motion are then computed from the gyro-centre Poisson brackets [3], given by

$$ B_{\parallel s}^* \dot{\mathbf{X}} = \frac{\mathbf{B}_s^*}{m_s} \partial_{p_{\parallel}} \bar{H} - \frac{\nabla \bar{H} \times \mathbf{B}}{q_s B}, \quad (10) $$

$$ m_s B_{\parallel s}^* \dot{p}_{\parallel} = -\mathbf{B}_s^* \cdot \nabla \bar{H} - q_s B_{\parallel s}^* \partial_t \mathcal{J}[A_{\parallel s}]. \quad (11) $$

After some algebra, one gets

$$ \dot{\mathbf{X}} = \left( p_{\parallel} - \frac{q_s}{m_s} \mathcal{J}[A_{\parallel H}] \right) \mathbf{b}_s^* + \mathbf{v}_E + \mathbf{v}_D + \mathbf{v}_{AS} + \mathbf{v}_{AH}, \quad (12) $$
where

$$ \mathbf{v}_E = \frac{\mathbf{B} \times \nabla \mathcal{J}[\Phi]}{B B_{\parallel s}^*}, \qquad (13) $$

$$ \mathbf{v}_D = \left[ m_s \left( p_{\parallel} + \frac{q_s}{m_s} \mathcal{J}[A_{\parallel s}] \right) \left( p_{\parallel} - \frac{q_s}{m_s} \mathcal{J}[A_{\parallel H}] \right) + \mu B \right] \times \frac{\mathbf{B} \times \nabla B}{q_s B^2 B_{\parallel s}^*}, \qquad (14) $$

$$ \mathbf{v}_{AS} = - \left( p_{\parallel} - \frac{q_s}{m_s} \mathcal{J}[A_{\parallel H}] \right) \frac{\mathbf{B} \times \nabla \mathcal{J}[A_{\parallel s}]}{B B_{\parallel s}^*}, \qquad (15) $$

$$ \mathbf{v}_{AH} = -p_{\parallel} \frac{\mathbf{B} \times \nabla \mathcal{J}[A_{\parallel H}]}{B B_{\parallel s}^*}, \qquad (16) $$

standing for the $E \times B$ velocity, magnetic drift velocity (predominantly), magnetic fluttering velocity due to $A_{\parallel s}$ and $A_{\parallel H}$, respectively. For $p_{\parallel}$, we get

$$ m_s \dot{p}_{\parallel} = -q_s \partial_t \mathcal{J}[A_{\parallel s}] - \left( \mathbf{b}_s^* - \frac{\mathbf{B} \times \nabla \mathcal{J}[A_{\parallel s}]}{B B_{\parallel s}^*} \right) \cdot \left( \mu \nabla B + q_s \nabla \mathcal{J}[\Phi] - q_s p_{\parallel} \nabla \mathcal{J}[A_{\parallel H}] \right) + \frac{m_s}{B} \left( p_{\parallel} + \frac{q_s}{m_s} \mathcal{J}[A_{\parallel s}] \right) (\mathbf{v}_E + \mathbf{v}_{AH}) \cdot \nabla B. \qquad (17) $$

The actual equations implemented numerically are written in normalised units. In GYSELA, we can define a reference species with mass $m_{\text{ref}} = m_p$ (the proton mass) and charge $q_{\text{ref}} = e$ (the elementary charge), temperature $T_{\text{ref}}$, and number density $n_{\text{ref}}$. The latter two quantities are usually chosen to be those of the main ion at mid-radius. On top of these, we also need to define a reference magnetic strength $B_0$, usually taken to be the field strength on the magnetic axis. Let $\omega_{ci} = q_{\text{ref}}B_0/m_{\text{ref}}$, $v_{\text{Tref}} = c_s = \sqrt{T_{\text{ref}}/m_{\text{ref}}}$, $\rho_{\text{ref}} = m_{\text{ref}}v_{\text{Tref}}/(q_{\text{ref}}B_0)$, $\beta_{\text{GYS}} = \mu_0 n_{\text{ref}} T_{\text{ref}}/B_0^2$. Other normalised quantities are given in Table A.2 in Appendix A. As an exception, we define
$\hat{p}_{\parallel} = \sqrt{m_s/m_{\text{ref}}}p_{\parallel}/v_{\text{Tref}}$, which is normalised to the thermal velocity species-wise, with $m_s$ the mass for species "s". The property of the system is predominantly determined by three dimensionless quantities that need to be prescribed: $\rho_{*} = \rho_{\text{ref}}/a$, where $a$ is the minor radius, $\beta_{\text{GYS}}$, and $v^{*}$, the neoclassical collisional frequency. The normalised equations are given in Appendix A.

## 2.2. The Poisson-Ampère equations

The gyrokinetic Poisson equation in the long-wavelength approximation (LWA) is given by

$$ -\sum_s \nabla \cdot \left( \frac{m_s n_{s0}}{B^2} \nabla_{\perp} \Phi \right) = \sum_s q_s n_s, \qquad (18) $$

where $\nabla_{\perp} = \nabla - \mathbf{b}\mathbf{b} \cdot \nabla \approx \nabla - \hat{\mathbf{e}}_{\varphi} \hat{\mathbf{e}}_{\varphi} \cdot \nabla$, where $\hat{\mathbf{e}}_{\varphi}$ is the unit vector in $\varphi$ direction. The particle density can be computed by

$$ n_s = 2\pi \int d\mu dp_{\parallel} \mathcal{J}^T[B_{\parallel s}^* \bar{F}_s], \qquad (19) $$

and $\mathcal{J}^T[\dots]$ the adjoint gyro-average operator. The density being used in the polarisation term, $n_{s0}$, can be approximated by the initial density, or updated after a number of time steps. For the Fourier solver in the original GYSELA literature [17], $B^2$ in the polarisation term is approximated by $B_0^2$, the reference field strength at the centre of the tokamak. This approximation is relaxed in the spline solver [28].

For smaller-scale modes, the Padé approximation can be used, with the Poisson equation given by

$$ -\sum_s \nabla \cdot \left( \frac{m_s n_{s0}}{B^2} \nabla_{\perp} \Phi \right) = [1 - (\nabla \cdot \rho_i^2 \nabla_{\perp})] \sum_s q_s n_s, \qquad (20) $$
where $\rho_i = \sqrt{m_i T_i}/(q_i B)$ is the Larmor radius of the main ion species.

The parallel Ampère equation, after taking the LWA approximation and noting $k_\perp a \gg 1$, can be written as

$$0 = \nabla \cdot (1 - \beta_{\text{tot}}) \nabla_{\perp} A_{\parallel H} - \sum_{s} \frac{\mu_0 q_s^2 n_{0,s}}{m_s} A_{\parallel H} \\ + \nabla_{\perp}^2 A_{\parallel S} + \mu_0 \sum_{s} \delta J_{\parallel s}, \qquad (21)$$

where $\mu_0$ is the vacuum magnetic permeability, and

$$\beta_{\text{tot}} = \mu_0 \sum_{s} \frac{n_{s0} T_{s0}}{B^2}. \qquad (22)$$

The perturbed current density (including the skin current) is defined as

$$\sum_{s} \delta \bar{J}_{\parallel, s} = 2\pi \sum_{s} q_s \int d\mu d p_{\parallel} \mathcal{J}^{T}[p_{\parallel} B_{\parallel s}^{*} \bar{F}_{s}] - J_{\parallel, 0}, \qquad (23)$$

with $J_{\parallel,0} = \mu_0^{-1} \mathbf{b} \cdot \nabla \times \mathbf{B}$ the equilibrium parallel current. One can initialise the electron to have a shifted-Maxwellian distribution function to match the equilibrium density $J_{\parallel,0}$, such that $\sum_s \delta \bar{J}_{\parallel, s} = 0$ at equilibrium. This is important for current-driven modes. For simulations without current-driven instabilities, e.g. ITG, we initialise the electrons as an un-shifted Maxwellian, meanwhile removing $J_{\parallel,0}$ from the right-hand side of (23).

We notice two main differences between (21) and the formula being used in ORB5 [14]. First, the skin current term proportional to $A_{\parallel H}$ is not double gyro-averaged: this difference is due to the different definition of $p_{\parallel}$ in (2). Second, we also included the $\beta_{\text{tot}} \nabla_{\perp} A_{\parallel H}$ coming from the polarisation current. One may refer to Appendix B for a derivation of (21) using the second-order Hamiltonian. Coincidentally, (21) is the same as that of ORB5
in the LWA limit, after expanding the double gyro-average operator in the ORB5 skin-current term.

In the low-beta limit, the $\beta_{tot}$ term is dropped, and the Ampère equation returns to its familiar form given by

$$ -\nabla_{\perp}^2 A_{\parallel H} + \sum_s \frac{\mu_0 q_s^2 n_{0,s}}{m_s} A_{\parallel H} = \mu_0 \sum_s \delta J_{\parallel s} + \nabla_{\perp}^2 A_{\parallel S}, \quad (24) $$

in which $A_{\parallel H}$ is solved given the right-hand side.

Up to now, we have not given the details of how $A_{\parallel S}$ evolves in time. This will be introduced in Section 3.1.

## 2.3. Local conservation laws

In this section, we derive flux-surface averaged conservation laws for which GYSELA will be tested against. The conservation laws are computed without $A_{\parallel H}$, since $A_{\parallel H}$ can be combined into $A_{\parallel S}$ and set to zero after the pullback transformation introduced in Section 3.1.

Let $\lambda(\mathbf{X}, p_{\parallel}, \mu, t)$ be the quantity one wants to measure, from the Vlasov equation (1), one gets that

$$ \begin{aligned} & \partial_t (B_{\parallel s}^* \lambda \bar{F}_s) + \nabla \cdot (B_{\parallel s}^* \dot{\mathbf{X}} \lambda \bar{F}_s) + \partial_{p_{\parallel}} (B_{\parallel s}^* \dot{p}_{\parallel} \lambda \bar{F}_s) \\ & = B_{\parallel s}^* \bar{F}_s \dot{\lambda} + B_{\parallel s}^* \lambda \mathcal{R}(\bar{F}_s). \end{aligned} \quad (25) $$

Suppose $\chi$ is a flux label (e.g. $\chi = r$ in circular geometry), integrating (25) over $\mu, p_{\parallel}, \theta, \varphi$ gives

$$ \partial_t \Lambda + \partial_\chi \Gamma_\Lambda = S_\Lambda + S_{\Lambda, \mathcal{R}}, \quad (26) $$

where $\Lambda = \langle \lambda \rangle$, $\Gamma_\Lambda = \langle \lambda \dot{X}^\chi \rangle$, $S_\Lambda = \langle \dot{\lambda} \rangle$, $S_{\lambda, \mathcal{R}} = \langle \lambda \mathcal{R} \rangle$, $\dot{X}^\chi$ the contravariant component of $\dot{\mathbf{X}}$ in $\chi$ direction, and

$$ \langle \cdots \rangle = 2\pi \int \cdots J B_{\parallel s}^* \bar{F}_s d\theta d\varphi dp_{\parallel} d\mu, \quad (27) $$
integrating the phase space and averaging over a flux surface, with $J$ being the Jacobian of the spatial coordinates.

### 2.3.1. Particle conservation

Inserting $\lambda = 1$ into (26) gives the particle balance equation, written as

$$ \partial_t N + \partial_\chi \Gamma_N = S_{N,\mathcal{R}}, \qquad (28) $$

with $N = \langle 1 \rangle$ being the gyrocentre density on a flux surface, and the particle radial flux given by

$$ \Gamma_N = \langle v_D^\chi + v_E^\chi + v_{AS}^\chi \rangle, \qquad (29) $$

with their physical meaning being magnetic drift fluxes (14), the $E \times B$ fluxes (13) and the magnetic fluttering fluxes (14), respectively. The last term $S_{N,\mathcal{R}}$ contains the contribution of particle sources, Krook operator, collisions, etc. In the absence of a particle source and Krook operator, one should have $S_{N,\mathcal{R}} = 0$ since the collision operator conserves the number of particles.

### 2.3.2. Energy conservation

Let $\lambda$ be the parallel energy $E_\parallel = \frac{1}{2}m_s p_\parallel^2$ and the perpendicular energy $E_\perp = \mu B$, respectively gives

$$ \partial_t \mathcal{E} + \partial_\chi (\Gamma_{\mathcal{E}\parallel} + \Gamma_{\mathcal{E}\perp}) = S_{\mathcal{E}} + S_{\mathcal{E},\mathcal{R}} \qquad (30) $$

with $\mathcal{E} = \langle E_\parallel + E_\perp \rangle$, $S_\mathcal{E} = S_{\mathcal{E}\parallel} + S_{\mathcal{E}\perp}$ and

$$ \Gamma_{\mathcal{E}\parallel/\perp} = \langle E_{\parallel/\perp} (v_D^\chi + v_E^\chi + v_{AS}^\chi) \rangle. \qquad (31) $$

Summing over $S_{\mathcal{E}\parallel} = \langle m_s \dot{p}_\parallel p_\parallel \rangle$ and $S_{\mathcal{E}\perp} = \langle \mu \dot{\mathbf{X}} \cdot \nabla B \rangle$, a lot of the terms cancel with each other. The remaining terms are

$$ S_{\mathcal{E}} = \langle q_s p_\parallel (-\partial_t \mathcal{J}[A_\parallel S] - \mathbf{b}_s^* \cdot \nabla \Phi) - q_s (\mathbf{v}_D + \mathbf{v}_{AS}) \cdot \nabla \Phi \rangle. \qquad (32) $$
The physical meaning of the first and second terms are just $J_\parallel E_\parallel$ and $J_\perp$. $E_\perp$ for each species, standing for the energy exchange between waves and particles in the parallel and perpendicular directions, respectively. The last term $S_{\epsilon, \mathcal{R}}$ contains the contribution of heating, Krook operator, collisions (energy transfer between species), etc.

### 2.3.3. Toroidal angular momentum

The toroidal canonical momentum is given by

$$P_{\varphi} = -q_s \psi_p + m_s p_{\parallel} b_{\varphi} + q_s \mathcal{J}[A_{\parallel S}] b_{\varphi}, \quad (33)$$

where $\psi_p$ is the poloidal flux and $b_\varphi$ the co-variant component of $b$ in the $\varphi$ direction. It can be shown by substituting in the equation of motion that

$$\frac{dP_{\varphi}}{dt} = -\partial_{\varphi}H + q_s p_{\parallel} \partial_{\varphi} \mathcal{J}[A_{\parallel S}] \quad (34)$$

Following Abiteboul [29], we define the momentum, Reynolds stress, and torque respectively as

$$\mathcal{L}_{\varphi} = \langle m_s p_{\parallel} b_{\varphi} + q_s \mathcal{J}[A_{\parallel S}] b_{\varphi} \rangle, \quad (35)$$

$$\Pi_{\varphi}^{\chi} = \langle (m_s p_{\parallel} b_{\varphi} + q_s \mathcal{J}[A_{\parallel S}] b_{\varphi})(v_D^{\chi} + v_E^{\chi} + v_{AS}^{\chi}) \rangle, \quad (36)$$

$$\mathcal{T}_{\varphi} = \langle \partial_{\varphi} H - q_s p_{\parallel} \partial_{\varphi} \mathcal{J}[A_{\parallel S}] \rangle. \quad (37)$$

The conservation equation should be written as

$$\partial_t \langle P_{\varphi} \rangle + \partial_{\chi} \Gamma_{P_{\varphi}} = -\mathcal{T}_{\varphi} + S_{P_{\varphi}, \mathcal{R}}. \quad (38)$$

One can use the relationship $\langle P_{\varphi} \rangle = \mathcal{L}_{\varphi} - q_s \psi_p N$, and $\Gamma_{P_{\varphi}} = \Pi_{\varphi}^{\chi} - q_s \psi_p \Gamma_N$, such that (38) becomes

$$\partial_t \mathcal{L}_{\varphi} + \partial_{\chi} \Pi_{\varphi}^{\chi} + \mathcal{T}_{\varphi} = q_s \Gamma_N \frac{d\psi_p}{d\chi} + S_{P_{\varphi}, \mathcal{R}} + q_s \psi_p S_{N, \mathcal{R}}, \quad (39)$$
where the particle conservation (28) is used to eliminate some of the terms. The term $S_{P,\mathcal{R}}$ contains the contribution of momentum source, Krook operator, collisions (drag between species), etc.

# 3. Numerical Methods

## 3.1. Predictor for the symplectic vector potential

So far $A_{\parallel S}$ is a free variable and the separation of $A_{\parallel}$ into $A_{\parallel S}$ and $A_{\parallel H}$ is arbitrary. A governing equation for $\partial_t A_{\parallel S}$ needs to be introduced to uniquely determine its time evolution. Such an equation must serve two purposes: to remove $\partial_t A_{\parallel S}$ from (17) so the numerical scheme can remain explicit in time, and to ensure that $A_{\parallel H}$ remains small to avoid the cancellation problem. In GYSELA, we mainly follow the choice by Mishchenko [23, 24], i.e. $A_{\parallel S}$ evolves mainly according to ideal MHD, with a new term to also capture the non-ideal contributions.

### 3.1.1. Original mixed-variable scheme

The original mixed-variable scheme evolves $A_{\parallel S}$ according to the ideal Ohm's law (hereby we name it IMHD predictor), given by

$$ \partial_t A_{\parallel S} = -\nabla_{\parallel} \Phi, \qquad (40) $$

while $A_{\parallel H}$ is solved by the Ampère equation (24). In GYSELA, the parallel gradient is computed by cubic splines in poloidal and toroidal directions. Along with the mixed-variable scheme, a pull-back transformation at the end of each time step is introduced, namely

$$ A_{\parallel S, \text{new}} = A_{\parallel S, \text{old}} + A_{\parallel H, \text{old}}, \quad A_{\parallel H, \text{new}} = 0, \qquad (41) $$
and a change of variable of $p_\parallel$

$$p_{\parallel, \text{new}} = p_{\parallel, \text{old}} - \frac{q_s}{m_s} \mathcal{J}[A_{\parallel H, \text{old}}], \qquad (42)$$

where the subscripts “old” and ”new” correspond to the quantities before and after the transformation, respectively. This gives the pull-back transform of $\bar{F}_s$, written as

$$\bar{F}_s(p_{\parallel, \text{new}}) = \bar{F}_s(p_{\parallel, \text{old}} - \frac{q_s}{m_s} \mathcal{J}[A_{\parallel H, \text{old}}]). \qquad (43)$$

### 3.1.2. Mixed-variable with non-ideal part

While the IMHD predictor keeps $A_{\parallel H}$ small for Alfvén waves, it has a draw-back for nearly-electrostatic modes (e.g. drift waves, GAM, etc). Indeed, these fields are almost static in Alfvén time scale. Solving (40) with a constant $\nabla_\parallel \Phi$ would leads to $A_{\parallel S} \sim -\Delta t \nabla_\parallel \Phi$, even though $A_\parallel$ should be close to zero. The consequence is that, within each time step, $A_{\parallel H}$ would need to compensate the accumulated $A_{\parallel S}$ to keep the total $A_\parallel$ nearly zero, giving $A_{\parallel H} \sim \Delta t \nabla_\parallel \Phi$.

To further reduce $A_{\parallel H}$ and therefore to preserve the symplecticity of the Strang splitting, we introduce an additional term in (40) to take into account the non-ideal contribution, written as

$$\partial_t A_{\parallel S} = -\nabla_\parallel \Phi + \partial_t A_{\parallel N}, \qquad (44)$$

where $A_{\parallel N}$ now takes the role of compensating the term $-\nabla_\parallel \Phi$ in the electrostatic limit, to keep $A_{\parallel S}$ near zero and ease the burden of $A_{\parallel H}$. In principle, one can derive the non-ideal Ohm’s law and a non-ideal predictor [30]. In this work, we instead extrapolate $A_{\parallel N}$ from the time history. We can choose $A_{\parallel N}$ to be a linear function of time within a time step, that is,

$$A_{\parallel N} = \Delta A_{\parallel N} t, \qquad (45)$$
such that $\partial_t A_{\parallel N} = \Delta A_{\parallel N}$. After each time step, we can adjust $\Delta A_{\parallel N}$ accordingly to

$$ \Delta A_{\parallel N, \text{new}} = \Delta A_{\parallel N, \text{old}} + \frac{A_{\parallel H, \text{old}}}{\Delta t}. \qquad (46) $$

This implies that the slope of $A_{\parallel N}$ now becomes the slope of the “true” non-ideal part of the past time step including the contribution from $A_{\parallel H}$. Now $A_{\parallel N}$ perfectly covers the non-ideal contribution should it be a linear function of time (e.g. when $E$ is electrostatic and $\Phi$ is changing slowly). A deviation from being a linear function would mean a quadratic dependency on time, which should now be captured by $A_{\parallel H}$ by solving the Ampère equation. Therefore, by introducing $A_{\parallel N}$, the size of $A_{\parallel H}$ reduces from $O(\Delta t)$ to $O(\Delta t^2)$, a much smaller quantity compared to the original mixed-variable scheme.

At $t = 0$, $A_{\parallel S}$ is solved using the symplectic version of the Ampère equation, i.e. solving (24) for $A_{\parallel S}$ assuming $A_{\parallel H} = 0$. The non-ideal part $A_{\parallel N}$ is also set to zero.

## 3.2. Fourier-based solver for the Ampère equation

GYSELA can currently support both circular geometry with no Shafranov shift, or the more realistic Culham geometry with shaping. In this work, we only report results using the circular geometry.

The Ampère equation is solved in the Fourier space of the two angles. In circular geometry, for each toroidal mode $n$, the operator $\nabla_{\perp}^2$ can be written as

$$ \nabla_{\perp}^2 \approx \partial_r^2 + \frac{1}{r} \partial_r - \frac{m^2}{r^2}, \qquad (47) $$

where $m$ is the poloidal mode number. The radial derivative is discretised
by second-order finite differencing. Note that an Ampère solver based on B-splines [31, 28] that can take arbitrary 2D flux surfaces is under development.

### 3.3. Advection of the distribution function

With both $A_{\parallel S}$ and $A_{\parallel H}$ being non-zero, the following relationship still holds:

$$ \partial_t (B_{\parallel s}^*) + \nabla \cdot (B_{\parallel s}^* \dot{\mathbf{X}}) + \partial_{p_{\parallel}} (B_{\parallel s}^* \dot{p}_{\parallel}) = 0, \quad (48) $$

which allows us to rewrite (1) in the advective form as

$$ \frac{D \bar{F}_s}{Dt} = \partial_t \bar{F}_s + \dot{\mathbf{X}} \cdot \nabla \bar{F}_s + \dot{p}_{\parallel} \partial_{p_{\parallel}} \bar{F}_s = \mathcal{R}(\bar{F}_s), \quad (49) $$

i.e. $\bar{F}_s$ is a constant of motion along gyro-centre trajectories without source and collisions. This forms the basic principle of the semi-Lagrangian scheme in GYSELA. The distribution function $\bar{F}_s$ is discretised on a 5D grid. The trajectories of the grid points are solved backward in time to find their foot at the previous time step. Since $\bar{F}_s$ is a constant along the trajectories, its value on the new grid can be computed by cubic spline interpolation at their foot on the old grid. GYSELA uses Strang split, that is, the advection in the 4D space is approximated by a combination of 1D and 2D advections and interpolations. This may break the conservation laws, which will need to be monitored and checked for each simulation.

The advection steps are recapped as follows.

Input: $\bar{F}_s(p_{\parallel}, t)$, $\Delta t$, $\Phi = (\Phi, A_{\parallel S}, \partial_t A_{\parallel N}, A_{\parallel H})$

1. Advection $\Delta t/2$ in $p_{\parallel}$.

2. Advection $\Delta t/2$ in $\varphi$ or along the field line.

3. Advection $\Delta t$ in $(r, \theta)$.
4. Advection $\Delta t/2$ in $\varphi$ or along the field line.

5. Advection $\Delta t/2$ in $p_\parallel$.

Output: $\bar{F}_s(p_\parallel, t + \Delta t)$

One can choose step 2 and 4 to be 1D advection along $\varphi$ (non-aligned advection) [17], or 2D advection along the field line (aligned advection) [32]. The latter interpolates along a field line, where the distribution changes slowly. It is equivalent to considering $\varphi$ as a piece-wise field-aligned coordinate [33], helping to reduce the need for a high grid resolution in $\varphi$. Importantly, all sub-steps in the Strang splitting should be second-order in time, for the entire scheme to be second-order. Steps 2 and 4 were solved with first-order accuracy in the original aligned scheme [32], but are promoted to second-order in the current work. It has been checked that failing to do so will lead to kinetic ballooning modes with an incorrect growth rate.

### 3.4. Time-stepping algorithm

The Vlasov equation (1), the Poisson equation (18) or (20), the Ampère equation (24), and the predictor for $A_\parallel S$ (44) form a closed system of equations and need to be solved self-consistently. EM GYSELA follows the predictor-corrector scheme of the original GYSELA. For the $i$-th time step, we perform the following calculations. The quantities at the beginning of the time step and at the end are denoted by superscript (i) and (i + 1), respectively.

1. Predictor step

(a) Compute $A_\parallel S^{(i+\frac{1}{4})} = A_\parallel S^{(i)} - \frac{\Delta t}{4} \nabla_\parallel \Phi^{(i)} + \frac{\Delta t}{4} \Delta A_\parallel S^{(i)}$.

(b) Advection with time step $\Delta t/2$ from $\bar{F}_s^{(i)}$, using $\Phi = (\Phi^{(i)}, A_\parallel S^{(i+\frac{1}{4})}, \Delta A_\parallel S^{(i)}, 0)$, obtain $\bar{F}_s^{(i+\frac{1}{2})}$.
(c) Solve the Poisson equation with $\bar{F}_s^{(i+\frac{1}{2})}$ to get $\Phi^{(i+\frac{1}{2})}$.

(d) Compute $A_{\parallel S}^{(i+\frac{1}{2})} = A_{\parallel S}^{(i)} - \frac{\Delta t}{4} \nabla_{\parallel} (\Phi^{(i)} + \Phi^{(i+\frac{1}{2})}) + \frac{\Delta t}{2} \Delta A_{\parallel N}^{(i)}$.

(e) Solve the mixed-variable Ampère equation with $\bar{F}_s^{(i+\frac{1}{2})}$ and $A_{\parallel S}^{(i+\frac{1}{2})}$ to get $A_{\parallel H}^{(i+\frac{1}{2})}$.

2. Corrector step

(a) Compute the collision operator $\mathcal{R}$ to update $\bar{F}_s^{(i)}$.

(b) Advection with time step $\Delta t$ from $\bar{F}_s^{(i)}$, using $\bar{\Phi} = (\Phi^{(i+\frac{1}{2})}, A_{\parallel S}^{(i+\frac{1}{2})}, \Delta A_{\parallel N}^{(i)}, A_{\parallel H}^{(i+\frac{1}{2})})$, obtain $\bar{F}_s^{(i+1)}$.

(c) Solve the Poisson equation with $\bar{F}_s^{(i+1)}$ to get $\Phi^{(i+1)}$.

(d) Compute $A_{\parallel S}^{(i+1)} = A_{\parallel S}^{(i)} - \frac{\Delta t}{2} \nabla_{\parallel} (\Phi^{(i)} + \Phi^{(i+1)}) + \Delta t \Delta A_{\parallel N}^{(i)}$.

(e) Solve the mixed-variable Ampère equation with $\bar{F}_s^{(i+1)}$ and $A_{\parallel S}^{(i+1)}$ to get $A_{\parallel H}^{(i+1)}$.

3. Pullback

(a) Accumulate $\Delta A_{\parallel N}^{(i+1)}$ and $A_{\parallel S}^{(i+1)}$ following (41) and (46).

(b) Pullback transformation into $p_{\parallel,\text{new}}$ and reset $A_{\parallel H}$ to zero according to (42) and (43).

4. Compute diagnostics

We note that step 1 (a) (d) and 2 (d) are found to be the most stable and accurate choices from numerical experiments. In fact, we postulate that $A_{\parallel S}$ should be solved with a higher order than $\Phi$ in the predictor step, since its accuracy strongly affects the accuracy of $A_{\parallel H}$, which then affects the stability sensitively.

3.5. CFL conditions

The advection in the Vlasov equation remains stable [27] as long as the trajectories of the nearby grid points do not cross each other, and the
non-symplectic components (e.g. $A_{||H}$) remain small, thanks to the semi-Lagrangian scheme. However, the electromagnetic field are evolved using the explicit predictor-corrector scheme, so CFL conditions still apply. A time advancing scheme with a careful choice of operator splitting might ensure stability [34], which could be a dedicated subject of future research. With the predictor-corrector scheme, while a theoretical proof of the CFL condition is difficult due to the complexity of the equations, we discover via numerical experiments that the maximum time step $\Delta t$ without triggering numerical instabilities is empirically given by

$$k_{||}v_A\Delta t < C_1, \quad k_{||}v_{Te}\Delta t < C_2, \qquad (50)$$

where $C_1$ and $C_2$ are constants close to unity, $v_A = B/\sqrt{\mu_0m_i n_i}$ is the Alfvén velocity and $v_{Te} = \sqrt{T_e/m_e}$ the electron thermal velocity. These conditions are associated with the Alfvén waves and the fast electron parallel motion, respectively. The latter originates from the pressure gradient term in the (hidden) parallel Ohm’s law, instead of the advection itself. Empirically, we found $C_1 \approx 1.5$ and $C_2 \approx 1$. In practice, the maximum $k_{||}$ and therefore the maximum time step is determined by the grid resolution or the filter. We first estimate the time step size from (50) and conduct a short test run without applying seed perturbations to $F$, such that any exponentially growing spurious modes are easy to detect. Then we adjust $\Delta t$ slightly until there are no spurious modes appearing.

### 3.6. Filtering

As mentioned in the previous section, the maximum $\Delta t$ is determined by the CFL conditions (50). In tokamak geometry $k_{||} \approx (m/q-n)/R_0$, where $m$
and $n$ are the poloidal and toroidal mode numbers, respectively. As a global code, GYSELA contains all the $m$ and $n$ modes supported by the grid: $|m| \le N_\theta/2$ and $|n| \le N_\varphi/2$ for non-aligned advection, and $|n - m/q| \le N_\varphi/2$ for aligned advection. To ensure the CFL condition of all modes is satisfied, the time step $\Delta t$ will sometimes be unreasonably small. Another consideration is the violation of the gyrokinetic ordering $k_\parallel\rho_i \ll k_\perp\rho_i \sim 1$ away from the diagonal $m = nq$, especially for modes with $m \approx 0$ and a very high $n$. As a consequence, a diagonal filter is needed in most EM simulations, which keeps only the modes with $|m - nq| \le \Delta m$. The choice of $\Delta m$ should be large enough to capture the physical spectrum. Without filtering, it is still possible to reduce the time step such that no numerical instability can grow, but would usually require $\omega_{ci}\Delta t \ll 1$.

It was found [35] that, to avoid numerical instabilities, one cannot directly apply the diagonal filter to the charge and current density on the right-hand side of the Poisson (18) and Ampère (24) equations, nor can one apply it only to the computed potentials. Rather, the filter should be embedded in the linear operator, such that the new linear equation is given by

$$[CF\mathcal{L}F^{-1}C] \cdot CF\phi = CF\rho, \qquad (51)$$

where $\phi$ is the field to be solved, $\rho$ the charge or current density on the right-hand side, $\mathcal{L}$ the original linear operator such as the Laplacian, $F$ the Fourier transform and $C$ the filter, i.e. a diagonal matrix containing zeros or ones specifying which Fourier components are kept. The equation can be solved in a reduced Fourier space $\hat{\phi} = CF\phi$ with the new operator $\mathcal{L}' = CF\mathcal{L}F^{-1}C$.

The implementation is straightforward for the Fourier-based solver in GYSELA: one can solve the equation directly in Fourier space and apply a
mode-wise Dirichlet boundary condition at certain radii where the mode is cut off due to the changing $q$ profile. The implementation in the GYSELA spline solver is more subtle and is currently under development. We will report the progress in our future publications.

### 3.7. Numerical dissipation

Interpolation in the semi-Lagrangian scheme introduces numerical dissipation. The dissipation rate $\gamma_{num}$ can be large for electrons due to their fast parallel motion. We need $\gamma_{num} \ll \gamma$, the physical growth rate, for a simulation to be physical.

One can compute the numerical dissipation by solving outside GYSELA a reduced parallel transport equation given by

$$ \partial_t f + \frac{v_{\parallel}}{qR_0} \partial_{\theta} f + \frac{v_{\parallel}}{R_0} \partial_{\varphi} f = 0, \qquad (52) $$

using the same scheme as GYSELA for a given $m, n$ mode and fixed $v_{\parallel}$ and $q$, then measure its rate of damping. Empirically, the dissipation rate reduces when the resolution in $\theta$ and $\varphi$ increases, and is the highest when the particle moves 1/3 of a grid for each time step. From numerical experiments, we discovered that for non-aligned advection, $\gamma_{num} \sim n^4$, while for aligned advection, $\gamma_{num} \sim k_{\parallel}^4$ (equivalent to a hyper-diffusion in the parallel direction), but stays low for $k_{\parallel} = 0$. These are considerations one needs to account for when choosing the grid resolution and mode of advection. Failing to do so will lead to a spurious growth rate at high toroidal mode number $n$, since the electrons receive much higher dissipation than the ions, equivalent to a drift-wave drive from finite resistivity.
In general, we recommend the aligned scheme for nonlinear kinetic electron simulation, not only for its low numerical dissipation rate for field-aligned modes, but also for its lower demand in the number of $\phi$ grid.

### 3.8. Initialisation

As a full-$F$ code, GYSELA usually experiences a strong relaxation phase when the initial distribution function is specified as a local Maxwellian. With the electrons being fully kinetic and the vector potential allowed to evolve, the plasma will attempt to relax into a self-consistent Grad-Shafranov equilibrium. Unlike electrostatic simulations ($A_\parallel = 0$) with adiabatic electrons and a fixed density profile, the electron density and temperature profiles are significantly flattened in EM simulations, mostly due to the much higher oscillation frequency of the Alfvén waves and the associated magnetic fluttering during the relaxation phase.

We have introduced a “quiet-start” procedure to step over the initial relaxation phase without causing too much profile relaxation. We choose to initialise the main ion as a canonical Maxwellian [36], where the orbit-averaged poloidal flux is given by

$$ \psi_{p,\text{avg}} = \psi_p - \frac{m_s}{q_s} (Rv_\parallel - R_0 v_{\parallel,\text{avg}}), \quad (53) $$

in which

$$ v_{\parallel,\text{avg}} = \operatorname{sgn}(v_\parallel) \sqrt{\frac{2}{m_s}(E - \mu B_0)H(E - \mu B_{\text{max}})}, \quad (54) $$

with $B_{\text{max}}$ being the maximum field strength on flux surface $\psi_p$ and $H$ the Heaviside function. The density and temperature profiles are specified as a function of $\bar{\psi}_p$. To maintain quasi-neutrality in the initial stage with a poloidally varying ion density profile, the electrons are initialised as a local
Maxwellian, but with the density everywhere matching the ion density. We
then start a simulation with only the $n = 0$ modes enabled. This allows us to
use a very low grid resolution, particularly in the toroidal direction (typically
with $N_\varphi = 8$), and take larger time steps (typically 5 to 10 times the $\Delta t$ used
in full simulations). During the quiet-start stage, collision is enabled, but the
strength of the source is reduced due to the much smaller collisional transport
compared to turbulent transport. After the initial relaxation damped away
and a neoclassical equilibrium is established, we interpolate the distribution
function onto a much finer grid, applying seed perturbations, and start our
turbulence simulation with all $n \neq 0$ modes enabled.

In this work, we have elected not to initialise the equilibrium current in
the electron distribution function and hence not to subtract the equilibrium
current from the right-hand side of (23) (no current-driven instabilities were
tested). In principle, one should also remove the neoclassical current in
nonlinear simulations so it is not double-counted. Starting from an analytical
estimation, one would need to run long neoclassical simulations, wait until
the neoclassical current reaches a steady state, then make sure it is subtracted
self-consistently from the background current after a few iterations. However,
we have not subtracted the neoclassical currents as our $\beta$ is relatively low
and the simulations are short compared to the time needed for neoclassical
currents to fully develop.
# 4. Results

## 4.1. Shear Alfvén waves

In our first test, we will measure the frequency of a shear Alfvén wave and compare it to the theoretical value in large aspect ratio tokamaks, given by

$$ \omega = \left(\frac{m}{q} - n\right) \frac{B_0}{R_0 \sqrt{\mu_0 m_i n_i}} = \left(\frac{m}{q} - n\right) \frac{\epsilon \rho_*}{\sqrt{\beta_{\text{ref}}/2}}, \quad (55) $$

where $m_i = m_{\text{ref}}$ and $n_i = n_{\text{ref}}$ are the mass and density of the main ion, respectively. Here $\beta_{\text{ref}} = 2\mu_0 n_{\text{ref}} T_{\text{ref}} / B_0^2 = 2\beta_{\text{GYS}}$. Motivated by Biancalani et al. [37], we choose $\rho_* = 1/90$, inverse aspect ratio $\epsilon = 0.1$, and the ion/electron mass ratio $m_i/m_e = 200$. All the profiles are flat in the simulation domain. The grid resolution is $(N_r, N_\theta, N_\varphi, N_{p\parallel}, N_\mu) = (64, 32, 32, 32, 32)$. A $m = 1, n = 1$ perturbation with amplitude $\delta f / F_0 = 10^{-6}$ is launched as initial perturbation for ions. We then perform a simulation with $\Delta t = 5\omega_{ci}^{-1}$ for a total of $t = 3000\omega_{ci}^{-1}$. The frequency is then measured from counting the peaks of the $m = n = 1$ mode of $\Phi$ on the surface $r = 0.5a$. We performed two sets of scans. First, we set $\beta_{\text{ref}} = 0.2\%$ and scan over $q$ as shown in Figure 1. Next, we fix $q = 2$ and scan over $\beta_{\text{ref}}$. In both cases, we find a good match between GYSELA and the analytical expression in (55).

Using the physical setup of $q = 2$ and $\beta_{\text{ref}} = 0.2\%$, we test the numerical convergence by decreasing the time step size $\Delta t$ from 8 to $0.25\omega_{ci}^{-1}$, running a total of $t = 40\omega_{ci}^{-1}$. We then measure the relative error, defined by $Err_\Phi = \|\Phi - \Phi_{\text{ref}}\| / \|\Phi_{\text{ref}}\|$, for the value of $\Phi$ on the $r = 0.5a$ surface, where $\|\cdot\|$ stands for the $L^2$ norm. The reference potential $\Phi_{\text{ref}}$ is taken from running
at a smaller time step, $\Delta t = 0.125\omega_{ci}^{-1}$, with the same scheme. The number of grid points is increased to $(N_r, N_\theta, N_\varphi, N_{p\parallel}, N_\mu) = (128, 128, 128, 64, 32)$ so that the error is not bounded by the resolution in space and velocity.

The relationship between Err$_\Phi$ and $\Delta t$ is plotted in Figure 2 (a), for the mixed-variable scheme with and without the non-ideal correction $A_{\parallel N}$. The figure shows that both schemes initially converge second-order in time, as expected by the predictor-corrector scheme and the Strang split. When the time step further decreases, the scheme without $A_{\parallel N}$ drops to a slightly lower order, while the scheme with $A_{\parallel N}$ continues to follow the trend of $O(\Delta t^2)$. Overall, the error of the former is lower than the latter. We note that a further increase in spatial resolution will make a better match to the $O(\Delta t^2)$ line for small $\Delta t$.

We also plot the maximum value of $|A_{\parallel H}|$ on the $r = 0.5a$ surface, normalised by $\max|A_{\parallel S}|$, in Figure 2 (b) as a function of $\Delta t$ for the two schemes. The value of $A_{\parallel H}$ is proportional to $\Delta t$ in the original mixed-variable scheme, while it decreases as $\Delta t^2$ with $A_{\parallel N}$ included and is several orders of magnitude lower. This leads to a better convergence property of the scheme with $A_{\parallel N}$ as shown in Figure 2 (a).

## 4.2. Linear benchmark in the Cyclone-based case

A standard benchmark [38] for electromagnetic gyrokinetic codes is the transition from ITG to KBM as $\beta$ increases using the Cyclone-based case (CBC) [39]. In CBC, the ion and electron density profile is taken to be

$$n_i = n_e = n_0 \exp \left[ -\kappa_n w_n \epsilon \tanh \left( \frac{r - r_0}{w_n a} \right) \right], \quad (56)$$
<figure><img src="image_2.png" /><figcaption>Figure 1: A comparison of the shear Alfvén frequency in GYSELA to theory in (55). Left: a scan over $q$. Right: a scan over $\beta_{\text{ref}}$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: (a) The relative $L^2$ error of $\Phi$ at $r/a = 0.5$ for the shear Alfvén wave test case as a function of time step $\Delta t$, for the mixed-variable scheme without and with the $A_{\parallel N}$ term. The desired order of the predictor-corrector scheme is $O(\Delta t^2)$. (b) The maximum value of $A_{\parallel H}$ at $r/a = 0.5$, normalised to the maximum value of $A_{\parallel S}$, as a function of $\Delta t$.</figcaption></figure>
<table><caption>Table 1: Parameters used in the Cyclone-based case (CBC).</caption><thead><tr><th>ρ*</th><th>ε</th><th><i>R</i>/<i>L</i><sub>n</sub></th><th><i>R</i>/<i>L</i><sub>Ti</sub>, <i>R</i>/<i>L</i><sub>Te</sub></th><th>w<sub>n</sub></th><th>w<sub>Ti</sub>, w<sub>Te</sub></th></tr></thead><tbody><tr><td>1/180</td><td>0.36</td><td>2.23</td><td>6.96</td><td>0.3</td><td>0.3</td></tr><tr><td><i>n</i><sub>0</sub>/<i>n</i><sub>ref</sub></td><td><i>T</i><sub>0,i</sub>/<i>T</i><sub>ref</sub>, <i>T</i><sub>0,e</sub>/<i>T</i><sub>ref</sub></td><td><i>r</i><sub>0</sub>/<i>a</i></td><td><i>m</i><sub>i</sub>/<i>m</i><sub>ref</sub></td><td><i>m</i><sub>e</sub>/<i>m</i><sub>ref</sub></td><td></td></tr><tr><td>1.0</td><td>1.0</td><td>0.5</td><td>1.0</td><td>1/1836.15</td><td></td></tr></tbody></table>

where $r_0$ is the radius with maximum gradient, $n_0$ the density at $r = r_0$, $\kappa_n$ linked to the density gradient length at $r = r_0$ via $R/L_n = \kappa_n$, and $w_n$ the typical width of the profile. Similarly, for temperature,

$$T_{i/e} = T_{0,i/e} \exp \left[ -\kappa_{T,i/e} w_{T,i/e} \epsilon \tanh \left( \frac{r-r_0}{w_{T,i/e}a} \right) \right], \quad (57)$$

with similar definitions of the constants, where the temperature gradient length at $r = r_0$ can be written as $R/L_{T,e/i} = \kappa_{T,e/i}$. The $q$ profile is chosen as

$$q(r) = 0.86 - 0.16 \left(\frac{r}{a}\right) + 2.52 \left(\frac{r}{a}\right)^2. \quad (58)$$

All the other parameters are given in Table 1.

The first test is a comparison of growth rate and frequency for $n = 19$ mode with various $\beta$. We run GYSELA with the LWA solver (18) and a grid resolution of ($N_r, N_\theta, N_\phi, N_{p_\|}, N_\mu$) = (256, 256, 8, 64, 32). The radial simulation domain is between $r = 0.2a$ and $0.8a$, giving a grid interval of $\Delta r = 0.42\rho_i$. A buffer region of $0.05a$ is placed on both sides where a small radial diffusion is applied to remove noise from the boundary. To save computation time, we use field-aligned advection and simulate 1/19 of a torus, hence a small $N_\phi = 8$ is allowed. The domain of simulation in parallel velocity space is between $\pm 6$ mid-radius thermal velocity ($\approx \pm 4$
thermal velocity on axis) for each species, and $0 \le \mu B_0 \le 18T_{\text{ref}}$ for magnetic moment. The result is found to be insensitive to the number of grid points in velocity space. A diagonal filter with $\Delta m = 10$ is used. The time step is chosen to be $\Delta t = 1\omega_{ci}^{-1} = 2 \times 10^{-3}(c_s/R_0)^{-1}$, estimated from the CFL condition of the largest $k_{\parallel} = \Delta m/(R_0 q_{\min})$ mode such that $k_{\parallel}v_{Te}\Delta t \approx 1$. The simulation is initiated with a perturbation in the ion distribution function of all the near-diagonal modes for $n = 19$ on the order of $\delta F/F = 1 \times 10^{-7}$. The simulations last until $t = 10000\omega_{ci}^{-1}$, except for the cases with a small growth rate in which longer runs with $t = 20000\omega_{ci}^{-1}$ were performed. The growth rate and frequency are computed by fitting the amplitude of $\Phi$ at $r = 0.5a, \theta = \varphi = 0$ to an exponentially growing cosine function.

<figure><img src="image_4.png" /><figcaption>Figure 3: Comparison of (a) the linear growth rate and (b) the real frequency of the $n = 19$ mode as a function of $\beta_{\text{ref}}$ in the CBC scenario, for GENE and GYSELA. The transition from ITG to KBM happens at $\beta_{\text{ref}} = 0.015$.</figcaption></figure>

The growth rate and frequency of the $n = 19$ mode are plotted in Figure 3. As $\beta$ increases, ITG is stabilised with the growth rate decreasing. At a certain critical $\beta$, the KBM is destabilised, with the growth rate increasing
<figure><img src="image_5.png" /><figcaption>Figure 4: The GYSELA mode structure of $\Phi$ and $A_\|$ on the cross-section $\varphi = 0$. (a) $\Phi, \beta_{\text{ref}} = 0.5\%$, (b) $A_\|, \beta_{\text{ref}} = 0.5\%$, (c) $\Phi, \beta_{\text{ref}} = 2.5\%$, (d) $A_\|, \beta_{\text{ref}} = 2.5\%$.</figcaption></figure>
again. The transition is also reflected in the discontinuity of the growth rate at the transition point. The GYSELA result shows a very good agreement with the global GENE code [38], including the transition threshold. The field $\Phi$ on the poloidal cross-section is given in Figure 4. showing a clear ballooning structure that agrees well with ORB5 and XGC [40].

The second linear test is performed at a fixed $\beta_{\text{ref}} = 1.01\%$, for a range of toroidal mode numbers $n$. We run GYSELA with the Padé solver (20) and a grid resolution of ($N_r, N_{\theta}, N_{\varphi}, N_{p_{\|}}, N_{\mu}$) = (512, 512, 16, 64, 32) for 1/5 of a torus, using again the field-aligned advection scheme, so $n = 5, 10, \dots, 60$ are allowed in the same simulation. The filter width $\Delta m$ is increased to 23 to incorporate higher $n$ modes, with the time step $\Delta t$ reduced to $0.5\omega_{ci}^{-1} = 10^{-3}(c_s/R_0)^{-1}$ due to the CFL condition, for a total of $20000\omega_{ci}^{-1}$. We then measure the growth rate and frequency of each mode number $n$. When there are multiple competing modes for the same $n$ growing at different radii, we pick the fastest growing one in the domain $0.45 \le r \le 0.55a$, and average the growth rate around its vicinity. The comparison between GYSELA and GENE is given in Figure 5, which shows a typical transition between ITG and the trapped electron mode/electron temperature gradient mode (TEM/ETG) as the growth rate increases, decreases, then increases again with $n$. The frequency jumps at the transition point $n = 35$ from positive (in the direction of ion diamagnetic velocity) to negative (electron diamagnetic velocity). The agreement between GYSELA and GENE is very good for $n \le 55$. Notably, GYSELA correctly captured the transition point. For $n = 60$ (not shown), the GYSELA growth rate ceases to increase with $n$, indicating a need to further boost the grid resolution and reduce the time
step.

<figure><img src="image_6.png" /><figcaption>Figure 5: Comparison of (a) the linear growth rate and (b) the real frequency of the $\beta_{\text{ref}} = 1.01\%$ case as a function of $n$ in the CBC scenario, for GENE and GYSELA. The transition from ITG to electron modes happens at $n = 35$. Approximately, $k_{\theta}\rho_{i} = 2nq/\rho_{*}$, so $n = 35$ corresponds to $k_{\theta}\rho_{i} \approx 0.54$.</figcaption></figure>

### 4.3. Nonlinear flux-driven electromagnetic ITG

To demonstrate GYSELA's capability of nonlinear flux-driven EM simulation, we start from the same equilibrium setup as Section 4.2. To save computational budget for testing, we selected a smaller system size $\rho_{*} = 1/100$ and a smaller ion-electron mass ratio $m_{i}/m_{e} = 100$. These are also the parameters chosen for GKNET [9] to study nonlinear ITG-KBM. The plasma $\beta_{\text{ref}}$ is set to 0.5% so the linear instability is in the ITG regime. We run GYSELA with the LWA solver, the align advection scheme, a grid resolution of $(N_{r}, N_{\theta}, N_{\phi}, N_{p_{\|}}, N_{\mu}) = (256, 256, 32, 128, 64)$, and $\Delta t = 1\omega_{ci}^{-1} \approx 0.0036(c_{s}/R_{0})^{-1}$. The toroidal domain is half a torus, so only modes with an even $n$ are allowed. The filter size is set to $\Delta m = 20$. The field solver updates its operator every 100 time steps to keep track of the new density profile to
be used in the polarisation and skin current terms. For a flux-driven sim-
ulation, collisions, heat exhaust and heat source are necessary. We enable
intra- and inter-species collisions with the effective collision rate $\nu^* = 0.1$
(banana regime). The heat exhaust is provided by a buffer region centred at
$r = 0.95a$ with a typical width $\Delta r = 0.05a$ and a radial diffusion coefficient
$D_r = 0.02$. Finally, we enable a heat source for both ions and electrons, with
typical temperature $T_{s,srce} = 1.1$ and strength $S_0 = 0.02$ and $0.01$ for ions
and electrons, respectively. The radial profile of the source is a Gaussian
function of width $0.1a$ centred at $r = 0$. Details of the source are given in
Appendix C. No particle source is added as the particle confinement time
is usually much longer compared to the energy confinement time.

The equilibrium was initiated by the quiet-start strategy for a total time
of $5000\omega_{ci}^{-1}$ using a time step $\Delta t = 5\omega_{ci}^{-1}$ and 8 toroidal points and a re-
duced heat source such that the profiles are not evolving significantly as
only neoclassical transport are included, while the other parameters are kept
unchanged. After that, a perturbation of $\delta F/F = 10^{-6}$ is added to ions
for all the near-diagonal modes. The simulation is then performed on the
ROME partition of the IRENE cluster with 64 nodes (128 AMD 7H12 CPUs
@ 2.6GHz per node) for 20000 iterations, costing approximately 68 hours of
wall time corresponding to a total of 0.56 million CPU hours.

The time trace of the flux-surface-averaged amplitude of $\Phi$ for different
toroidal mode numbers $n$ is plotted in Figure 6, showing linear growth then
the saturation at $t \approx 7000\omega_{ci}^{-1}$. The most unstable mode has $n = 12$ in the
linear stage, similar to GKNET [9] and TRIMEG-GKX [41]. The growth
rate of high-$n$ modes are twice of the low-$n$ modes in the late linear stage,
due to nonlinear excitation of those modes. The two-dimensional plot of $\Phi$
and $A_\parallel$ at $t = 4000\omega_{ci}^{-1}$ (linear stage), $7000\omega_{ci}^{-1}$ (early nonlinear stage) and
$20000\omega_{ci}^{-1}$ (steady-state) are shown in Figure 7.

The simulation enters the nonlinear stage at around $t = 7000$, with the
generation of zonal flows that tear apart the radially elongated structures
and break them into small eddies. This can also be shown by the oscil-
lation of the $n = 0$ mode at saturation in Figure 6. After that, turbu-
lence gradually spreads from mid-radius to the entire domain, establishing
steady-state density and temperature profiles, as shown in Figure 8. The
temperature gradient length relaxes to $R/L_{Ti} \sim 4.5$ with a flat profile in
the steady state. Figure 9 shows the averaged heat and particle flux over
the radial domain $0.3 \le r \le 0.7a$ as a function of time. The fluxes in the
steady-state ($t \ge 10000\omega_{ci}^{-1}$) are calculated to be $Q_i/Q_{GB} = 4.66 \pm 1.72$,
$Q_e/Q_{GB} = 2.71 \pm 1.16$, and $\Gamma/\Gamma_{GB} = 0.41 \pm 0.33$, where $Q_{GB} = n_{\text{ref}}T_{\text{ref}}c_s\rho_*^2$
and $\Gamma_{GB} = n_{\text{ref}}c_s\rho_*^2$ are the heat and particle fluxes in gyro-Bohm units, re-
spectively. The ratio of ion and electron heat flux is roughly 2:1, since the ion
heat source is twice as strong as the one for the electrons. The ion heat diffu-
sivity $\chi_i/\chi_{GB} = Q_i/Q_{GB} \times T_{\text{ref}}L_{Ti}/(T_i a) \approx 2.88$ is lower than that of the flux
tube $\delta f$ code GYRO in CBC condition [42] ($\chi_i/\chi_{GB} = 6.15$ at $R/L_T = 6.9$)
by a factor of 2 due to the flattening of the temperature profiles. It is higher
than the respective local $\delta f$ diffusivity $\chi_i/\chi_{GB} = 1.0$ at $R/L_T = 4.5$ due to
global effects near marginal stability ($R/L_T \sim 4.0$) [43].

Finally, we ran an extra simulation for $1000\omega_{ci}^{-1}$ with no heat source to
test the conservation laws. The ion and electron particle conservation are
shown in Figure 10. The conservation error is typically below 1%, except at
<figure><img src="image_7.png" /><figcaption>Figure 6: The amplitude for different Fourier modes of $\Phi$ at $r = 0.5a$ as a function of time.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: The $n \neq 0$ component of $\Phi$ and $A_\parallel$ on the $r-\theta$ plane at $t = 4000\omega_{ci}^{-1}$ (linear stage), $7000\omega_{ci}^{-1}$ (saturation), and $20000\omega_{ci}^{-1}$ (steady-state), respectively. Top row: $\Phi$; Bottom row: $A_\parallel$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: The evolution of the ion temperature gradient length $R/L_{T_i}$, the electron temperature gradient length $R/L_{T_e}$, and the density gradient length $R/L_n$ over time.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9: The ion heat flux, electron heat flux, and particle flux in gyro-Bohm unit averaged over $r = 0.3a$ to $r = 0.7a$ as a function of time.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10: Local particle conservation test for ions (left) and electrons (right), averaged over $20000\omega_{ci}^{-1} \le t \le 21000\omega_{ci}^{-1}$. The relative error is below 1% except at the edge $r/a \ge 0.85$, where the diffusive buffer region is located.</figcaption></figure>

the edge where a diffusive boundary condition is applied. Similarly, we plot the energy conservation and angular momentum conservation summing over species in Figure 11, again with a typical error again below 1%. The reason for summing over species is to average out the contribution from inter-species collisions, which should otherwise be accounted for separately.

## 5. Conclusion

We have developed and validated the electromagnetic (EM) extension of the global full-F gyrokinetic code GYSELA, enabling self-consistent simulations of electromagnetic turbulence and transport in magnetised plasmas. Comprehensive verification was carried out through a series of benchmark tests. The shear Alfvén wave simulations confirmed the expected dispersion relation and demonstrated the second-order convergence of the time integration scheme. Linear benchmarks in the Cyclone Base Case (CBC)
<figure><img src="image_12.png" /><figcaption>Figure 11: Local energy (left) and local toroidal angular momentum (right) conservation test, averaged over $20000\omega_{ci}^{-1} \le t \le 21000\omega_{ci}^{-1}$, summed over species. The relative error is below 1% except at the edge $r/a \ge 0.85$, where the diffusive buffer region is located.</figcaption></figure>

showed excellent agreement with established global gyrokinetic codes such as GENE, including the correct capture of the ITG-KBM transition and the mode structure. Finally, nonlinear flux-driven electromagnetic ITG simulations demonstrated the robustness of the implementation, achieving stable steady states with good particle, energy, and momentum conservation.

These developments establish the EM version of GYSELA as a reliable and versatile framework for investigating multi-scale electromagnetic plasma dynamics. Its capability for global, flux-driven simulations opens new opportunities to study the impact of electromagnetic effects on turbulent transport in high-beta plasmas. Moreover, the present work lays the foundation for incorporating fast-particle physics in future studies, enabling investigations of interactions between energetic particles, Alfvén eigenmodes, and background turbulence.

The GYSELA code is currently being upgraded to Gysela-X++ [44], which will feature X-point geometry and GPU acceleration. The advances
<table><caption>Table A.2: Normalised quantities.</caption><thead><tr><th>̂t</th><th>̂X</th><th>̂∇</th><th>̂Z<sub>s</sub></th><th>̂A<sub>s</sub></th><th>̂n<sub>s</sub></th><th>̂T<sub>s</sub></th></tr></thead><tbody><tr><td>ω<sub>ci</sub>t</td><td>X/ρ<sub>ref</sub></td><td>ρ<sub>ref</sub>∇</td><td>q<sub>s</sub>/q<sub>ref</sub></td><td>m<sub>s</sub>/m<sub>ref</sub></td><td>n<sub>s</sub>/n<sub>ref</sub></td><td>T<sub>s</sub>/T<sub>ref</sub></td></tr><tr><td>̂B</td><td>̂J</td><td>̂A<sub>∥</sub></td><td>̂Φ</td><td>̂v</td><td>̂μ</td><td>̂F</td></tr><tr><td>B/B<sub>0</sub></td><td>Jρ<sub>ref</sub>/B<sub>0</sub></td><td>A<sub>∥</sub>/(B<sub>0</sub>ρ<sub>ref</sub>)</td><td>q<sub>ref</sub>Φ/T<sub>ref</sub></td><td>v/v<sub>Tref</sub></td><td>μB<sub>0</sub>/T<sub>ref</sub></td><td>n<sub>ref</sub>F/v<sub>Tref</sub><sup>3</sup></td></tr></tbody></table>

presented in this work will serve as a key reference for the electromagnetic scheme in the future Gysela-X++.

## Acknowledgments

We would like to thank Alessandro Biancalani, Zhixin Lu and Roman Hatzky for useful discussions. This work is partly funded by Ministry of Education (MOE) AcRF Tier 1 grants RS02/23 and RG156/23, and National Research Foundation Singapore (NRF) core funding “Fusion Science for Clean Energy”, under the framework of Singapore Alliance with France for Fusion Energy (SAFE). The computational work for this article was partially performed on resources of the National Supercomputing Centre (NSCC), Singapore (Project 12003525). This project was provided with computer and storage resources by GENCI at TGCC thanks to the grant 2025-A0180502224 on the supercomputer Joliot Curie’s ROME and SKL partitions.
Appendix A. Normalised equations implemented in GYSELA

Using these normalised quantities, the contra-variant components of the
velocities are given by

$$
\begin{align*}
\hat{v}_E^i &= \frac{\hat{\mathbf{B}} \times \hat{\nabla} \mathcal{J}[\hat{\Phi}]}{\hat{B} \hat{B}_{\|s}^*} = \frac{1}{\hat{B}_{\|s}^*} [\mathcal{J}[\hat{\Phi}], \hat{X}^i], \\
\hat{v}_D^i &= \left[ \left( \hat{p}_{\|} + \frac{\hat{Z}_s}{\sqrt{\hat{A}_s}} \mathcal{J}[\hat{A}_{\|s}] \right) \left( \hat{p}_{\|} - \frac{\hat{Z}_s}{\sqrt{\hat{A}_s}} \mathcal{J}[\hat{A}_{\|H}] \right) + \hat{\mu} \hat{B} \right] \frac{1}{\hat{Z}_s \hat{B} \hat{B}_{\|s}^*} [\hat{B}, \hat{X}^i], \\
\hat{v}_{AS}^i &= - \frac{1}{\sqrt{\hat{A}_s}} \left( \hat{p}_{\|} - \frac{\hat{Z}_s}{\sqrt{\hat{A}_s}} \mathcal{J}[\hat{A}_{\|H}] \right) \frac{1}{\hat{B}_{\|s}^*} [\mathcal{J}[\hat{A}_{\|S}], \hat{X}^i], \\
\hat{v}_{AH}^i &= - \frac{\hat{p}_{\|}}{\sqrt{\hat{A}_s}} \frac{1}{\hat{B}_{\|s}^*} [\mathcal{J}[\hat{A}_{\|H}], \hat{X}^i], \tag{A.1}
\end{align*}
$$

where the Poisson bracket is given by $[F, G] = \mathbf{b} \cdot \hat{\nabla} F \times \hat{\nabla} G$. For $\hat{p}_{\|}$,

$$
\begin{align}
\sqrt{\hat{A}_s \hat{p}_{\|}} &= -\hat{Z}_s \partial_t \mathcal{J}[\hat{A}_{\|S}] - \left( \mathbf{b}_s^* - \frac{\hat{\mathbf{B}} \times \hat{\nabla} \mathcal{J}[\hat{A}_{\|S}]}{\hat{B} \hat{B}_{\|s}^*} \right) \nonumber \\
&\quad \cdot \left( \hat{\mu} \hat{\nabla} \hat{B} + \hat{Z}_s \hat{\nabla} \mathcal{J}[\hat{\Phi}] - \frac{\hat{Z}_s}{\sqrt{\hat{A}_s}} \hat{p}_{\|} \hat{\nabla} \mathcal{J}[\hat{A}_{\|H}] \right) \nonumber \\
&\quad + \frac{\sqrt{\hat{A}_s}}{\hat{B}} \left( \hat{p}_{\|} + \frac{\hat{Z}_s}{\sqrt{\hat{A}_s}} \mathcal{J}[\hat{A}_{\|S}] \right) (\hat{\mathbf{v}}_E + \hat{\mathbf{v}}_{AH}) \cdot \hat{\nabla} \hat{B}. \tag{A.2}
\end{align}
$$

The normalised number density and parallel current are computed from
the normalised distribution given by

$$
\begin{align}
\hat{n}_s &= 2\pi \int d\hat{\mu} d\hat{p}_{\|} \mathcal{J}^T[\hat{B}_{\|s}^* \hat{F}_s], \nonumber \\
\delta \hat{J}_{\|,s} &= 2\pi \frac{\beta_{ref} \hat{Z}_s}{\sqrt{\hat{A}_s}} \int d\hat{\mu} d\hat{p}_{\|} \mathcal{J}^T[\hat{p}_{\|} \hat{B}_{\|s}^* \hat{F}_s] - \hat{J}_{\|,eq}. \tag{A.3}
\end{align}
$$

The Poisson equation with LWA approximation now writes

$$
\sum_s \hat{Z}_s \hat{n}_s = - \sum_s \hat{\nabla} \cdot \frac{\hat{A}_s \hat{n}_{s0}}{\hat{B}^2} \hat{\nabla}_{\perp} \hat{\Phi}, \quad (A.4)
$$
while the parallel Ampère equation is given by

$$0 = (\hat{\nabla}_{\perp}^2 A_{\parallel H} + \hat{\nabla}_{\perp}^2 A_{\parallel S}) + \sum_s \delta \hat{J}_{\parallel,s} - \beta_{ref} \sum_s \frac{\hat{Z}_s^2 \hat{n}_{0,s}}{\hat{A}_s} A_{\parallel H}. \quad (A.5)$$

## Appendix B. Derivation of the mixed-variable gyrokinetics

The full procedure of producing the gyrokinetic Vlasov-Maxwell equations can be found from Brizard and Hahm (B&H) [3] Section IV and V, which cover the case of $A_{\parallel}$ in both the symplectic and Hamiltonian form. More specifically, a parameter $\beta$ is introduced, such that $A_{\parallel S} = \beta A_{\parallel}$ and $A_{\parallel H} = (1 - \beta)A_{\parallel}$. In B&H, only the cases with $\beta = 0$ or 1 are considered. But it is straightforward to generalise the procedure to cover arbitrary values of $\beta$. In this appendix, we will follow this procedure to derive the gyrocentre Hamiltonian up to the second order, then the Poisson-Ampère equations by using the action principle. In this section, we will use an overhead bar for the gyrocentre variables $\bar{Z} = (\bar{X}, \bar{p}_{\parallel}, \bar{\mu}, \bar{\zeta}, t)$ and no bar for the guiding-centre variables $\bar{Z} = (\bar{X}, v_{\parallel}, \mu, \zeta, t)$ to avoid confusions. Here, $\zeta$ is the gyro-angle. The difference between the particle location and the guiding centre location is marked by $\rho = \mathbf{x} - \mathbf{X}$.

### Appendix B.1. Derivation of the gyrocentre Hamiltonian

Adding $\Phi$ and $A_{\parallel}$, the 1-form in the guiding centre coordinates read

$$\gamma = \Gamma_0 + \Gamma_1 - (H_0 + H_1)dt, \quad (B.1)$$

where

$$\Gamma_0 = \left( \mathbf{A} + \frac{m_s v_{\parallel}}{q_s} \mathbf{b} \right) \cdot d(\mathbf{X} + \rho) + \frac{m_s}{q_s} \mu d\zeta, \quad (B.2)$$
and

$$H_0 = \frac{m_s}{2} v_{\parallel}^2 + \mu B. \qquad (B.3)$$

The perturbed symplectic part of the 1-form in guiding centre coordinates is given by (Eq. (153) in B&H)

$$\begin{aligned} \Gamma_1 &= q_s A_{\parallel}(\mathbf{X} + \rho, t)\mathbf{b} \cdot d(\mathbf{X} + \rho) \\ &= q_s A_{\parallel}(\mathbf{X} + \rho, t)\mathbf{b} \cdot \left(d\mathbf{X} + \frac{\partial \rho}{\partial \mu}d\mu + \frac{\partial \rho}{\partial \zeta}d\zeta\right). \end{aligned} \qquad (B.4)$$

The $d\mu$ and $d\zeta$ terms vanish since $\mathbf{b} \cdot \rho = 0$. This gives the components $\Gamma_{1a} = q_s A_{\parallel}(\mathbf{X} + \rho, t)b_a$, if $a$ is an index of the spatial dimensions, and $\Gamma_{1a} = 0$ otherwise. The first order perturbed Hamiltonian is just $H_1 = q_s \Phi(\mathbf{X} + \rho)$. The effective first order Hamilationian is therefore given by (Eq. (125) in B&H)

$$K_1 = H_1 - \Gamma_{1a} \dot{Z}_0^a = q_s (\Phi - v_{\parallel} A_{\parallel}), \qquad (B.5)$$

where $\dot{Z}_0^a$ is the guiding centre equilibrium motion. In this context, only the parallel components of the spatial dimensions are useful, i.e. $\dot{Z}_0^a = v_{\parallel} b^a$ if $a$ is indexing spatial dimensions.

As stated by the mixed-variable scheme, one retains only $A_{\parallel S}$ in the gyrocentre equivalent of $\Gamma_1$, giving (Eq. (156) B&H)

$$\bar{\Gamma}_1 = q_s \mathcal{J}[A_{\parallel S}(\bar{\mathbf{X}}, t)]\mathbf{b} \cdot d\bar{\mathbf{X}}, \qquad (B.6)$$

so that $\bar{\Gamma}_{1a} = q_s \mathcal{J}[A_{\parallel}]b_a$ for spatial dimensions, and 0 otherwise. The first-order gyrocentre Hamiltonian is derived from Eq. (126) of B&H, given by

$$\bar{H}_1 = \mathcal{J}[K_1] + \bar{\Gamma}_{1a} \dot{Z}_0^a = q_s (\mathcal{J}[\Phi] - \bar{p}_{\parallel} \mathcal{J}[A_{\parallel H}]), \qquad (B.7)$$
where the guiding centre variables are replaced by gyrocentre ones.

The guiding centre perturbed field is given by (considering an arbitrary $\beta$ in Eq. (170) of B&H),

$$ \delta\psi_{gc} = \Phi - p_{\parallel}(A_{\parallel S} + A_{\parallel H}) + p_{\parallel}\mathcal{J}[A_{\parallel S}], \qquad (B.8) $$

B&H defines the gyro-angle dependent part of $\delta\psi_{gc}$ as

$$ q_s \delta\tilde{\psi}_{gc} = q_s \delta\psi_{gc} - q_s \mathcal{J}[\delta\psi_{gc}] = \mathcal{H}_1 - \mathcal{J}[\mathcal{H}_1], \qquad (B.9) $$

where $\mathcal{H}_1 = q_s\Phi - q_s p_{\parallel}(A_{\parallel S} + A_{\parallel H})$. The first order generating function is given by (Eq. (172) B&H)

$$ S_1 = \frac{q_s}{\Omega_0} \delta\tilde{\Psi}_{gc} = \frac{q_s}{\Omega_s} \int^{\bar{\zeta}} \delta\tilde{\psi}_{gc} d\bar{\zeta} \qquad (B.10) $$

where $\Omega_s$ is the gyro-frequency.

Once these quantities are defined, one can readily substitute them into the equation of second-order Hamiltonian (Eq. (129) B&H). After some algebra (details can be found in [45] Section H.5), one arrives at

$$ \begin{aligned} \bar{H}_2 = & \frac{q_s^2}{2m_s} \{ \mathcal{J}[(A_{\parallel H} + A_{\parallel S})^2] - \mathcal{J}[A_{\parallel S}]^2 - 2\mathcal{J}[A_{\parallel S}] \mathcal{J}[A_{\parallel H}] \} \\ & - \frac{1}{2B} \partial_{\mu} [\mathcal{J}[\mathcal{H}_1^2] - (\mathcal{J}[\mathcal{H}_1])^2], \end{aligned} \qquad (B.11) $$

## Appendix B.2. The Poisson-Ampère equations

The total action functional is given by

$$ S = S_{em} + \sum_{s} S_{fp,s}, \qquad (B.12) $$

where

$$ S_{em} = \frac{1}{2} \int dt dx \left( \epsilon_0 |\nabla \Phi|^2 - \mu_0^{-1} |\nabla_{\perp} A_{\parallel S} + \nabla_{\perp} A_{\parallel H}|^2 \right), \qquad (B.13) $$
$$ S_{fp,s} = - \int dt d\bar{Z} \bar{F}_{1s} \bar{H}_1 + \bar{F}_s \bar{H}_2, \qquad (B.14) $$

where the former is expressed in real space and the latter in gyrocentre phase space $\bar{Z}$ of species s. Now $\bar{F}_{1s}$ is the first-order “perturbed” distribution function, such that $\bar{F}_{1s}\bar{H}_1$ is second-order. For the $\bar{F}_s\bar{H}_2$ term, one should in principle replace $\bar{F}_s$ by $\bar{F}_{s,0}$, the unperturbed zeroth-order distribution, such that $S_{fp,s}$ is consistent with second order. One can argue that the difference between the slow guiding centre and gyrocentre is second order, and is already accounted for by the second-order Hamiltonian. One can just replace gyrocentre coordinates $\bar{Z}$ by guiding centre coordinates $Z$, so

$$ S_{fp,s} = S_{skin,s} + S_{part,s}, \qquad (B.15) $$

with

$$ S_{skin,s} = - \frac{q_s^2}{2m_s} \int dt dZ \{ \mathcal{J}[(A_{\parallel H} + A_{\parallel S})^2] - \mathcal{J}[A_{\parallel S}]^2 \\ - 2\mathcal{J}[A_{\parallel S}] \mathcal{J}[A_{\parallel H}] \} F_s, \qquad (B.16) $$

$$ S_{part,s} = -q_s \int dt dZ (F_s - F_{s,eq}) (\mathcal{J}[\Phi] - p_{\parallel} \mathcal{J}[A_{\parallel H}]) \\ - \frac{1}{2B} \int dt dZ \partial_{\mu} F_s [\mathcal{J}[\mathcal{H}_1^2] - (\mathcal{J}[\mathcal{H}_1])^2], \qquad (B.17) $$

where an integration by part is applied to move the $\mu$ derivative to $F_s$. The next step is to change back to real coordinate and note that $\int dt dZ F \mathcal{J}[G] = \int dt dZ G \mathcal{J}^T[B_{\parallel s}^* F]$, where $dZ = dxdv'$ and $dv' = d\mu dp_{\parallel} d\zeta$, so

$$ \begin{aligned} S_{skin,s} = & -\frac{q_s^2}{2m_s} \int dt d\mathbf{z} (A_{\parallel H} + A_{\parallel S})^2 \mathcal{J}^T[B_{\parallel s}^* F_s] \\ & + \frac{q_s^2}{m_s} \int dt d\mathbf{z} A_{\parallel H} \mathcal{J}^T[B_{\parallel s}^* F_s \mathcal{J}[A_{\parallel S}]], \\ & + \frac{q_s^2}{2m_s} \int dt d\mathbf{z} A_{\parallel S} \mathcal{J}^T[B_{\parallel s}^* F_s \mathcal{J}[A_{\parallel S}]], \end{aligned} \qquad (B.18) $$
$$
\begin{align}
S_{part,s} = & -q_s \int dt dz \Phi \mathcal{J}^T [B_{\parallel s}^* \delta F_s] \nonumber \\
& + q_s \int dt dz A_{\parallel H} \mathcal{J}^T [B_{\parallel s}^* p_{\parallel} \delta F_s] \nonumber \\
& - \frac{1}{2B} \int dt dz (\mathcal{H}_1^2 \mathcal{J}^T [\partial_{\mu} F_s] - \mathcal{H}_1 \mathcal{J}^T [\partial_{\mu} F_s (\mathcal{J}[\mathcal{H}_1])]). \tag{B.19}
\end{align}
$$

Finally we need to perform variation with respect to Φ and A_∥H.

$$
\frac{\delta S_{em}}{\delta \Phi} = -\epsilon_0 \nabla^2 \Phi \qquad (B.20)
$$

$$
\frac{\delta S_{em}}{\delta A_{\parallel H}} = \mu_0^{-1} \delta A_{\parallel H} (\nabla_{\perp}^2 A_{\parallel H} + \nabla_{\perp}^2 A_{\parallel S}), \quad (B.21)
$$

$$
\frac{\delta S_{\text{skin},s}}{\delta A_{\parallel H}} = - \frac{q_s^2}{m_s} \int d\mathbf{v}' \{ A_{\parallel} \mathcal{J}^T [B_{\parallel s}^* F_s] - \mathcal{J}^T [B_{\parallel s}^* F_s \mathcal{J}[A_{\parallel s}]] \}, \quad (B.22)
$$

$$
\begin{equation}
\begin{split}
\frac{\delta S_{part}}{\delta \Phi} = {}& -q_s \int d\mathbf{v}' \mathcal{J}^T [B_{\parallel s}^* \delta F_s] - \frac{q_s^2}{B} \int d\mathbf{v}' (\Phi - p_{\parallel} A_{\parallel}) \mathcal{J}^T [B_{\parallel s}^* \partial_{\mu} F_s] \\
& + \frac{q_s^2}{B} \int d\mathbf{v}' \mathcal{J}^T [B_{\parallel s}^* \partial_{\mu} F_s (\mathcal{J}[\Phi] - p_{\parallel} \mathcal{J}[A_{\parallel}])],
\end{split}
\tag{B.23}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
\frac{\delta S_{part}}{\delta A_{\parallel H}} = {}& q_s \int d\mathbf{v}' \mathcal{J}^T [p_{\parallel} B_{\parallel s}^* \delta F_s] \\
& + \frac{q_s^2}{B} \int d\mathbf{v}' (\Phi - p_{\parallel} A_{\parallel}) p_{\parallel} \mathcal{J}^T [B_{\parallel s}^* \partial_{\mu} F_s] \\
& - \frac{q_s^2}{B} \int d\mathbf{v}' \mathcal{J}^T [B_{\parallel s}^* \partial_{\mu} F_s (\mathcal{J}[\Phi] - p_{\parallel} \mathcal{J}[A_{\parallel}])].
\end{aligned}
\tag{B.24}
\end{equation}
$$

Then, equating $\delta S$ to zero gives the gyrokinetic Poisson equation

$$
\begin{align}
0 ={}& -\epsilon_0 \nabla^2 \Phi - \sum_s q_s \delta n_s - \sum_s \int d\mathbf{v}' (\Phi - p_{\parallel} A_{\parallel}) \mathcal{J}^T [\partial_{\mu} F_s] \nonumber \\
& + \sum_s \frac{q_s^2}{B} \int d\mathbf{v}' \mathcal{J}^T [B_{\parallel s}^* \partial_{\mu} F_s (\mathcal{J}[\Phi] - p_{\parallel} \mathcal{J}[A_{\parallel}])], \tag{B.25}
\end{align}
$$
and the gyrokinetics Ampère equation

$$
\begin{equation}
\begin{split}
0 ={}& \mu_0^{-1} (\nabla_{\perp}^2 A_{\parallel H} + \nabla_{\perp}^2 A_{\parallel S}) + \sum_s \delta J_{\parallel s} \\
     & - \sum_s \frac{q_s^2}{m_s} (A_{\parallel H} + A_{\parallel S}) \int d\mathbf{v}' \mathcal{J}^T [B_{\parallel s}^* F_s] \\
     & + \sum_s \frac{q_s^2}{m_s} \int d\mathbf{v}' \mathcal{J}^T [B_{\parallel s}^* F_s \mathcal{J}[A_{\parallel S}]] \\
     & + \sum_s \frac{q_s^2}{B} \int d\mathbf{v}' (\Phi - p_{\parallel} A_{\parallel}) p_{\parallel} \mathcal{J}^T [B_{\parallel s}^* \partial_{\mu} F_s] \\
     & - \sum_s \frac{q_s^2}{B} \int d\mathbf{v}' p_{\parallel} \mathcal{J}^T [B_{\parallel s}^* \partial_{\mu} F_s (\mathcal{J}[\Phi] - p_{\parallel} \mathcal{J}[A_{\parallel}])],
\end{split}
\tag{B.26}
\end{equation}
$$

with density and current defined as

$$
\delta n_s = \int d\mathbf{v}' \mathcal{J}^T [B_{\parallel s}^* \delta F_s], \qquad (B.27)
$$

$$
\delta \bar{J}_{\parallel, s} = q_s \int d\mathbf{v}' \mathcal{J}^T [p_{\parallel} B_{\parallel s}^* \delta F_s]. \quad (B.28)
$$

These equations are consistent with those derived by using a pull-back
transformation of the distribution function [46]. They are exact - no ap-
proximation are made so far within the nonlinear gyrokinetics framework of
B&H. Note that we did not approximate F_s in the skin current term and
polarisation term by anything yet for completeness.

Appendix B.3. Approximations

Some approximations are needed to proceed.

Appendix B.3.1. Plasma approximation

One can drop the $\epsilon_0$ term and use the so-called plasma approximation.
Then the Poisson equation becomes the quasi-neutrality equation. This is
allowed for perturbations with $k\lambda_D \ll 1$.
Appendix B.3.2. Boussinesq approximation

In computing $\partial_\mu F_s$ in the polarisation term, one can replace it with a shifted Maxwellian distribution with no perturbation, so $\partial_\mu F_{s0} = -BF_{s0}/T_{s0}$. This is just to say $\bar{F}_s\bar{H}_2 \approx \bar{F}_{s0}\bar{H}_2$, consistent with second order. However, there is no restriction on what $\bar{F}_{s,0}$ can be. As we are running a full-F code, it is more useful to think $\bar{F}_{s,0}$ to evolve as a slow function of time.

Another point is with energetic particles which do not have a Maxwellian distribution function. Also, the gyroaverage operator might not be trivial for $F_{s0}$. Nevertheless, the energetic particles (EPs) are usually quite small in density so even if there are some errors it might not be catastrophic. We can discuss EPs separately.

With this approximation,

$$ \sum_s q_s \delta n_s = \sum_s q_s^2 \int d\mathbf{v}' B_{\|s}^* \left\{ \frac{F_{s0}}{T_{s0}} (\Phi - p_{\|} A_{\|}) - \mathcal{J}^T \left[ \frac{F_{s0}}{T_{s0}} (\mathcal{J}[\Phi] - p_{\|} \mathcal{J}[A_{\|}]) \right] \right\}. \quad (B.29) $$

Now if we elect to use Boussinesq approximation, we must do so also for the skin-current term, such that the ordering remains consistent.

$$ \begin{aligned} 0 = & \mu_0^{-1}(\nabla_{\perp}^2 A_{\|H} + \nabla_{\perp}^2 A_{\|S}) + \sum_s \delta J_{\|s} - \sum_s \frac{q_s^2 n_{0,s}}{m_s} A_{\|H} \\ & - \sum_s \frac{q_s^2}{m_s} \int d\mathbf{v}' B_{\|s}^* (F_s A_{\|S} - \mathcal{J}^T [F_s \mathcal{J}[A_{\|S}]] \\ & - \sum_s q_s^2 \int d\mathbf{v}' B_{\|s}^* p_{\|} \left\{ \frac{F_{s0}}{T_{s0}} (\Phi - p_{\|} A_{\|}) \right. \\ & \left. - \mathcal{J}^T \left[ \frac{F_{s0}}{T_{s0}} (\mathcal{J}[\Phi] - p_{\|} \mathcal{J}[A_{\|}]) \right] \right\}. \end{aligned} \quad (B.30) $$
Appendix B.3.3. Long-wave-length approximation

We need to expand the gyroaverage operator and its conjugate for practical numerical computation. One of the usual approximation is the “long-wave-length” approximation, which gives

$$ \mathcal{J}^T[F\mathcal{J}[H_1]] \approx F\mathcal{H}_1 + \nabla \cdot \left( \frac{m_s\mu}{q_s^2 B} F \nabla_{\perp} \mathcal{H}_1 \right). \quad (B.31) $$

Using this approximation, the Poisson equation looks like

$$ \sum_s q_s \delta n_s = -\nabla \cdot \sum_s \int d\mathbf{v}' B_{\parallel s}^* \frac{m_s \mu}{B} \frac{F_{s0}}{T_{s0}} \nabla_{\perp} (\Phi - p_{\parallel} A_{\parallel}), \quad (B.32) $$

and the Ampère equation becomes

$$ \begin{aligned} 0 ={}& \mu_0^{-1}(\nabla_{\perp}^2 A_{\parallel H} + \nabla_{\perp}^2 A_{\parallel S}) + \sum_s \delta J_{\parallel s} - \sum_s \frac{q_s^2 n_{0,s}}{m_s} A_{\parallel H} \\ & + \nabla \cdot \sum_s \int d\mathbf{v}' B_{\parallel s}^* \left( \frac{\mu F_s}{B} \nabla_{\perp} A_{\parallel S} \right) \\ & + \nabla \cdot \sum_s m_s \int d\mathbf{v}' B_{\parallel s}^* \frac{F_{s0}\mu}{T_{s0}B} \nabla_{\perp} (\Phi - p_{\parallel} A_{\parallel}). \end{aligned} \quad (B.33) $$

Appendix B.3.4. Un-shifted Maxwellian

To eliminate the cross-terms, one can assume an un-shifted Maxwellian distribution for $F_{s0}$. This will make $F_{s0}$ an even function in $s_0$, and thus, removing all integral odd in $p_{\parallel}$. We also need to note that for an unshifted Maxwellian, the moment of $\mu B$ will give $n_{s0}T_{s0}$, while the moment of $\mu B p_{\parallel}^2$ gives $n_{s0}T_{s0}^2/m_s$. Therefore, the Poisson equation becomes (18), which is the usual Poisson for LWA. The Ampère equation becomes (24).
# Appendix C. Energy source used in nonlinear simulation

We use a pure heat source in Section 4.3, with its form given by [17, 47]

$$ S_{\text{heat}} = \left[ \bar{v}_{G\parallel,s}^2 - \frac{1}{2} - \frac{J_{\parallel B}}{2 - J_{\parallel B}^2} (2 - \bar{\mu}_s)(2\bar{v}_{G\parallel,s} - J_{\parallel B}) \right] \times \frac{2S_0}{(2\pi T_{s,srce}/m_s)^{3/2} T_{s,srce}} S_r(r) e^{-\bar{v}_{G\parallel,s}^2 - \bar{\mu}_s}, \quad (C.1) $$

where $\bar{\mu}_s = \mu_s B/T_{s,srce}$, $\bar{v}_{G\parallel,s} = v_\parallel/\sqrt{2T_{s,srce}/m_s}$, and $J_{\parallel B} = \sqrt{2m_s T_{s,srce}}/(q_s B^2)\mu_0 J_\parallel$.

The radial envelope $S_r(r)$ is a Gaussian function centred at $r = 0$ with a given width, normalised to unity when integrated over the radius.

# Declaration of generative AI and AI-assisted technologies in the manuscript preparation process

During the preparation of this work the author(s) used ChatGPT-5 in order to polish the language. After using this tool/service, the author(s) reviewed and edited the content as needed and take(s) full responsibility for the content of the published article.

# References

[1] J. Citrin, P. Mantica, Overview of tokamak turbulence stabilization by fast ions, Plasma Physics and Controlled Fusion 65 (3) (2023). doi:10.1088/1361-6587/acab2b.

[2] Y.-S. Na, T. Hahm, P. Diamond, A. Di Siena, J. Garcia, Z. Lin, How fast ions mitigate turbulence and enhance confinement in tokamak fusion plasmas, Nature Reviews Physics (2025) 1–13.
[3] a. Brizard, T. Hahm, Foundations of nonlinear gyrokinetic theory, Reviews of Modern Physics 79 (2) (2007) 421–468.
doi:10.1103/RevModPhys.79.421.
URL http://link.aps.org/doi/10.1103/RevModPhys.79.421

[4] X. Garbet, Y. Idomura, L. Villard, T. H. Watanabe, Gyrokinetic simulations of turbulent transport, Nuclear Fusion 50 (4) (2010).
doi:10.1088/0029-5515/50/4/043002.

[5] R. Kleiber, M. Borchardt, R. Hatzky, A. Könies, H. Leyh, A. Mishchenko, J. Riemann, C. Slaby, J. García-Regaña, E. Sánchez, M. Cole, EUTERPE: A global gyrokinetic code for stellarator geometry, Computer Physics Communications 295 (2024) 109013.
doi:10.1016/j.cpc.2023.109013.
URL https://linkinghub.elsevier.com/retrieve/pii/S0010465523003582

[6] T. Görler, X. Lapillonne, S. Brunner, T. Dannert, F. Jenko, F. Merz, D. Told, The global version of the gyrokinetic turbulence code GENE, Journal of Computational Physics 230 (18) (2011) 7053–7071.
doi:10.1016/j.jcp.2011.05.034.
URL https://linkinghub.elsevier.com/retrieve/pii/S0021999111003457

[7] D. Michels, A. Stegmeir, P. Ulbl, D. Jarema, F. Jenko, GENE-X: A full-f gyrokinetic turbulence code based on the flux-coordinate independent approach, Computer Physics Communications 264 (2021) 107986. doi:10.1016/j.cpc.2021.107986.
URL https://linkinghub.elsevier.com/retrieve/pii/S0010465521000989
[8] M. Maurer, A. Bañón Navarro, T. Dannert, M. Restelli, F. Hindenlang, T. Görler, D. Told, D. Jarema, G. Merlo, F. Jenko, GENE-3D: A global gyrokinetic turbulence code for stellarators, Journal of Computational Physics 420 (2020) 109694, publisher: Elsevier Inc. doi:10.1016/j.jcp.2020.109694.

URL https://doi.org/10.1016/j.jcp.2020.109694

[9] A. Ishizawa, K. Imadera, Y. Nakamura, Y. Kishimoto, Global gyrokinetic simulation of turbulence driven by kinetic ballooning mode, Physics of Plasmas 26 (8) (2019) 082301. doi:10.1063/1.5100308.

URL https://pubs.aip.org/pop/article/26/8/082301/631822/Global-gyrokinetic-si

[10] J. Seo, J.-M. Kwon, S. Yi, L. Qi, A new hybrid simulation model for tokamak plasma turbulence, Computer Physics Communications 258 (2021) 107626. doi:10.1016/j.cpc.2020.107626.

URL https://linkinghub.elsevier.com/retrieve/pii/S0010465520303003

[11] A. G. Peeters, Y. Camenen, F. J. Casson, W. A. Hornsby, A. P. Snodin, D. Strintzi, G. Szepesi, The nonlinear gyro-kinetic flux tube code GKW, Computer Physics Communications 180 (12) (2009) 2650-2672, publisher: Elsevier B.V. doi:10.1016/j.cpc.2009.07.001.

URL http://dx.doi.org/10.1016/j.cpc.2009.07.001

[12] Z. Lin, T. S. Hahm, W. Lee, W. M. Tang, R. B. White, Turbulent transport reduction by zonal flows: Massively parallel simulations, Science 281 (5384) (1998) 1835-1837.

[13] L. Ye, Y. Xu, X. Xiao, Z. Dai, S. Wang, A gyrokinetic continuum code
based on the numerical lie transform (nlt) method, Journal of Computational Physics 316 (2016) 180–192.

[14] E. Lanti, N. Ohana, N. Tronko, T. Hayward-Schneider, A. Bottino, B. F. McMillan, A. Mishchenko, A. Scheinberg, A. Biancalani, P. Angelino, S. Brunner, J. Dominski, P. Donnel, C. Gheller, R. Hatzky, A. Jocksch, S. Jolliet, Z. X. Lu, J. P. Martin Collar, I. Novikau, E. Sonnendrücker, T. Vernay, L. Villard, ORB5: A global electromagnetic gyrokinetic code using the PIC approach in toroidal geometry, Computer Physics Communications 251, arXiv: 1905.01906 (2020). doi:10.1016/j.cpc.2019.107072.

[15] Z. Lu, G. Meng, M. Hoelzl, P. Lauber, The development of an implicit full $f$ method for electromagnetic particle simulations of Alfvén waves and energetic particle physics, Journal of Computational Physics 440 (2021) 110384. doi:10.1016/j.jcp.2021.110384.

URL https://linkinghub.elsevier.com/retrieve/pii/S0021999121002795

[16] R. Hager, S. Ku, A. Y. Sharma, C. S. Chang, R. M. Churchill, A. Scheinberg, Electromagnetic total-$f$ algorithm for gyrokinetic particle-in-cell simulations of boundary plasma in XGC, Physics of Plasmas 29 (11) (2022) 112308. doi:10.1063/5.0097855.

URL https://pubs.aip.org/pop/article/29/11/112308/2843979/Electromagnetic-tot

[17] V. Grandgirard, J. Abiteboul, J. Bigot, T. Cartier-Michaud, N. Crouseilles, G. Dif-Pradalier, C. Ehrlacher, D. Esteve, X. Garbet, P. Ghendrih, G. Latu, M. Mehrenberger, C. Norscini, C. Passeron, F. Rozar,
Y. Sarazin, E. Sonnendrücker, A. Strugarek, D. Zarzoso, A 5D gyrokinetic full-f global semi-Lagrangian code for flux-driven ion turbulence simulations, Computer Physics Communications 207 (2016) 35–68, publisher: Elsevier B.V. doi:10.1016/j.cpc.2016.05.007.

URL http://dx.doi.org/10.1016/j.cpc.2016.05.007

[18] G. Dif-Pradalier, P. H. Diamond, V. Grandgirard, Y. Sarazin, J. Abiteboul, X. Garbet, P. Ghendrih, A. Strugarek, S. Ku, C. S. Chang, On the validity of the local diffusive paradigm in turbulent plasma transport, Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 82 (2) (2010) 1–4. doi:10.1103/PhysRevE.82.025401.

[19] G. Dif-Pradalier, G. Hornung, P. Ghendrih, Y. Sarazin, F. Clairet, L. Vermare, P. H. Diamond, J. Abiteboul, T. Cartier-Michaud, C. Ehrlacher, D. Estève, X. Garbet, V. Grandgirard, O. D. Gürcan, P. Hennequin, Y. Kosuga, G. Latu, P. Maget, P. Morel, C. Norscini, R. Sabot, A. Storelli, Finding the Elusive E × B Staircase in Magnetized Plasmas, Physical Review Letters 114 (8) (2015) 085004. doi:10.1103/PhysRevLett.114.085004.

URL https://link.aps.org/doi/10.1103/PhysRevLett.114.085004

[20] C. Gillot, Model reduction for tokamak plasma turbulence: beyond fluid and quasi-linear descriptions, Ph.D. thesis, Aix Marseille Université (2020).

[21] Y. Chen, S. Parker, Gyrokinetic turbulence simulations with kinetic electrons, Physics of Plasmas 8 (5) (2001) 2095–2100.
doi:10.1063/1.1351828.

URL https://pubs.aip.org/pop/article/8/5/2095/858650/Gyrokinetic-turbulence-s

[22] W. W. Lee, Gyrokinetic particle simulation model, Journal of Computational Physics 72 (1) (1987) 243–269.

[23] A. Mishchenko, A. Könies, R. Kleiber, M. Cole, Pullback transformation in gyrokinetic electromagnetic simulations, Physics of Plasmas 21 (9) (2014).

[24] A. Mishchenko, M. Cole, R. Kleiber, A. Könies, New variables for gyrokinetic electromagnetic simulations, Physics of Plasmas 21 (5) (2014).

[25] R. Kleiber, R. Hatzky, A. Könies, A. Mishchenko, E. Sonnendrücker, An explicit large time step particle-in-cell scheme for nonlinear gyrokinetic simulations in the electromagnetic regime, Physics of Plasmas 23 (3) (2016) 032501. doi:10.1063/1.4942788.

URL https://pubs.aip.org/pop/article/23/3/032501/1017481/An-explicit-large-time

[26] A. Mishchenko, A. Bottino, A. Biancalani, R. Hatzky, T. Hayward-Schneider, N. Ohana, E. Lanti, S. Brunner, L. Villard, M. Borchardt, R. Kleiber, A. Könies, Pullback scheme implementation in ORB5, Computer Physics Communications 238 (2019) 194–202.
doi:10.1016/j.cpc.2018.12.002.

URL https://linkinghub.elsevier.com/retrieve/pii/S0010465518304181

[27] E. Sonnendrücker, J. Roche, P. Bertrand, A. Ghizzo, The Semi-Lagrangian Method for the Numerical Resolution of the Vlasov Equation, Journal of Computational Physics 149 (2) (1999) 201–220.
doi:10.1006/jcph.1998.6148.

URL https://linkinghub.elsevier.com/retrieve/pii/S0021999198961484

[28] P. Donnel, K. Obrejan, Y. Sarazin, R. Bigué, E. Bourne, G. Brochard, L. de Gianni, G. Dif-Pradalier, X. Garbet, V. Grandgirard, et al., An adaptive quasi-neutrality solver for full-f flux-driven gyrokinetic simulations of tokamak plasmas in presence of poloidal asymmetries (2025).

[29] J. Abiteboul, X. Garbet, V. Grandgirard, S. J. Allfrey, P. Ghendrih, G. Latu, Y. Sarazin, A. Strugarek, Conservation equations and calculation of mean flows in gyrokinetics, Physics of Plasmas 18 (8) (2011) 082503. doi:10.1063/1.3620407.

URL https://pubs.aip.org/pop/article/18/8/082503/316799/Conservation-equation

[30] Z. Lu, G. Meng, R. Hatzky, E. Sonnendrücker, A. Mishchenko, M. Hoelzl, Generalized mixed variable-pullback scheme with non-ideal ohm's law for electromagnetic gyrokinetic simulations, arXiv preprint arXiv:2509.19553 (2025).

[31] E. Bourne, P. Leleux, K. Kormann, C. Kruse, V. Grandgirard, Y. Güçlü, M. J. Kühn, U. Rüde, E. Sonnendrücker, E. Zoni, Solver comparison for Poisson-like equations on tokamak geometries, Journal of Computational Physics 488 (2023) 112249. doi:10.1016/j.jcp.2023.112249.

URL https://linkinghub.elsevier.com/retrieve/pii/S0021999123003443

[32] G. Latu, M. Mehrenberger, Y. Güçlü, M. Ottaviani, E. Sonnendrücker, Field-Aligned Interpolation for Semi-Lagrangian Gyrokinetic Simulations, Journal of Scientific Computing 74 (3) (2018) 1601–1650.
doi:10.1007/s10915-017-0509-5.

URL http://link.springer.com/10.1007/s10915-017-0509-5

[33] F. Hariri, M. Ottaviani, A flux-coordinate independent field-aligned approach to plasma turbulence simulations, Computer Physics Communications 184 (11) (2013) 2419-2429.

[34] Y. Li, P. J. Morrison, S. Possanner, E. Sonnendrücker, On the poisson brackets of hybrid plasma models with kinetic ions and massless electrons, arXiv preprint arXiv:2509.11160 (2025).

[35] S. Jolliet, B. McMillan, L. Villard, T. Vernay, P. Angelino, T. Tran, S. Brunner, A. Bottino, Y. Idomura, Parallel filtering in global gyrokinetic simulations, Journal of Computational Physics 231 (3) (2012) 745-758. doi:10.1016/j.jcp.2011.01.029.

URL https://linkinghub.elsevier.com/retrieve/pii/S0021999111000544

[36] P. Angelino, A. Bottino, R. Hatzky, S. Jolliet, O. Sauter, T. M. Tran, L. Villard, On the definition of a kinetic equilibrium in global gyrokinetic simulations, Physics of Plasmas 13 (5) (2006) 052304. doi:10.1063/1.2193947.

URL https://pubs.aip.org/pop/article/13/5/052304/1032218/On-the-definition-of

[37] a. Biancalani, a. Bottino, P. Lauber, D. Zarzoso, Numerical validation of the electromagnetic gyrokinetic code NEMORB on global axisymmetric modes, Nuclear Fusion 54 (10) (2014) 104004. doi:10.1088/0029-5515/54/10/104004.

URL http://stacks.iop.org/0029-5515/54/i=10/a=104004?key=crossref.de735287789
[38] T. Görler, N. Tronko, W. A. Hornsby, A. Bottino, R. Kleiber, C. Norscini, V. Grandgirard, F. Jenko, E. Sonnendrücker, Intercode comparison of gyrokinetic global electromagnetic modes, Physics of Plasmas 23 (7) (2016) 072503. doi:10.1063/1.4954915.

URL https://pubs.aip.org/pop/article/23/7/072503/594661/Intercode-comparison-

[39] A. M. Dimits, G. Bateman, M. A. Beer, B. I. Cohen, W. Dorland, G. W. Hammett, C. Kim, J. E. Kinsey, M. Kotschenreuther, A. H. Kritz, L. L. Lao, J. Mandrekas, W. M. Nevins, S. E. Parker, A. J. Redd, D. E. Shumaker, R. Sydora, J. Weiland, Comparisons and physics basis of tokamak transport models and turbulence simulations, Physics of Plasmas 7 (3) (2000) 969–983. doi:10.1063/1.873896.

[40] M. D. J. Cole, A. Mishchenko, A. Bottino, C. S. Chang, Tokamak ITG-KBM transition benchmarking with the mixed variables/pullback transformation electromagnetic gyrokinetic scheme, Physics of Plasmas 28 (3) (2021) 034501. doi:10.1063/5.0030937.

URL https://pubs.aip.org/pop/article/28/3/034501/835105/Tokamak-ITG-KBM-trans

[41] Z. Lu, G. Meng, R. Hatzky, P. Lauber, M. Hoelzl, Trimeg-gkx: an electromagnetic gyrokinetic particle code with a piecewise field-aligned finite element method for micro-and macro-instability studies in tokamak core plasmas, arXiv preprint arXiv:2504.21837 (2025).

[42] J. Kinsey, A. Electrons, D. D. Gradients, H. Dilution, D. H. Dilutions, The gyro nonlinear gyrokinetic simulation database, Tech. rep., Technical report, General Atomics, San Diego, CA (2007).
[43] C. Gillot, G. Dif-Pradalier, Y. Sarazin, C. Bourdelle, A. B. Navarro, Y. Camenen, J. Citrin, A. Di Siena, X. Garbet, P. Ghendrih, et al., The problem of capturing marginality in model reductions of turbulence, Plasma Physics and Controlled Fusion 65 (5) (2023) 055012.

[44] E. Bourne, V. Grandgirard, Y. Asahi, J. Bigot, P. Donnel, A. Hoffmann, A. Kara, P. Krah, B. Legouix, E. Malaboeuf, et al., Gyselalib++: A portable c++ library for semi-lagrangian kinetic and gyrokinetic simulations, Journal of Open Source Software 10 (113) (2025) 8582.

[45] X. Garbet, M. Lesur, Gyrokinetics, lecture (Feb. 2023).
URL https://hal.science/hal-03974985

[46] A. Di Siena, Implementation and investigation of the impact of different background distributions in gyrokinetic plasma turbulence studies, Ph.D. thesis, Ulm University (Germany). Fakultät für Naturwissenschaften (2019).

[47] Y. Sarazin, V. Grandgirard, J. Abiteboul, S. Allfrey, X. Garbet, P. Ghendrih, G. Latu, A. Strugarek, G. Dif-Pradalier, Large scale dynamics in flux driven gyrokinetic turbulence, Nuclear Fusion 50 (5) (2010) 054004.