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. (Reference Frei, Jorge and Ricci2020), 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 Reference Hazeltine and Meiss2003)

(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. Reference Frei, Ball, Hoffmann, Jorge, Ricci and Stenger2021). 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. Reference Frei, Jorge and Ricci2020), 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. (Reference Frei, Jorge and Ricci2020) 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 Reference Dorland and Hammett1993). 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. Reference Frei, Hoffmann and Ricci2022b). 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. Reference Lapillonne, Brunner, Dannert, Jolliet, Marinoni, Villard, Görler, Jenko and Merz2009)

(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 Reference Ball and Brunner2021). 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. Reference Beer, Cowley and Hammett1995). 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. Reference Beer, Cowley and Hammett1995; Lapillonne et al. Reference Lapillonne, Brunner, Dannert, Jolliet, Marinoni, Villard, Görler, Jenko and Merz2009)

(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 Reference Connor, Hastie and Taylor1978)

(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. Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000). Despite its known inconsistencies (Lapillonne et al. Reference Lapillonne, Brunner, Dannert, Jolliet, Marinoni, Villard, Görler, Jenko and Merz2009), 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. Reference Jenko, Dorland, Kotschenreuther and Rogers2000) 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 Reference Jorge, Ricci and Loureiro2017; Mandell, Dorland & Landreman Reference Mandell, Dorland and Landreman2018; Jorge et al. Reference Jorge, Frei and Ricci2019; Frei et al. Reference Frei, Jorge and Ricci2020), 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 Reference Gradshteyn and Ryzhik2014)

(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 Reference Belli and Candy2011). However, the use of the Hermite–Laguerre basis, which has a long history in plasma physics (see, e.g. Grant & Feix Reference Grant and Feix1967; Hirshman & Sigmar Reference Hirshman and Sigmar1976; Madsen Reference Madsen2013; Schekochihin et al. Reference Schekochihin, Parker, Highcock and Dellar2016; Jorge et al. Reference Jorge, Ricci and Loureiro2017; Mandell et al. Reference Mandell, Dorland and Landreman2018), provides a direct relation to the fluid quantities that are evolved by Braginskii-like fluid models (Zeiler et al. Reference Zeiler, Drake and Rogers1997). 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 Reference Gradshteyn and Ryzhik2014)

(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. Reference Mandell, Dorland and Landreman2018, Reference Mandell, Dorland, Abel, Gaur, Kim, Martin and Qian2022), 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. (Reference Frei, Ball, Hoffmann, Jorge, Ricci and Stenger2021, Reference Frei, Hoffmann and Ricci2022b, Reference Frei, Ernst and Riccia). In contrast to the GX code (Mandell et al. Reference Mandell, Dorland, Abel, Gaur, Kim, Martin and Qian2022) that implements a Dougherty collision operator, we consider here the linearized GK Coulomb (Li & Ernt Reference Li and Ernst2011; Pan & Ernst Reference Pan and Ernst2019; Frei et al. Reference Frei, Ball, Hoffmann, Jorge, Ricci and Stenger2021), the Sugama (Sugama et al. Reference Sugama, Watanabe and Nunami2009), the IS (Sugama et al. Reference Sugama, Matsuoka, Satake, Nunami and Watanabe2019) collision operators, and like-species Dougherty (Dougherty Reference Dougherty1964) 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 Reference Helander and Sigmar2002)

(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. Reference Jorge, Ricci and Loureiro2017) or a semi-collisional closure Loureiro, Schekochihin & Zocco Reference Loureiro, Schekochihin and Zocco2013), 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 Reference Lee1987; Lin et al. Reference Lin, Nishimura, Xiao, Holod, Zhang and Chen2007).

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. Reference Mandell, Dorland, Abel, Gaur, Kim, Martin and Qian2022) has been developed for this purpose and successfully achieved this goal.

