# An adaptive quasi-neutrality solver for full-F flux-driven gyrokinetic simulations of tokamak plasmas in presence of poloidal asymmetries

Peter Donnel, Kevin Obrejan, Yanick Sarazin, Roméo Bigué, Emily Bourne, Guillaume Brochard, Ludovica de Gianni, Guilhem Dif-Pradalier, Xavier Garbet, Virginie Grandgirard, et al.

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

## To cite this version:

Peter Donnel, Kevin Obrejan, Yanick Sarazin, Roméo Bigué, Emily Bourne, et al. An adaptive quasi-neutrality solver for full-F flux-driven gyrokinetic simulations of tokamak plasmas in presence of poloidal asymmetries. Open Plasma Science, In press, Volume 1, 10.46298/ops.15867. hal-05111064v3

HAL Id: hal-05111064

https://hal.science/hal-05111064v3

Submitted on 13 Jan 2026

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>
**REGULAR ARTICLE**

**History**

Received June 12, 2025

Revised November 27, 2025

Accepted December 9, 2025

Published January 5, 2026

**Identifiers**

DOI 10.46298/ops.15867

HAL hal-05111064

ArXiv -

Zenodo -

**Supplementary Material**

**Licence**

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

©The Authors

# An adaptive quasi-neutrality solver for full-F flux-driven gyrokinetic simulations of tokamak plasmas in presence of poloidal asymmetries

Peter Donnel<sup>1,*</sup>, Kevin Obrejan<sup>1</sup>, Yanick Sarazin<sup>1</sup>, Roméo Bigué<sup>1,2</sup>, Emily Bourne<sup>3</sup>, Guillaume Brochard<sup>1</sup>, Ludovica De Gianni<sup>1</sup>, Guilhem Dif-Pradalier<sup>1</sup>, Xavier Garbet<sup>1,4</sup>, Virginie Grandgirard<sup>1</sup>, Philipp Krah<sup>1</sup>, Yann Munschy<sup>1</sup>, and Matthieu Protais<sup>1</sup>

**Abstract**

Gyrokinetic codes are used to simulate transport in tokamak plasmas. In these codes, the distribution functions evolve self consistently with an electromagnetic field. To compute the temporal evolution of the electrostatic potential, a quasi-neutrality equation is solved. In some gyrokinetic codes, the quasi-neutrality solver assumes that the background densities and temperatures are constant in time and on flux surfaces. This assumption, which implicitly uses the so-called $\delta F$ approach, can break up, in particular at the edge of the plasma which can display large and time evolving poloidal asymmetries.

In this paper, a numerical solver of the quasi-neutrality equation accounting for time evolving poloidal asymmetries is presented. This solver is compatible with all electron models (adiabatic, kinetic or hybrid) and written for the long wavelength or the Padé approximations for the polarisation term. The impact of such an improvement is carefully reported on different types of simulations, illustrating when the $\delta F$ approach for quasi-neutrality is valid and when it fails. A procedure to limit the numerical cost of updating the background profiles in the quasi-neutrality solver is also presented.

**Keywords**— quasi-neutrality, gyrokinetic, tokamak

<footnote n="1"><sup>1</sup>CEA, IRFM, F-13108 Saint Paul-lez-Durance, France.</footnote>

<footnote n="2"><sup>2</sup>PIIM Laboratory, UMR 6633 CNRS-University of Provence, Marseille.</footnote>

<footnote n="3"><sup>3</sup>SCITAS, EPFL, CH-1015 Lausanne, Switzerland.</footnote>

<footnote n="4"><sup>4</sup>School of Physical and Mathematical Sciences, Nanyang Technological University, 637371 Singapore.</footnote>

<footnote>* Corresponding author: peter.donnel@cea.fr</footnote>

<footnote>Cite as: Donnel *et al.*, An adaptive quasi-neutrality solver for full-F flux-driven gyrokinetic simulations of tokamak plasmas in presence of poloidal asymmetries, *Open Plasma Science* 1, 5 (2025), doi: 10.46298/ops.15867</footnote>
# Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Models for the quasi-neutrality equation</td><td>3</td></tr><tr><td>2.1</td><td>Treatment of the polarisation term</td><td>4</td></tr><tr><td>2.2</td><td>Different electron models</td><td>5</td></tr><tr><td>2.3</td><td>Impact of the magnetic geometry</td><td>6</td></tr><tr><td>2.4</td><td>Modification of the quasi-neutrality in the presence of a limiter</td><td>6</td></tr><tr><td>3</td><td>Numerical resolution of the quasi-neutrality equation</td><td>7</td></tr><tr><td>3.1</td><td>Time evolution of the parameters of the QN equation</td><td>8</td></tr><tr><td>3.2</td><td>Treatment of the zonal mode in the adiabatic electron model</td><td>8</td></tr><tr><td>3.2.1</td><td>Fixed-point method in the long wavelength approximation</td><td>9</td></tr><tr><td>3.2.2</td><td>Fixed-point method in the Padé approximation</td><td>9</td></tr><tr><td>4</td><td>Numerical tests</td><td>10</td></tr><tr><td>4.1</td><td>Geodesic Acoustic Modes (GAMs) dynamics</td><td>10</td></tr><tr><td>4.1.1</td><td>Adiabatic electrons</td><td>10</td></tr><tr><td>4.1.2</td><td>Full kinetic electrons</td><td>13</td></tr><tr><td>4.2</td><td>Impact of updating coefficients in the presence of a limiter and kinetic electrons</td><td>15</td></tr><tr><td>4.2.1</td><td>Simulation with a uniform boundary condition</td><td>15</td></tr><tr><td>4.2.2</td><td>Simulation with a limiter</td><td>18</td></tr><tr><td colspan="2">Acknowledgments</td><td>21</td></tr><tr><td>A</td><td>Convergence of the fixed-point method for the zonal mode in the adiabatic electron model</td><td>21</td></tr><tr><td>A.1</td><td>Diagonalisation of the differential operator</td><td>21</td></tr><tr><td>A.2</td><td>Convergence properties</td><td>22</td></tr><tr><td colspan="2">References</td><td>25</td></tr></table>

# 1 Introduction

Tokamak plasmas are subject to turbulent transport. The best approach to simulate this turbulence is to use gyrokinetic codes [GIVW10] that solve the temporal evolution of distribution functions for the species of the plasma. These distribution functions evolve self-consistently with the electric potential (and the potential vector for electromagnetic simulations). To compute the electric potential, one needs to solve the quasi-neutrality equation which corresponds to the limit of the Poisson equation when considering spatial scales larger than the Debye length.

The polarisation term, which corresponds to the difference between the particle and gyrocenter densities, is almost always computed assuming the distribution function is close to a Maxwellian distribution function, i.e., $F_s = F_{Ms} + \delta F_s$ where $F_{Ms}$ corresponds to a Maxwellian distribution. In $\delta F$ gyrokinetic codes, this Maxwellian background is supposed to be fixed in time and the associated density and temperature are flux surface functions. This is a consistent choice as this background is not evolved at all in the course of the simulations. On the other hand, this assumption is a priori false in a full-F code which evolves the full distribution function. Nevertheless, this hypothesis is sometimes made even in full-F codes because of the numerical
cost to update the coefficients of the quasi-neutrality solver. This assumption of a constant polarization co-efficient is quite good for core simulations as long as the total simulation time is short compared with the particle confinement time. It explains why this assumption is often found in gyrokinetic codes that are used to study such plasmas, like GKNET [Obr17, IKI24], GT5D, and up to recently GYSELA [GAB$^{+}$16] and ORB5 [LOT$^{+}$20]. On the other hand, large poloidal asymmetries have been observed experimentally [TCL$^{+}$14] and numerically [CTL$^{+}$15, CCC$^{+}$17] in the edge. This is a consequence of the poloidally localized plasma/wall interaction, which is independent of the exact configuration (limiter or divertor). It is therefore necessary for the code to be able to handle such large poloidal asymmetries of the background distribution function when the simulation domain includes the Scrape-Off Layer (SOL). This is for instance the case in XGC that uses an adaptive background scheme [HKS$^{+}$22] or GENE-X [MSU$^{+}$21] (private communication). ORB5 has recently been upgraded to update the polarization density [MVB$^{+}$24] while still assuming this density is a flux function.

In this paper, we report on the new capability of the GYSELA code [GAB$^{+}$16] to treat time evolving poloidal asymmetries in the quasi-neutrality equation. This background distribution function is still assumed to be axisymmetric to lighten the numerical cost of the quasi-neutrality solver. This assumption is valid as long as the equilibrium magnetic field is axisymmetric but might fail for instance in presence of ripple, resonant magnetic perturbations or for stellarator configurations.

The outline of this paper is as follows. In section 2, the different quasi-neutrality models implemented in GYSELA are presented. They depend on the chosen electron response (adiabatic, hybrid or kinetic), the approximation of the polarisation density (long wavelength or Padé approximation) as well as the presence or not of a limiter. Section 3 describes how the quasi-neutrality equation is numerically solved in GYSELA, with a particular focus on how to treat the adiabatic electron response. In section 4, the quasi-neutrality solver is successfully verified on the dynamics of geodesic acoustic modes. Then the impact of updating the coefficients of the quasi-neutrality is shown on different typical simulations, highlighting the necessity to take into account poloidal asymmetries when the SOL is simulated with a kinetic electron response.

## 2 Models for the quasi-neutrality equation

The gyrokinetic Poisson equation reads [BH07, GIVW10]:

$$ -\nabla^2 \phi = \frac{1}{\epsilon_0} \sum_s e_s n_s = \frac{1}{\epsilon_0} \sum_s e_s (\bar{n}_s + n_{pol,s}), \quad (1) $$

