Moment-based approach to the flux-tube linear gyrokinetic model

Abstract

This work reports on the development and numerical implementation of the linear electromagnetic gyrokinetic (GK) model in a tokamak flux-tube geometry using a moment approach based on the expansion of the perturbed distribution function on a velocity-space Hermite–Laguerre polynomials basis. A hierarchy of equations of the expansion coefficients, referred to as the gyro-moments (GMs), is derived. We verify the numerical implementation of the GM hierarchy in the collisionless limit by performing a comparison with the continuum GK code GENE, recovering the linear properties of the ion temperature gradient, trapped electron, kinetic ballooning and microtearing modes, as well as the collisionless damping of zonal flows. An analysis of the distribution functions and ballooning eigenmode structures is performed. The present investigation reveals the ability of the GM approach to describe fine velocity-space-scale structures appearing near the trapped and passing boundary and kinetic effects associated with parallel and perpendicular particle drifts. In addition, the effects of collisions are studied using advanced collision operators, including the GK Coulomb collision operator. The main findings are that the number of GMs necessary for convergence decreases with plasma collisionality and is lower for pressure gradient-driven modes, such as in H-mode pedestal regions, compared with instabilities driven by trapped particles and magnetic gradient drifts often found in the core. The accuracy of approximations often used to model collisions (relative to the GK Coulomb operator) is studied in the case of trapped electron modes, showing differences between collision operator models that increase with collisionality and electron temperature gradient, consistent with the results of Pan et al. (Phys. Rev. E, vol. 103, 2021, L051202). Such differences are not observed in other edge microinstabilities, such as microtearing modes. The importance of a proper collision operator model is also confirmed by analysing the collisional damping of geodesic acoustic modes and zonal flows.

1. Introduction

Linear and nonlinear gyrokinetic (GK) simulations are the tools of reference for describing low-frequency (compared with the ion gyrofrequency, $\varOmega _i$) and small scales (of the order of the ion gyroradius, $\rho _i$) electromagnetic microinstabilities occurring in the core and edge regions of fusion devices (Told et al. 2008; Holland et al. 2011; Navarro et al. 2015; Kotschenreuther et al. 2017; Neiser et al. 2019). On the other hand, the use of GK in the turbulent simulation of the boundary region, which includes both the edge and the scrape-off-layer (SOL), remains challenging, despite the recent development of edge particle and continuum GK codes (Churchill et al. 2017; Mandell et al. 2020; Michels et al. 2021). Gyrokinetic simulations of the boundary are currently restricted by (i) their considerable computational cost, (ii) the presence of large-scale fluctuations, which are not present in the core, and (iii) the challenge of describing the high-collisionality regime using proper collision operator models, such as the Fokker–Planck Landau collision operator (Landau 1936), referred to as the Coulomb operator in this work. Turbulence in the SOL region is most often simulated by models based on drift-reduced Braginskii-like fluid equations, which evolve the lowest-order particle fluid moments (density, temperature and velocity) (Zeiler, Drake & Rogers 1997). Braginskii-like fluid simulations of the SOL turbulence have shown their ability to model the SOL in a complex magnetic field topology (see, e.g. Stegmeir et al. 2019; Giacomin, Stenger & Ricci 2020; Bufferand et al. 2021), in good agreement with experimental results (see, e.g. De Oliviera et al. 2022; Galassi et al. 2022). The validity of Braginskii-like models relies on the high-collisionality assumption, quantified by the smallness of the ratio of the particle mean-free path to the parallel scale length, $\lambda _{mfp} / L_\parallel \ll 1$. This scaling might not be appropriate to describe the entire collisionality range of the SOL and, more generally, in the boundary region. In particular, the high plasma temperature at the top of the pedestal and local transient events (such as edge localized modes) can significantly lower the plasma collisionality, even in the SOL, calling for a kinetic description of the boundary region. Aiming to bridge the gap between fluid and GK simulations, a moment approach to the GK model based on a Hermite–Laguerre decomposition of the full gyrocentre distribution function (full-F) was recently introduced in Frei, Jorge & Ricci (2020). This model, which we refer to as the gyro-moment (GM) approach, is derived in a generalized GK ordering appropriate to the boundary region and is valid for an arbitrary level of collisionality since it implements the full GK Coulomb collision operator (Jorge, Frei & Ricci 2019). The ability of the GM approach to describe drift waves (Jorge, Ricci & Loureiro 2018) and ion-scale instabilities (Frei, Hoffmann & Ricci 2022b) efficiently has been demonstrated at an arbitrary level of collisionality using the GK Coulomb collision operator and other advanced collision operator models (Frei et al. 2021; Frei, Ernst & Ricci 2022a). However, these investigations are limited to electrostatic and local linear studies neglecting, for instance, electromagnetic and trapped particle effects, excluding therefore instabilities such as the trapped electron modes (TEM), recognized as one of the main drives of electron heat transport in the boundary region (Rafiq et al. 2009; Schmitz et al. 2012), as well as the kinetic ballooning modes (KBM), which can limit, for instance, the maximal achievable pressure gradient in H-mode pedestals (Snyder et al. 2009; Wan et al. 2012).

The present work aims to extend previous GM investigations (Jorge et al. 2018, 2019; Frei et al. 2022b) to a tokamak flux-tube configuration. More precisely, the GK model we consider in this work, based on the $\delta f$ and linearized version of Frei et al. (2020), includes ion and electron species, trapped and passing particles, finite electromagnetic effects and collisions using advanced collision operators, such as the GK Coulomb (Li & Ernst 2011; Jorge et al. 2019; Pan & Ernst 2019), Sugama (Sugama, Watanabe & Nunami 2009), and improved Sugama (IS) (Sugama et al. 2019) collision operators (Jorge et al. 2019; Frei et al. 2021, 2022a). We remark that the GK model considered here is similar to the one implemented in the GPU-native code GX (Mandell et al. 2022), which is designed for fusion reactor optimization based on turbulent calculations and includes a Dougherty collision operator (Dougherty 1964). The linearized GM hierarchy equation that we develop allows us to investigate the linear properties of the ion temperature mode (ITG) with adiabatic and kinetic electrons, the TEM, the KBM, the microtearing mode (MTMs) and the dynamics of zonal flows (ZFs) including geodesic acoustic modes (GAMs) and ZF damping in regimes relevant to the boundary region, from the low-collisionality banana to the high-collisionality Pfirsch–Schlüter regime. Our numerical results are tested and verified in the collisionless limit with the state-of-the-art continuum GK code GENE (Jenko et al. 2000; Görler et al. 2011). More precisely, we compare the linear growth rates and mode frequencies, and investigate the velocity-space and the ballooning eigenmode structures. In particular, a careful investigation of the velocity-space structures of the distribution functions allows us to assess the convergence properties of the GM approach and identify the optimal number of GMs that need to be retained in the simulations. In addition, the present comparison provides physical insights into the performance of the GM approach to describe important microinstabilities. Finding excellent agreement with GENE in all the cases explored in the present work, we demonstrate that the GM approach can accurately capture kinetic physics such as, e.g. resonances due to parallel and perpendicular drifts of passing particles, trapped particles, magnetic gradient drift resonance and small-scale velocity-space features near the passing and trapped boundary. Furthermore, it is found that the number of GMs necessary to achieve convergence in the collisionless limit is often of the same order as the number of velocity-space grid points used in GENE, based on finite difference schemes. As expected, the number of GMs is significantly reduced as the level of collisionality increases (Frei et al. 2021). More interestingly, this is also true at low collisionality in the case of instabilities developing in steep pressure gradient conditions such as the ones appearing in H-mode operations. We remark that the comparisons presented here can also be extended to other GK codes using different discretization techniques in velocity space, such as the GS2 (Dorland 2000) and GKW (Peeters et al. 2009) codes. In addition to a comparison with the GENE code, we also perform a convergence study of the GM approach in the collisionless limit with a general electromagnetic dispersion relation of the GK model that we analytically derive.

In the high-collisionality Pfirsch–Schlüter regime, the regularization of the velocity-space distribution functions and the availability of advanced collision operator models expressed in terms of GMs allow us to derive reduced-fluid models as an asymptotic limit of the GM hierarchy equation. This illustrates the multi-fidelity aspect of the GM approach. A collision operator model comparison is carried out in this work by considering instabilities relevant to the edge region. More precisely, deviations in the TEM linear growth rates (up to $15\,\%$) between the GK Coulomb and other collision operators at collisionalities relevant to edge H-mode conditions are found, consistent with Pan, Ernst & Crandall (2020); Pan, Ernst & Hatch (2021). The amplitude of these deviations depends on the pressure gradients that drive the instability, such as the electron pressure gradient, and are absent for other edge instabilities such as MTMs, in contrast to the case of MTMs shown in Pan et al. (2021). In all cases, the IS operator model provides the smallest deviations with respect to the GK Coulomb. Finally, the impact of collisions on GAM dynamics and ZF damping is studied. It is shown that, in general, energy diffusion, conservation laws and finite Larmor radius (FLR) terms in the collision operator models cannot be ignored when predicting their correct, long-time evolution, consistent with Pan et al. (2020, 2021). In view of the importance of turbulent transport and its self-consistent interaction with ZFs in the boundary region, the present study highlights that a systematic assessment of the physics fidelity of collision operators is necessary for a detailed and correct description of the turbulent plasma dynamics in the boundary region and that the GM approach is an ideal tool to carry out such investigations.

The rest of this paper is structured as follows. In (2), we present the flux-tube linear GK model that we project onto the Hermite–Laguerre basis yielding the GM hierarchy equation, and discuss its numerical implementation. In § 3, we investigate the description within the GM approach of kinetic effects associated with drifts of passing particles. Section 4 presents a comprehensive collisionless study of microinstabilities and ZF dynamics with a detailed comparison against the GENE code. Collisional effects are introduced in § 5 where the high-collisional limit of the GM hierarchy is derived and the collisionality dependence of edge instabilities is revealed. In § 6, we use the GM approach to investigate microinstabilities at steep pressure gradients, typically found in low-collisionality H-mode conditions. Finally, a discussion and an outlook are presented in § 7. Appendix A reports on convergence studies of the GM approach using an electromagnetic GK dispersion relation.

2. Flux-tube GM model

The flux-tube approach allows for the simulation of plasma turbulence in a computational domain that extends along a magnetic field line and over a narrow region. The flux-tube configuration is motivated by the smallness of the ratio of the typical perpendicular turbulent scale length, which is of the order the ion Larmor radius $\rho _i$ (for ion-scale turbulence), to the perpendicular equilibrium scale $L_\perp$, $\rho _i / L_\perp \ll 1$, and by the anisotropic nature of turbulence along and perpendicular to the equilibrium magnetic field lines (Beer, Cowley & Hammett 1995; Xanthopoulos & Jenko 2006). While the flux-tube approach can be justified in the core region, the flux-tube model allows us to assess the use of the GM approach in the study of microinstabilities, which are relevant to the boundary region.

The presentation section is structured as follows. In § 2.1, we present the linearized GK model. The development of this model in a flux-tube geometry is reported in § 2.2. The GM approach based on a Hermite–Laguerre decomposition of the perturbed distribution functions is introduced in § 2.3. The collision operators used in this work are listed in § 2.4, and, finally, the numerical implementation of the GM hierarchy equation is discussed in (2.5).

2.1. The GK model

We consider the linearized electromagnetic GK Boltzmann equation in the presence of an equilibrium magnetic field, as well as density and temperature gradients. The flux-tube assumption of separation between the turbulent (of the order of $\rho _i$) and the equilibrium (of the order of $L_\perp$) scales allows us to neglect the equilibrium profiles as considered constant across the computational domain. In the following, we use the gyrocentre phase-space coordinates $\boldsymbol {Z} = (\boldsymbol {R}, \mu, v_\parallel,\theta )$, where $\boldsymbol {R} = \boldsymbol {r} - \boldsymbol{\rho} _a$ is the gyrocentre position, with $\boldsymbol{r}$ the particle position and $\boldsymbol{\rho} _a(\boldsymbol{R}, \mu, \theta ) = \boldsymbol{b} \times \boldsymbol{v}/\varOmega _a$ its gyroradius. Here, $\boldsymbol {b} = \boldsymbol{B} / B$, $\varOmega _a = q_a B / m_a$, $a$ is the particle species, $\mu = m_a v_\perp ^2/[2 B(\boldsymbol {R})]$ is the magnetic moment, $v_\parallel = \boldsymbol{b} \boldsymbol {\cdot } \boldsymbol{v}$ is the component of the velocity parallel to the equilibrium magnetic field and, finally, $\theta$ is the gyroangle. Contrary to Frei et al. (2020), we assume that the gyrocentre distribution function, $F_a = F_a (\boldsymbol {R}, \mu, v_\parallel,t)$, is a perturbed Maxwellian, i.e. $F_{a} = F_{Ma} + {g}_a$, with ${g}_a = {g}_a(\boldsymbol {R} ,\mu, v_\parallel, t)$ the perturbation with respect to the local Maxwellian distribution function $F_{Ma} = N {\rm e}^{- s_{\parallel a}^2 - x_a} / ({\rm \pi} ^{3/2} v_{Ta}^3)$, where ${g}_{a}/ F_{Ma} \ll 1$, $N = N_i(\boldsymbol {R}) = N_e(\boldsymbol {R})$ the background gyrocentre density (assuming $q_i = + e$ for simplicity), $s_{ \parallel a} = v_\parallel / v_{Ta}(\boldsymbol {R})$, $x_a= \mu B(\boldsymbol {R}) / T_{a}(\boldsymbol {R})$ and $v_{T_a}^2(\boldsymbol {R}) = 2 T_{a}(\boldsymbol {R})/m_a$. Under these assumptions, the linearized electromagnetic GK Boltzmann equation for the Fourier modes ${g}_a(\boldsymbol{k}_\perp, \ell, \mu, v_\parallel, t)$ (with $\ell$ the arc-length coordinate along a magnetic field line) is (Hazeltine & Meiss 2003)

(2.1)\begin{equation} \frac{\partial}{\partial t } {g}_a + {\rm i} \omega_{Ba} h_a + v_\parallel \boldsymbol{\nabla}_\parallel h_a - \frac{\mu}{m_a} (\boldsymbol{b} \boldsymbol{\cdot} \boldsymbol{\nabla} B)\frac{\partial }{\partial v_\parallel} h_a - {\rm i} \omega_{Ta}^* \frac{ e \chi_a }{T_e}F_{Ma} = \mathcal{C}_a, \end{equation}

where we introduce the gyro-averaged electromagnetic field, $\chi _a = {\rm J}_0(b_a \sqrt {x_a}) ( \phi - v_\parallel \psi )$, with the perturbed electrostatic potential, $\phi = \phi (\boldsymbol{k}_\perp, \ell, t)$ and the component parallel to $\boldsymbol{B}$ of the perturbed magnetic vector potential, $\psi = \psi (\boldsymbol{k}_\perp, \ell, t)$, defined such that the transverse component of the perturbed magnetic field is $\delta \boldsymbol{B}_\perp \simeq \boldsymbol {\nabla }_\perp \psi \times \boldsymbol{b}$. The perpendicular wavevector is defined as $\boldsymbol{k}_\perp = \boldsymbol{k} - (\boldsymbol{b} \boldsymbol {\cdot } \boldsymbol{k}) \boldsymbol{b}$ and $\ell$ is the arc length describing the direction along $\boldsymbol{B}$, such that the parallel gradient is $\boldsymbol {\nabla }_\parallel = \boldsymbol {b} \boldsymbol {\cdot } \boldsymbol {\nabla } = \partial _\ell$. In addition, we introduce the magnetic drift frequency $\omega _{Ba} = \boldsymbol{v}_{Da} \boldsymbol {\cdot } \boldsymbol{k}$, with $\boldsymbol{v}_{Da} = \mu \boldsymbol {b} \times \boldsymbol {\nabla } \ln B /q_a + v_\parallel ^2 / \varOmega _a \boldsymbol {b} \times \boldsymbol{\kappa}$ being the combination of the $\boldsymbol {\nabla } B$ and curvature drifts, and the diamagnetic frequency $\omega _{Ta}^* = [\omega _N + \omega _{T_a}( x_a + s_{\parallel a}^2 -3/2) ]$, with $\omega _N = T_e \boldsymbol {b} \times \boldsymbol {\nabla } \ln N \boldsymbol {\cdot } \boldsymbol{k}/ (eB)$ and $\omega _{T_a} = T_e \boldsymbol{b} \times \boldsymbol {\nabla } \ln T_a \boldsymbol {\cdot } \boldsymbol{k}/(eB)$. We remark that, using the magnetohydrodynamics (MHD) equilibrium condition, $\boldsymbol{J} \times \boldsymbol{B} = \boldsymbol {\nabla } P$ (with $P = \sum _{a}N_a T_a$ the total equilibrium pressure), and Ampere's law, $\boldsymbol {\nabla } \times \boldsymbol {B} = 4 {\rm \pi}\boldsymbol{J}$, the magnetic curvature can be expressed as $\boldsymbol{\kappa} = \boldsymbol{b} \boldsymbol {\cdot } (\boldsymbol {\nabla } \boldsymbol{b}) = \boldsymbol {\nabla }_\perp \ln B + (4 {\rm \pi}\boldsymbol {\nabla } P )/ B^2$, such that the magnetic drift frequency, $\omega _{Ba}$, becomes $\omega _{Ba} = v_{Ta}^2 ( x_a+ 2 s_{\parallel a}^2) R_B/ (2 \varOmega _a) + v_{Ta}^2 s_{\parallel a}^2 / \varOmega _a \boldsymbol{b} \times (4 {\rm \pi}\boldsymbol {\nabla } P) / B^2 \boldsymbol {\cdot } \boldsymbol{k}$, where $R_B = ( \boldsymbol{b} \times \boldsymbol {\nabla } \ln B ) \boldsymbol {\cdot } \boldsymbol{k}$. FLR effects give rise to the zeroth-order Bessel function, ${\rm J}_0(b_a \sqrt {x_a})$, where the argument $b_a = k_\perp v_{T_a} /\varOmega _a$ is the normalized perpendicular wavevector, with $k_\perp = |\boldsymbol{k}_\perp |$. The non-adiabatic part of the perturbed gyrocentre distribution function ${g}_a$ that appears in (2.1), $h_a = h_a (\boldsymbol{k}_\perp, \ell, \mu, v_\parallel,t)$, is defined by