A comparison between the continuum GK GENE code (Jenko et al. Reference Jenko, Dorland, Kotschenreuther and Rogers2000; Görler et al. Reference Görler, Lapillonne, Brunner, Dannert, Jenko, Aghdam, Marcus, McMillan, Merz and Sauter2011) 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.1. Parallel streaming and recurrence phenomena

Passing particles are known to generate fine filament-like structures in $v_\parallel$ (Idomura et al. Reference Idomura, Ida, Kano, Aiba and Tokuda2008), 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. Reference Hammett, Beer, Dorland, Cowley and Smith1993). 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 Reference Rosenbluth and Hinton1998). 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 Reference Gradshteyn and Ryzhik2014), 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. (Reference Frei, Hoffmann and Ricci2022b) and be evaluated using numerical algorithms (Gürcan Reference Gürcan2014)), 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. Reference Frei, Hoffmann and Ricci2022b).

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. Reference Sugama and Watanabe2006), 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. Reference Frei, Hoffmann and Ricci2022b).

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.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. Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000). The cyclone base case is widely used to verify GK codes (Merlo et al. Reference Merlo, Sauter, Brunner, Burckel, Camenen, Casson, Dorland, Fable, Görler and Jenko2016; Tronko et al. Reference Tronko, Bottino, Görler, Sonnendrücker, Told and Villard2017). 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. Reference Frei, Hoffmann and Ricci2022b) 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. Reference Mandell, Dorland, Abel, Gaur, Kim, Martin and Qian2022), 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. Reference Frei, Hoffmann and Ricci2022b).

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. Reference Mandell, Dorland, Abel, Gaur, Kim, Martin and Qian2022).

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 Reference Lee1987; Lin et al. Reference Lin, Nishimura, Xiao, Holod, Zhang and Chen2007), that can limit the explicit time stepping scheme. While using a reduced ion mass can limit the non-adiabatic electron response (Dominski et al. Reference Dominski, Brunner, Görler, Jenko, Told and Villard2015), 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 Reference Hallatschek and Dorland2005), 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. Reference Dominski, Brunner, Görler, Jenko, Told and Villard2015; Ajay, Brunner & Ball Reference Ajay, Brunner and Ball2021). 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. Reference Merlo, Sauter, Brunner, Burckel, Camenen, Casson, Dorland, Fable, Görler and Jenko2016), 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. Reference Frei, Hoffmann and Ricci2022b). 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 Reference Weiland and Hirose1992; Citrin et al. Reference Citrin, Garcia, Görler, Jenko, Mantica, Told, Bourdelle, Hatch, Hogeweij and Johnson2014), 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. Reference Connor, Hastie and Taylor1978; Tang, Connor & Hastie Reference Tang, Connor and Hastie1980; Aleynikova & Zocco Reference Aleynikova and Zocco2017).

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 Reference Belli and Candy2010). 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. Reference Dickinson, Roach, Saarelma, Scannell, Kirk and Wilson2012). More precisely, MTMs are usually driven unstable by a combination of finite electron temperature and collisionality (even small) in the core region (Catto & Rosenbluth Reference Catto and Rosenbluth1981). MTMs also exist in the edge region in the collisionless limit, driven unstable by the electron magnetic drift resonance effects (Applegate et al. Reference Applegate, Roach, Connor, Cowley, Dorland, Hastie and Joiner2007; Dickinson et al. Reference Dickinson, Roach, Saarelma, Scannell, Kirk and Wilson2013).

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. Reference Dickinson, Roach, Saarelma, Scannell, Kirk and Wilson2013), 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. Reference Doerk, Jenko, Görler, Told, Pueschel and Hatch2012; Dickinson et al. Reference Dickinson, Roach, Saarelma, Scannell, Kirk and Wilson2013) (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. Reference Frei, Hoffmann and Ricci2022b). 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. Reference Winsor, Johnson and Dawson1968). 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 (Reference Hinton and Rosenbluth1999), as well as with the damping rate and frequency, $\gamma _G$ and $\omega _G$, given in Sugama et al. (Reference Sugama and Watanabe2006). 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. Reference Frei, Ball, Hoffmann, Jorge, Ricci and Stenger2021, Reference Frei, Hoffmann and Ricci2022b).

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 Reference Hinton and Rosenbluth1999), 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. (Reference Sugama and Watanabe2006) 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. Reference Diamond, Itoh, Itoh and Hahm2005). Rosenbluth & Hinton (Reference Rosenbluth and Hinton1998) 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 (Reference Xiao and Catto2006) 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. Reference Merlo, Sauter, Brunner, Burckel, Camenen, Casson, Dorland, Fable, Görler and Jenko2016), in contrast to gyrofluid models (see, e.g. Beer & Hammett Reference Beer and Hammett1996). 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 Reference Rosenbluth and Hinton1998). 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.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. Reference Jorge, Ricci and Loureiro2017; Frei et al. Reference Frei, Hoffmann and Ricci2022b), where $\epsilon _\nu \ll 1$ being the ratio of the electron mean-free path to the typical parallel scale length (Chapman & Cowling Reference Chapman and Cowling1941). 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. Reference Frei, Jorge and Ricci2020). 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 Reference Beer and Hammett1996; Staebler, Kinsey & Waltz Reference Staebler, Kinsey and Waltz2005) 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. (Reference Frei, Ernst and Ricci2022a) 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. Reference Frei, Ball, Hoffmann, Jorge, Ricci and Stenger2021).

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. Reference Ernst, Lang, Nevins, Hoffman, Chen, Dorland and Parker2009). 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. Reference Pan, Ernst and Crandall2020). 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. (Reference Pan, Ernst and Crandall2020, Reference Pan, Ernst and Hatch2021) 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. (Reference Pan, Ernst and Crandall2020, Reference Pan, Ernst and Hatch2021), 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. (Reference Pan, Ernst and Crandall2020, Reference Pan, Ernst and Hatch2021), 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. (Reference Pan, Ernst and Hatch2021) 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. Reference Pan, Ernst and Crandall2020, Reference Pan, Ernst and Hatch2021). 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. Reference Pan, Ernst and Hatch2021), 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. (Reference Pan, Ernst and Hatch2021), 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. Reference Sugama, Matsuoka, Satake, Nunami and Watanabe2019) 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 Reference Hazeltine and Strauss1976; Catto & Rosenbluth Reference Catto and Rosenbluth1981), MTMs found in the edge region display a different collisionality dependence. Indeed, edge GK simulations of MTMs (Doerk et al. Reference Doerk, Jenko, Görler, Told, Pueschel and Hatch2012; Dickinson et al. Reference Dickinson, Roach, Saarelma, Scannell, Kirk and Wilson2013) 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. (Reference Pan, Ernst and Hatch2021) 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. Reference Lebedev, Yushmanov, Diamond, Novakovskii and Smolyakov1996; Novakovskii et al. Reference Novakovskii, Liu, Sagdeev and Rosenbluth1997; Gao Reference Gao2013), 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. (Reference Novakovskii, Liu, Sagdeev and Rosenbluth1997) and Gao (Reference Gao2013), 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 (Reference Gao2013) 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 Reference Gao2013) 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 (Reference Li and Gao2015). 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 (Reference Hinton and Rosenbluth1999) 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 (Reference Xiao, Catto and Molvig2007) extends the analytical neoclassical prediction of Hinton & Rosenbluth (Reference Hinton and Rosenbluth1999) 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. Reference Frei, Ball, Hoffmann, Jorge, Ricci and Stenger2021). 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. (Reference Frei, Ball, Hoffmann, Jorge, Ricci and Stenger2021), 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. (Reference Pan, Ernst and Crandall2020, Reference Pan, Ernst and Hatch2021), 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$.