where $\phi$ is the electric potential, $e_s$ the charge of the species $s$. We split the particle density $n_s$ of species $s$ into the polarisation density $n_{pol,s}$ and the gyrocenter density $\bar{n}_s$ defined as:

$$ \bar{n}_s(\mathbf{x}, t) = \int_0^\infty d\mu \int_{-\infty}^\infty dv_\parallel \mathcal{J}_{v,s} \bar{f}_s(\mathbf{x}, v_\parallel, \mu, t), \quad (2) $$

where $\bar{f}_s$ is the gyrocenter distribution function of species $s$ and $\mathcal{J}_{v,s} = 2\pi \frac{B_{\parallel,s}^*}{m_s}$ is the velocity Jacobian with $B_{\parallel,s}^* = B + \frac{m_s}{e_s} v_{\parallel} \mathbf{b} \cdot (\nabla \times \mathbf{b})$ which corresponds to the volume element in the guiding-center velocity space.

When considering spatial scales larger than the Debye length, deviations to the quasi-neutrality are screened by the plasma. Hence, the gyrokinetic Poisson equation reduces to the quasi-neutrality equation in this limit

$$ -\sum_s e_s n_{pol,s} = \sum_s e_s \bar{n}_s \quad (3) $$
In the rest of this paper, we consider that the initial distribution functions are chosen such that the polarisation densities are equal to zero, i.e., $\sum_s e_s \bar{n}_s (t=0) = 0$ and $\phi (t=0) = 0$. Therefore, one has:

$$-\sum_s e_s \delta n_{pol,s} = \sum_s e_s \delta \bar{n}_s, \qquad (4)$$

where the sign $\delta$ is used for deviations with respect to the initial state.

## 2.1 Treatment of the polarisation term

The polarisation term reads at leading order [DMB$^{+}$17]:

$$n_{pol,s}(\mathbf{x}, t) = \int J_Z d\mathbf{Z} \delta(\mathbf{X} + \mathbf{\rho}_s - \mathbf{x}) \frac{e_s}{B(\mathbf{x})} \tilde{\phi}(\mathbf{X}, \mu, t) \frac{\partial \bar{f}_s}{\partial \mu}(\mathbf{X}, v_{\parallel}, \mu, t), \qquad (5)$$

where $\mathbf{Z}$ represents the 6D phase space, $J_Z$ is the associated Jacobian, $\mathbf{X}$ is the gyro-center position, $\mathbf{\rho}_s$ the Larmor radius, $B$ is the background magnetic field supposed constant in time in this paper, and $\tilde{\phi} = \phi - \langle \phi \rangle_c$ with $\langle \phi \rangle_c$ being the gyro-averaged electrostatic potential. Note that compared to [DMB$^{+}$17], the time evolution of the distribution function is kept to be consistent with a full-f approach.

There are two approximate versions of the polarization term in the limit where the thermal Larmor radius is small compared to the typical transverse wavelengths of turbulence, i.e., $k_{\perp} \rho_s \ll 1$. The simplest one, often referred to as the long wavelength approximation reads:

$$e_s \delta n_{pol,s}^{LWA}(\mathbf{X}, t) = \nabla_{\perp} \cdot (\alpha_s \nabla_{\perp} \phi), \qquad (6)$$

with

$$\alpha_s = \frac{m_s \langle \bar{n}_s \rangle_{\varphi}}{\langle B \rangle_{\varphi}^2}. \qquad (7)$$

Two important remarks need to be raised at this point. On the one hand, the polarization term is proportional to the mass of the species considered. For this reason, the polarization term is neglected for electrons. The second remark is that only the axisymmetric components of the gyrocenter density and of the magnetic field are kept in the expression of $\alpha_s$. Therefore, $\alpha_s$ is a 2D function and not a 3D function. This is a key approximation to reduce the numerical cost of the quasi-neutrality solver. Of course this approximation breaks down if the magnetic field is not axisymmetric or if large toroidal variations of the density exist.

A better approximate version of the polarization term is called the Padé approximation and reads:

$$e_s \delta n_{pol,s}^{Pade}(\mathbf{X}, t) = [1 - \nabla_{\perp} \cdot \rho_s^2 \nabla_{\perp}]^{-1} [\nabla_{\perp} \cdot (\alpha_s \nabla_{\perp} \phi)]. \qquad (8)$$

The numerical evaluation of the Padé approximation of the polarization term is not straightforward. The problem is solved by multiplying the quasi-neutrality equation (4) by $\prod_i [1 - \nabla_{\perp} \cdot \rho_i^2 \nabla_{\perp}]$ as done in [LDB$^{+}$16]. Note that $i$ is referring to ions only as the Larmor radius of electrons is neglected.

$$-\sum_i \prod_{j \neq i} [1 - \nabla_{\perp} \cdot \rho_j^2 \nabla_{\perp}] [\nabla_{\perp} \cdot (\alpha_i \nabla_{\perp} \phi)] = \prod_i [1 - \nabla_{\perp} \cdot \rho_i^2 \nabla_{\perp}] \left( \sum_s e_s \delta \bar{n}_s \right), \qquad (9)$$

where $i$ and $j$ are referring to ion species. Remembering that we work in the limit $k_{\perp} \rho_i \ll 1$, we have, by keeping only the correction terms up to the second order:

$$-\nabla_{\perp} \cdot (\alpha \nabla_{\perp} \phi) = [1 - \nabla_{\perp} \cdot \kappa \nabla_{\perp}] \left( \sum_s e_s \delta \bar{n}_s \right) = \mathcal{P} \left( \sum_s e_s \delta \bar{n}_s \right), \qquad (10)$$
where $\alpha = \sum_i \alpha_i$, $\kappa = \sum_i \rho_i^2 = \sum_i \frac{m_i \langle T_i \rangle_{\phi}}{e^2 \langle B \rangle_{\phi}^2}$ and the Padé operator $\mathcal{P} = [1 - \nabla_{\perp} \cdot \kappa \nabla_{\perp}]$.

Let us note that higher order, non-linear terms exist in the polarization term. Their impact have been assessed in the long wavelength approximation for the case of Ion Temperature Gradient (ITG) turbulence [Ido12]. The corrections being small, these non-linear terms are neglected in this work.

## 2.2 Different electron models

Three electron models are considered in this work: the adiabatic electron (AE), the fully kinetic electron (FKE) and the hybrid kinetic electron (HKE) models. An important remark is that we use the drift-kinetic approximation for electrons. It means that we neglect the size of the electron Larmor radius. This choice is consistent with the fact that the electron polarisation term has been neglected in the previous section.

The AE model assumes that electrons are always at equilibrium with electrostatic potential fluctuations, hence there is no need to compute the evolution of the electron distribution function, limiting the numerical cost. The quasi-neutrality equation with AE reads:

$$ -\nabla_{\perp} \cdot (\alpha \nabla_{\perp} \phi^{\text{AE}}) + \mathcal{P} [\beta^{\text{AE}} (\phi^{\text{AE}} - \langle \phi^{\text{AE}} \rangle_{\text{FS}})] = \mathcal{P} \left( \sum_i e_i \delta \bar{n}_i \right), \quad (11) $$

with $\beta^{\text{AE}} = e^2 \frac{\langle \bar{n}_e \rangle_{\text{FS}}}{\langle T_e \rangle_{\text{FS}}}$ . Another advantage of the AE model is that, since the electron dynamics is not resolved, it allows to take large time steps constrained by the ion dynamics only, reducing furthermore the numerical cost of simulations. But the AE model discards all the kinetic physics associated with electrons, limiting its validity domain.

The FKE model does not make approximations on the electron response and is the simplest one when considering the quasi-neutrality problem. Indeed for the FKE model, the quasi-neutrality equation reads:

$$ -\nabla_{\perp} \cdot (\alpha \nabla_{\perp} \phi^{\text{FKE}}) = \mathcal{P} \left( \sum_i e_i \delta \bar{n}_i - e \delta \bar{n}_e \right). \quad (12) $$

A drawback of the FKE model is that it often requires a small time step and the inclusion of electromagnetic effects, to avoid the so-called $\omega_H$ mode which is a high frequency mode corresponding to the electrostatic limit of the kinetic Alfven wave [Ido16]. For these reasons, FKE simulations are the most expensive numerically.

To limit the numerical cost of simulations while keeping a part of the electron physics, especially to include the trapped electron modes instability, a third model considering only trapped electrons to have a kinetic response while passing electrons have an adiabatic response (for non-axisymmetric modes). Two solutions were proposed for the treatment of axisymmetric modes. In the historical one [Ido16], all axisymmetric modes are computed using the FKE model, but only the flux surface average $\langle \phi^{\text{FKE}} \rangle_{\text{FS}}$ is kept in the end, so that

$$ \phi_{old}^{\text{HKE}} = \phi_{n \neq 0}^{\text{TKE}} + \langle \phi^{\text{FKE}} \rangle_{\text{FS}}, $$

where $\phi^{\text{TKE}}$ is computed by assuming passing electrons as adiabatic. This allows to ensure ambipolarity and conservation of toroidal momentum. But the dynamics of Geodesic Acoustic Modes (GAMs) is significantly affected by this model. This model has been upgraded in [LMB$^{+18}$] to improve the treatment of the GAM dynamic by keeping also the contribution of the modes $m \neq 0, n = 0$. These are computed by assuming an adiabatic response for passing electrons, hence

$$ \phi_{new}^{\text{HKE}} = \phi^{\text{TKE}} - \langle \phi^{\text{TKE}} \rangle_{\text{FS}} + \langle \phi^{\text{FKE}} \rangle_{\text{FS}}. $$