(2.2)\begin{equation} h_a = {g}_a + \frac{q_a}{T_{a}} F_{Ma} \chi_a. \end{equation}

On the right-hand side of (2.1), the effect of collisions is described by the linearized collision operator $\mathcal {C}_a = \sum _{b} \mathcal {C}_{ab}(\boldsymbol{k}_\perp, \ell, \mu, v_\parallel )$ (Frei et al. 2021). The GK Boltzmann equation, (2.1), is closed by the GK quasi-neutrality condition

(2.3)\begin{equation} \sum_a \frac{ q_a^2}{T_a} \left( 1 - \varGamma_0(a_a) \right) \phi = \sum_a q_a \frac{1}{N_a} 2 {\rm \pi}\int \,{\rm d} \mu \,{\rm d} v_\parallel \frac{B}{m_a} {\rm J}_0(b_a \sqrt{x_a}) {g}_a, \end{equation}

that provides the self-consistent electrostatic potential (Frei et al. 2020), where $a_a= b_a^2 / 2$ and $\varGamma _0(x) = {\rm I}_0(x) {\rm e}^{- x}$, with ${\rm I}_0$ the modified Bessel function of order zero, and by the GK Ampere's law

(2.4)\begin{equation} \left(\frac{ k_\perp^2}{4 {\rm \pi}}+ \sum_a \frac{q_a^2 N_a}{m_a} \varGamma_0(a_a) \right) \psi = \sum_a q_a 2 {\rm \pi}\int \,{\rm d} \mu \,{\rm d} v_\parallel \frac{B}{m_a} {\rm J}_0(b_a \sqrt{x_a} ) v_\parallel {g}_a , \end{equation}

that provides the Fourier component of the perturbed magnetic vector potential $\psi$. We remark that the linear GK model in (2.1), (2.3) and (2.4) can be obtained from the full-F model presented in Frei et al. (2020) by neglecting nonlinearities and the terms in the guiding-centre transformation arising from the large amplitude and long wavelength components of the fluctuating electromagnetic fields.

In the present work, the adiabatic electron approximation is also considered. In this case, electron inertia is neglected, such that the parallel electric field balances the parallel pressure gradient, and therefore the electron density follows the perturbed electrostatic potential $\phi$. Imposing that the perturbed electron density vanishes on average on a flux surface, the GK quasi-neutrality condition, (2.3), can be simplified

(2.5)\begin{equation} \frac{q_i^2 }{T_i} \left( 1 - \varGamma(a_i) \right) \phi + \frac{e^2}{T_e} \left( \phi - \left\langle \phi \right\rangle_{fs} \right) = \frac{q_i}{N_i} \int \,{\rm d} \mu \,{\rm d} v_\parallel \,{\rm d} \theta \frac{B}{m_i} {\rm J}_0(b_i \sqrt{x_i}) {g}_i, \end{equation}

where $\left \langle \dots \right \rangle _{fs}$ denotes the flux-surface average operator (Dorland & Hammett 1993). The adiabatic electron approximation allows us to remove the fast electron dynamics that limits, for instance, the time step in turbulent simulations, thus allowing the study of ion-driven instabilities such as the ITG (Frei et al. 2022b). However, retaining the electron dynamics is essential in describing electromagnetic effects and instabilities driven unstable by trapped electrons.

2.2. Field-aligned coordinate system and flux-tube model

Taking advantage of the highly anisotropic turbulence along and across the magnetic field lines, we define a coordinate system with one coordinate aligned with the magnetic field line. To this aim, we introduce the Clebsch-type, field-aligned coordinate system $(x,y,z)$ and write the equilibrium magnetic field $\boldsymbol{B}$ as

(2.6)\begin{equation} \boldsymbol{B} = B_0 \boldsymbol{\nabla} x \times \boldsymbol{\nabla} y, \end{equation}

where $B_0$ is the reference magnetic field strength. Given (2.6), the coordinates $(x,y)$ generate a plane perpendicular to the magnetic field since $\boldsymbol {B} \boldsymbol {\cdot } \boldsymbol {\nabla } x = \boldsymbol {B} \boldsymbol {\cdot } \boldsymbol {\nabla } y =0$. On the other hand, the coordinate $z$ is used to describe the direction along the equilibrium magnetic field line. Among the Clebsch coordinates, we choose to consider (Lapillonne et al. 2009)

(2.7a–c)\begin{equation} x = X( \psi_p - \psi_p(0)), \quad y = Y (q(\psi_p) \chi - \phi_t), \quad z = \chi, \end{equation}

where $\psi _p$ is the poloidal flux label, $\psi _p(0)$ is the value of $\psi _p$ at the centre of the flux tube, $- {\rm \pi}\le \chi \le + {\rm \pi}$ is the straight-field line angle chosen to describe the parallel direction, $q(\psi _p)$ is the local safety factor, and $\phi _t$ the geometrical toroidal angle. Therefore, the coordinate $x$ is a radial magnetic flux-surface label while $y$ labels the magnetic field lines on a flux surface (binormal coordinate), with $X$ and $Y$ being normalization constants chosen such that $x$ and $y$ have the unit of length. The Jacobian of the coordinates system is $\mathcal {J}_{xyz} = (\boldsymbol {\nabla } x \boldsymbol {\cdot } \boldsymbol {\nabla } y \times \boldsymbol {\nabla } z)^{-1}$.

In the flux-tube model, the $x$ and $y$ directions are treated in Fourier space by assuming periodic boundary conditions along them (Ball & Brunner 2021). We thus introduce the perpendicular wavenumber vector $\boldsymbol{k}_\perp = k_x \boldsymbol {\nabla } x + k_y \boldsymbol {\nabla } y$, $k_x$ and $k_y$ being the radial and binormal wavenumbers, respectively. A real valued fluctuating quantity $A(x,y,z)$ is therefore expressed as

(2.8)\begin{equation} A(x,y,z) = \sum_{k_{x}, k_y} \mathcal{A}(k_x, k_y,z) \exp({{\rm i} k_x x + {\rm i} k_y y}), \end{equation}

with $\mathcal {A}(k_x , k_y, z)$ the Fourier components of $A$. The periodic boundary condition in $x$ is justified in the local approximation, whereby constant radial equilibrium gradients are considered, while the safety factor $q(\psi _p)$ is linearized around the centre of the flux-tube domain located at $x = 0$, i.e. we write $q(\psi _p) \simeq q [1 + x s / (X \psi _p(0)) ]$ and introduce the magnetic shear $s = (\psi _p(0)/q) \,{\rm d}q /{\rm d}\psi _p$, with $q = q(\psi _p(0))$ the safety factor at the centre of the flux tube (Beer et al. 1995). The periodic boundary condition in $y$ stems from the $2 {\rm \pi}$ periodicity in the geometrical toroidal angle $\phi _t$ (see (2.7a–c)). The periodicity in the straight-field line angle $\chi$ imposes the boundary conditions along $z$ (Beer et al. 1995; Lapillonne et al. 2009)

(2.9)\begin{equation} \mathcal{A}(k_x , k_y, z = {\rm \pi}) = \mathcal{A}(k_x + 2 {\rm \pi}s k_y, k_y, z ={-} {\rm \pi}). \end{equation}

The ballooning eigenmode function of the fluctuating quantity $\mathcal {A}$, denoted by $\mathcal {A}_B$, can be constructed by coupling the $(k_x,z)$ linear modes through the ballooning transformation (Connor, Hastie & Taylor 1978)

(2.10)\begin{equation} \mathcal{A}_B(\chi) = \mathcal{A}( k_x + n_{k_x} 2 {\rm \pi}s k_y, k_y, z ), \end{equation}

where $- \infty \le \chi = z + 2 {\rm \pi}n_{k_x} \le \infty$ (with $-{\rm \pi} \leq z \leq {\rm \pi}$) is the extended ballooning angle.

We note that the norm of the perpendicular wavenumber $\boldsymbol{k}_\perp$, that enters in, e.g. the Bessel function ${\rm J}_0$ appearing in (2.1), is expressed by

(2.11)\begin{equation} k_\perp{=} \sqrt{K_x k_x + g^{xy} k_x k_y + g^{yy} k_y^2}, \end{equation}

where we introduce the effective radial wavenumber $K_x = \boldsymbol {\nabla } x \boldsymbol {\cdot } \boldsymbol{k}_\perp = g^{xx} k_x + g^{xy} k_y$ and the geometrical coefficients given by the metric tensor elements $g^{xx} = \boldsymbol {\nabla } x \boldsymbol {\cdot } \boldsymbol {\nabla } x$, $g^{xy} = \boldsymbol {\nabla } x \boldsymbol {\cdot } \boldsymbol {\nabla } y$, $g^{yy} = \boldsymbol {\nabla } y \boldsymbol {\cdot } \boldsymbol {\nabla } y$ (similar definitions are used for $g^{yz}$, $g^{xz}$ and $g^{zz}$). We remark that $g^{xy} \neq 0$ since the $x$ and $y$ coordinates are not orthogonal.

Using the fact that the equilibrium density and temperature varies only along $x$ (i.e. $\boldsymbol {\nabla } N = \boldsymbol {\nabla } x \partial _x N$ and $\boldsymbol {\nabla } T_a = \boldsymbol {\nabla } x \partial _x T_a$), the linearized GK Boltzmann equation, (2.1), describing the time evolution of ${g}_a = {g}_a (k_{x},k_y, z, \mu, v_\parallel )$, reads in the $(x,y,z)$ coordinate system, as

(2.12)\begin{gather} \frac{\partial}{\partial t } {g}_a + \frac{v_{Ta}}{\mathcal{J}_{xyz}} \frac{s_{{\parallel} a}}{\hat{B}} \frac{\partial}{\partial z} h_a + {\rm i} \omega_{Ba}h_a - \frac{x_a v_{Ta}}{2} \frac{1}{\mathcal{J}_{xyz} \hat{B} } \frac{\partial }{\partial z} \ln B \frac{\partial }{\partial {s_{{\parallel} a}}} h_a \nonumber\\ + {\rm i} \omega_{Ta}^* \frac{e \chi_a}{T_e} F_{aM} = \mathcal{C}_{a}, \end{gather}

where $\hat {B}^2 = B^2 / B_0^2 =g^{xx} g^{yy} - g^{xy} g^{xy}$, and the frequencies

(2.13)\begin{equation} \omega_{Ba} = \frac{v_{Ta}^2}{2 \varOmega_a} \left( x_a+ 2 s_{{\parallel} a}^2 \right) C_{x,y}(B) - \frac{v_{Ta}^2}{2 \varOmega_a} s_{{\parallel} a}^2 \frac{\hat{B}}{L_\perp} \frac{\alpha}{q^2 } , \end{equation}

and

(2.14)\begin{equation} \omega_{Ta}^* = \frac{1}{L_\perp}\left[ R_N + R_{Ta} \left( x_a + s_{ {\parallel} a}^2 -\frac{3}{2}\right) \right] \frac{T_e k_y}{e B}. \end{equation}

Here, the normalized density and temperature gradients, $R_N = - L_\perp \partial _x \ln N$ and $R_{Ta} = - L_\perp \partial _x \ln T_{a}$, respectively, and the MHD parameter $\alpha = q^2 \beta _e \sum _{a} \tau _a ( R_N + R_{Ta})$. The flux-tube approach allows us to approximate the density and temperature gradient lengths by their local values evaluated at $x =0$, $L_N$ and $L_{T_a}$, respectively, such that $\partial _x \ln N_a = - 1/L_N$ and $\partial _x \ln T_a = - 1/L_{T_a}$. The curvature operator, $C_{x,y}(B)$ in (2.13), is defined by

(2.15)\begin{equation} C_{x,y}(B) = \mathcal{C}_x(\ln B) k_x + \mathcal{C}_y( \ln B) k_y, \end{equation}

where $\mathcal {C}_x(A) = (\varGamma _1\partial _y A +\varGamma _2\partial _z A ) / \hat {B}$, $\mathcal {C}_y(A)= (\varGamma _3\partial _z A - \varGamma _1 \partial _x A)/ \hat {B}$ (with $\varGamma _1 = g^{xy} g^{yx}-g^{xx} g^{yy}$, $\varGamma _2 = g^{xz} g^{yx}-g^{xx} g^{yz}$ and $\varGamma _3 = g^{xz} g^{yy}-g^{xy} g^{yz}$).

In the present numerical implementation, we consider concentric and circular flux surfaces modelled by the $s-\alpha$ model (Dimits et al. 2000). Despite its known inconsistencies (Lapillonne et al. 2009), the $s-\alpha$ model provides an efficient and easy-to-implement model that can be used to validate simulation codes when the details of the magnetic geometry are not important. In the $s-\alpha$ model, the normalized amplitude of the magnetic field is given by $\hat {B} = B / B_0 = 1 / (1 + \epsilon \cos z)$ where $\epsilon$ is the inverse aspect ratio assumed to be small, $\epsilon \ll 1$. It follows that $\mathcal {J}_{xyz} \hat {B} = q R_0$ (with $R_0$ the major radius of the tokamak device) and the non-zero metric elements are $g^{xx} =1$, $g^{xy} = s z$, $g^{yy} = 1 + z^2 s^2$. We choose the reference equilibrium length $L_\perp$ to be the major radius of the tokamak device, i.e. we set $L_\perp = R_0$. The parallel derivative of the magnetic field strength $B$ and the curvature operator $C_{x,y}(B)$ are therefore expressed by

(2.16)\begin{gather} \frac{\partial }{\partial z} \ln B = \epsilon \sin z, \end{gather}

(2.17)\begin{gather}C_{x,y}(B ) ={-} \frac{\hat{B}}{R_0}(\sin z K_x + \cos z k_y ), \end{gather}

with $K_x = k_x + s z k_y$. Given the expressions of the metric elements, the perpendicular wavenumber $k_\perp$, defined in (2.11), becomes

(2.18)\begin{equation} k_\perp{=}\sqrt{ k_x K_x + s z k_x k_y + (1 + s^2 z^2) k_y^2}. \end{equation}

The linearized electromagnetic GK Boltzmann equation, given in (2.1), coupled with the GK field equations, (2.3) and (2.4), constitute a closed set of partial differential equations. Within a continuum numerical approach, this set of equations is discretized using a two-dimensional velocity-space grid where the velocity-space derivatives and integrals contained in (2.1) and in the collision operator $\mathcal {C}_{ab}$ are evaluated numerically. For instance, the widely used GK continuum code GENE (Jenko et al. 2000) uses a uniform grid in the $(v_\parallel, \mu )$ coordinates in its local and linear flux-tube implementation. Using a different approach, we develop the GK model into a set of fluid-like equations by expanding the distribution function on a polynomial basis in the velocity-space coordinates $(v_\parallel, \mu )$.

2.3. Gyro-moment expansion

We use a GM approach based on a Hermite–Laguerre expansion of the perturbed distribution function ${g}_a$ to solve the electromagnetic linearized GK equation given in (2.12). More precisely, the perturbed gyrocentre distribution function, ${g}_a$, is expanded onto a Hermite–Laguerre polynomial basis (Jorge, Ricci & Loureiro 2017; Mandell, Dorland & Landreman 2018; Jorge et al. 2019; Frei et al. 2020), such that

(2.19)\begin{equation} {g}_a = \sum_{p = 0}^\infty \sum_{j = 0}^\infty N_a^{pj} \frac{H_p(s_{{\parallel} a}) L_j(x_a)}{\sqrt{2^p p!}} F_{Ma}. \end{equation}

In (2.19), we introduce the physicist's Hermite and Laguerre polynomials, $H_p$ and $L_j$, that can be defined via their Rodrigues’ formulas (Gradshteyn & Ryzhik 2014)

(2.20a)\begin{gather} H_p(x) = ({-}1)^p {\rm e}^{x^2} \frac{{\rm d}^p}{{{\rm d}\,x}^p} \left( {\rm e}^{- x^2} \right), \end{gather}

(2.20b)\begin{gather}L_j(x) = \frac{{\rm e}^{x}}{j!} \frac{{\rm d}^j }{{{\rm d}\,x}^j} \left ( {\rm e}^{- x } x^j\right), \end{gather}

and we note their orthogonality relations

(2.21a)\begin{gather} \int_{- \infty}^\infty \,{{\rm d}x} H_p(x) H_{p'}(x) {\rm e}^{- x^2} = 2^p p! \sqrt{\rm \pi} \delta_{p}^{p'}, \end{gather}

(2.21b)\begin{gather}\int_0^\infty \,{{\rm d}\,x} L_j(x) L_{j'}(x) {\rm e}^{{-}x} = \delta_j^{j'}. \end{gather}

Using the orthogonality relations, the Hermite–Laguerre velocity moments of ${g}_a$, i.e. the GMs $N_a^{pj}$, are defined by

(2.22)\begin{equation} N_a^{pj}(k_x, k_y, z) = \frac{1}{N} 2 {\rm \pi}\int \,{\rm d} \mu \,{\rm d} v_{{\parallel}} \frac{B}{m_a} {g}_a \frac{H_p(s_{{\parallel} a}) L_j(x_a)}{\sqrt{2^p p!}}, \end{equation}

with $N = \int \,{\rm d} \mu \,{\rm d} v_{\parallel } \,{\rm d} \theta B F_{Ma} / m_a$ the background gyrocentre density. We remark that any polynomial basis could, in principle, be used to expand the perturbed distribution function ${g}_a$. For instance, a polynomial basis of interest for high-collisional plasmas, based on Legendre and associated Laguerre polynomials in the pitch-angle and speed coordinates $\xi = v_\parallel / v$ and $v$ (or energy $v^2$) respectively, can be used (Belli & Candy 2011). However, the use of the Hermite–Laguerre basis, which has a long history in plasma physics (see, e.g. Grant & Feix 1967; Hirshman & Sigmar 1976; Madsen 2013; Schekochihin et al. 2016; Jorge et al. 2017; Mandell et al. 2018), provides a direct relation to the fluid quantities that are evolved by Braginskii-like fluid models (Zeiler et al. 1997). For instance, $N_a^{10}$ is associated with the normalized parallel velocity, $u_{a \parallel }$, while $N_a^{20}$ and $N_a^{01}$ to the parallel and perpendicular temperatures, $T_{\parallel a}$ and $T_{\perp a}$.

The Bessel function ${\rm J}_0$ (appearing in both (2.1) and (2.3) and arising from FLR effects) and, more generally ${\rm J}_m$, with $m > 0$, can be conveniently expanded onto associated Laguerre polynomials, $L^m_n(x) = (-1)^m \,{\rm d}^m L_{n + m}(x) / {{\rm d}\,x}^m$, as (Gradshteyn & Ryzhik 2014)

(2.23)\begin{equation} {\rm J}_m(b_a \sqrt{x_a}) = \left(\frac{b_a \sqrt{x_a}}{2}\right)^{m}\sum_{n=0}^\infty \frac{n!\mathcal{K}_{n}(b_a )}{(n + m)!} L^m_n(x_a), \end{equation}

where we introduce the velocity-independent expansion coefficients

(2.24)\begin{equation} \mathcal{K}_{n}(b_a) = \frac{1}{n!}\left( \frac{b_a}{2}\right)^{2n } {\rm e}^{- b_a^2 /4}. \end{equation}

To simplify our notation, in the rest of the paper we normalize the time $t$ to $R_0 / c_{s}$ (with $c_{s}^2 = T_e / m_i$ the ion sound speed), the perpendicular wavenumbers $k_\perp$, $k_x$ and $k_y$ to $\rho _s = c_{s} / \varOmega _i$ the ion sound gyroradius (with $\varOmega _i = q_i B_0 / m_i$ the ion gyrofrequency defined with the reference magnetic field $B_0$), the particle mass $m_a$ to $m_i$, the particle charge $q_a$ to the electron charge $e$, the temperature $T_a$ to the electron equilibrium temperature $T_e$, the electrostatic potential $\phi$ to $T_{e} / e$, and the magnetic vector potential $\psi$ to $\rho _s B_0$.

We now project the linearized GK Boltzmann equation onto the Hermite–Laguerre basis by multiplying (2.1) by $B H_p L_j / \sqrt {2^p p!}$ and integrating over the velocity space. This yields the linearized GM hierarchy equation defined by

(2.25)\begin{align} & \frac{\partial }{\partial t} N_a^{pj} + \frac{ L_\perp}{\mathcal{J}_{xyz}} \frac{1}{\hat{B}} \frac{\sqrt{ \tau_a}}{\sigma_a} \left\{ \left( \sqrt{p+1} \frac{\partial}{\partial z} n_a^{p+1j}+ \sqrt{p}\frac{\partial}{\partial z} n_a^{p-1j} \right) \right. \nonumber\\ & \quad \left. - \frac{\partial}{\partial z} \ln B \left( (j+1)\sqrt{p+1} n_a^{p+1j} - j \sqrt{p} n_a^{p-1j} - j \sqrt{p+1} n_a^{p+1j-1} + \sqrt{p} (j +1) n_a^{p-1j+1} \right) \right\} \nonumber\\ & \quad + \left( \frac{{\rm i} \tau_a L_\perp}{q_a \hat{B}} C_{x,y}(B) + \frac{{\rm i} \tau_a }{q_a} \frac{({-}1) \alpha}{q^2} k_y \right) \left( \sqrt{(p+1)(p+2)}n_a^{p+2j} +(2p+1) n_a^{pj} \right. \nonumber\\ & \quad \left. + \sqrt{p(p-1)} n_a^{p-2j}- j n_a^{pj-1} - (j+1) n_a^{pj+1} \right) + \frac{{\rm i} \tau_a L_\perp}{q_a \hat{B}} C_{x,y}(B) ( 2 j +1) n_a^{pj} \nonumber\\ & \quad + {\rm i} \left[ \mathcal{K}_{j} \delta_p^0 R_N + R_{T_a} \left(\frac{1}{\sqrt{2}}\mathcal{K}_{j} \delta_p^2 + \delta_p^0 \left(2 j \mathcal{K}_{j} - j \mathcal{K}_{j-1 } - (j+1)\mathcal{K}_{j+1}\right) \right) \right] k_y \phi \nonumber\\ & \quad - {\rm i} \frac{\sqrt{2 \tau_a}}{\sigma_a}\left[ \frac{\mathcal{K}_{j} \delta_p^1 }{\sqrt{2}} R_N + R_{T_a} \left(\frac{\sqrt{3}}{2}\mathcal{K}_{j} \delta_p^3 + \frac{\delta_p^1}{\sqrt{2}} \left((2 j +1)\mathcal{K}_{j} - j \mathcal{K}_{j-1 } - (j+1)\mathcal{K}_{j+1}\right) \right) \right] k_y \psi \nonumber\\ & \qquad = \mathcal{C}_{a}^{pj}, \end{align}

with $\sigma _a = \sqrt {m_a / m_i}$ and $\tau _a = T_a /T_e$. In (2.25), we define $\mathcal {C}_{a}^{pj} = \sum _{b} \mathcal {C}_{ab}^{pj}$ with $\mathcal {C}_{ab}^{pj} =\mathcal {C}_{ab}^{pj}(k_x, k_y, z)$ the Hermite–Laguerre expansion of the linearized collision operator between species $a$ and $b$

(2.26)\begin{equation} \mathcal{C}_{ab}^{pj} = 2 {\rm \pi}\int \,{\rm d} \mu \,{\rm d} v_\parallel \frac{B}{m_a} \frac{H_p(s_{{\parallel} a}) L_j(x_a)}{ \sqrt{2^p p!}} \mathcal{C}_{ab}. \end{equation}

We remark that, in the case of GK collision operators, the linearized collision operator, $\mathcal {C}_{ab}^{pj}$, depends on $k_x$, $k_y$ and $z$ through the modulus of the perpendicular wavenumber $k_\perp$ (see (2.18)). On the other hand, $\mathcal {C}_{ab}^{pj}$ becomes independent of $k_\perp$, if drift-kinetic (DK) collision operators are used. In (2.25), we also introduce the non-adiabatic GMs $n_a^{pj}$, that are obtained by projecting (2.2) onto the Hermite–Laguerre basis, yielding

(2.27)\begin{equation} n_a^{pj} = N^{pj}_a + \frac{q_a}{ \tau_a } \mathcal{K}_{j} \left( \phi \delta_p^0 - \frac{ \sqrt{\tau_a} }{ \sigma_a} \delta_p^1 \psi \right). \end{equation}

Finally, the GK quasineutrality condition and the GK Ampere's law, (2.3) and (2.4), are normalized and expressed in terms of GMs as follows:

(2.28)\begin{equation} \sum_{a} \frac{q_a^2}{\tau_a}\left( 1 - \sum_{n=0}^{\infty} \mathcal{K}_{n}^2 \right) \phi = \sum_a q_a \sum_{n=0}^{\infty} \mathcal{K}_{n} N_a^{0n}, \end{equation}

and

(2.29)\begin{equation} \left( 2 k_\perp^2 + \beta_e \sum_a \frac{q_a^2}{\sigma_a^2} \sum_{n=0}^{\infty} \mathcal{K}_{n}^2 \right) \psi = \beta_e \sum_a q_a \frac{\sqrt{\tau_a}}{\sigma_a} \sum_{n=0}^{\infty} \mathcal{K}_{n} N_a^{1n}, \end{equation}

respectively, where $\beta _e = 8 {\rm \pi}N T_e / B_0^2$ is the electron plasma beta. On the other hand, assuming adiabatic electrons, the GK quasi-neutrality equation, (2.5), becomes

(2.30)\begin{equation} \left[ 1 + \frac{q_i^2}{\tau_i}\left( 1 - \sum_{n=0}^\infty \mathcal{K}_{n}^2 \right) \right] \phi - \left\langle \phi \right\rangle_{fs} = q_i \sum_{n=0}^\infty \mathcal{K}_{n} N_i^{0n}, \end{equation}

where the flux-surface-averaged operator of a function $f$ is expressed as $\left \langle f \right \rangle _{fs}= \int \,{{\rm d} y} \int \,{\rm d} z \mathcal {J}_{xyz} f / \int \,{\rm d}z \int \,{{\rm d} y} \mathcal {J}_{xyz}$. We remark that the argument $b_a = \sigma _a \sqrt {2 \tau _a} k_\perp / \hat {B}$ of the kernel functions, $\mathcal {K}_{j} = \mathcal {K}_{j}(b_a)$ defined in (2.24), depends on geometrical quantities, through $k_\perp$ given in (2.11), and on the magnetic field strength $B$, through its $\rho _a$ dependence. We remark that a similar Hermite–Laguerre approach in the $\delta f$ limit of the GK model has been recently formulated and implemented in the GX code (Mandell et al. 2018, 2022), showing a promising numerical efficiency to simulate the collisionless core region to optimize future reactor designs.

2.4. Linearized collision operator models

To model the effects of collisions $\mathcal {C}_{ab}^{pj}$ on the right-hand side of (2.25), we use the GM expansion of advanced collision operator models previously derived and benchmarked in Frei et al. (2021, 2022b, a). In contrast to the GX code (Mandell et al. 2022) that implements a Dougherty collision operator, we consider here the linearized GK Coulomb (Li & Ernt 2011; Pan & Ernst 2019; Frei et al. 2021), the Sugama (Sugama et al. 2009), the IS (Sugama et al. 2019) collision operators, and like-species Dougherty (Dougherty 1964) collision operators.

Collisional effects are described by means of the ion–ion collision frequency (normalized to the ion transit time $R_0 / c_{s}$)

(2.31)\begin{equation} \nu_{ii} = \frac{4 \sqrt{\rm \pi}}{3} \frac{R_0 N e^4 \ln \varLambda}{ c_{s} m_i^{1/2} T_i^{3/2} } , \end{equation}

with $\ln \varLambda$ the Coulomb logarithm. The normalized electron–ion collision frequency is then

(2.32)\begin{equation} \nu_{ei} = \frac{\nu_{ii}} { \sqrt{m_e / m_i }} \left( \frac{T_i}{T_e} \right)^{3/2}. \end{equation}

The electron and ion neoclassical collisionalities, $\nu _{e}^*$ and $\nu _i^*$, respectively, are then expressed by (Helander & Sigmar 2002)

(2.33a,b)\begin{equation} \nu_e^* = \frac{\sqrt{2} q}{\epsilon^{3/2}} \frac{T_i^{3/2} }{T_e^{3/2}} \nu_{ii}, \quad \nu_i^* = \frac{q}{\sqrt{2} \epsilon^{3/2}} \left( \frac{T_e}{T_i} \right)^{1/2} \nu_{ii}, \end{equation}

being the collisionless banana regime achieved when $\nu _e^* \lesssim 1$ and the high-collisional Pfirsch–Schlüter regime when $\nu _e^* \gtrsim 1 / \epsilon ^{3/2}$ for the electrons.

2.5. Numerical implementation

To solve numerically the linearized GM hierarchy equation, (2.25), we evolve a finite number of GMs, $(p,j) \leq (P,J)$. Throughout the present work, we consider the same $(P,J)$ for both electrons and ions. In addition, we use a simple closure by truncation by imposing $N_a^{pj} = 0$ for $(p,j) > (P,J)$. While rigorous asymptotic closures can be used (e.g. a high-collisional closure (Jorge et al. 2017) or a semi-collisional closure Loureiro, Schekochihin & Zocco 2013), the closure by truncation appears to be sufficiently accurate for the purposes of the present linear study.

For the spatial discretization, we use a single $k_y$ mode in an axisymmetric equilibrium and evolve a finite number, $2 N_{k_x}+1$, of $k_x$ modes (the $k_x$ modes are coupled through the parallel boundary condition at finite shear according to (2.9)). The values of the $k_x$ modes allowed in the system are imposed by (2.9) and are labelled by $k_{x,n} = \delta k_x \pm n_{k_x} 2 {\rm \pi}s k_y$ with $n_{k_x} =0, 1, \ldots, N_{k_x}$, where $\delta k_x = - z_0 k_y s$. However, for simplicity, we centre the grid of radial modes around the $k_x =0$ mode and neglect the effects of the finite ballooning angle $z_0$ by setting $\delta k_x =0$, if not specified otherwise. The $z$ direction, $- {\rm \pi}< z \leq {\rm \pi}$, is discretized using $N_z$ grid points that are uniformly distributed and the parallel derivatives, appearing in (2.25), are evaluated using a fourth-order centred finite difference scheme. Hyperdiffusion in $z$, proportional to $\sim \eta _z \partial ^4_z$, is added on the right-hand side of (2.25) to avoid artificial numerical oscillations. Since a finite number of $k_x$ modes are evolved, boundary conditions for the $n_{k_x} = \pm N_{k_x}$ modes are needed for $n_a^{pj}$. While different choices of boundary conditions exist, we consider

(2.34)\begin{equation} n_a^{pj}( - N_{k_x} 2 {\rm \pi}s k_y, k_y, - {\rm \pi}) = n_a^{pj}( + N_{k_x} 2 {\rm \pi}s k_y, k_y, {\rm \pi}), \end{equation}

for all $(p,j) \leq (P,J)$. For comparison, we remark that homogeneous Dirichlet boundary conditions are used in GENE. However, by increasing $N_{k_x}$ and $N_z$, our tests show that our results are not affected by the boundary conditions we impose along $z$.

An explicit fourth-order Runge–Kutta scheme is used to perform the time integration of (2.25). We denote with $\Delta t$ the time step and $t_n$ the discrete time values. We remark that the largest possible time step, $\Delta t$, when the electron dynamics is included, is limited by the presence of the high-frequency wave $\omega _H$ (Lee 1987; Lin et al. 2007).

In the present work, the complex frequency of the linear modes, $\omega = \omega _r + {\rm i} \gamma$ (where $\omega _r$ is the real mode frequency and $\gamma$ is the mode growth rate), is computed by using the weighted average

(2.35)\begin{equation} \omega^n(k_y) = \frac{\sum_{k_x, z} \omega_l^n(k_x , k_y, z) W(k_x , k_y,, z)}{\sum_{k_x, z} W(k_x , k_y,, z)}, \end{equation}

of the local complex frequency $\omega _l^n(k_x , k_y,,z) = \ln [ \phi ^{n}(k_x , k_y, z) / \phi ^{n-1}(k_x , k_y,, z)] / \Delta t$ (where $\phi ^n$ is the perturbed electrostatic potential at time $t = t_n$). Choosing $W( k_x, k_y, z) = \phi ^{n-1}( k_x, k_y, z)$, we evolve (2.25) until

(2.36)\begin{equation} \frac{\sum_{k_x, z} | \omega_l^n( k_x, k_y, z) - \omega^n (k_y)|^2 W( k_x, k_y, z)}{\sum_{k_x, z} W( k_x, k_y, z)} < \delta, \end{equation}

being $\delta = 10^{-4}$ for all the linear computations presented here. We note that we initialize the evolution of the GM hierarchy by imposing a perturbed density of constant amplitude along $z$ for all $k_x$ modes. Finally, we remark that a Cartesian message passing interface (MPI) domain decomposition along $k_x$, $z$ and the Hermite index $p$ is used. While the present parallelization is sufficient for the applications presented in this work, we believe that better parallelization strategies can be applied to achieve high computing performances and good scalability, also in comparison with present GK codes. For instance, the GPU-native GX code (Mandell et al. 2022) has been developed for this purpose and successfully achieved this goal.

A comparison between the continuum GK GENE code (Jenko et al. 2000; Görler et al. 2011) and the GM approach is presented in § 4. In the GENE code, the velocity space is discretized by uniformly distributed grid points between the normalized intervals $s_{\parallel } = v_\parallel / v_{Ta} \in [ - s_{\parallel M} , + s_{\parallel M} ]$ and $x = \mu B /T_a \in [0, x_M ]$ (typically $s_{\parallel M} = 3$ and $x_M = 9$ in our calculations) with a fixed number of grid points in each direction that we denote by $N_{v_\parallel }$ and $N_{\mu }$, respectively. The velocity-space derivatives and integrals are then approximated using finite difference methods. Hence, the numerical approximation of the distribution function, ${g}_a$, is given through the value of ${g}_a$ on a set of discrete grid points. On the other hand, within the GM approach, the numerical approximation of ${g}_a$ is given by the Hermite–Laguerre expansion coefficients, $N_a^{pj}$, such that the distribution function is reconstructed thanks to the truncated expansion in (2.19), given $P$ and $J$.

3. Representation of passing particle drifts in the GM approach