In GYSELA, we use this improved HKE model.
The HKE algorithm is done in three steps. The first step is to solve the FKE potential via eq. (12). Then compute its flux surface average which is put on the right hand side of eq. (13), corresponding to the case where the passing electrons are treated adiabatically while trapped electrons are treated kinetically:

$$-\nabla_{\perp} \cdot (\alpha \nabla_{\perp} \phi^{\text{TKE}}) + \mathcal{P} (\beta^{\text{TKE}} \phi^{\text{TKE}}) = \mathcal{P} \left( \sum_{i} e_i \delta \bar{n}_i - e \delta \bar{n}_{e, \text{trap.}} + \beta^{\text{TKE}} \langle \phi^{\text{FKE}} \rangle_{\text{FS}} \right), \quad (13)$$

with $\beta^{\text{TKE}} = e^2 \langle \bar{n}_{e, \text{pas.}} \rangle_{\text{FS}} / \langle T_e \rangle_{\text{FS}}$ where $\bar{n}_{e, \text{pas.}}$ (respectively $\bar{n}_{e, \text{trap.}}$) is the density of passing (respectively trapped) electrons.

Once eq. (13) solved, the final step for the HKE model consists in replacing the flux surface average of this equation by the one computed by the FKE model.

$$\phi^{\text{HKE}} = \phi^{\text{TKE}} - \langle \phi^{\text{TKE}} \rangle_{\text{FS}} + \langle \phi^{\text{FKE}} \rangle_{\text{FS}} \quad (14)$$

Note that by injecting eq. (12) in eq. (13), and defining $\delta\phi = \phi^{\text{TKE}} - \phi^{\text{FKE}}$ one gets:

$$-\nabla_{\perp} \cdot (\alpha \nabla_{\perp} \delta\phi) + \mathcal{P} (\beta^{\text{TKE}} \delta\phi) = \mathcal{P} [e \delta \bar{n}_{e, \text{pas.}} - \beta^{\text{TKE}} (\phi^{\text{FKE}} - \langle \phi^{\text{FKE}} \rangle_{\text{FS}})] . \quad (15)$$

The interpretation of eq. (15) is that the difference between $\phi^{\text{TKE}}$ and $\phi^{\text{FKE}}$ comes from the non-adiabatic response of passing electrons as expected.

## 2.3 Impact of the magnetic geometry

While the general expression of the quasi-neutrality equation is unchanged when passing from a circular geometry to an arbitrary geometry, the differential operators are modified. In this work, we assimilate the parallel direction with the toroidal direction, i.e., $e_{\parallel} = e_{\varphi}$. With this assumption, the differential operators read:

$$\nabla_{\perp} \phi = \frac{\partial \phi}{\partial x^j} g^{jk} e_k = \left( \frac{\partial \phi}{\partial r} g^{rr} + \frac{\partial \phi}{\partial \theta} g^{\theta r} \right) e_r + \left( \frac{\partial \phi}{\partial r} g^{r \theta} + \frac{\partial \phi}{\partial \theta} g^{\theta \theta} \right) e_{\theta} \quad (16)$$

and

$$\nabla \cdot \mathbf{X}_{\perp} = \frac{1}{|\mathcal{J}|} \frac{\partial}{\partial x^j} (|\mathcal{J}|X^j) = \frac{1}{|\mathcal{J}|} \frac{\partial}{\partial r} (|\mathcal{J}|X^r) + \frac{1}{|\mathcal{J}|} \frac{\partial}{\partial \theta} (|\mathcal{J}|X^{\theta}), \quad (17)$$

where $g_{ij}$ is the 3D metric tensor and $|\mathcal{J}| = \sqrt{|\det(g_{ij})|}$ the norm of the associated Jacobian. Using these expressions, it can be shown that:

$$\nabla_{\perp} \cdot (F \nabla_{\perp} G) = g^{rr} \frac{\partial F}{\partial r} \frac{\partial G}{\partial r} + g^{r\theta} \left( \frac{\partial F}{\partial r} \frac{\partial G}{\partial \theta} + \frac{\partial F}{\partial \theta} \frac{\partial G}{\partial r} \right) + g^{\theta\theta} \frac{\partial F}{\partial \theta} \frac{\partial G}{\partial \theta} \\ + \frac{F}{|\mathcal{J}|} \left( \frac{\partial g_{\theta\theta}}{\partial r} \frac{\partial G}{\partial r} - \frac{\partial g_{r\theta}}{\partial r} \frac{\partial G}{\partial \theta} - \frac{\partial g_{r\theta}}{\partial \theta} \frac{\partial G}{\partial r} + \frac{\partial g_{rr}}{\partial \theta} \frac{\partial G}{\partial \theta} \right) \\ + F \left( g^{rr} \frac{\partial^2 G}{\partial r^2} + 2g^{r\theta} \frac{\partial^2 G}{\partial r \partial \theta} + g^{\theta\theta} \frac{\partial^2 G}{\partial \theta^2} \right). \quad (18)$$

## 2.4 Modification of the quasi-neutrality in the presence of a limiter

In the presence of a limiter, the AE quasi-neutrality equation is modified to take into account the fact that in the scrape-off layer, the potential should scale as $\phi = \Lambda T_e$ where $\Lambda$ is defined as:

$$\Lambda = -\frac{1}{2} \ln \left[ 2\pi \frac{m_e}{m_i} \left( 1 + \frac{T_i}{T_e} \right) \right]. \quad (19)$$
The modified AE equation has been presented in [DPGS$^{+22}$]. It reads:

$$ -\nabla_{\perp} \cdot (\alpha \nabla_{\perp} \phi^{\text{AE}}) + \mathcal{P} [\beta^{\text{AE}} (\phi^{\text{AE}} - (1 - \mathcal{M}^{\text{SOL}}) \langle \phi^{\text{AE}} \rangle_{\text{FS}})] = \\ \mathcal{P} \left( \sum_{i} e_i \delta \bar{n}_i + \beta^{\text{AE}} [\Lambda (\mathcal{M}^{\text{SOL}} - \mathcal{M}^{\text{mat}}) (\langle T_e \rangle_{\text{FS}} - T_e^{\text{b.c.}}) + (\mathcal{M}^{\text{mat}} - \mathcal{M}^{\text{wall}}) \phi^{\text{bias}}] \right). \quad (20) $$

The different masks for a circular cross-section are represented in Fig. 1. $\mathcal{M}^{\text{SOL}}$ is a poloidally symmetric function that smoothly but rapidly transition from zero in the closed field lines area to one in the open field lines area. $\mathcal{M}^{\text{wall}}$ is another poloidally symmetric mask with the same features as $\mathcal{M}^{\text{SOL}}$ but with a transition at a higher minor radius. Finally, $\mathcal{M}^{\text{mat}}$ is a poloidally localised function that represents all solid boundaries, i.e., the limiter and the wall. Different shapes are available for these masks as shown in Fig. 1. The physical interpretation of all areas delimited by the different masks is summarized in Fig. 1. Here $T_e^{\text{b.c.}}$ refers to the electron temperature at the outer boundary of the simulation domain. $\phi^{\text{bias}}$ is added to be able to model a biased limiter.

When the electron distribution function is solved in the scrape-off layer, one needs to use the FKE model. Indeed, it has been realized that only the full kinetic electron model allows fulfilling the parallel force balance in the scrape-off layer (see section 5.6.1 of [Mun24]). For this reason, for the HKE model in the presence of a limiter, one needs to modify the computation of $\phi^{\text{TKE}}$, i.e., eq. (13) which becomes:

$$ -\nabla_{\perp} \cdot (\alpha \nabla_{\perp} \phi^{\text{TKE}}) + \mathcal{P} [(1 - \mathcal{M}^{\text{SOL}}) \beta^{\text{TKE}} \phi^{\text{TKE}}] = \mathcal{P} \left( \sum_{i} e_i \delta \bar{n}_i \right) \\ - e \mathcal{P} [(1 - \mathcal{M}^{\text{SOL}}) \delta \bar{n}_{e, \text{trap.}} + \mathcal{M}^{\text{SOL}} \delta \bar{n}_e] \\ + \mathcal{P} [(1 - \mathcal{M}^{\text{SOL}}) \beta^{\text{TKE}} \langle \phi^{\text{FKE}} \rangle_{\text{FS}}]. \quad (21) $$

Note that with this model, the presence of the limiter is not taken into account in the quasi-neutrality equation. This modified HKE model allows running simulations with a simplified limiter which uses a Krook term to penalize both ion and electron distribution functions (see section 4.2). Our aim is to eventually change this equation to take into account 'sheath' type boundary conditions. This is left for a future work.

## 3 Numerical resolution of the quasi-neutrality equation

For its main computation, the spline-based solver uses a linear algebra library implemented in Selalib [BLK$^{+23}$] which has been upgraded to handle general 2D equations of the form:

$$ -\nabla_{\perp} \cdot (F_1(r, \theta) \nabla_{\perp} \phi) - \nabla_{\perp} \cdot (V(r, \theta) \phi) + F_2(r, \theta) \phi = \text{RHS}(r, \theta), \quad (22) $$

where the $F_i$ are smooth scalar fields and $V$ is a vector field on the poloidal plane. The solver assumes $\phi = 0$ on the edge of the domain. The differential operator described in section 2.3 is used.

The linear solver uses a finite element-like approach to solve the equation in a weak form using a conjugate gradient method. To optimize the numerical efficiency, the solution of the previous time step is used as an initial guess for the current time step. An important remark is that the solver assumes that the domain on which the solver is applied includes the pole ($r = 0$). For more details on this solver, see [BLK$^{+23}$].