To interpret the investigations of microinstabilities in § 4, we first study analytically and numerically the GM approach description of kinetic effects associated with the parallel streaming and perpendicular drifts of passing particles. Particle resonances driven by these drifts play an important role, e.g. in geodesic acoustic mode (GAM) oscillations, in the ZF dynamics, and more generally, in the collisionless mechanisms of microinstabilities (Winsor, Johnson & Dawson 1968; Rosenbluth & Hinton 1998). In addition, the parallel streaming of passing particles and the finite orbit width (FOW) effects associated with magnetic gradient drifts can create fine-scale velocity-space structures in the distribution function (Idomura et al. 2008). It was recently reported (Frei et al. 2022b) that magnetic gradient drifts broaden the GM spectrum (both Hermite and Laguerre moments), while the parallel streaming of passing particles usually leads to the requirement of a larger number of Hermite than Laguerre GMs. Due to their importance, in particular at low collisionality (e.g. in the banana regime), we identify situations where a large number of GMs is necessary to resolve fine velocity-space structures. To investigate the representations of kinetic effects using the GM approach and if not stated otherwise, we consider the shearless limit ($s=0$), the safety factor $q = 1.4$, and the inverse aspect ratio $\epsilon = 0.1$. In addition, we focus on passing ions with adiabatic electrons and, therefore, omit the species label $a$ in this section for simplicity.

In the remainder of the present section, we study the parallel streaming of passing particles and illustrate the associated recurrence phenomena in § 3.1. A comparison with the GENE code confirms the ability of the GM method in the description of fine $v_\parallel$ structures. FOW effects driven by the perpendicular magnetic drifts are assessed in § 3.2.

3.1. Parallel streaming and recurrence phenomena

Passing particles are known to generate fine filament-like structures in $v_\parallel$ (Idomura et al. 2008), on scales that decrease linearly with time. To illustrate the appearance of these fine-scale structures and their effect on the GMs, we consider a simple one-dimensional model for the distribution function ${g} = {g}(\ell,v_\parallel,t)$ that describes the streaming of particles along the magnetic field lines (Hammett et al. 1993). Expressed in physical units, this reads

(3.1)\begin{equation} \frac{\partial }{\partial t} {g} + v_\parallel \partial_\ell {g} = 0, \end{equation}

with the initial condition ${g}(\ell,v_\parallel,0) = h(v_\parallel ) \cos (k_\parallel \ell )$, being $h(v_\parallel )$ a continuous function of $v_\parallel$ and $\ell$ the curvilinear coordinate along the magnetic field lines. The solution of (3.1), ${g}(\ell,v_\parallel,t) = h(v_\parallel ) \cos [k_\parallel ( \ell - v_\parallel t)]$, shows an effective wavenumber in velocity space $k_{v_\parallel } = k_\parallel t$ that increases linearly with time. Therefore, finer and finer-scale structures in $v_\parallel$ appear progressively. To understand the properties of the GM approach to solve (3.1), we introduce the Hermite moments, $N^p=\int d v_\parallel {g} H_p(s_\parallel ) {\rm e}^{- s_\parallel ^2}/ \sqrt {{\rm \pi} 2^p p! }$. Assuming $h(v_\parallel ) = h_0$ constant, the analytical expressions of $N^p$, satisfying the moment hierarchy equation, $\partial _t N^p + v_T (\sqrt {p+1} \partial _\ell N^{p+1} + \sqrt {p} \partial _\ell N^{p-1})/\sqrt {2} =0$ associated with (3.1), can be obtained by projecting the analytical solutions of ${g}$. One finds

(3.2)\begin{equation} N^{p} = \begin{cases} h_0 \cos(k_\parallel \ell) \dfrac{({-}1)^{p/2} 2^{p/2}}{\sqrt{2^{p}p!}} \left( \dfrac{ \omega_t t}{\sqrt{2}} \right)^{p} {\rm e}^{ - (\omega_t t)^2/4 }, & \quad p =2n \\ h_0 \cos(k_\parallel \ell) \dfrac{({-}1)^{(p-1)/2} 2^{p/2}}{\sqrt{2^{p}p!}} \left( \dfrac{\omega_t t}{\sqrt{2}} \right)^{p} {\rm e}^{ - (\omega_t t)^2/4 }, & \quad p = 2n +1, \end{cases} \end{equation}

where we introduce the transit frequency, $\omega _t = k_\parallel v_T$. The filamentation in $v_\parallel$ yields the propagation of a wavepacket in the Hermite spectrum to higher values of $p$ as time increases, with the maximum of the spectrum occurring at $\omega _t t = \sqrt { 2p}$. The increase of the effective wavenumber in velocity space, $k_{v_\parallel }$, with time challenges both the continuum numerical algorithms and the GM approach. In fact, $\lambda _{v_\parallel } = 2 {\rm \pi}/ k_{v_\parallel }$ typically sets the minimal distance between the grid points $\Delta v_\parallel$ in $v_\parallel$. Similarly, the minimal number $P$ of Hermite polynomials necessary for convergence increases with $k_{v_\parallel }$. An approximate expression of $k_{v_\parallel }$, that can be represented by an Hermite polynomial of order $p$, can be derived by noticing that the distance between the roots of the Hermite polynomials is of the order of ${\rm \pi} v_T /\sqrt {2 p}$, yielding $k_{v_\parallel } \simeq 2 \sqrt {2 p} / v_T \sim \sqrt {p} / v_T$.

As a consequence of the finite velocity-space resolution, a recurrence phenomenon occurs, which limits the validity of the numerical solutions. The recurrence manifests as a time-periodic perturbations, which have a purely numerical origin. Recurrence is observed both in the continuum method and in the GM approach. The recurrence time, $T_R$, is the time necessary for the structures in the distribution function to develop on a scale comparable to the numerical resolution, i.e. $k_\parallel T_R \sim k_{v_\parallel }^{{\rm max}}$. Within a continuum approach, $T_R$ is estimated as $T_R \simeq 2 {\rm \pi}q R_0 / \Delta v_\parallel$ (considering $k_\parallel \simeq 1 / q R_0$ typical of an interchange mode), while one has

(3.3)\begin{equation} T_R \simeq 2 \sqrt{2 P} \frac{q R_0}{ v_T }, \end{equation}

within the GM approach. Therefore, in continuum GK codes, the recurrence time is expected to scale linearly with the number of grid points $N_{v_\parallel }$, while $T_R$ scales less favourably in the GM approach as $\sqrt {P}$, according to (3.3).

To illustrate the recurrence phenomenon, we consider the time evolution of the flux-surface averaged electrostatic potential, $\left \langle \phi \right \rangle _{fs}$, in the absence of density and temperature gradients, at long radial wavelength and with a small and finite collisionality ($\nu _{ii} \simeq 10^{-4}$). The electrostatic potential, $\left \langle \phi \right \rangle _{fs}$, evolves into oscillations, associated with GAM) (the collisionless dynamics of GAMs is investigated in § 4.5) that are ultimately damped. We perform the simulations for different values of $P$ (with $J = 16$) and repeat the same simulations with GENE, varying the number of grid points $N_{v_\parallel }$ (with $N_\mu = 16$). The results are shown in figure 1. The $T_R$ estimates for both cases agree with the analytical scalings. We note that the amplitudes of the oscillations due to the recurrence are smaller in the GM approach than in GENE because of the finite collisionality used in the GM calculations, which damps higher-order GMs. Finally, we remark that the analytical estimate of the collisionless ZF residual $\varpi$, defined in (4.1) is in agreement with the simulation results (see § 4.5).

Figure 1. Recurrence effects observed in the GM approach for increasing values of $P$ with $J = 16$ (a) and in GENE for increasing values of $N_{v_\parallel }$ with $N_\mu = 16$ (b). The normalized (in units of $R_0 / c_s$) recurrence times are estimated with $T_R \simeq \sqrt {2} {\rm \pi}q N_{v_\parallel }$ for GENE and $T_R \simeq 2 q \sqrt { P}$ for the GM simulations (see (3.3)) and are shown by the dashed coloured lines. The black dashed line represents the collisionless ZF residual $\varpi$ given in (4.1) (Rosenbluth & Hinton 1998). We note that the numerical hyperdiffusion along $z$ is set to zero in all cases. Here, the parameters are $\epsilon = 0.1$, $q = 1.4$ and $k_x = 0.05$.

Finally, we consider the perturbed ion distribution function during the GAM oscillations at $t \omega _G \simeq 10$ (with $\omega _G \sim q v_T / R_0$ the typical GAM frequency) and compare the ion perturbed distribution functions at the outboard midplane, $z = 0$, obtained from GENE and the GM approach in figure 2. For GENE simulations, we use $N_{v_\parallel } = 1024$ and $N_\mu = 16$, which yield $\lambda _{v_\parallel }^{{\rm min}} \simeq 0.003 v_{T}$. For the GM approach, we use $(P,J) = (256,16)$, therefore setting $\lambda _{ v_\parallel }^{{\rm min}} = {\rm \pi}v_{T} / \sqrt {2 P} \simeq 0.14 v_T$. We observe that at $t \omega _G \simeq 10$, the GM hierarchy is able to capture the main features of the $v_\parallel$ filamentation due to the parallel streaming of passing particles. Finally, we remark that the fine-scale structures in $v_\parallel$, such as the ones displayed in figure 2, are expected to be smeared out in nonlinear simulations due to resonant interactions, such as phase-mixing.

Figure 2. Modulus of the normalized (to the maximum) ion distribution function at the outboard midplane obtained with the GM approach with $(P,J) = (256, 16)$ (a) and using GENE with $( N_{v_\parallel }, N_\mu ) = (1024, 16)$ for reference (b) during the GAM oscillations shown in figure 1 at time $t \omega _G = 10$, which is before the recurrence time $T_R$ in both cases. The dashed blue line is the particle trapping boundary. The parameters are the same as in figure 1.

3.2. Effects of perpendicular magnetic drifts

Similarly to the parallel streaming of passing particles, the perpendicular drifts associated with the magnetic gradient and curvature frequency, $\omega _{Ba}$, drive resonance phenomena. Here, we consider the resonance driven by FOW effects also associated with $\omega _{Ba}$ and, more precisely, with the radial component of the perpendicular magnetic gradient drifts, $\boldsymbol{v}_{Da} \boldsymbol {\cdot } \boldsymbol {\nabla } x$, appearing in (2.1).

To analytically investigate the representation of FOW effects in the GM approach, we consider the collisionless time evolution of a radial perturbation, such that $\boldsymbol{k} = k_x \boldsymbol {\nabla } x$, in the absence of density and temperature gradients ($\omega _{Ta}^* = 0$) and neglect terms in (2.25) related to the parallel variation of $B$ (i.e. $\boldsymbol {b} \boldsymbol {\cdot } \boldsymbol {\nabla } B = 0$). Therefore, we focus on passing particles using concentric, circular, flux surface in the small inverse aspect ratio limit. In the electrostatic limit, multiplying the GK Boltzmann equation, (2.1), by the phase-factor ${\rm e}^{ {\rm i} \mathcal {Q} \cos z}$ with $\mathcal {Q} = \epsilon k_x \rho _p [ v_\parallel / v_T + \mu B v_T /(2 v_\parallel T)]$, $\rho _p = v_T / \varOmega _p$ being the poloidal gyroradius and $\varOmega _p = e B_p /m$ the poloidal gyrofrequency, yields an equation for the non-adiabatic response $h$

(3.4)\begin{equation} \left( \frac{\partial }{\partial t} + \frac{v_\parallel}{q R_0}\frac{\partial }{\partial z} \right) {\rm e}^{{\rm i} \mathcal{Q} \cos z} h = \frac{\partial }{\partial t} \left( {\rm e}^{{\rm i} \mathcal{Q} \cos z} \frac{e {\rm J}_0 \phi}{T} F_{M} \right). \end{equation}

We remark that the factor $\mathcal {Q}$, proportional to $\rho _p k_x$, is associated with FOW effects due to the radial drifts, $\boldsymbol {\nabla } x \boldsymbol {\cdot } \boldsymbol {v}_{Da}$, of passing particles.

In order to obtain the first insight on the impact of the FOW effects on the GM spectrum, we solve (3.4) by introducing the Fourier decomposition $h = \sum _l h_{l} {\rm e}^{{\rm i}l z -{\rm i} \omega t }$ and $e \phi / T = \sum _m \phi _{m} {\rm e}^{{\rm i}m z - {\rm i} \omega t}$. With the help of the Jacobi–Anger identity, ${\rm e}^{{\rm i} \mathcal {Q} \cos z} = \sum _n i^n {\rm J}_n(\mathcal {Q}) {\rm e}^{{\rm i} n z}$ (Gradshteyn & Ryzhik 2014), and evaluating the convolutions arising from the products of $z$-dependent quantities, such as ${\rm e}^{{\rm i} \mathcal {Q} \cos z} h$ and ${\rm e}^{{\rm i} \mathcal {Q} \cos z} \phi$, (3.4) can be solved for $h_m$, obtaining

(3.5)\begin{equation} h_m = \sum_{l,l'} i^{l' -l} {\rm J}_l(\mathcal{Q}) {\rm J}_{l'}(\mathcal{Q}) \frac{\omega}{\omega - v_\parallel( m +l)/(q R_0)} {\rm J}_0(b \sqrt{x})\phi_{m + l -l'} F_{M}. \end{equation}

Projecting $g_m = \int \,{\rm d} z g {\rm e}^{- {\rm i} m z - {\rm i} \omega t}$ with $h_m$ expressed by using (3.5) onto the Hermite–Laguerre basis yields the collisionless expression of the Fourier component of the GM of $g_m$, i.e. $N^{pj}_m = \int \,{\rm d} z N^{pj} {\rm e}^{ - {\rm i} z m}$, given by

(3.6)\begin{equation} N^{pj}_m ={-} \mathcal{K}_{j}(b) \delta_p^0 \phi_m + \sum_{l,l'} i^{l' -l} \phi_{m + l -l'} \frac{{\rm I}_{ll'm}^{pj}}{\sqrt{2^p p!}}, \end{equation}

having defined the resonant velocity-space integral

(3.7)\begin{equation} I_{ll'm}^{pj} = \frac{1}{\sqrt{\rm \pi}} \int_{- \infty}^\infty \,{\rm d} s_\parallel \int_0^\infty \,{{\rm d}\,x} {\rm J}_l(\mathcal{Q}) {\rm J}_{l'}(\mathcal{Q}) \frac{\omega {\rm e}^{- s_\parallel^2 - x}}{\omega - v_\parallel( m +l)/(q R_0)} H_p(s_\parallel) L_j(x) {\rm J}_0(b\sqrt{x}). \end{equation}