In GYSELA, eq. (22) is solved independently for each poloidal cross section as fluctuations of the gyrocenter densities used in the RHS of the QN equation are a priori 3D functions. On the other hand, we use toroidally (respectively flux surface) averaged densities and temperatures for computing $\alpha, \kappa$ (respectively $\beta$) on the left hand side of the QN equation. Regarding $\alpha$ and $\kappa$, such averages are consistent with the present ordering in
<figure><img src="image_4.png" /><figcaption>Figure 1: Representation of the different masks present in the code.</figcaption></figure>

$\rho_*$ at which the equations are solved: only slowly evolving quantities (adiabatic) have to be considered in the expression of the polarisation density - apart from the electric potential itself. As far as $\beta$ is concerned, its form ensures that the adiabatic response of the electrons is vanishing when averaged on flux surfaces. These features have also the advantage to lighten the memory footprint of the quasi-neutrality solver as the matrix system to be solved is constructed and stored only once instead of $N_\phi$ times, where $N_\phi$ is the number of grid points in the toroidal direction.

## 3.1 Time evolution of the parameters of the QN equation

Let us define $\tau^{-1} = \partial \ln \|\langle \bar{X} \rangle_\varphi\| / \partial t$, then $\tau$ represents the typical evolution time for the evolution of the profiles. This is typically a transport time scale which is way larger than the time step used in gyrokinetic codes $\Delta t \ll \tau$. This time scale separation allows one to compute the matrix system to be solved in the QN equation not at each time step but instead every $N_{\text{update}}$ time steps, reducing dramatically the numerical cost as the construction of the matrix system is quite expensive. One has to make sure that $N_{\text{update}} \ll \tau / \Delta t$. In practice, for the case described for instance in section 4.2.1, updating the coefficient every 10 time steps adds an overhead of 10% to the computation time of the quasi-neutrality solver compared to never updating them. Since the quasi-neutrality only amounts to a small fraction of the overall simulation time, this approach is not very costly. In addition, the (re-)construction of the matrix scales as $O(N)$ while the solving of the system scales as $O(N^{3/2})$ where $N = N_r \cdot N_\theta$ is the number of grid points in the poloidal plane. So the relative cost of this update of the coefficients decreases with the size of the simulated system.

## 3.2 Treatment of the zonal mode in the adiabatic electron model

For the FKE model eq. (12), the spline solver eq. (22) can directly be used to compute the electric potential without difficulty. For the HKE model, the flux surface average of the potential needs to be computed to get the adiabatic response of passing electrons. Fortunately, the zonal mode is obtained via the FKE model eq. (12) so that when computing $\phi^{\text{TKE}}$ in eq. (13) (or eq. (21) in the presence of a limiter), the flux surface
average appears in the RHS of the equation and is already computed. So the procedure proposed in section 2.2 allows one to use directly the spline solver eq. (22).

The AE model eq. (11) (or eq. (20) in the presence of a limiter) is more difficult to treat with the spline solver as the flux surface average appears on the LHS of the equation. To compute the electric potential within this model, a fixed-point method is proposed. There is no proof of convergence of this method in the general case, but a proof can be derived in the specific case of a circular cross section, in the large aspect ratio limit and with constant coefficients α and β. This proof is given in appendix A. In the rest of this subsection, the fixed-point method algorithm is presented for the long wavelength approximation first and then generalised to the Padé version.

### 3.2.1 Fixed-point method in the long wavelength approximation

The quasi-neutrality equation in the presence of adiabatic electrons contains an extra term compared with the Selalib solver eq. (22). To circumvent this issue while also making use of the linear solver, a fixed-point scheme combining two equations is used to compute the toroidal average of the potential $\hat{\phi} = \int \phi^{AE} d\varphi/(2\pi)$ from which the flux surface averaged potential $\langle \phi^{AE} \rangle_{FS} = \langle \hat{\phi} \rangle_{FS}$ is computed and moved to the RHS of eq. (11) (or eq. (20) in the presence of a limiter). The equation for $\hat{\phi}$ is obtained by computing the toroidal average of eq. (20):

$$ -\nabla_{\perp} \cdot (\alpha \nabla_{\perp} \hat{\phi}) + [\beta^{AE} (\hat{\phi} - (1 - \mathcal{M}^{\text{SOL}}) \langle \hat{\phi} \rangle_{\text{FS}})] = \hat{\rho}, \quad (23) $$

where

$$ \hat{\rho} = \sum_{i} e_{i} \langle \delta \bar{n}_{i} \rangle_{\varphi} + \beta^{\text{AE}} \left[ \Lambda (\mathcal{M}^{\text{SOL}} - \mathcal{M}^{\text{mat}}) (\langle T_{e} \rangle_{\text{FS}} - T_{e}^{\text{b.c.}}) + (\mathcal{M}^{\text{mat}} - \mathcal{M}^{\text{wall}}) \phi^{\text{bias}} \right]. \quad (24) $$

The two forms of the equation which are used in the fixed-point method are:

$$ -\nabla_{\perp} \cdot (\alpha \nabla_{\perp} \hat{\phi}^{k+1}) + \beta^{\text{AE}} \mathcal{M}^{\text{SOL}} \hat{\phi}^{k+1} = \hat{\rho} - \beta^{\text{AE}} (\hat{\phi}^{k} - (1 - \mathcal{M}^{\text{SOL}}) \langle \hat{\phi}^{k} \rangle_{\text{FS}}), \quad (25) $$

which is efficient at converging low radial modes, and

$$ -\nabla_{\perp} \cdot (\alpha \nabla_{\perp} \hat{\phi}^{k+1}) + \beta^{\text{AE}} \hat{\phi}^{k+1} = \hat{\rho} + \beta^{\text{AE}} (1 - \mathcal{M}^{\text{SOL}}) \langle \hat{\phi}^{k} \rangle_{\text{FS}}, \quad (26) $$

which is efficient for high ones. Here, the notation $\hat{\phi}^k$ refers to the previous solution (coming from either initialisation or the previous iteration) and $\hat{\phi}^{k+1}$ is the updated solution.

The Selalib solver can be used to solve both eq. (25) and eq. (26) as the flux surface average term is always on the RHS and not part of the matrix which needs to be solved for. The electrostatic potential from the previous timestep is re-used as initial guess $\phi^0$ to initialise the procedure. After solving eq. (25) once, we iteratively solve eq. (26) until the convergence criterion

$$ \| \langle \phi^{k+1} \rangle_{\text{FS}} - \langle \phi^k \rangle_{\text{FS}} \| \le \text{tol} \cdot \| \langle \phi^{\text{new}} \rangle_{\text{FS}} \| \quad (27) $$

is fulfilled. In practice, this combined scheme only takes a few iterations to converge.

### 3.2.2 Fixed-point method in the Padé approximation

The method is similar to the one developed for the long wavelength approximation case but with more terms. The equation for $\hat{\phi}$ is obtained by computing the toroidal average of eq. (20):

$$ \begin{aligned} & -\nabla_{\perp} \cdot [(\alpha + \kappa\beta) \nabla_{\perp} \hat{\phi}] - \nabla_{\perp} \cdot (\kappa \hat{\phi} \nabla_{\perp} \beta) + \beta \hat{\phi} \\ &= \hat{\rho} - \nabla_{\perp} \cdot [\kappa \nabla_{\perp} \hat{\rho}] - \nabla_{\perp} \cdot [\kappa \nabla_{\perp} (\beta (1 - \mathcal{M}^{\text{SOL}}) \langle \hat{\phi} \rangle_{\text{FS}})] + \beta (1 - \mathcal{M}^{\text{SOL}}) \langle \hat{\phi} \rangle_{\text{FS}}. \end{aligned} \quad (28) $$
Since the general solver eq. (22) is compatible with the additional terms, a similar set of fixed-point methods can be used to pre-compute the terms in $\langle \phi \rangle_{\text{FS}}$. In this case, eq. (25) is replaced by:

$$ -\nabla_{\perp} \cdot \left[ (\alpha + \kappa\beta\mathcal{M}^{\text{SOL}}) \nabla_{\perp} \hat{\phi}^{n+1} \right] - \nabla_{\perp} \cdot \left\{ \kappa\hat{\phi}^{n+1}\nabla_{\perp} [\beta(1-\gamma)] \right\} + \beta(1-\gamma)\hat{\phi}^{n+1} \\ = \hat{\rho} - \nabla_{\perp} \cdot [\kappa\nabla_{\perp}\hat{\rho}] - \beta(\hat{\phi}^n - (1-\mathcal{M}^{\text{SOL}})\langle \hat{\phi}^n \rangle_{\text{FS}}) \quad (29) $$

and eq. (26) by:

$$ -\nabla_{\perp} \cdot \left[ (\alpha + \kappa\beta) \nabla_{\perp} \hat{\phi}^{n+1} \right] - \nabla_{\perp} \cdot \left[ \kappa \hat{\phi}^{n+1} \nabla_{\perp} \beta \right] + \beta \hat{\phi}^{n+1} = \hat{\rho} - \nabla_{\perp} \cdot [\kappa \nabla_{\perp} \hat{\rho}] + \beta (1 - \mathcal{M}^{\text{SOL}}) \langle \hat{\phi}^n \rangle_{\text{FS}}. \quad (30) $$

# 4 Numerical tests

## 4.1 Geodesic Acoustic Modes (GAMs) dynamics

GAMs are axisymmetric and oscillating modes of the electric potential in tokamak plasmas [CSI21]. Checking their properties (frequency, damping rate) is a standard test for gyrokinetic codes [BBE$^{+}$17], as theoretical predictions and benchmarks are well documented. We perform only a subset of the tests that exist in the literature to check that the Poisson solver is well implemented.

For all the simulations in this section, the aspect ratio is $A = 10$ and $\rho_{\star} = \rho_{i}/a = 1/160$ with the thermal Larmor radius defined as $\rho_{i} = \frac{v_{Ti}}{\omega_{ci}}$ where $\omega_{ci}$ is the cyclotron frequency. The simulation are done with flat density, temperature and safety factor profiles and without collisions. There is an initial perturbation of the form $\phi_{00}(t=0) = 1 - \cos(k_{r}r)$ with $k_{r} = 2\pi/a$, which corresponds to $k_{r}\rho_{i} = 0.039$ which is the same value as the one used in [BBE$^{+}$17, GGE$^{+}$19]. Note that there is an apparent factor $\sqrt{2}$ when comparing with [BBE$^{+}$17] which comes from the definition of the Larmor radius. We can therefore extract some results from [BBE$^{+}$17]. This is in particular the case for results coming from the codes ORB5 [LOT$^{+}$20] and GENE [GLB$^{+}$11]. For the tests performed in this section, we use 'Culham' equilibria [CCH$^{+}$88] without Shafranov shift. Note that this choice is different from the one used in ORB5 and GENE which both use numerical equilibria generated by CHEASE [LBS96]. It could explain the small but finite difference of results observed when performing the elongation scan Fig. 4, 5. The numerical discretization for the numerical tests in this section is $N_{r} = 64$, $N_{\theta} = 128$, $N_{\phi} = 8$, $N_{v_{\parallel}} = 512$ and $N_{\mu} = 64$. The velocity grid is uniform with $v_{\parallel,max} = 7$ and $\mu_{max} = 12$. The radial grid span the entire plasma volume $r \in [0, a]$.

### 4.1.1 Adiabatic electrons

There exists a vast literature on GAMs with adiabatic electrons. We will perform two tests. The first one is a safety factor scan for a circular concentric geometry. The second test is an elongation scan. The two tests are done in the case of a broad radial structure $k_{r}\rho_{i} \ll 1$.

For the scan in safety factor, the reference theoretical prediction is the one derived by Sugama [SW06, SW08]. The prediction of the GAM frequency reads:

$$ \frac{\omega_{\text{GAM}}^{\text{Sugama}}}{q\omega_{Ti}} = \sqrt{f_{T}\left[1 + \frac{f_{S1}}{q^{2}f_{T}^{2}}\right]}, \quad (31) $$

with $\omega_{Ti} = \frac{\sqrt{2v_{Ti}}}{qR_{0}}$ the thermal transit frequency of ions, $v_{Ti} = \sqrt{\frac{T_{i}}{m_{i}}}$ the thermal velocity of ions, $f_{T} = \frac{7}{4} + \tau_{e}$,
$f_{S1} = \frac{23}{8} + 2\tau_e + \frac{1}{2}\tau_e^2$ and $\tau_e = \frac{T_E}{T_i}$. The prediction for the damping rate reads:

$$ \frac{\gamma_i^{\text{Sugama}}}{q\omega_{Ti}} = -\frac{\sqrt{\pi}}{2} q^3 f_T \left[ \exp(-\hat{w}^2) (\hat{w}^2 + 2\tau_e + 1) + \frac{(qk_r\rho_i)^2}{2} \exp\left(-\frac{\hat{w}^2}{4}\right) \left(\frac{\hat{w}^4}{128} + f_{S2}\hat{w}^2 + f_{S3}\right) \right] \quad (32) $$

where $f_{S2} = \frac{1+\tau_e}{16}$, $f_{S3} = \frac{3}{8} + \frac{7\tau_e}{16} + \frac{5\tau_e}{32}$ and $\hat{w} = \omega_{\text{GAM}}/\omega_{ti}$. Note that there is an apparent factor $\sqrt{2}$ when comparing with [SW06, SW08] which comes from the different definition of the Larmor radius.

For the scan in elongation, the theoretical prediction of reference is the one derived by Gao [GPW$^{+}$09]. This prediction takes into account corrections due to the shear of the elongation $s_κ = \frac{r}{κ} \frac{\partial κ}{\partial r}$. For all the equilibria in this paper $s_κ < 2 \cdot 10^{-3}$. Therefore we take $s_κ = 0$ in the theoretical prediction for the sake of simplicity. The same choice was also done in [BBE$^{+}$17]. With this assumption, the GAM frequency is given by:

$$ \frac{\omega_{\text{GAM}}^{\text{Gao}}}{q\omega_{Ti}} = \sqrt{\frac{2f_T}{\kappa^2 + 1}} \left(1 + \frac{\kappa^2 + 1}{4} \frac{f_{S1}}{q^2 f_T^2}\right) \quad (33) $$

and the damping rate is given by:

$$ \frac{\gamma_i^{\text{Gao}}}{q\omega_{Ti}} = -\frac{\sqrt{\pi}}{2} \frac{\hat{w}^6}{qf_T} \exp(-\hat{w}^2) \quad (34) $$

Note that there is a typo in the original paper [GPW$^{+}$09] which was corrected in [BBE$^{+}$17]. An important point is that the damping rate prediction by Gao eq. (34) does not take into account corrections due to the radial structure of the mode while Sugama's prediction eq. (32) goes up to second order via the term proportional to $(qk_r\rho_i)^2$. Therefore it is expected that Sugama's prediction is more accurate for large safety factor when considering circular cross sections.

For the simulations with adiabatic electrons presented in this subsection, we use a time step $\Delta t = 25\omega_{ci}^{-1}$. In Fig. 2, the GAM frequency as a function of the safety factor given by GYSELA is compared to two theoretical predictions eq. (31),(33). A good match is found, especially with Gao's prediction at low safety factor that seems to fit better numerical results. In Fig. 3, the GAM damping rate as a function of the safety factor given by GYSELA is compared to the two theoretical predictions eq. (32),(34). Again, a good match is found with Sugama's prediction for safety factors below 3. Above this value, the numerical results display a higher damping rate as the one predicted theoretically. But this was already reported in Fig.1 of [BBE$^{+}$17] for GYSELA, ORB5 and GENE.

In Fig. 4, the GAM frequency as a function of the elongation is compared between GYSELA results, GENE and ORB5 results extracted from [BBE$^{+}$17] and the theoretical prediction by Gao, eq. (33). The trend is retrieved with GYSELA. The mismatch between the three codes increases when the elongation increases. We think this is the consequence of using a different magnetic geometry. Indeed ORB5 and GENE use magnetic equilibria generated by CHEASE, whereas GYSELA uses the 'Culham' model to generate its magnetic equilibria. Both models are expected to diverge when shaping increases. In Fig. 5, the GAM damping rate as a function of the elongation is compared between GYSELA results, GENE and ORB5 results extracted from [BBE$^{+}$17] and the theoretical prediction by Gao, eq. (34). Again the trend is found in good qualitative agreement between the three gyrokinetic codes, but the mismatch between the codes increases with shaping for the same reason.
<figure><img src="image_5.png" /><figcaption>Figure 2: GAM frequency as a function of the safety factor.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 3: GAM damping rate as a function of the safety factor.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 4: GAM frequency as a function of the elongation.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 5: GAM damping rate as a function of the elongation.</figcaption></figure>

4.1.2 Full kinetic electrons

On top of the ion contribution reported in the previous subsection, there is a non-negligible contribution of electrons to the damping rate of GAMs due to a resonance between the bounce/transit frequency of barely
trapped/passing electrons and the GAM real frequency. This contribution takes the form [EGG$^{+}$18]:

$$ \frac{\gamma_e}{q\omega_{Ti}} \approx -(0.315 + 0.30\epsilon) \left(1 + \frac{2 + \tau_e}{2\hat{w}^2}\right)^2 q\sqrt{\frac{\tau_e m_e}{2 m_i}} \mathcal{D}(\sigma_*, \epsilon) , \quad (35) $$

where $\sigma_*=\hat{\omega}\sqrt{2m_e/(m_i\epsilon\tau_e)}$ and $\mathcal{D}(\sigma_*, \epsilon)$ is a priori an unknown function. To evaluate the function $\mathcal{D}(\sigma_*, \epsilon)$ we perform a mass ratio scan at constant safety factor $q=3.5$ and temperature ratio $\tau_e=1.0$, a scan in temperature ratio $\tau_e$ at constant safety factor $q=2.0$ and mass ratio $m_i/m_e=1600$, a scan in safety factor at constant temperature ratio $\tau_e=1.0$ and mass ratio $m_i/m_e=1600$. The results of these scans are presented on Fig. 6. To evaluate the function $\mathcal{D}$, we use the following expression:

$$ \mathcal{D}_{GYS} = \frac{\gamma_{GYS} - \gamma_i^{\text{Sugama}}}{q\omega_{Ti}} \left[ (0.315 + 0.30\epsilon) \left(1 + \frac{2 + \tau_e}{2\hat{w}_{GYS}^2}\right)^2 q\sqrt{\frac{\tau_e m_e}{2 m_i}} \right]^{-1} . \quad (36) $$

<figure><img src="image_9.png" /><figcaption>Figure 6: Evaluation of $\mathcal{D}_{GYS}$, based on eq. (36). The data come from a mass ratio scan at $q=3.5$ and $\tau_e=1$ (red dots), a temperature ratio scan at constant $q=2.0$ and $m_i/m_e=1600$ (green dots) and a safety factor scan at constant $\tau_e=1$ and $m_i/m_e=1600$ (orange dots).</figcaption></figure>