While a closed analytical expression of the resonant integral $I_{ll'm}^{pj}$, given in (3.7), can be obtained in terms of generalized plasma dispersion relations by following Frei et al. (2022b) and be evaluated using numerical algorithms (Gürcan 2014)), this is rather complex and outside the scope of the present work. Instead, we focus here on physical insights on FOW effects that can be obtained directly by the inspection of the analytical form of the integral $I_{ll'm}^{pj}$. We first observe that FLR (of the order of $b$) and FOW (of the order of $\epsilon k_x \rho _p \sim q b$) effects can be neglected in $I_{ll'm}^{pj}$ in the long radial wavelength limit $k_x \ll 1$, since ${\rm J}_0(b \sqrt {x}) \sim 1$, ${\rm J}_l(\mathcal {Q}) \sim 1$ for $l =0$, and ${\rm J}_\ell (\mathcal {Q}) \sim 0$ for $l \neq 0$. In the same limit, the resonant term contributes to the GMs throughout the $j=0$ term because of the Laguerre orthogonality relation given in (2.21b). On the other hand, when $k_x \rho _p \sim 1$ (but $k_x \rho _s \ll 1$), FOW effects drive $j > 0$ GMs because of the $\mu$ dependence of $\mathcal {Q}$ in the arguments of ${\rm J}_l(\mathcal {Q})$ and the presence of Laguerre polynomials $L_j$ with $j >0$, that couples the Fourier harmonic $l$. As $k_x \rho _p \gtrsim 1$ and $k_x \rho _s \sim 1$, FLR effects drive GMs also through the $x$ dependence of ${\rm J}_0 ( b \sqrt {x})$ (Frei et al. 2022b).

We numerically illustrate the effects of resonance driven by FOW and FLR effects by evolving (3.4), i.e. by solving the GM hierarchy in (2.25) neglecting the background gradients ($R_N = R_{Ta} =0$) and the parallel gradient of the magnetic field $B$ ($\partial _z \ln B =0$), but retaining the parallel streaming of passing particles. In figure 3, we plot the modulus of the GM spectrum averaged over $z$, defined by

(3.8)\begin{equation} \left\langle \left| N_a^{pj} \right| \right\rangle_z = \frac{ \int \,{\rm d} z \mathcal{J}_{xyz} \left| N_a^{pj} \right|}{\int \,{\rm d} z \mathcal{J}_{xyz} }, \end{equation}

obtained numerically during the GAM oscillations, which are an eigensolution of (3.4) (Sugama et al. 2006), at time $t \omega _G \simeq 2$ (see § 4.5) for different values of $k_x$. We evolve $(P,J) = (64,24)$ GMs. As $k_x$ increases, the GM spectrum broadens in both $p$ and $j$ directions since high-order GMs are driven by FOW and FLR effects. While the FOW contributes with the parallel streaming in the Hermite GMs because of the $s_\parallel$ dependence in $y$ associated with the curvature drift, the increased broadening in Laguerre direction with $k_x$ is associated with the FLR and $\boldsymbol {\nabla } B$ drift yielding the $x$ dependence in $y$. We remark that the same broadening mechanism of the GM spectrum was identified in the case of toroidal ITG (Frei et al. 2022b).

Figure 3. Normalized (to the maximum value) GM spectrum for $k_x = 0.05$ (a), $k_x =0.5$ (b) and $k_x =1$ (c) during the GAM oscillation at a time $t \omega _G \simeq 2$. The GM spectrum is represented on a logarithmic scale and artificially saturated for visualization purposes. Here, we consider $q = 1.4$, $\epsilon = 0.1$.

4. Collisionless microinstability and comparison with GENE

We now turn to the investigation of the collisionless properties of microinstabilities using the GM approach. In particular, we focus on the linear study of the ITG, TEM, KBM and MTM and consider also the dynamics of GAM and ZFs. We perform a detailed comparison with the continuum GK code GENE, which uses a finite difference method in velocity space and the same velocity-space coordinates as the GM approach. The linear growth rates, real mode frequencies, ballooning eigenmode structures, and the associated velocity-space structures are compared with GENE results as a function of the number $(P,J)$ of GMs. We find that the GM approach is in excellent agreement with GENE, and that convergence is most often achieved with a number of GMs of the same order as the number of grid points used in GENE, i.e. $P \sim N_{v_\parallel }$ and $J \sim N_\mu$, despite the presence of strong kinetic features (see § 3). Interestingly, we find that a small number of GMs is needed for convergence for pressure gradients driven mode (such as the KBM), while it is increased when sharp gradients in the distribution functions appear (e.g. in the TEM). The present section provides a verification of the GM approach, which is shown to be able to represent the collisionless limit of the essential microinstabilities that are responsible for the anomalous turbulent transport in the boundary of fusion devices.

The present section considers tests of increasing complexity. In § 4.1, we first perform the ITG cyclone base case test with adiabatic electrons (Dimits et al. 2000). Then, in § 4.2, we illustrate the transition from the ITG mode to the TEM by introducing kinetic electrons in our model, focusing on the electrostatic limit. Electromagnetic effects are then considered, studying the KBMs in § 4.3 and the MTMs in § 4.4. Finally, we study the collisionless GAM and ZF dynamics in § 4.5. In Appendix A, as a further collisionless study, we focus on the local and strong ballooning limit of the flux-tube model, allowing us to derive analytically an electromagnetic GK dispersion relation, which we compare with the solution of the GM approach in the same limit.

4.1. Cyclone base case with adiabatic electrons

As a first linear collisionless test, we consider the electrostatic ITG cyclone base case scenario with adiabatic electrons (Dimits et al. 2000). The cyclone base case is widely used to verify GK codes (Merlo et al. 2016; Tronko et al. 2017). In the cyclone base case scenario, the safety factor, magnetic shear and inverse aspect ratio are fixed at $q = 1.4$, $s = 0.8$, and $\epsilon = 0.18$, respectively. Additionally, we set the MHD parameter $\alpha = 0$ also for the rest of the present work, if not mentioned otherwise. Physical dissipation in the GMs is introduced by using the GK Dougherty collision operator (Frei et al. 2022b) with a small but finite value of collisionality ($\nu _{ei} = \nu _{ii} = 10^{-4}$). The ion density and temperature gradients are $R_N = R / L_N = 2.22$ and $R_{T_i} =R / L_{T_i} = 6.9$, corresponding to a value of $\eta = L_N / L_{T_i} \simeq 3$, which is above the ITG mode linear threshold. We choose $N_{k_x} = 5$ and $N_z = 24$. In addition to GENE, we compare our results with the GX code (Mandell et al. 2022), which uses a similar polynomial decomposition as the one used in this work. If not indicated, we use a high velocity-space resolution of $(N_{v_\parallel }, N_\mu ) = (128,24)$ in GENE as a reference.

The ITG growth rate, $\gamma$ (normalized to $c_s /R_0$), is plotted in figure 4 as a function of the binormal wavenumber $k_y$ (normalized to the ion sound Larmor radius $\rho _s$) for different temperature gradients $R_{T_i}$. Different number of GMs, $(P,J)$, are considered also for the GX code. First, we remark that our results coincide with GX for all values of $(P,J)$. In addition, both spectral velocity-space codes agree well with the GENE code when $(P,J) \gtrsim (32,16)$. Second, we note that the GM approach provides a better estimate of the ITG growth rate at long wavelength, even when low values of $(P ,J)$ are used, showing that FOW and FLR effects require a large number of Laguerre GMs for their description. This is needed for the gyro-averaging, as one can infer from (2.23) (Frei et al. 2022b).

Figure 4. The ITG growth rate $\gamma$ and real mode frequency $\omega _r$ as a function of the binormal wavenumber $k_y$ for various ion temperature gradients $R_{T_i}$. Different numbers $(P,J)$ of GMs are considered, and the results are compared with the continuum GK code GENE (red lines) and pseudo-spectral code GX (light coloured lines) (Mandell et al. 2022).

Finally, we perform the ballooning transformation, given in (2.10), to compare the ballooning eigenmode function $\phi _B$, as obtained from the GM approach and from GENE. These are plotted in figure 5. We observe that the functions $\phi _B$ are in good agreement, peaking at the outboard midplane position. The inspection of the normalized GM spectrum, defined in (3.8) and also shown in figure 5, reveals that the velocity space is indeed well resolved with $(P,J) = (32,16)$. Finally, we observe that convergence is achieved when $P > J$, a situation typically found in all cases discussed in the present paper.

Figure 5. Real part (blue lines), imaginary part (red lines) and modulus (black lines) of the ballooning eigenmode function $\phi _B(\chi )$ normalized to $\phi _B(0)$ (a), obtained using the GM (solid lines) and GENE (dashed lines). Normalized GM spectrum for the $k_x =0$ and $k_x = \pm 2 {\rm \pi}s k_y$ modes is plotted on (b,c). The logarithmic scale is artificially saturated. Here, $R_{T_i} = 6$, ${k_y = 0.3}$ and adiabatic electrons are considered.

4.2. Ion temperature gradient and TEM

We now introduce the trapped and passing electron dynamics allowing us to investigate the transition between the ITG and TEM. The electron dynamics introduces fast waves such as the high-frequency wave, $\omega _H^2 = (k_\parallel ^2 / k_\perp ^2) (m_i / m_e) \varOmega _i^2$ (Lee 1987; Lin et al. 2007), that can limit the explicit time stepping scheme. While using a reduced ion mass can limit the non-adiabatic electron response (Dominski et al. 2015), we consider $\sqrt {m_i/ m_e} \simeq 19.24$ for numerical reasons. The reduced ion mass used in this work is sufficient to capture the main features of the non-adiabatic electron response to investigate the GM convergence. Due to the localized and fine radial structures associated with the non-adiabatic electron response (Hallatschek & Dorland 2005), we evolve a larger number of radial modes (i.e. $N_{k_x} = 11$), and increase the number of parallel grid points to $N_z = 24$ to properly resolve the tails. We use the same resolution in GENE. Electromagnetic effects are neglected in this section.

The growth rate and real mode frequency of the most unstable mode are shown in figure 6 as a function of the binormal wavenumber $k_y$, using the same parameters as in figure 4 and considering a finite electron temperature gradient, $R / L_{T_e} = R / L_{T_i} = 6.96$. The GM approach agrees with GENE at high velocity-space resolution for all wavelengths, when roughly the same number of GMs as number of grid points, i.e. $(P,J) \sim (N_{v_\parallel }, N_\mu ) = (32,16)$, are used. A transition from ITG to TEM is identified near $k_y \simeq 0.5$ when the mode propagation changes from the ion ($\omega _r > 0$) to electron ($\omega _r < 0$) diamagnetic direction.

Figure 6. The ITG and TEM growth rate $\gamma$ (a) and real mode frequency $\omega _r$ (b) as a function of the binormal wavenumber $k_y$ for different values of $(P,J)$ (circle makers). GENE simulations are shown by the cross markers for different resolutions $(N_{v_\parallel }, N_\mu )$. The dashed line in the right panel corresponds to the ion diamagnetic direction for $\omega _r > 0$ and to the electron diamagnetic direction for $\omega _r < 0$.

The effects of the electron dynamics are illustrated by investigating the modulus of the electrostatic ballooning eigenmode function $\phi _B$ (see (2.10)). We consider the same parameters as in figure 6 and ${k_y = 0.3}$ at different ballooning angles, $z_0 = - \delta k_x / s k_y$, and show the results in figure 7 using $(P,J) = (32, 16)$ and GENE (we also show the GENE results with a realistic mass ratio for deuterium plasmas). First, we observe that extended tails in the mode envelope of $\phi _B$ are present and are associated with the non-adiabatic response of passing electrons (Dominski et al. 2015; Ajay, Brunner & Ball 2021). Second, while the mode at $\delta k_x = 0$ and $\delta k_x = 0.1$ is identified as ITG, a transition to TEM is observed at $\delta k_x \gtrsim 0.2$ at $k_y \gtrsim 0.3$, in contrast to the ITG–TEM transition occurring at $k_y \gtrsim 0.5$ with $\delta k_x = 0$ in figure 6. An excellent agreement is observed with GENE at the outboard midplane ($\chi = 0$), where the most unstable part of the mode is localized, while the small differences that appear in the tails, near $\chi / {\rm \pi}\gtrsim 2$, in the case of the TEM (${\delta k_x = 0.2}$) are attributed to numerical reasons (Merlo et al. 2016), as confirmed by increasing the number of grid points, $N_z$, and the number of radial modes, $N_{k_x}$. On the other hand, the value of the parallel diffusion used has little effects on the results.

Figure 7. Modulus of the electrostatic ballooning eigenmode function $\phi _B(\chi )$ (normalized to $\phi _B(0)$) obtained using the GM approach with $(P,J) = (32,16)$ (dashed black lines) and using GENE (solid red lines with $\sqrt {m_i/ m_e} \simeq 19.24$ and solid blue lines with a realistic mass ratio for deuterium plasmas, i.e. $\sqrt {m_i/ m_e} \simeq 60.59$) for increasing values of $\delta k_x$ (from a–c). We consider an ITG mode ($\delta k_x =0$ and $\delta k_x =0.1$) and a TEM ($\delta k_x =0.2$). The $\chi$ range considered for the numerical solution is truncated for visual reasons. Here, the same parameters as figure 6 are used, except ${k_y = 0.3}$.

To investigate the presence of velocity-space structures driven by, e.g. trapped particles, we compare in figure 8 the modulus of the deviation of the electron distribution function, ${g}_e$, from a Maxwellian, which is proportional to the non-adiabatic distribution function $h_e$ (see (2.2)), as obtained using GENE and the GM approach with $(P,J) = (32,16)$. We focus on the case of the ITG mode (at ${k_y = 0.3}$) and of the TEM (at $k_y = 1.3$) at the outboard midplane ($z =0$ and $k_x =0$). While a good qualitative agreement is found in the ITG case, larger deviations are observed in the TEM case in particular near $s_{\parallel e} = v_{\parallel } / v_{Te}= 0$ and along the trapped and passing boundary (shown by the dashed blue lines) where a strong gradient is observed in the GENE case. The deviations between GENE and the GM approach are also visualized on the right panels of figure 8, where the distribution functions ${g}_e$ are plotted as a function of $x_e$ at $s_{\parallel e} = 0$. While $(P,J) = (32,16)$ is in good agreement with GENE for the ITG case, differences remain at $x_e \gtrsim 2.5$ between GENE and the GMs for the TEM case, despite the convergence in the growth rate with $(P,J) = (32,16)$ (see figure 6). These deviations are associated with the finite number of GMs used in our simulations. In fact, the effects of unresolved GMs can be investigated by considering the normalized electron GM spectrum, $|N_e^{pj}|$, associated with the distribution displayed in figure 8 and plotted in figure 9. As observed, the GM spectrum fills the whole space and decays only by two orders of magnitude in the Hermite direction going from $p=0$ to $p=32$, highlighting the presence of fine structures along $v_\parallel$ in both ITG and TEM. Also, we notice that the decay in the Laguerre direction $j$ is faster in the ITG than in the TEM case, explaining the different levels of deviation observed in the right panel of figure 8. The effects of the magnetic gradient drifts, associated with the ${\rm i} \omega _{Ba}$ term in (2.1), can also be identified by the band-like structures in the GM spectrum of both cases (Frei et al. 2022b). However, despite the presence of underresolved velocity-space structures by the GM approach, convergence of the growth rate is achieved in figure 6 with $(P,J) \sim (32, 16)$.

Figure 8. Deviation of the distribution from a Maxwellian, $|{g}_e| - F_M$, at the outboard midplane for to the ITG mode at ${k_y = 0.3}$ (a–c) and of the TEM at $k_y = 1.3$ (d–f), obtained using GENE (a,d) and the GM approach with $(P,J) = (32,16)$ (b,e). The trapped and passing boundary is shown by the dashed blue lines. The modulus of distribution function ${g}_e$ along $s_{\parallel e} = 0$ is also shown (c, f) for different values of $(P,J)$ and GENE. The same parameters as in figure 7 are used.

Figure 9. Modulus of the normalized electron GM spectrum associated with the ITG (a) and with the TEM mode (b) plotted on a logarithmic scale (colour bars are artificially saturated at $10^{-5}$). The same parameters as in figure 8 are used.

Finally, we focus on the case of a TEM developing at long perpendicular wavelengths. This instability appears when the ion temperature gradient is below the ITG linear threshold. More precisely, we evaluate the growth rate and real mode frequency of the most unstable mode as the normalized ion temperature gradient, $R_{T_i}$, is varied at fixed binormal wavenumber and density and electron temperature gradients, i.e. $k_y = 0.25$, $R_N = 3$ and $R_{T_e} = 4.5$. The results are shown in figure 10, where the TEM mode ($\omega _r < 0$) is observed for $R_{T_i} < 4$ and the ITG mode is the most unstable mode when $R_{T_i} \gtrsim 4$ ($\omega _r > 0$). While convergence is achieved with $(P, J) = (32,16)$ for the ITG mode (when $R_{T_i} \gtrsim 4$), a larger number of GMs is required for the TEM at weaker $R_{T_i}$, i.e. $(P,J) = (128,24)$. The number of GM needed for convergence is therefore even larger than the TEMs appearing at larger $k_y$ (see figure 6). We remark that achieving convergence in GENE requires approximately $(N_{v_\parallel }, N_\mu ) \gtrsim (64, 16)$. We notice that the real mode frequency, $\omega _r$, is less sensitive to the resolution in velocity space. The lack of convergence of the GM approach in the case of TEM at $k_y = 0.25$ is explained by the presence of sharp velocity-space gradients that occur near the trapped and passing boundary, a feature stronger than the one developing at $k_y = 1.3$ (see figure 8).

Figure 10. The ITG and TEM growth rate $\gamma$ (a) and frequency $\omega _r$ (b) as a function of the ion normalized temperature gradient, $R_{T_i}$, for $k_y = 0.25$ and different values of $(P,J)$. GENE results are shown by the cross markers.

4.3. Kinetic ballooning modes

We now turn to collisionless microinstabilities appearing when electromagnetic effects are considered. While electromagnetic effects are known to be most often stabilizing (Weiland & Hirose 1992; Citrin et al. 2014), they can trigger the KBM if the electron plasma beta, $\beta _e = 8 {\rm \pi}N T_e / B_0^2$, is above a certain threshold (Connor et al. 1978; Tang, Connor & Hastie 1980; Aleynikova & Zocco 2017).

The KBM mode is an ideal MHD mode resulting from the interplay between pressure gradients, magnetic curvature, and field line bending, modified by kinetic effects. This mode typically develops at long parallel wavelengths and perpendicular wavelengths of the order of the ion gyroradius, $k_y \rho _i \lesssim 1$ (Belli & Candy 2010). To study the KBM, we consider the parameters $R_N = 3$, $R_{T_e} = 4.5$, $R_{T_i} = 8$ and $k_y = 0.25$, solving the GM hierarchy equation, (2.25), coupled to the GK Ampere's law expressed in terms of GMs given in (2.29) in addition to the GK quasineutrality condition in (2.28). A scan over $\beta _e$ is performed for various $(P,J)$. The results are displayed in figure 11 and are compared with GENE at different velocity-space resolutions. We first observe a discontinuous jump in the mode frequency, $\omega _r$, near $\beta _e \simeq \beta _e^c = 0.012$, corresponding to the transition between the KBM and ITG modes, which are stabilized by electromagnetic effects. We remark that the value of $\beta _e^c$ in figure 10 is less than $5\,\%$ smaller with respect to the linear threshold derived from fluid MHD theory, i.e. $\beta _e^{{\rm MHD}}$, where the kinetic effects are neglected. Second, while the GM approach requires a number of GMs of the same order as the number of grid points used in GENE in the case of the ITG mode, i.e. $(P,J) \gtrsim (32,16)$, the KBM mode is well described by fewer GMs, i.e. $(P,J) \gtrsim (16,8)$, a number of GMs smaller than the number of grid points necessary in GENE to achieve convergence.

Figure 11. The ITG and KBM growth rate $\gamma$ (a) and real mode frequency $\omega _r$ (b) as a function of $\beta _e$ for different values of $(P,J)$ (circle markers) compared with the GENE results (cross markers) for different values of $(N_{v_\parallel }, N_\mu )$. The ideal MHD threshold of $\beta _e^{{\rm MHD}} = 0.6 s/[q_0^2 (2 R_N + R_{T_e} + R_{T_i})] \simeq 0.0132$ is shown by the vertical dotted-dashed lines.

The low resolution of the GM approach in the case of KBM can be explained by the fact that the KBM presents fewer fine-scale structures of the distribution function compared with the ITG and TEM cases (see figure 12). Also, we observe that the GM spectrum is well resolved, contrary the ITG and TEM cases shown in figure 9. The case of the KBM mode in figure 12 exemplifies the small number of GMs often required for pressure gradient-driven modes, with kinetic effects playing a minor role.

Figure 12. Modulus of ${g}_e$ (normalized to its maximum) at the outboard midplane in the case of the KBM for $\beta _e = 0.03$ (see figure 11) obtained using GENE (a) and using $(P,J) = (32,16)$ GMs (b), with the corresponding modulus of the normalized electron GM spectrum (c).

Finally, we investigate the ballooning eigenmode function associated with the perturbed magnetic vector potential, $\psi$. We plot the ballooning eigenmode function $\psi _B$ (see (2.10)) for the KBM mode developing at $\beta _e = 0.03$, with $(P,J) = (32,16)$, and compare it with GENE in the left panel of figure 13. The KBM mode is characterized by the ballooning-parity, such that $\psi _B$ is anti-symmetric around the outboard midplane located at $\chi = 0$, i.e. $\psi _B( -\chi ) = - \psi _B(\chi )$, while the electrostatic potential eigenmode function, $\phi _B$, is symmetric (but not shown). A good agreement in the perturbed magnetic potential $\psi$ is observed between the GM approach and GENE.

Figure 13. Real (blue) and imaginary (red) parts of the ballooning eigenmode function $\psi _B$ (normalized to the electrostatic potential $\phi _B(0)$) in the case of KBM mode when $\beta _e = 0.03$ (a) and in the case of MTM at ${k_y = 0.3}$ (b) obtained using GENE (dotted lines) and the GM approach with $(P,J) = (32, 16)$ (solid lines). The same parameters as in figures 11 and 14 are used respectively. The $\chi$ range is truncated for visual reasons.

4.4. Microtearing modes

As a final collisionless microinstability investigated using the GM approach, we consider the MTMs, which are driven unstable at finite $\beta _e$ values if the electron temperature gradient is above a linear threshold (Dickinson et al. 2012). More precisely, MTMs are usually driven unstable by a combination of finite electron temperature and collisionality (even small) in the core region (Catto & Rosenbluth 1981). MTMs also exist in the edge region in the collisionless limit, driven unstable by the electron magnetic drift resonance effects (Applegate et al. 2007; Dickinson et al. 2013).

Here, we focus on MTMs appearing in edge conditions because of the role of electron magnetic drift resonance effects that often require a larger number of GMs (see figure 9) and the fact that it persists at a vanishing value of collisionality, in contrast to core MTMs. We consider a safety factor $q =4$, a magnetic shear $s = 2.4$, gradients of density and electron temperature $R_N = 3$ and $R_{Te} =8$, respectively, and an electron plasma beta of $\beta _e = 0.02$, above the linear thresholds for the MTM onset. While the ion kinetic response is ignored in previous linear MTM studies (see, e.g. Dickinson et al. 2013), we include it but neglect gradients in the ion temperature, i.e. $R_{Ti} = 0$. In contrast to the core MTMs that are extended along the parallel direction, the ballooning MTM eigenmode structure is considerably less elongated at the higher safety factor and larger shear of the edge. Therefore, we use $N_{k_x} = 11$ and $N_z = 64$.

A scan over the binormal wavenumber, $k_y$, is shown in figure 14 for different numbers of GMs and with results of GENE. First, we remark that a good agreement is found with GENE when $(P,J) \gtrsim (32,16)$. Second, the MTM growth rate peaks near ${k_y = 0.3}$, while the real mode frequency increases in magnitude linearly with the electron diamagnetic frequency, i.e. $\omega _r \sim \omega _e^*$. Third, a larger number of GMs is required to achieve convergence compared with the KBM case and that number increases with $k_y$, which is a consequence of the role of the electron magnetic drift motion (proportional to ${\rm i} \omega _{Be}$ in (2.1)) in the collisionless destabilization mechanism of MTMs (Doerk et al. 2012; Dickinson et al. 2013) (see § 3.2). In contrast to KBMs, MTMs are characterized and identified by a tearing parity where $\psi _B$ is even around the outboard midplane position, i.e. $\psi _B(-\chi ) = \psi _B(\chi )$, while $\phi _B$ is odd. The ballooning eigenmode function, $\psi _B$, in the case of the MTM at ${k_y = 0.3}$ is shown on the right panel of figure 13, revealing its tearing parity and in excellent good agreement with GENE.

Figure 14. The MTM growth rate $\gamma$ (a) and real mode frequency $\omega _r$ (b) as a function of $k_y$ for different values of $(P,J)$ (circle markers) with the GENE results (cross markers) for different values of $(N_{v_\parallel }, N_\mu )$.

The role of the electron magnetic drift motions in the MTM destabilization mechanism is visualized by considering the electron distribution function and its GM spectrum, both displayed in figure 15. While a good agreement between the electron distribution functions obtained using GENE and the GM approach is observed, the effects of electron magnetic drifts can be identified by the presence of band-like structures that extends in the Laguerre direction in the GM spectrum (Frei et al. 2022b). This explains the broad GM spectrum observed in the MTM simulations compared with the KBM case displayed in figure 12.

Figure 15. Modulus of ${g}_e$, (normalized to its maximum) for the MTM at ${k_y = 0.3}$ obtained using GENE (a) and with $(P,J) = (32,16)$ (b) with the modulus of the normalized electron GM spectrum $|N_e^{pj}|$ (c).

4.5. Collisionless GAM dynamics and ZF damping

As a final collisionless test, we consider the time evolution of an initial seeded and radially dependent density perturbation without equilibrium pressure gradients and with adiabatic electrons. The initial density perturbation creates a perturbed poloidal flow rapidly evolving into poloidally non-symmetric and radially localized oscillations, associated with GAMs (Winsor et al. 1968). GAMs are damped by collisionless processes, such as parallel streaming and FOW effects due to passing particles (see § 3).

To investigate the collisionless GAM dynamics, we consider $q = 1.4$, $\epsilon = 0.1$ and $s=0$. We simulate the time evolution of the flux-surface-averaged electrostatic potential, $\left \langle \phi \right \rangle _{fs}$, by considering an initial perturbed density with a radial wavenumber $k_x = 0.01$. Because of the fine velocity-space structures associated with GAMs (see § 3.1), we use a large number of GMs, i.e. $(P,J) = (800,16)$ and a small but finite collisionality to limit the effects of the recurrence avoiding the use of artificial velocity-space hyperdiffusion (collisions do not significantly affect the GAM dynamics in the banana regime, $\nu _i^* \lesssim 1$ (see § 5.3). We compare our numerical results with the analytical time prediction derived in Hinton & Rosenbluth (1999), as well as with the damping rate and frequency, $\gamma _G$ and $\omega _G$, given in Sugama et al. (2006). The results are plotted in figure 16 where a GENE simulation is also shown for comparison. The GAM oscillations are in good agreement with the analytical predictions, as well as with GENE simulations. The GAM damping $\gamma _G$ and frequency $\omega _G$, computed numerically by fitting the time trace of figure 16 with the model $\phi _z(t) / \phi _z(0) - \varpi \simeq A \cos (\omega _G t) \exp (- \gamma _G t)$ (with $A$ a fitting constant), are compared with GENE as a function of the parallel velocity resolutions (i.e. as a function of $P$ and $N_{v_\parallel }$) at various low collisionality in the banana regime. A good agreement is observed for the GAM damping in the banana regime with the GENE results. Finally, we remark that the convergence of the GM approach improves with collisionality, consistent with previous studies (Frei et al. 2021, 2022b).

Figure 16. (a) Comparison of the time evolution of $\left \langle \phi \right \rangle _{fs}(t) / \left \langle \phi \right \rangle _{fs} (0)$ between GENE with $(N_{v_\parallel }, N_\mu ) = (128, 24)$ (red solid line with markers) and the GM approach with $(P ,J) = (800,16)$ (cyan solid line) in the banana regime ($\nu _i^* = 0.003$). The collisionless analytical time evolution (black dotted) is obtained from the Hinton–Rosenbluth analytical results (Hinton & Rosenbluth 1999), i.e. $\left \langle \phi \right \rangle _{fs}(t) / \left \langle \phi \right \rangle _{fs}(0) \simeq (1 - \varpi ) \exp ( - \gamma _G t) \cos (\omega _G t) + \varpi$, with $\gamma _G$ and $\omega _G$ obtained from Sugama et al. (2006) and the collisionless residual $\varpi$ defined in (4.1) (solid black line). (b) Convergence of $\gamma _G$ as a function of the number of parallel grid points $N_{v_\parallel }$ ($N_\mu = 24$) for GENE (dashed lines) and as a function of $P$ ($J=18$) for the GMs (solid lines) at different banana collisionalities. Here, $q = 1.4$, $\epsilon = 0.1$ and $k_x = 0.01$.

Following the damping of the GAM oscillations, a non-vanishing residual is observed, known as the ZF residual. The ZFs are axisymmetric and primarily poloidal flows that play an important role in saturating turbulence (Diamond et al. 2005). Rosenbluth & Hinton (1998) show that the ratio of the flux-surface-averaged electrostatic potential, $\left \langle \phi \right \rangle _{fs}(t)$, to its initial value, $\left \langle \phi \right \rangle _{fs}(0)$, converges to a non-vanishing residual level approximated by

(4.1)\begin{equation} \frac{ \left\langle \phi \right\rangle_{fs}(\infty)}{ \left\langle \phi \right\rangle_{fs}(0)} \to \varpi = \frac{1}{1 + q^2 \varTheta /\epsilon^2 }, \end{equation}

where the numerical factor $\varTheta =1.635 \varepsilon ^{3 / 2}+0.5 \varepsilon ^{2}+0.36 \varepsilon ^{5 / 2}$ is derived in Xiao & Catto (2006) including higher-order terms in the small inverse aspect ratio $\epsilon$. The analytical prediction of the collisionless ZF residual, given in (4.1), is obtained by assuming concentric and circular flux surfaces in the $\epsilon \ll 1$ limit and a perpendicular wavelength longer than the ion gyro-radius, $k_x \ll 1$. Equation (4.1) is confirmed by a number of GK codes (Merlo et al. 2016), in contrast to gyrofluid models (see, e.g. Beer & Hammett 1996). Indeed, gyrofluid models evolve a considerably smaller number of moments than the calculations shown in figure 17 and they use closures based on consideration of the properties of linear instabilities, yielding an artificial damping of the ZF residual (Rosenbluth & Hinton 1998). In order to compare our numerical results with (4.1), we average the simulated ZF residual over a time window that extends from a time $t$ to a time $t + \tau$ (with $t \gg 1 / \gamma _G$ and $\tau \sim 20$). We show the time-averaged ZF residual of $\left \langle \phi \right \rangle _{fs}(\infty ) / \left \langle \phi _z \right \rangle _{fs} (0)$ as a function of $\epsilon$ in figure 17 obtained from the GM approach with $(P,J) = (128,16)$. We observe that the time-averaged collisionless ZF residual agrees well with the analytical prediction $\varpi$ given in (4.1). This confirms that the GM approach can correctly reproduce the collisionless ZF damping process even with a simple closure by truncation, in contrast to previous gyrofluid models.

Figure 17. Time-averaged collisionless ZF residual as a function of the inverse aspect ratio, $\epsilon$, obtained with $(P,J) = (128,16)$ GMs (red markers). The solid black line is the analytical prediction $\varpi$ given in (4.1). The same parameters as in figure 16 are used.

5. High-collisional limit and collisional effects on microinstabilities

While collisional effects are often neglected in the core, they can no longer be ignored near the separatrix and in the SOL region because of the low plasma temperatures. For instance, $\nu _e^* \sim 0.03$ is expected at the top of the pedestal in ITER and $\nu _e^* \gtrsim 50$ near the separatrix. The high collisionality significantly affects the linear properties of edge microinstabilities, in particular the TEMs and MTMs which have been identified to play a major role in the pedestal region (Fulton et al. 2014; Hatch et al. 2016; Garcia et al. 2022).

We study the collisional dependence of TEMs and MTMs using the GM approach in this section. In particular, we consider advanced collision operator models, such as the Coulomb, the Sugama and the IS collision operators (Frei et al. 2021, 2022a). Our results confirm that the IS operator approaches the Coulomb operator better than the Sugama operator in the high-collisional Pfirsch–Schlüter regime (Frei et al. 2022a), while the Sugama operator often underestimates the linear growth rates when FLR terms in the collision operator cannot be ignored. Consistently with previous works (Pan et al. 2020, 2021) and using the GM approach, we show that the presence of FLR collisional terms yields a stabilization of the TEM and MTM modes at high collisionality and that the accuracy (relative to the Coulomb operator) of collision operator models depends on physical parameters such as, e.g. the electron temperature gradient. In addition, we show that a high-collisional reduced GM model is able to capture the main trend of the TEM and MTM linear growth rates in the Pfirsch–Schlüter regime. Finally, because the GAMs and ZFs are often observed in the edge region, we also assess the effect of collisions and the choice of collision operators on their dynamics (Pan et al. 2020, 2021).

The present section is structured as follows. In § 5.1, we first derive the high-collisional limit of the GM flux-tube model, yielding a reduced high-collisional $6$GM model. Second, we investigate the collisionality dependence of TEMs and of the MTMs in typical edge parameters, from the banana (e.g. top of H-mode pedestals) to the Pfirsch–Schlüter collisionality regimes (e.g. the bottom of pedestal and SOL) in § 5.2. Finally, we study the collisional effects on the GAM dynamics and on the ZF damping in §§ 5.3 and 5.4, respectively.

5.1. Linear high-collisional limit

To consider the high-collisional limit, we assume $N_a^{00} \sim N_a^{10} \sim N_a^{01} \sim N_a^{20}$ and $N_a^{30} \sim N_a^{11} \sim \epsilon _\nu N_a^{00}$ (Jorge et al. 2017; Frei et al. 2022b), where $\epsilon _\nu \ll 1$ being the ratio of the electron mean-free path to the typical parallel scale length (Chapman & Cowling 1941). We also neglect all higher-order GMs with $p + 2j > 3$. Evaluating the GM hierarchy equation, (2.25), with $(p,j) = (0,0)$, $(1,0)$, $(2,0)$ and $(0,1)$, we obtain the evolution equations for the lowest-order GMs associated with the perturbed gyrocentre density $N_a$, parallel velocity $u_{ \parallel a}$, parallel and perpendicular temperatures $T_{\parallel a}$ and $T_{\perp a}$, respectively. Finally, considering $(p,j) = (3,0)$ and $(1,1)$, we obtain the evolution equations for the parallel and perpendicular heat fluxes, $Q_\parallel$ and $Q_\perp$. The evolution equations are derived using the relations between the GMs and the fluctuations of the gyrocentre fluid quantities, $N_a = N_a^{00}$, $u_{\parallel a} = v_{Ta} N_a^{10} / \sqrt {2}$, $T_{\parallel a} / T_a= \sqrt {2} N_a^{20} + N_a$ and $T_{\perp a} / T_a= N_a - N_a^{01}$ (Frei et al. 2020). Assuming the MHD parameter $\alpha =0$, these equations are given in physical units by

(5.1a)\begin{align} & \frac{\partial }{\partial t} N_a + \boldsymbol{\nabla}_\parallel u_{{\parallel} a}^\psi - u_{{\parallel} a}^\psi \boldsymbol{\nabla}_\parallel \ln B + \frac{{\rm i} R_B }{q_a B} \left( T_{{\parallel} a} + T_{{\perp} a} + q_a ( 2 \mathcal{K}_{0} - \mathcal{K}_{1}) \phi \right) \nonumber\\ & \quad + {\rm i} \left( \mathcal{K}_{0} \omega_N - \omega_{T_a} \mathcal{K}_{1} \right) \frac{e \phi}{T_e} = 0, \end{align}

(5.1b)\begin{align} & m_a \frac{\partial }{\partial t} u_{{\parallel} a} + \boldsymbol{\nabla}_\parallel T_{{\parallel} a} + q_a \boldsymbol{\nabla}_\parallel \left(\mathcal{K}_{0} \phi \right) - (T_{{\parallel} a} - T_{{\perp} a} + q_a \mathcal{K}_{1} \phi) \boldsymbol{\nabla}_\parallel \ln B \nonumber\\ & \quad + \frac{{\rm i} m_a R_B}{2 \varOmega_a} \left( Q_{{\parallel} a} + 4 v_{Ta}^2 u_{{\parallel} a}^\psi - Q_{{\perp} a} + \frac{2 T_a \varOmega_a}{m_a} \mathcal{K}_{1} \frac{\psi}{B}\right) \nonumber\\ & \quad - {\rm i} \frac{\sqrt{2} e T_a }{m_a T_e} \left( \frac{\mathcal{K}_{0}}{\sqrt{2}} \omega_N + \omega_{T_a} ( \mathcal{K}_{0} - \mathcal{K}_{1}) \right) \psi = \mathcal{C}_{a}^{10}, \end{align}

(5.1c)\begin{align}& \frac{1}{T_a} \frac{\partial }{\partial t} T_{{\parallel} a} + \boldsymbol{\nabla}_{{\parallel}} \left( \frac{Q_{{\parallel} a}}{v_{Ta}^2} + 3 u_{{\parallel} a}^\psi\right) - \boldsymbol{\nabla}_\parallel \ln B \left(\frac{Q_{{\parallel} a}}{v_{Ta}^2} + \frac{2 Q_{{\perp} a}}{v_{Ta}^2} + u_{{\parallel} a}^\psi - \frac{2 q_a}{m_a} \mathcal{K}_{1} \psi \right) \nonumber\\ & \quad + \frac{{\rm i} v_{Ta}^2 R_B}{2 \varOmega_a T_a} \left( 7 T_{{\parallel} a} - 4 N_a T_a + T_{{\perp} a } + q_a \phi ( 4 \mathcal{K}_{0} - \mathcal{K}_{1}) \right) \nonumber\\ & \quad + {\rm i} \left( \mathcal{K}_{0} ( \omega_N +\omega_{T_a} ) - \omega_{T_a} \mathcal{K}_{1} \right) \frac{ e\phi}{T_e}= \sqrt{2} \mathcal{C}_{a}^{20}, \end{align}

(5.1d)\begin{align}& \frac{1}{T_a} \frac{\partial }{\partial t} T_{{\perp} a} + \boldsymbol{\nabla}_\parallel \left( u_{{\parallel} a}^\psi - \frac{Q_{{\perp} a}}{v_{Ta}^2} + \frac{q_a}{T_a}\mathcal{K}_{1} \psi \right) - 2 \boldsymbol{\nabla}_\parallel \ln B \left( u_{{\parallel} a}^\psi - \frac{Q_{{\perp} a}}{v_{Ta}^2} + \frac{q_a}{T_a}\mathcal{K}_{1} \psi \right) \nonumber\\ & \quad + \frac{{\rm i} v_{Ta}^2 R_B}{2 \varOmega_a T_a}\left( T_{{\parallel} a} + 5 T_{{\perp} a} - 3 N_a T_a + q_a \phi \left( 2 \mathcal{K}_{2} + 3 \mathcal{K}_{0} - 5 \mathcal{K}_{1}\right) \right) \nonumber\\ & \quad + {\rm i} \left( \mathcal{K}_{0} ( \omega_N + \omega_{Ta}) - \mathcal{K}_{1} ( \omega_N + 3 \omega_{Ta}) + 2 \mathcal{K}_{2} \omega_{Ta} \right) \frac{e \phi}{T_e} ={-} \mathcal{C}_{a}^{01}, \end{align}

where we introduce $u_{\parallel a}^\psi = u_{\parallel a} - q_a \mathcal {K}_{0} \psi / m_a$. Similarly, we derive for the parallel and perpendicular heat fluxes, $Q_{\parallel a} = \sqrt {3 }v_{Ta}^3 N_a^{30}$ and $Q_{\perp a} = v_{Ta}^3 N_a^{11} / \sqrt {2}$

(5.2a)\begin{align} & \frac{1}{ \sqrt{3} v_{T_a}^3}\frac{\partial }{\partial t} Q_{{\parallel} a} + \frac{\sqrt{3} }{2} v_{Ta} \boldsymbol{\nabla}_\parallel \left( \frac{T_{{\parallel} a}}{T_a} - N_a \right) + \frac{{\rm i} v_{Ta}^2}{2 \varOmega_a} R_B \left( \frac{8 Q_{{\parallel} a}}{ \sqrt{3} v_{Ta}^3} + \frac{2 \sqrt{3} u_{{\parallel} a}^{\psi} }{v_{Ta}} \right) \nonumber\\ & \quad - {\rm i} v_{Ta} \frac{\sqrt{3}}{2 }\omega_{T_a} \mathcal{K}_{0} \psi = \mathcal{C}_{a}^{30}, \end{align}

(5.2b)\begin{align} & \frac{ \sqrt{2}}{ v_{Ta}^3 } \frac{\partial }{\partial t} Q_{ {\perp} a} + \frac{v_{Ta} }{\sqrt{2}} \boldsymbol{\nabla}_\parallel \left( N_a -\frac{T_{{\perp} a}}{T_a}\right) + \frac{v_{Ta}}{\sqrt{2}} \frac{q_a}{T_a} \boldsymbol{\nabla}_\parallel \left( \mathcal{K}_{1} \phi\right) \nonumber\\ & \quad + \frac{v_{Ta} }{\sqrt{2}} \left( \frac{T_{{\parallel} a}}{T_a} - \frac{T_{{\perp} a}}{T_a} + \frac{q_a}{T_a}\mathcal{K}_{1} \phi \right)\boldsymbol{\nabla}_\parallel \ln B \nonumber\\ & \quad + \frac{{\rm i} v_{Ta}^2}{2 \varOmega_a} R_B \left( \frac{6 \sqrt{2} Q_{{\perp} a}}{v_{Ta}^3} - 3 \sqrt{2} \frac{q_a }{T_a} \mathcal{K}_{1} \psi + \frac{T_{{\perp} a}}{T_a} - N_a - \frac{q_a }{ T_a} \mathcal{K}_{1} \phi \right) \nonumber\\ & \quad - {\rm i} v_{Ta}\left( \mathcal{K}_{1} \omega_N + \omega_{T_a} \left(3\mathcal{K}_{1} - \mathcal{K}_{0} - 2 \mathcal{K}_{2}\right) \right) \frac{e \psi}{ \sqrt{2} T_e} = \mathcal{C}_{a}^{11}, \end{align}

where the GMs, $N_a^{pj}$, with $p + 2j > 3$ are neglected. The evolution equations of the lowest-order gyrocentre fluid quantities, (5.1) and (5.2), are closed by the GK quasineutrality condition and GK Ampere's, given (2.28) and (2.29), where the higher-order GMs that appear in these equations are neglected. Equations (5.1) and (5.2) constitute a set of linearized fluid-like equations that evolve self-consistently the $6$ lowest-order GMs per particle species. We refer to this as the high-collisional $6$GM model. We remark that, while the $6$GM simply neglects all higher-order GMs by using a closure by truncation, previous gyrofluid models (see, e.g. Beer & Hammett 1996; Staebler, Kinsey & Waltz 2005) are based on ad hoc collisionless closure relations for these GMs derived by mimicking the linear collisionless response. However, these previous models either neglect or use simplified collision operator models, in contrast to the 6GM model presented here. In fact, closed analytical expressions for the $\mathcal {C}_{a}^{ps}$ terms appearing on the right-hand sides of (5.1) and (5.2) can be used in the case of the DK Coulomb collision operator reported in the appendix of Frei et al. (2022a) to model collisions accurately in the $6$GM model. While other collision operator models can also be considered, the use of the DK Coulomb operator guarantees a relatively simple (yet accurate) description of collisional effects, in contrast to the GK Coulomb operator which relies on the evaluation of a large number of sums depending on the value of $k_\perp$ (Frei et al. 2021).

5.2. Collisional effects on TEM and MTM microinstabilities

We first consider the collisional effects on a density gradient-driven TEM appearing with safety factor $q = 3$, magnetic shear $s = 0.8$, and inverse aspect ratio $\epsilon = 0.3$. A finite density gradient of $R_N = 4$, weaker than typical density gradients found in the middle of H-mode pedestals, is used and different values of electron temperature gradient are considered. The ITG drive is neglected for simplicity in this section by considering $R_{T_i} = 0$. We also use $T_i /T_e =1$ and introduce electromagnetic effects with $\beta _e = 10^{-4}$ (below the KBM linear threshold). Given these parameters, a density gradient-driven TEM is identified in the collisionless limit with a peak growth rate located near $k_y = 0.5$, propagating in the ion diamagnetic direction, i.e. $\omega _r > 0$. We study the effect of collisions on this density gradient-driven TEM at $k_y = 0.5$.

Since, typically, $\nu _{ei} R_0 / c_s \gtrsim 1$ at the top and bottom of H-mode pedestals, while $\nu _{ei} R_0 / c_s \ll 1$ in the core, we scan the electron collisionality, $\nu _e^*$, over several orders of magnitude and compute the TEM growth rate, $\gamma$, and the real mode frequency, $\omega _r$, using the DK and GK Coulomb, Sugama, and IS operators. To perform our numerical investigations, we use $(P,J) = (16,8)$, which is sufficient to guarantee convergence over the full collisionality range considered here.

The results of our analysis are shown in figure 18 in the cases of a pure density gradient driven TEM (i.e. $\eta _e = R_{T_e} / R_N =0$) and in the case of a TEM driven by equal density and electron temperature gradients (i.e. $\eta _e = 1$). We also plot the predictions of the high-collisional $6$GM model, derived in § 5.1, for comparison. First, we observe that the TEM growth rate decreases with $\nu _e^*$ in the banana regime ($\nu _e^* \lesssim 1$), while it increases with $\nu _e^*$ in the Pfirsch–Schlüter regime ($\nu _e^* \gtrsim 1$), in all cases. In addition, collisions tend to increase the TEM real mode frequency in all cases. It is noticeable that the purely density-driven TEM mode ($\eta _e = 0$) propagates in the ion diamagnetic direction ($\omega _r > 0$) and has a negative frequency when $\eta _e = 1$ (Ernst et al. 2009). Second, it is remarkable that the GK operators damp more strongly the TEM than the DK operators and that the presence of FLR collisional terms has a smaller effect on $\omega _r$ (Pan et al. 2020). In addition, we notice that the $6$GM (which ignores the FLR collisional term) overestimates the TEM growth rate and real mode frequency when $\nu _e^* \gtrsim 1$, but still captures the correct trend of the growth rate compared with the DK Coulomb. The agreement of the $6$GM model with the full GM hierarchy improves at a collisionality much larger than the ones considered in figure 18, i.e. when $\nu _e^* \gtrsim 50$, but not shown here.

Figure 18. The TEM growth rate (a,c) and real mode frequency (b,d) as a function of the electron collisionality, $\nu _e^*$, using the DK and GK Coulomb, Sugama and IS collision operators with $(P,J) = (16, 8)$, for $\eta _e = 0$ (a,b) and $\eta _e = 1$ (c,d). The results from the high-collisional $6$GM model are plotted for comparison (black cross markers). Here, $k_y = 0.5$.

As found in Pan et al. (2020, 2021) for the case of the GK Sugama relative to the GK Coulomb operator, it is noticeable that, despite the small differences observed between the Coulomb, Sugama, and IS operators in the case of purely density gradient-driven TEM ($\eta _e =0$), the presence of finite electron temperature gradient produces a non-negligible underestimate (up to $15\,\%$) of the TEM growth rate by the (DK and GK) Sugama and IS operators compared with respect to the (DK and GK) Coulomb operator.

We also notice that the IS operator approaches the predictions of the GK Coulomb when $\eta _e = 1$ and $\nu _e^* \gtrsim 1$ better than the Sugama one. The study of the TEM growth rate in Pan et al. (2020, 2021), as confirmed here, suggests that the accuracy of collision operator models (and the presence of FLR terms) compared with the Coulomb operator depends on the physical parameters considered, for example, these deviations increase with collisionality.

Following Pan et al. (2020, 2021), who analyse the dependence on the electron temperature gradient and collisionality, we first scan the TEM growth rate and frequency as a function of $\eta _e$ and $\nu _e^*$ using the GK Coulomb collision operator and repeat the calculations with the DK Coulomb, GK Sugama and GK IS operators. Then, the relative deviations of the TEM growth rate, $\sigma (\gamma ) = |\gamma - \gamma _C| / \gamma _C$ (with $\gamma _C$ the growth rate obtained using the GK Coulomb) is computed for all the different operators and the results are displayed in figure 19, reproducing figure 3(b) of Pan et al. (2021) with the GM approach and the addition of the GK IS operator. First, we observe that the effects of FLR collisional damping are clearly visible due to the deviations (up to $20\,\%$) appearing for $\nu _e^* \gtrsim 1$ when the DK Coulomb operator is used (Pan et al. 2020, 2021). Second, the deviations between the GK Sugama and GK IS from GK Coulomb are strongly dependent on the electron temperature gradient. Consistent with (Pan et al. 2021), for all collisionalities, $\sigma (\gamma )$ peaks near $\eta _e \sim 1.2$ and increases with collisionality reaching a maximum value of the order of $15\,\%$ for the GK Sugama and a value of $8\,\%$ for the GK IS. As explained in Pan et al. (2021), the influence of the electron temperature gradients on the accuracy of the Sugama and IS operators originates from the approximation in their field component, which are formulated as a truncated expansion of the $v^2$ moments of the distribution function and driven by finite $R_{Te}$ (see (2.25) with $p=0$ and $p=2$), explaining the qualitative dependence seen in figure 19. In addition, we remark that the GK IS performs better than the GK Sugama. This can also be explained by the fact that IS operator (Sugama et al. 2019) contains correction terms that are proportional to the difference between $v^2$ moments of the Sugama and Coulomb operators. The importance of these terms increases with $R_{Te}$.

Figure 19. Relative deviations of the TEM growth rate with respect to the case of the GK Coulomb, $\sigma (\gamma )$, when the DK Coulomb (a), GK Sugama (b) and GK IS (c) are used. The solid white line is the transition from ion to electron diamagnetic directions. Same parameters as in figure 18.

Finally, we investigate the collisional dependence of MTMs. Contrary to the MTM linear investigations in the core region that report the peak of the growth rate occurring near $\nu _{ei} /\omega _r \sim 1$ (with $\omega _r$ is the real MTM mode frequency) and vanish in the collisionless limit (Hazeltine & Strauss 1976; Catto & Rosenbluth 1981), MTMs found in the edge region display a different collisionality dependence. Indeed, edge GK simulations of MTMs (Doerk et al. 2012; Dickinson et al. 2013) suggest that the MTM growth rate does not vanish in the collisionless limit and remains nearly constant in the weak collisionality regime, $\nu _{ei} / \omega _r \ll 1$, while collisions have a stabilizing effect in the high-collisional limit, $\nu _{ei } / \omega _r \gg 1$. Hence, we scan the MTM growth rate and real mode frequency at $k_y = 0.5$ as a function of the electron collisionality, $\nu _{e}^*$, with the same parameters of the MTM described in § 4.4 and using the Coulomb, Sugama and IS operators. The results are shown in figure 20, where the high-collisional $6$GM model result is plotted as well for comparison. First, we remark that, in agreement with previous collisional MTM investigations, the growth rate is stabilized by collisions and flattens out for $\nu _{ei} / \omega _r \ll 1$. Interestingly, it is found that the choice of the GK operator does not significantly affect the MTM growth rate for $\nu _e^* \lesssim 1$, yielding a larger growth rate than the DK operators, while the latter have a stabilizing effect on the MTM followed by an increase of the real mode frequency $\omega _r$, not present in the GK operators. We also notice the good agreement between the $6$GM model and the DK Coulomb at high collisionality. Finally, in contrast to the TEM case (see figure 19) and the results of Pan et al. (2021) that consider a different MTM case based on JET pedestal parameters, the difference between the different collision operator models does not show a strong dependence on the electron temperature gradient in the case of the MTM considered here.

Figure 20. The MTM growth rate (a) and real mode frequency (b) as a function of the electron collisionality, $\nu _e^*$, using the DK and GK Coulomb, Sugama and IS collision operators with $(P,J) = (16, 8)$. Here, the parameters are the same as in figure 14 with $k_y = 0.5$.

5.3. Collisional effects on GAM dynamics

We now investigate the role of collisions in the GAM dynamics present in the edge region using the same assumptions as in § 4.5, i.e. adiabatic electrons). Hence, only the ion–ion collisions are considered in this section. Only a few theoretical works investigate the effect of collisions on the GAM dynamics (Lebedev et al. 1996; Novakovskii et al. 1997; Gao 2013), despite the fact that collisional effects can affect qualitatively and quantitatively the GAM damping and frequency when $\nu _{ii} \gtrsim 1$. Differences are observed between the collision operator models (see, e.g. Novakovskii et al. (1997) and Gao (2013), which consider a Hirschman–Sigmar–Clarke operator and a Krook operator, respectively), and it is usually found that collisionality decreases the GAM frequency, $\omega _G$, while it has a more complex effect on the GAM damping, $\gamma _G$. More precisely, the GAM damping is essentially proportional to the collisionality when $\nu _{ii} \lesssim 1$, i.e. $\gamma _G \sim \nu _{ii}$. On the other hand, the GAM damping is reduced, and collisional effects on the GAM frequency become important when $\nu _{ii} \gtrsim 1$.

To investigate the effect of collisions and collision operator models on the GAM dynamics, we consider the collisional dispersion relation derived by Gao (2013) in the limit of adiabatic electrons and long radial wavelengths, where ion–ion collisional effects are modelled with a particle conserving Krook operator

(5.3)\begin{equation} \mathcal{C}_i ={-} \nu_{ii} \left[ {\rm J}_0 h_i - \frac{F_{Mi}}{N} \int \,{\rm d} \boldsymbol{v} {\rm J}_0 h_i \right]. \end{equation}

We remark that the Krook operator in (5.3) conserves particles, but does not conserve momentum and energy. In our normalized units, the GAM dispersion relation derived by (Gao 2013) is

(5.4)\begin{gather} \frac{\xi - {\rm i} \hat{\nu}}{ \xi } \frac{1}{q^2} + \left[ \frac{1}{2} - \frac{1}{2 \xi^2} + \left( \xi^2 + 1 + \frac{1}{2 \xi^2} \right) \left( 1 + \xi Z(\xi)\right) \right] \nonumber\\ - \frac{1}{4 \xi^3} \left( \xi + {\rm i} \hat{\nu} \right) \frac{[1 - (2 \xi^2 +1)(1 + \xi Z(\xi))]^2}{2 + ({\rm i}\hat{\nu} + \xi)Z (\xi) } =0, \end{gather}

with $\xi = q (\omega _G + {\rm i} \gamma _G + {\rm i} \nu _{ii}) / \sqrt {2}$, $\hat {\nu } =q \nu _{ii} / \sqrt {2}$ and $Z(\xi ) = \int \,{{\rm d}\,x} {\rm e}^{-x^2} / (x - \xi ) / 2 {\rm \pi}$ the plasma dispersion function. We compare the analytical result in (5.4) with the GM approach simulations using the same operator in figure 21. To this aim, we project the Krook collision operator, (5.3), onto the Hermite–Laguerre basis in the DK limit, yielding

(5.5)\begin{equation} \mathcal{C}_i^{pj} ={-} \nu_{ii} \left( N_i^{pj} - \delta_{p}^0 \delta_j^0 N_i^{00} \right), \end{equation}

and compute $\gamma _G$ and $\omega _G$ as a function of $\nu _{ii}$ for different values of the safety factor $q$. To highlight the effect of collision operator models, the calculations are also performed using the DK Coulomb and DK Dougherty collision operators, which conserve momentum and energy. We first remark that convergence is achieved with $(P,J) = (24,8)$, a smaller number of GMs than in the collisionless case (see figure 16). Second, we notice the GAM damping and frequency, $\gamma _G$ and $\omega _G$, obtained from the numerical simulations using the Krook operator, (5.3), and the analytical prediction in (5.4) agree. Third, while all the collision operators present the same qualitative behaviour with collisionality in the predictions of $\gamma _G$ and $\omega _G$, significant quantitative differences can be observed. In fact, while $\gamma _G$ increases with $\nu _{ii}$ for $\nu _{ii} \lesssim 1$, such that $\gamma _G \sim \nu _{ii}$ for all operators, and eventually decreases for $\nu _{ii} \gtrsim 1$, the Krook operator overestimates the GAM damping and underestimates the GAM frequency. These deviations from the other collision operators are due to the lack of conservation properties of the Krook operator. Similar observations on the comparison between the Krook operator and other collision operator models (including an energy and momentum conserving Krook, a pitch-angle scattering, and the Hirschman–Sigmar–Clarke collision operators) are reported in Li & Gao (2015). We remark that the DK Dougherty collision operator yields a stronger GAM damping than the DK Coulomb operator. Not shown are the results from the Sugama and IS operators that yield results similar to the DK Coulomb, with a better agreement achieved by the IS operator at high collisionality.

Figure 21. GAM damping, $\gamma _G$ and frequency, $\omega _G$ as a function of the collisionality, $\nu _{ii}$, obtained from the dispersion relation (5.4) (black markers) and by using the Krook (red markers), the DK Coulomb (blue markers) and the DK Dougherty (green markers) collision operators. Different values of the safety factor are considered ($q=3$ with solid lines and $q =5$ with dashed lines), with $\epsilon = 0.1$.

5.4. Collisional ZF damping

The collisional damping of ZFs was first addressed in Hinton & Rosenbluth (1999) in the banana regime for radial wavelengths much larger than the ion gyroradius. Their work demonstrates that the long-time evolution of $\left \langle \phi \right \rangle _{fs}$ follows a slow exponential decay that converges to a finite value that scales as $B_p^2 / B^2$ (with $B_p$ the modulus of the poloidal magnetic field). More recently, by using a momentum conserving pitch-angle scattering operator for long radial wavelengths, Xiao, Catto & Molvig (2007) extends the analytical neoclassical prediction of Hinton & Rosenbluth (1999) to arbitrary finite collisionality and demonstrates that the long-time ZF residual follows

(5.6)\begin{equation} \frac{\left\langle\phi \right\rangle_{fs}(\infty)}{\left\langle\phi \right\rangle_{fs}(0)} \to \varsigma = \frac{\beta}{1 + \beta}, \end{equation}

where $\beta = \epsilon ^2 /q^2$. We compare the analytical prediction in (5.6) with the GM approach considering the Coulomb, the Sugama as well as the pitch-angle scattering operator, and the Dougherty collision operators, two operators not present in our previous ZF collisional damping tests (see, e.g. Frei et al. 2021). The presence of collisions allows us to evolve a smaller number of GMs than in the collisionless case to achieve convergences, i.e. $(P,J) =(24,12)$ (see figure 17).

Figure 22 shows the time evolution of $\left \langle \phi \right \rangle _{fs}$ for three increasing radial wavenumbers, $k_x = 0.05, 0.1$ and $0.2$, with a collisionality level in the Pfirsch–Schlüter regime, i.e. $\nu _i^* = 3.13$. The DK operators are used for $k_x = 0.05$, while the GK operators are considered for the larger values of $k_x$. Despite the small (but finite) values of radial wavenumbers, FOW effects are important at these parameters because the associated radial wavelengths are of the order of the poloidal ion gyroradius $\rho _p$, i.e. $k_x \rho _p \lesssim 1$ (see § 3). We first observe that the long-time ZF residual agrees with (5.6) for all operators when $k_x = 0.05$. Second, the effect of energy diffusion (absent in the pitch-angle scattering operator but present in the other operators) enhances the collisional ZF damping. Third, the presence of FLR terms in the collision operators yields a stronger ZF damping. This can be deduced by comparing the deviation between the GK Coulomb and the DK Coulomb operator in the $k_x = 0.1$ and $k_x = 0.2$ cases. We also notice the effects of FLR terms associated with the ion polarization term, which reduces the ZF residual, as it can be seen by comparing the analytical prediction of (5.6) and the DK Coulomb operator. Fourth, as previously observed in Frei et al. (2021), the GK Sugama collision operator provides a better approximation of the GK Coulomb than the other operators, while the GK Dougherty produces the strongest ZF damping. In Pan et al. (2020, 2021), the ZF damping is shown to be weaker with the GK Coulomb operator than with the Sugama operator by using GENE. Finally, we remark that the oscillations appearing at early times when the pitch-angle operator is used (absent in all other operators) demonstrate that energy diffusion is important in the collisional damping of high-order GMs. Indeed, these oscillations, which do not affect the long-time ZF residual, are absent in the operators that implement energy diffusion and also disappear with the pitch-angle operator when the number of GMs is increased.

Figure 22. Collisional ZF damping for increasing radial wavelengths $k_x = 0.05$ (a), $k_x = 0.1$ (b) and $k_x = 0.2$ (c) when $\nu _i^* = 3.13$. The DK collision operators are used when $k_x = 0.05$, while the GK collision operators are considered for $k_x = 0.1$ and $k_x =0.2$. The collisionless and collisional residuals, $\varpi$ (see figure 17) and $\varsigma$ respectively, are plotted with the black dashed and blue dashed lines. In the $k_x = 0.1$ and $k_x =0.2$ cases, the results using the DK Coulomb (blue dotted) are also shown for comparisons. Here, $q = 1.4$ and $\epsilon = 0.1$.

6. Microinstabilities in steep pressure gradient conditions

The presence of steep pressure gradients in the edge pedestals, when $R_0 / L_N \sim R_{T_{e,i}} \gtrsim 10$, leads to microinstabilities that can significantly differ from the ones usually encountered in the edge of L-mode discharges or in the core (Fulton et al. 2014; Xie & Xiao 2015; Xie & Li 2016; Han et al. 2017; Kotschenreuther et al. 2017; Xie, Xiao & Lin 2017b; Pueschel et al. 2019). In weak equilibrium gradient conditions, microinstabilities are often characterized by a conventional ballooning eigenmode function, with the electrostatic potential featuring an even mode parity around the outboard midplane position ($\chi = 0$) and peaking at the same location with a well-defined mode propagation direction. On the other hand, numerical studies (Fulton et al. 2014; Xie & Xiao 2015) reveal the existence of modes with unconventional parallel mode structures peaking at $\chi \neq 0$ when the gradients are increased to values relevant to the H-mode pedestals, i.e. $R_N \sim R_{T_{e,i}} \gtrsim 10$. In addition, transition in the mode parity can occur, often related to discontinuous jumps in the mode frequency and to changes in the mode propagation direction (e.g. from the ion to the electron diamagnetic direction or vice versa). The presence of these unconventional modes can potentially influence the level of particle and heat turbulent transport in the H-mode pedestal (Fulton et al. 2014; Xie et al. 2017b; Pueschel et al. 2019), and can possibility affect the commonly used mode identification criteria (Dickinson et al. 2012; Xie, Lu & Li 2018; Pueschel et al. 2019).

In the present study, we follow the nomenclature used in previous investigations (see, e.g. Xie et al. 2017b; Pueschel et al. 2019). We characterize the unstable modes by introducing a label, $\ell \ge 0$, associated with the structure of the ballooning eigenmode function and, in particular, the mode parity and number of peaks in the parallel direction. For instance, the $\ell = 0$ mode defines the conventional mode structure with even parity and peaking at the outboard midplane (with no secondary peak). On the other hand, the $\ell > 0$ modes are characterized by multiple peaks present at different parallel locations. Even values of $\ell$ denote even parity modes, and vice versa.

The transition from the $\ell =0$ modes to $\ell > 0$ can be identified by discontinuous jumps in the mode frequency $\omega _r$ and by the appearance of multiple peaks in the ballooning eigenmode function. We verify our results obtained using the GM approach with the direct GENE eigensolver, because of the presence of subdominant unstable modes with similar growth rates and related to the sensitivity of the initial value solver used in this work to the initial conditions (Xie et al. 2017a). For our investigation, we consider the parameters $q = 2.7$, $s =0.5$, and $\epsilon = 0.18$ in the low-collisionality banana regime with $\beta _e = 10^{-4}$. Since the $\ell > 0$ modes usually have large parallel wavenumbers (see below), we use $N_{k_x} = 10$, $N_z =32$ points and $(P,J) =(24,8)$ GMs. We consider the unstable modes occurring at a binormal wavenumber $k_y = 0.25$, which corresponds to the peak growth rate at the parameters used in this section.

To illustrate the appearance of the $\ell > 0$ modes, we plot the growth rate, $\gamma$, and real mode frequency, $\omega _r$, as a function of the normalized density gradient $R_N$ in figure 23, as obtained by using the GM approach and the GENE direct eigensolver in the case of $\eta _{e,i} = 1$ (i.e. $R_{Te}$ and $R_{T_i}$ equivalent to the density gradient $R_N$) and $\eta _{e,i} =0$ (i.e. absence of temperature gradients). A discontinuous jump in the real frequency $\omega _r$ is observed in all cases, and the ballooning eigenmode functions, obtained with the GM approach below and above the identified density gradient threshold $R_N \simeq 50$, are analysed in figure 24 in the case of $\eta _i = \eta _e = 1$. When $R_N \lesssim 50$, the most unstable mode displays a conventional, $\ell =0$, ballooning mode structure. On the other hand, the most unstable mode for $R_N \gtrsim 50$ is characterized by an unconventional mode structure that peaks at $\chi = {\rm \pi}/2$ and $\chi = 3{\rm \pi} /2$, justifying the $\ell = 2$ label for this mode. This is in good agreement with the eigenvalue spectrum obtained with GENE. We remark that the $\ell =0$ and $\ell =2$ modes are both characterized by a ballooning parity. However, a steeper gradient is required to drive the $\ell =2$ mode unstable, since it has a larger parallel wavenumber, $k_\parallel \sim \ell / q R_0$ (see figure 24) . Therefore, it is more sensitive to the stabilization effects of Landau damping than the $\ell = 0$ mode. Finally, we notice that the $\ell = 0$ mode persists when $\eta _i = \eta _e =0$, while it disappears when the electrons are adiabatic, we identify it as a TEM. Similarly, we identify the $\ell = 2$ mode as a TEM. Therefore, our results confirm that the mode identification based on the sign of the real mode frequency is ambiguous at steep gradients (Ernst et al. 2009). Indeed, the most unstable mode when $R_N \lesssim 50$ changes continuously from the ion ($\omega _r > 0$) to the electron ($\omega _r < 0$) diamagnetic directions (see figure 23).

Figure 23. Real mode frequency, $\omega _r$, and growth rate, $\gamma$, are shown by the blue and red markers, respectively as a function of the normalized density gradient, $R_N$, obtained by the GM approach (coloured markers) in the case of $\eta _e = \eta _i = 1$ (a) and $\eta _e = \eta _i = 0$ (b). The results from the GENE direct eigensolver are plotted by the black markers. The dominant $\ell = 0$ mode, characterized by $\omega _r > 0$ when $R_N \lesssim 50$, transits to the $\ell = 2$ mode with $\omega _r < 0$ when $R_N \gtrsim 60$ in all cases.

Figure 24. Real (blue lines) and imaginary (red lines) parts of the ballooning eigenmode functions of the electrostatic potential $\phi _B$ (a) and of the magnetic vector potential $\psi _B$ (b) corresponding to the $\ell = 0$ mode when $R_N = 20$ (dashed lines) and to the $\ell = 2$ mode when $R_N = 80$ (solid lines), identified in figure 23 for $\eta _e = \eta _i =1$. The ballooning eigenmode functions, $\phi _B$ and $\psi _B$, are normalized to $\phi _B(0)$.

We finally investigate the GM spectrum of the $\ell =0$ and $\ell =2$ modes. A convergence study reveals that the number of Hermite GMs, $P$, is reduced when increasing pressure gradients, such that convergence is achieved when $P \gtrsim 30$ for $R_N \sim 10$, while $P \gtrsim 10$ is sufficient above $R_N \sim 50$, with a small number of Laguerre GMs, i.e. $J \sim 3$ for all cases. This shows that, in general, the number of GMs decreases with $R_N$. This can be understood from the fact that the $\ell > 0$ modes found in the H-mode pedestals are expected to be less sensitive to magnetic gradient drift resonance effects than instabilities usually found in the core (Connor, Hastie & Helander 2006). Since magnetic gradient drifts and FOW effects, proportional to ${\rm i} \omega _{Ba}$ in (2.1), are responsible for broadening the collisionless GM spectrum (see § 3), we expect that a small number of GMs is required to describe the $\ell > 0$ modes appearing at steep pressure gradients since modes, for which the parallel dynamics is essential, have a collisionless GM spectrum considerably less extended than the modes driven by magnetic gradient effects (Frei et al. 2022b). As a confirmation, we plot in figure 25 the collisionless normalized electron and ion GM spectrum of the $\ell = 0$ and $\ell = 2$ TEM modes when $R_N = 20$ and $80$, respectively. We note the fast decay of the spectrum in the Hermite direction in the case of $R_N = 80$ compared with $R_N = 20$. In addition, in the former case, band structures can be identified, which are driven by the resonance effects associated with the ${\rm i} \omega _{Ba}$ term (Frei et al. 2022b). Finally, we observe that the electron GM spectrum is much broader than the ion GM, demonstrating the role of electron dynamics. The inspection of the collisionless GM spectrum suggests that the GM approach enables the description of H-mode pedestals with a relatively low velocity-space resolution even at low collisionality compared with core conditions (see § 4).

Figure 25. Electron (a,c) and ion (b,d) GM spectrum of the $\ell = 0$ TEM mode when $R_N=20$ (a,b) and of the $\ell = 2$ TEM when $R_N = 80$ (c,d). Here, $\eta _{e,i} = 1$.

7. Conclusion

This work presents the first linear flux-tube GK simulations carried out by using the GM approach at arbitrary collisionality. The approach is based on the projection of the perturbed gyrocentre distributions onto a Hermite–Laguerre basis. Building on previous studies using the same approach but performed in the local limit, kinetic effects of trapped and passing particles and electromagnetic effects are retained for the first time. A comprehensive linear study of microinstabilities, which includes the ITG, TEM, KBM, MTM as well as GAM dynamics and ZF damping, is performed with detailed comparisons with the continuum GK code GENE in the collisionless limit.

We successfully compare the linear growth rates and mode frequencies, velocity-space structures of the distribution functions, and eigenmode structures with GENE at low collisionality. The amplitude of the ZF residual is also verified against analytical predictions showing the ability of the GM approach to overcome the limitations of previous gyrofluid models. These investigations assess the convergence properties of the GM approach and identify the optimal number of GMs in the presence of strong kinetic effects that feature sharp velocity-space structures due to resonances associated with the drift of passing particles and the presence of trapped particles. We show that the GM approach agrees with GENE when the considered number of GMs, $(P,J)$, roughly equals the number of grid points, $(N_{v_\parallel },N_{\mu } )$, used to discretize the velocity space in GENE. Indeed, we find that $P \sim N_{v_\parallel }$ and $J \sim N_{\mu }$ are necessary to achieve convergence in most cases when parameters relevant to the core region are used, such as low collisionality and weak pressure gradients. On the other hand, we demonstrate that the necessary number of GMs decreases with collisionality and a reduced number of GMs is sufficient, even in the low-collisionality regime, to achieve convergence in the case of modes such as KBM and modes destabilized in steep pressure gradients regions found, e.g. in H-mode pedestals. This allows us to speculate that the GM approach features convergence properties well adapted to perform future nonlinear simulations of the plasma boundary.

Taking advantage of the formulation of advanced collision operators, including the Coulomb, Sugama and, more recently, the IS collision operators within the GM approach, we investigate the TEM and MTMs (two important edge microinstabilities) exploring a collisionality from the low-collisionality banana to the high-collisionality Pfirsch–Schlüter regimes. We demonstrate that the FLR terms in the collision operators are essential since they reduce the level of collisionality where a significant stabilization of the TEM and a suppression of the MTM is observed. In addition, comparing the predictions of the different collision operator models with the GK Coulomb allows for the assessment of the accuracy of other collision operator models. In all cases, non-negligible deviations with the GK Coulomb are observed at collisionalities relevant to H-mode pedestals. While these deviations increase with collisionality in all cases, the most significant ones are found at finite electron temperature gradients, in particular, in the case of the TEM. Indeed, the GK Sugama operator underestimates the linear growth rate up to $15\,\%$ and the GK IS operator up to $8\,\%$. Finally, the impact of collisions on the GAM dynamics and ZF collisional damping show that the analytical details of collision operator models (e.g. conservation laws and energy diffusion) are essential to correctly predict their long-time evolution. In general, the present results demonstrate that a careful analysis of the collisional dependence of microinstabilities and, more generally, of the impact of the choice of collision operator model is necessary to carry out accurate collisional simulations of the plasma dynamics in the boundary region.

While the analysis presented in this work is limited to linear cases, the extension of the GM method to the nonlinear turbulent regime using advanced collision operators is underway (Hoffmann, Frei & Ricci 2023). We also remark that significant progress has been recently made in the development of collisionless nonlinear flux-tube simulations using a similar approach (Mandell et al. 2022). We also note that, although the numerical implementation of the GM hierarchy presented here is restricted to the flux-tube configuration and relies on the linearized GK $\delta f$ approach, the present study paves the way to future nonlinear simulations of the boundary region based on the GM approach, including a realistic geometry and full-F conditions. Ultimately, we expect that the GM method will enable comprehensive simulations with a reduced computational cost than high-fidelity GK simulations when applied to, e.g. the Pfirsch–Schlüter regime and low-collisionality H-mode pedestal conditions. At the same time, the GM approach provides an improved fluid description over the reduced Braginskii-like fluid model in the low-collisionality limit. This work represents a first step towards future full-F simulations of nonlinear turbulence using the GM approach.