Based on the various scans performed to obtain Fig. 6, we propose a polynomial fit of the form:

$$ \mathcal{D}_{GYS}^{\text{fit}}(\sigma_*) = 1.61\sigma_*^5 - 10.7\sigma_*^4 + 26.3\sigma_*^3 - 29.4\sigma_*^2 + 14.4\sigma_* - 0.388 \quad . \quad (37) $$

This fit is reported on Fig. 6. As one can see, this fit is not perfect, especially for the scan in $\tau_e$. This could come from the fact that the scan in $\tau_e$ has been performed at an intermediate safety factor value $q=2.0$, for which ion and electron contributions to GAM damping are of the same order of magnitude, thus reducing the accuracy of the method proposed to compute $\mathcal{D}$ with the formula eq. (36). We also removed the points with the lowest safety factor ($q<2.0$) in the safety factor scan for the same reasons. Despite this, we will use eq. (37) in eq. (35) for the safety factor scan represented in Fig. 7. Note that the average value of $\mathcal{D}$ in the
range of $\sigma_*$ tested is around 1.9. This value seems roughly twice larger than the one given in [GGE$^{+19}$]. But there is a factor 2 error in Eq. 19 of [GGE$^{+19}$] so that the two studies are consistent.

The safety factor scan are represented in Fig. 7. For this scan, a constant mass ratio $m_i/m_e = 1600$ and temperature ratio $\tau_e = 1$ were used. A good match is found with the theoretical prediction when summing the ion contribution, given by eq. (32), and the electron contribution, given by eq. (35), to compute the theoretical GAM damping rate prediction.

Note that the frequency of the GAM is expected to remain unchanged when adding electrons as reported numerically [BBE$^{+17}$, GGE$^{+19}$]. This property is recovered with this new solver (not shown here).

<figure><img src="image_10.png" /><figcaption>Figure 7: GAM damping rate as a function of the safety factor for the FKE model.</figcaption></figure>

All the tests performed in this section prove that the quasi-neutrality solver is correctly implemented.

## 4.2 Impact of updating coefficients in the presence of a limiter and kinetic electrons

In this section, we will show how the update of the coefficients $\alpha$ and $\beta$ in the quasi-neutrality equation impacts the results of the simulation. We will concentrate on two pairs of simulations with the hybrid electron model. Indeed this model is the most difficult in terms of quasi-neutrality equations and is rich enough to illustrate our purpose.

### 4.2.1 Simulation with a uniform boundary condition

In this section, we compare two simulations using the hybrid electron model eq. (12),(13),(14). A uniform Krook term is applied in the outer part of the simulation $\rho \ge 1$. To reduce the cost of the simulation, we use an artificially large mass ratio of $m_i/m_e = 100$. The thermal ion Larmor radius at mid radius $r/a = 0.5$ normalized to the minor radius of the tokamak is $\rho_*= 1/150$. The ion collisionality at mid radius is $\nu_* = 0.1$. We simulate a shaped plasma (see Fig. 10) to highlight the capacity of the code to simulate such magnetic geometries. The results are qualitatively the same with a circular geometry. The numerical resolution is
$N_r \times N_\theta \times N_\phi \times N_{v\parallel} \times N_\mu = 256 \times 256 \times 64 \times 128 \times 64$. We simulate half a torus and the radial domain ranges from the magnetic axis to $r_{max}/a = 1.2$. The time step is $\Delta t = 5\omega_{c0}^{-1}$, where $\omega_{c0}$ is the ion cyclotron frequency on the magnetic axis. The only difference between the two simulations is that the first one uses constant coefficients in the quasi-neutrality equation whereas the other one uses coefficients which are updated every 10 time steps, i.e., $\Delta t_{update} = 50\omega_{c0}^{-1}$. To reduce the length of the simulation, we perform a decaying turbulence simulation, i.e., without sources. Due to this choice, we expect a relatively rapid relaxation of temperature profiles. Despite being a bit symptomatic, this configuration can also represent a flux driven simulation for which the initial profiles are far from the steady state ones. The total simulation time is $t_{final} = 5 \cdot 10^4 \omega_{c0}^{-1}$, which is roughly five times what is required to reach the nonlinear regime in the simulation. This time is long enough to observe a significant evolution of profiles.

The ion density profiles at the beginning and the end of the simulations are reported on Fig. 8. The electron density profiles are qualitatively the same as expected from the quasi-neutrality condition. We see on this figure that the density did not evolve much and no significantly poloidal asymmetry is displayed. Therefore we expect that the polarisation coefficient $\alpha$ in the quasi-neutrality equation is more or less the same in both simulations.

<figure><img src="image_11.png" /><figcaption>Figure 8: Ion density profile in an hybrid electron simulation without limiter.</figcaption></figure>

The electron temperature profiles at the beginning and the end of the simulations are plotted on Fig. 9. As previously, no significant poloidal asymmetry is observed. On the other hand, the flux surface average profile evolved significantly due to the lack of any heat source that could have compensated for the energy transport. We therefore expect a significant evolution of the coefficient $\beta^{TKE}$ in Eq. (13) which corresponds to the adiabatic response of passing electrons. Note that the profiles at the end of the simulation are more or less the same in both simulations.

We now turn our attention to the non-axisymmetric components of the electric potential $\phi(R, Z, \varphi = 0) - \langle\phi(R, Z)\rangle_\varphi$ where $\langle\phi(R, Z)\rangle_\varphi$ is the toroidal average of the electric potential. This quantity at the end of the two simulations is plotted in Fig. 10. The turbulence intensity is 40% lower for the simulation with updated coefficients compared to the simulation without update. Due to this large difference, a smaller heat flux is expected when updating the coefficient (in this specific case only). Therefore, if we continue the simulation, we expect a difference in temperature profiles between the two simulations that will increase with time.

This pair of simulation illustrates why it is mandatory to update the coefficients in the quasi-neutrality equation for relatively long flux-driven simulations. As the background profiles are free to evolve, they can deviate significantly from initial profiles. The update of the quasi-neutrality coefficients then leads to signif-
<figure><img src="image_12.png" /><figcaption>Figure 9: Electron temperature profile in an hybrid electron simulation without limiter.</figcaption></figure>

icant changes in the level of turbulence predicted by the simulation, hence on the resulting transport.

<figure><img src="image_13.png" /><figcaption>Figure 10: Electric potential fluctuations in an hybrid electron simulation without limiter.</figcaption></figure>
### 4.2.2 Simulation with a limiter

In the previous section, we showed that updating the quasi-neutrality coefficients is mandatory for relative long flux driven simulations in the case of poloidally symmetric boundary conditions. In this section, we show that this update is far more critical when poloidally localized boundary conditions are used.

In the pair of simulations presented in this section an infinite penalization for both ions and electrons is used in a poloidally localized region depicted on Fig. 1, which models a limiter. To simplify the handling of the limiter, we chose a circular outer boundary. The target distribution function on which the distribution function is forced to relax in the limiter is a centered Maxwellian of temperature equal to the initial outer edge temperature, and the actual density. Therefore, this penalization acts as a momentum (in absolute value) and an energy sink but allows the density to evolve freely. Concerning the quasi-neutrality equation, we use the model proposed in section 2.4, i.e., an hybrid electron model in the core ($r < a$) and a full kinetic electron model in the Scrape-off layer ($r > a$). The actual set of equation used for computing the electric potential are eq. (12),(21),(14).

The main physical parameters are the same as in the previous section, i.e., $m_i/m_e = 100$, $\rho_*= 1/150$ and $\nu_* = 0.1$. The numerical resolution is $N_r \times N_\theta \times N_\varphi \times N_{v_\|} \times N_\mu = 512 \times 512 \times 64 \times 128 \times 64$. We simulate a full torus and the radial domain ranges from the magnetic axis to $r_{max}/a = 1.2$. Note that the radial and poloidal resolution are doubled compared to the simulation with homogeneous boundary condition. This is due to the necessity to discretize correctly the limiter. Indeed a pair of simulations with $N_r \times N_\theta \times N_\varphi \times N_{v_\|} \times N_\mu = 256 \times 256 \times 64 \times 128 \times 32$ was performed. Qualitatively, both simulations provide the same output, but the low resolution simulations display some features typical of a lack of spatial resolution close to the limiter. The time step is $\Delta t = 0.5\omega_{c0}^{-1}$, where $\omega_{c0}$ is the ion cyclotron frequency on the magnetic axis. This is an order of magnitude lower compared with a simulation without a limiter. The time step was reduced to avoid a numerical instability. We checked numerically that the larger the mass ratio $m_i/m_e$, the smaller the time step needs to be for a stable simulation. At this stage, this instability is not perfectly understood. Our hypothesis is that this instability is reminiscent of the $\omega_H$ mode, known to be detrimental for the time step when using the full kinetic electron model in electrostatic simulations. We will further investigate this instability in a future work. The only difference between the two simulations discussed in this section is the fact that the quasi-neutrality coefficients are updated (every 100 time steps, i.e., $\Delta t_{update} = 50\omega_{c0}^{-1}$) or not.

In Fig. 11, the axisymmetric ion density is displayed for different poloidal locations. The case without updating the quasi-neutrality coefficients, shown in Fig. 11a, displays a huge poloidal asymmetry at the transition between open and closed field lines (i.e., $r=a$). This asymmetry grows with time until the simulation crashes. On the other hand, the simulation with updated coefficients displays only a modest poloidal asymmetry at the same position. This pair of simulations illustrates the reason why updating the quasi-neutrality coefficients (with poloidal dependencies) is mandatory when simulating plasma with a poloidally localized boundary condition. Despite the ion density (with update) displaying only modest deviations to the initial profile, the update of quasi-neutrality coefficients indeed allows the system to reach an equilibrium in the scrape-off layer.

To highlight the fact that the simulation with updated quasi-neutrality coefficients is stable, we continued the simulation until the beginning of the nonlinear phase. In Fig. 12, the axisymmetric ion density is displayed for different poloidal locations at the beginning of the non-linear phase. Moderate poloidal asymmetries are visible in the core, but no large poloidal asymmetries are visible at the transition between open and closed field lines. In Fig. 13, the fluctuation of the electric potential at the same time are displayed. Turbulence is clearly visible in the core and the beginning of an instability is also appearing in the inner side of the scrape-off layer.
<figure><img src="image_14.png" /><figcaption>**Figure 11:** axisymmetric component of the electron density at a relatively early time of the simulation.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 12:** ion density</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 13: Electric potential fluctuations</figcaption></figure>

Conclusion

In this article we have presented a new and versatile quasi-neutrality solver for the global flux driven code GYSELA. It allows one to run simulations with adiabatic, full kinetic or hybrid electrons. We have proposed a new hybrid electron model allowing one to simulate at the same time core and scrape-off layer physics at a relatively moderate numerical cost. The solver can use either the long wavelength or the Padé approximation to compute the polarisation term. It can also handle shaped geometries.

The solver was successfully tested in a set of simulations of the GAM dynamics (frequency and damping rate) both with adiabatic electrons and full kinetic electrons (the two electron models having theoretical predictions to compare with).

Finally, we illustrated some configurations in which the update of the quasi-neutrality coefficients is mandatory. For simulations with poloidally symmetrical boundary conditions, a low update rate is enough as density and temperature profiles display moderate poloidal asymmetries and the radial profiles evolve slowly. For simulations with poloidally localized boundary conditions, the update is far more critical as it appears to stabilise the simulation in the scrape-off layer region.

This work can easily be extended to electromagnetic simulations as the Ampère equation actually takes a similar form as the quasi-neutrality equation. The electromagnetic version of this work will be presented in a forthcoming publication.
Acknowledgments

The authors thank J. Dominski and P. Ulbl for providing information on the quasi-neutrality treatment in XGC and GENE-X respectively. This project was provided with computer and storage resources by GENCI at TGCC and CINES thanks to the grant 2024-A0160502224 and 2025-A0180502224 on the supercomputer Joliot Curie's the SKL and ROME partitions as well as on the supercomputer Adastra's the GENOA partition. This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 - EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them. This work was supported by the EUROfusion Theory and Advanced Simulation Coordination (E-TASC) initiative under the TSVV-01 (Theory, Simulation, Verification and Validation) 'Physics of the L-H Transition and Pedestals' project.

A Convergence of the fixed-point method for the zonal mode in the adiabatic electron model

We do not provide a formal proof for the convergence of either (long wavelength or Padé approximation) fixed-point methods in the general case. Their formulation is based on considerations in a simplified case, where convergence can be proven, and was then tweaked for faster convergence in actual simulations. We consider the simpler equation

$$-\nabla_{\perp} \cdot (\alpha \nabla_{\perp} \phi) + \beta (\phi - \langle \phi \rangle_{\text{FS}}) = \rho, \quad (38)$$

where $\alpha$ and $\beta$ are strictly positive constant numbers and the metric is that of a circular geometry with a large aspect ratio so that

$$\nabla_{\perp} \cdot (\alpha \nabla_{\perp} \phi) = \alpha \left( \partial_{rr}^2 + \frac{1}{r} \partial_r + \partial_{\theta\theta}^2 \right) \phi \quad (39)$$

and $\langle \phi \rangle_{\text{FS}}$ is merely an averaging over the poloidal angle.

Denoting $\mathcal{L}_{\beta}(\phi) = -\nabla_{\perp} \cdot (\alpha \nabla_{\perp} \phi) + \beta \phi$, we can rewrite eq. (38) as either

$$\mathcal{L}_{\beta}(\phi) - \beta \langle \phi \rangle_{\text{FS}} = \rho \quad (40)$$

or

$$\mathcal{L}_0(\phi) + \beta (\phi - \langle \phi \rangle_{\text{FS}}) = \rho, \quad (41)$$

which are the basis for the fixed-point schemes.

A.1 Diagonalisation of the differential operator

A fixed-point scheme of the form $\phi^{k+1} = F(\phi^k)$ will converge if $F$ is a contraction mapping. When $F$ is an affine operator (composition of a linear operator and a translation), this means that all the eigenvalues of the linear part of $F$ must be strictly smaller than 1 in absolute value, and the convergence will be faster when those eigenvalues are closer to zero.

We thus start by determining the spectrum of the differential operator $\mathcal{L}_{\beta}$, i.e., finding $(\lambda, \phi_{\lambda})$ pairs such that $\mathcal{L}_{\beta}(\phi_{\lambda}) = \lambda \phi_{\lambda}$. Formally, this only corresponds to the point spectrum, while the diagonalisation of $\mathcal{L}_{\beta}$ would also require determining its continuous and residual spectra. However, the discretisation of space in the code reduces all this to a finite dimension case where only the point spectrum remains.
In our simplified setting, we can use a Fourier transform in $\theta$ to reduce the relation $\mathcal{L}_{\beta}(\phi) = \lambda\phi$ to a set of 1D differential equations:

$$ \phi_m'' + \frac{1}{r}\phi_m' + \left(\frac{\lambda - \beta}{\alpha} - \frac{m^2}{r^2}\right) \phi_m = 0, \quad (42) $$

where the prime denotes a radial derivative and $m$ is the poloidal mode number. The solutions of these equations are the Bessel functions $J_m$ and $Y_m$. Excluding $Y_m$ which diverges at $r = 0$, the solutions are thus of the form $J_m(\omega \cdot r)$ with $\omega = \sqrt{(\lambda - \beta)/\alpha}$. The Dirichlet boundary condition at $r = a$ implies $J_m(\omega a) = 0$, i.e.,

$$ \lambda = \beta + \alpha \left( \frac{Z_p^m}{a} \right)^2, \quad (43) $$

where $Z_p^m$ is the $p$-th zero of $J_m$, which is strictly positive. Thus low values of $p$ correspond to low radial modes, which oscillate only a few times over the radial domain while higher values of $p$ oscillates many times.

## A.2 Convergence properties

In addition to the eigenvalues of $\mathcal{L}$, the flux-surface averaging operator will also contribute by filtering some of the poloidal modes as $\langle \phi_m \rangle_{\text{FS}} = \delta_m \phi_m$ in our simplified setting. The first fixed-point scheme, based on eq. (40) is obtained by inverting the operator

$$ \phi^{k+1} = \mathcal{L}_{\beta}^{-1} \left( \rho + \beta \left\langle \phi^k \right\rangle_{\text{FS}} \right) \quad (44) $$

whose linear spectrum will be zero for $m \neq 0$ (those modes are simply set to a fixed value) and

$$ \frac{\beta}{\beta + \alpha (Z_p^0/a)^2} \quad (45) $$

for $m = 0$. Since these eigenvalues are always between 0 and 1, the scheme will converge, although convergence might be slow for low values of $p$ for which $\alpha (Z_p^0/a)^2 \sim (Z_p^0 \rho_*)^2 \ll \beta$ where $\rho_*$ is the ion gyroradius of reference normalised by the minor radius.

Similarly, the second scheme is based on eq. (41) after inverting the operator:

$$ \phi^{k+1} = \mathcal{L}_0^{-1} \left( \rho + \beta \left( \left\langle \phi^k \right\rangle_{\text{FS}} - \phi^k \right) \right) \quad (46) $$

whose spectrum is zero for $m = 0$ but

$$ \frac{\beta}{\alpha (Z_p^m/a)^2} = \frac{\beta}{\alpha} \left( \frac{a}{Z_p^m} \right)^2 \sim \frac{1}{(Z_p^m \rho_*)^2} \quad (47) $$

for $m \neq 0$, i.e., potentially diverging for low modes. For this reason, this second scheme is only used for a single iteration, which allows priming the modes in $m = 0$ before applying the first fixed-point scheme eq.(44).
References

[BBE<sup>+17</sup>] A. BIANCALANI, A. BOTTINO, C. EHRLACHER, V. GRANDGIRARD, G. MERLO, I. NOVIKAU, Z. QIU, E. SONNENDRÜCKER, X. GARBET, T. GÖRLER, S. LEERINK, F. PALERMO, and D. ZARZOSO, Cross-code gyrokinetic verification and benchmark on the linear collisionless dynamics of the geodesic acoustic mode, *Physics of Plasmas* **24** no. 6 (2017), 062512. doi:10.1063/1.4985571.

[BLK<sup>+23</sup>] E. BOURNE, P. LELEUX, K. KORMANN, C. KRUSE, V. GRANDGIRARD, Y. GÜÇLÜ, M. J. KÜHN, U. RÜDE, E. SONNENDRÜCKER, and 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.

[BH07] A. J. BRIZARD and T. S. HAHM, Foundations of nonlinear gyrokinetic theory, *Rev. Mod. Phys.* **79** no. 2 (2007), 421–468. doi:10.1103/RevModPhys.79.421.

[CTL<sup>+15</sup>] R. M. CHURCHILL, C. THEILER, B. LIPSCHULTZ, I. H. HUTCHINSON, M. L. REINKE, D. WHYTE, J. W. HUGHES, P. CATTO, M. LANDREMAN, D. ERNST, C. S. CHANG, R. HAGER, A. HUBBARD, P. ENNEVER, J. R. WALK, and A. C.-M. TEAM, Poloidal asymmetries in edge transport barriers, *Physics of Plasmas* **22** no. 5 (2015), 056104. doi:10.1063/1.4918353.

[CCC<sup>+17</sup>] R. CHURCHILL, J. CANIK, C. CHANG, R. HAGER, A. LEONARD, R. MAINGI, R. NAZIKIAN, and D. STOTLER, Kinetic simulations of scrape-off layer physics in the DIII-D tokamak, *Nuclear Materials and Energy* **12** (2017), 978–983, Proceedings of the 22nd International Conference on Plasma Surface Interactions 2016, 22nd PSI. doi:10.1016/j.nme.2016.12.013.

[CCH<sup>+88</sup>] J. W. CONNOR, S. C. COWLEY, R. J. HASTIE, T. C. HENDER, A. HOOD, and T. J. MARTIN, Tearing modes in toroidal geometry, *The Physics of Fluids* **31** no. 3 (1988), 577–590. doi:10.1063/1.866840.

[CSI21] G. CONWAY, A. SMOLYAKOV, and T. IDO, Geodesic acoustic modes in magnetic confinement devices, *Nuclear Fusion* **62** no. 1 (2021), 013001. doi:10.1088/1741-4326/ac0dd1.

[DPGS<sup>+22</sup>] G. DIF-PRADALIER, P. GHENDRIH, Y. SARAZIN, E. CASCHERA, F. CLAIRET, Y. CAMENEN, P. DONNEL, X. GARBET, V. GRANDGIRARD, Y. MUNSCHY, L. VERMARE, and F. WIDMER, Transport barrier onset and edge turbulence shortfall in fusion plasmas, *Communications Physics* **5** (2022). doi:10.1038/s42005-022-01004-z.

[DMB<sup>+17</sup>] J. DOMINSKI, B. F. MCMILLAN, S. BRUNNER, G. MERLO, T.-M. TRAN, and L. VILLARD, An arbitrary wavelength solver for global gyrokinetic simulations. Application to the study of fine radial structures on microturbulence due to non-adiabatic passing electron dynamics, *Physics of Plasmas* **24** no. 2 (2017), 022308. doi:10.1063/1.4976120.

[EGG<sup>+18</sup>] C. EHRLACHER, X. GARBET, V. GRANDGIRARD, Y. SARAZIN, P. DONNEL, E. CASCHERA, P. GHENDRIH, and D. ZARZOSO, Contribution of kinetic electrons to GAM damping, *Journal of Physics: Conference Series* **1125** no. 1 (2018), 012010. doi:10.1088/1742-6596/1125/1/012010.

[GPW<sup>+09</sup>] Z. GAO, L. PENG, P. WANG, J. DONG, and H. SANUKI, Plasma elongation effects on temperature gradient driven instabilities and geodesic acoustic modes, *Nuclear Fusion* **49** no. 4 (2009), 045014. doi:10.1088/0029-5515/49/4/045014.

[GIVW10] X. GARBET, Y. IDOMURA, L. VILLARD, and T. WATANABE, Gyrokinetic simulations of turbulent transport, *Nuclear Fusion* **50** no. 4 (2010), 043002. doi:10.1088/0029-5515/50/4/043002.
[GAB<sup>+</sup>16] 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, and D. ZARZOSO, A 5D gyrokinetic full-f global semi-lagrangian code for flux-driven ion turbulence simulations, *Computer Physics Communications* **207** (2016), 35–68. doi:10.1016/j.cpc.2016.05.007.

[GGE<sup>+</sup>19] V. GRANDGIRARD, X. GARBET, C. EHRLACHER, A. BIANCALANI, A. BOTTINO, I. NOVIKAU, Y. ASAHI, E. CASCHERA, G. DIF-PRADALIER, P. DONNEL, P. GHENDRIH, C. GILLOT, G. LATU, C. PASSERON, Y. SARAZIN, and D. ZARZOSO, Linear collisionless dynamics of the GAM with kinetic electrons: Comparison simulations/theory, *Physics of Plasmas* **26** no. 12 (2019), 122304. doi:10.1063/1.5113679.

[GLB<sup>+</sup>11] T. GÖRLER, X. LAPILLONNE, S. BRUNNER, T. DANNERT, F. JENKO, F. MERZ, and D. TOLD, The global version of the gyrokinetic turbulence code GENE, *Journal of Computational Physics* **230** no. 18 (2011), 7053–7071. doi:10.1016/j.jcp.2011.05.034.

[HKS<sup>+</sup>22] R. HAGER, S. KU, A. Y. SHARMA, C. S. CHANG, R. M. CHURCHILL, and A. SCHEINBERG, Electromagnetic total-f algorithm for gyrokinetic particle-in-cell simulations of boundary plasma in XGC, *Physics of Plasmas* **29** no. 11 (2022), 112308. doi:10.1063/5.0097855.

[Ido16] Y. IDOMURA, A new hybrid kinetic electron model for full-f gyrokinetic simulations, *Journal of Computational Physics* **313** (2016), 511–531. doi:10.1016/j.jcp.2016.02.057.

[Ido12] Y. IDOMURA, Accuracy of momentum transport calculations in full-f gyrokinetic simulations, *Computational Science and Discovery* **5** no. 1 (2012), 014018. doi:10.1088/1749-4699/5/1/014018.

[IKI24] K. IMADERA, Y. KISHIMOTO, and A. ISHIZAWA, Turbulent particle pinch in gyrokinetic flux-driven ITG/TEM turbulence, *Nuclear Fusion* **64** no. 8 (2024), 086006. doi:10.1088/1741-4326/ad52a4.

[LDB<sup>+</sup>16] E. LANTI, J. DOMINSKI, S. BRUNNER, B. MCMILLAN, and L. VILLARD, Padé approximation of the adiabatic electron contribution to the gyrokinetic quasi-neutrality equation in the ORB5 code, *Journal of Physics: Conference Series* **775** (2016), 012006. doi:10.1088/1742-6596/775/1/012006.

[LMB<sup>+</sup>18] E. LANTI, B. F. MCMILLAN, S. BRUNNER, N. OHANA, and L. VILLARD, Gradient- and flux-driven global gyrokinetic simulations of ITG and TEM turbulence with an improved hybrid kinetic electron model, *Journal of Physics: Conference Series* **1125** no. 1 (2018), 012014. doi:10.1088/1742-6596/1125/1/012014.

[LOT<sup>+</sup>20] E. LANTI, N. OHANA, N. TRONKO, T. HAYWARD-SCHNEIDER, A. BOTTINO, B. MCMILLAN, A. MISHCHENKO, A. SCHEINBERG, A. BIANCALANI, P. ANGELINO, S. BRUNNER, J. DOMINSKI, P. DONNEL, C. GHELLER, R. HATZKY, A. JOCKSCH, S. JOLLIET, Z. LU, J. MARTIN COLLAR, I. NOVIKAU, E. SONNENDRÜCKER, T. VERNAY, and L. VILLARD, ORB5: A global electromagnetic gyrokinetic code using the PIC approach in toroidal geometry, *Computer Physics Communications* **251** (2020), 107072. doi:10.1016/j.cpc.2019.107072.

[LBS96] H. LÜTJENS, A. BONDESON, and O. SAUTER, The CHEASE code for toroidal MHD equilibria, *Computer Physics Communications* **97** no. 3 (1996), 219–260. doi:10.1016/0010-4655(96)00046-X.

[MSU<sup>+</sup>21] D. MICHELS, A. STEGMEIR, P. ULBL, D. JAREMA, and 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.
[Mun24] Y. MUNSCHY, *Kinetic and Gyrokinetic physics of plasma-wall interaction in tokamaks*, PhD thesis, Aix-Marseille Universite, November 2024.

[MVB<sup>+</sup>24] M. MURUGAPPAN, L. VILLARD, S. BRUNNER, G. DI GIANNATALE, B. F. MCMILLAN, and A. BOTTINO, Gyrokinetic flux-driven simulations in mixed TEM/ITG regime using a delta-f PIC scheme with evolving background, *Physics of Plasmas* **31** no. 11 (2024), 113901. doi:10.1063/5.0233727.

[Obr17] K. OBREJAN, *Study of magnetic shaping effects on plasma flows and micro-instabilities in tokamak plasmas using the full-f gyrokinetic code based on a real space field solver*, PhD thesis, Kyoto University, 2017.

[SW06] H. SUGAMA and T.-H. WATANABE, Collisionless damping of geodesic acoustic modes, *Journal of Plasma Physics* **72** no. 6 (2006), 825–828. doi:10.1017/S0022377806004958.

[SW08] H. SUGAMA and T.-H. WATANABE, Erratum: ‘Collisionless damping of geodesic acoustic modes’
[J. Plasma Physics (2006) 72, 825], *Journal of Plasma Physics* **74** no. 1 (2008), 139–140.
doi:10.1017/S002237780700668X.

[TCL<sup>+</sup>14] C. THEILER, R. CHURCHILL, B. LIPSCHULTZ, M. LANDREMAN, D. ERNST, J. HUGHES, P. CATTO, F. PARRA, I. HUTCHINSON, M. REINKE, A. HUBBARD, E. MARMAR, J. TERRY, J. WALK, and THE ALCATOR C-MOD TEAM, Inboard and outboard radial electric field wells in the H- and I-mode pedestal of Alcator C-Mod and poloidal variations of impurity temperature, *Nuclear Fusion* **54** no. 8 (2014), 083017. doi:10.1088/0029-5515/54/8/083017.