Hostname: page-component-848d4c4894-pjpqr Total loading time: 0 Render date: 2024-07-07T03:37:33.582Z Has data issue: false hasContentIssue false

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

Published online by Cambridge University Press:  10 August 2023

B.J. Frei*
Affiliation:
Ecole Polytechnique Fédérale de Lausanne (EPFL), Swiss Plasma Center, CH-1015 Lausanne, Switzerland
A.C.D. Hoffmann
Affiliation:
Ecole Polytechnique Fédérale de Lausanne (EPFL), Swiss Plasma Center, CH-1015 Lausanne, Switzerland
P. Ricci
Affiliation:
Ecole Polytechnique Fédérale de Lausanne (EPFL), Swiss Plasma Center, CH-1015 Lausanne, Switzerland
S. Brunner
Affiliation:
Ecole Polytechnique Fédérale de Lausanne (EPFL), Swiss Plasma Center, CH-1015 Lausanne, Switzerland
Z. Tecchioll
Affiliation:
Ecole Polytechnique Fédérale de Lausanne (EPFL), Swiss Plasma Center, CH-1015 Lausanne, Switzerland
*
Email address for correspondence: baptiste.frei@epfl.ch
Rights & Permissions [Opens in a new window]

Abstract

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

Type
Research Article
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright
Copyright © The Author(s), 2023. Published by Cambridge University Press

1. Introduction

Linear and nonlinear gyrokinetic (GK) simulations are the tools of reference for describing low-frequency (compared with the ion gyrofrequency, $\varOmega _i$) and small scales (of the order of the ion gyroradius, $\rho _i$) electromagnetic microinstabilities occurring in the core and edge regions of fusion devices (Told et al. Reference Told, Jenko, Xanthopoulos, Horton and Wolfrum2008; Holland et al. Reference Holland, Schmitz, Rhodes, Peebles, Hillesheim, Wang, Zeng, Doyle, Smith and Prater2011; Navarro et al. Reference Navarro, Happel, Görler, Jenko, Abiteboul, Bustos, Doerk and Told2015; Kotschenreuther et al. Reference Kotschenreuther, Hatch, Mahajan, Valanju, Zheng and Liu2017; Neiser et al. Reference Neiser, Jenko, Carter, Schmitz, Told, Merlo, Bañón Navarro, Crandall, McKee and Yan2019). On the other hand, the use of GK in the turbulent simulation of the boundary region, which includes both the edge and the scrape-off-layer (SOL), remains challenging, despite the recent development of edge particle and continuum GK codes (Churchill et al. Reference Churchill, Chang, Ku and Dominski2017; Mandell et al. Reference Mandell, Hakim, Hammett and Francisquez2020; Michels et al. Reference Michels, Stegmeir, Ulbl, Jarema and Jenko2021). Gyrokinetic simulations of the boundary are currently restricted by (i) their considerable computational cost, (ii) the presence of large-scale fluctuations, which are not present in the core, and (iii) the challenge of describing the high-collisionality regime using proper collision operator models, such as the Fokker–Planck Landau collision operator (Landau Reference Landau1936), referred to as the Coulomb operator in this work. Turbulence in the SOL region is most often simulated by models based on drift-reduced Braginskii-like fluid equations, which evolve the lowest-order particle fluid moments (density, temperature and velocity) (Zeiler, Drake & Rogers Reference Zeiler, Drake and Rogers1997). Braginskii-like fluid simulations of the SOL turbulence have shown their ability to model the SOL in a complex magnetic field topology (see, e.g. Stegmeir et al. Reference Stegmeir, Ross, Body, Francisquez, Zholobenko, Coster, Maj, Manz, Jenko and Rogers2019; Giacomin, Stenger & Ricci Reference Giacomin, Stenger and Ricci2020; Bufferand et al. Reference Bufferand, Bucalossi, Ciraolo, Falchetto, Gallo, Ghendrih, Rivals, Tamain, Yang and Giorgiani2021), in good agreement with experimental results (see, e.g. De Oliviera et al. Reference De Oliviera, Body, Galassi, Theiler, Laribi, Tamain, Stegmeir, Giacomin, Zholobenko and Ricci2022; Galassi et al. Reference Galassi, Theiler, Body, Manke, Micheletti, Omotani, Wiesenberger, Baquero-Ruiz, Furno and Giacomin2022). The validity of Braginskii-like models relies on the high-collisionality assumption, quantified by the smallness of the ratio of the particle mean-free path to the parallel scale length, $\lambda _{mfp} / L_\parallel \ll 1$. This scaling might not be appropriate to describe the entire collisionality range of the SOL and, more generally, in the boundary region. In particular, the high plasma temperature at the top of the pedestal and local transient events (such as edge localized modes) can significantly lower the plasma collisionality, even in the SOL, calling for a kinetic description of the boundary region. Aiming to bridge the gap between fluid and GK simulations, a moment approach to the GK model based on a Hermite–Laguerre decomposition of the full gyrocentre distribution function (full-F) was recently introduced in Frei, Jorge & Ricci (Reference Frei, Jorge and Ricci2020). This model, which we refer to as the gyro-moment (GM) approach, is derived in a generalized GK ordering appropriate to the boundary region and is valid for an arbitrary level of collisionality since it implements the full GK Coulomb collision operator (Jorge, Frei & Ricci Reference Jorge, Frei and Ricci2019). The ability of the GM approach to describe drift waves (Jorge, Ricci & Loureiro Reference Jorge, Ricci and Loureiro2018) and ion-scale instabilities (Frei, Hoffmann & Ricci Reference Frei, Hoffmann and Ricci2022b) efficiently has been demonstrated at an arbitrary level of collisionality using the GK Coulomb collision operator and other advanced collision operator models (Frei et al. Reference Frei, Ball, Hoffmann, Jorge, Ricci and Stenger2021; Frei, Ernst & Ricci Reference Frei, Ernst and Ricci2022a). However, these investigations are limited to electrostatic and local linear studies neglecting, for instance, electromagnetic and trapped particle effects, excluding therefore instabilities such as the trapped electron modes (TEM), recognized as one of the main drives of electron heat transport in the boundary region (Rafiq et al. Reference Rafiq, Pankin, Bateman, Kritz and Halpern2009; Schmitz et al. Reference Schmitz, Holland, Rhodes, Wang, Zeng, White, Hillesheim, Peebles, Smith and Prater2012), as well as the kinetic ballooning modes (KBM), which can limit, for instance, the maximal achievable pressure gradient in H-mode pedestals (Snyder et al. Reference Snyder, Groebner, Leonard, Osborne and Wilson2009; Wan et al. Reference Wan, Parker, Chen, Yan, Groebner and Snyder2012).

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

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

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

2. Flux-tube GM model

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

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

2.1. The GK model

We consider the linearized electromagnetic GK Boltzmann equation in the presence of an equilibrium magnetic field, as well as density and temperature gradients. The flux-tube assumption of separation between the turbulent (of the order of $\rho _i$) and the equilibrium (of the order of $L_\perp$) scales allows us to neglect the equilibrium profiles as considered constant across the computational domain. In the following, we use the gyrocentre phase-space coordinates $\boldsymbol {Z} = (\boldsymbol {R}, \mu, v_\parallel,\theta )$, where $\boldsymbol {R} = \boldsymbol {r} - \boldsymbol{\rho} _a$ is the gyrocentre position, with $\boldsymbol{r}$ the particle position and $\boldsymbol{\rho} _a(\boldsymbol{R}, \mu, \theta ) = \boldsymbol{b} \times \boldsymbol{v}/\varOmega _a$ its gyroradius. Here, $\boldsymbol {b} = \boldsymbol{B} / B$, $\varOmega _a = q_a B / m_a$, $a$ is the particle species, $\mu = m_a v_\perp ^2/[2 B(\boldsymbol {R})]$ is the magnetic moment, $v_\parallel = \boldsymbol{b} \boldsymbol {\cdot } \boldsymbol{v}$ is the component of the velocity parallel to the equilibrium magnetic field and, finally, $\theta$ is the gyroangle. Contrary to Frei et al. (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.7ac)). 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. Representation of passing particle drifts in the GM approach

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

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

3.1. Parallel streaming and recurrence phenomena

Passing particles are known to generate fine filament-like structures in $v_\parallel$ (Idomura et al. 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. Collisionless microinstability and comparison with GENE

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

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

4.1. Cyclone base case with adiabatic electrons

As a first linear collisionless test, we consider the electrostatic ITG cyclone base case scenario with adiabatic electrons (Dimits et al. 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 ac). 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}$ (ac) and of the TEM at $k_y = 1.3$ (df), 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 (cf) 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. High-collisional limit and collisional effects on microinstabilities

While collisional effects are often neglected in the core, they can no longer be ignored near the separatrix and in the SOL region because of the low plasma temperatures. For instance, $\nu _e^* \sim 0.03$ is expected at the top of the pedestal in ITER and $\nu _e^* \gtrsim 50$ near the separatrix. The high collisionality significantly affects the linear properties of edge microinstabilities, in particular the TEMs and MTMs which have been identified to play a major role in the pedestal region (Fulton et al. Reference Fulton, Lin, Holod and Xiao2014; Hatch et al. Reference Hatch, Kotschenreuther, Mahajan, Valanju, Jenko, Told, Görler and Saarelma2016; Garcia et al. Reference Garcia, De La Luna, Sertoli, Casson, Mazzi, Štancar, Szepesi, Frigione, Garzotti and Rimini2022).

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

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

5.1. Linear high-collisional limit

To consider the high-collisional limit, we assume $N_a^{00} \sim N_a^{10} \sim N_a^{01} \sim N_a^{20}$ and $N_a^{30} \sim N_a^{11} \sim \epsilon _\nu N_a^{00}$ (Jorge et al. 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$.

6. Microinstabilities in steep pressure gradient conditions

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

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

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

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

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

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

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

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

7. Conclusion

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

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

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

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

Acknowledgements

The authors acknowledge helpful discussions with J. Ball, P. Donnel and R. Jorge.

Editor William Dorland thanks the referees for their advice in evaluating this article.

Funding

This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 – EUROfusion). Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them. The simulations presented herein were carried out in part on the CINECA Marconi supercomputer under the TSVVT421 project and in part at CSCS (Swiss National Supercomputing Center). This work was supported in part by the Swiss National Science Foundation.

Declaration of interests

The authors report no conflict of interest.

Appendix A. Collisionless, local and strong ballooning limit of the flux-tube model

In this appendix, we perform a collisionless, local and strong ballooning limit analysis of the GM approach. To this aim, we derive an electromagnetic GK dispersion relation by solving explicitly the GK model introduced in § 2.1. We treat the electron kinetically and make no ordering assumption neither on the amplitude of perpendicular wavenumber nor on the magnitude of the magnetic drift frequency ${\rm i} \omega _{Ba}$. The dispersion relation we obtain allows us to perform a local convergence analysis as a function of the number of GMs $(P,J)$ in the presence of non-adiabatic electrons and electromagnetic effects. We note that the local analysis performed in this section neglects the contributions from the trapped particles and, therefore, ignores modes driven unstable by trapped particle effects, such as TEM. Nevertheless, we remark that the contribution from the trapped particles can be included in the analysis by solving their bounced averaged kinetic equation. We derive the electromagnetic GK dispersion relation in § A.1 and study the convergence properties of the GM approach in the case of ITG and KBM in § A.2.

A.1. Local electromagnetic GK dispersion relation

We evaluate (2.1) at the outboard midplane location (i.e. $z = 0$ and $k_x =0$). As a consequence, the parallel gradient of the magnetic field strength vanishes ($\boldsymbol{b} \boldsymbol {\cdot } \boldsymbol {\nabla } B = 0$), and the contribution from the trapped particles is ignored. The local approximation allows us to introduce the parallel wavenumber $k_\parallel \simeq 1 / q \partial _z$ and the perpendicular wavenumber $k_\perp$, defined in (2.18), reduces to $k_\perp = k_y$. Therefore, the parallel and perpendicular wavenumbers, $k_\parallel$ and $k_\perp$, are treated as scalar values and input parameters in the local limit.

Neglecting collisions appearing on the right-hand side of (2.1) and Fourier transforming in time, an explicit expression for the perturbed gyrocentre distribution function ${g}_a$ can be obtained, i.e.

(A1)\begin{equation} {g}_a = \sum_{j=1}^3 \left( {g}_{a \phi}^{(j)} \phi + {g}_{a \psi}^{(j)}\psi \right), \end{equation}

where the electrostatic, ${g}_{a \phi }^{(j)}$, and electromagnetic ${g}_{a \psi }^{(j)}$, components of ${g}_a$ are defined by

(A2a)\begin{gather} {g}_{a \phi}^{(1)} ={-} \frac{q_a}{\tau_a} F_{Ma} {\rm J}_0(b_a \sqrt{x_a}), \end{gather}
(A2b)\begin{gather}{g}_{a \phi}^{(2)} = \frac{q_a}{\tau_a} \frac{\omega {\rm J}_0(b_a \sqrt{x_a}) F_{aM} } { \omega - \omega_{Ba}- z_{{\parallel} a} s_{{\parallel} a}/ \sigma_a}, \end{gather}
(A2c)\begin{gather}{g}_{a \phi}^{(3)} ={-} \frac{\omega_{Ta}^* {\rm J}_0(b_a \sqrt{x_a}) F_{Ma} } { \omega - \omega_{Ba}- z_{{\parallel} a} s_{{\parallel} a}/ \sigma_a}, \end{gather}

and

(A3a)\begin{gather} {g}_{a \psi}^{(1)} = \frac{\sqrt{2 }}{\sigma_a} \frac{q_a}{ \sqrt{\tau_a}} F_{Ma} s_{{\parallel} a} {\rm J}_0(b_a \sqrt{x_a}), \end{gather}
(A3b)\begin{gather}{g}_{a \psi}^{(2)} ={-} \frac{\sqrt{2 }}{\sigma_a}\frac{q_a}{\sqrt{\tau_a}} \frac{\omega s_{{\parallel} a} {\rm J}_0(b_a \sqrt{x_a}) F_{aM} } { \omega - \omega_{Ba}- z_{{\parallel} a} s_{{\parallel} a}/ \sigma_a}, \end{gather}
(A3c)\begin{gather}{g}_{a \psi}^{(3)} = \frac{\sqrt{2 \tau_a}}{\sigma_a} \frac{\omega_{Ta}^* s_{{\parallel} a} {\rm J}_0(b_a \sqrt{x_a}) F_{Ma} } { \omega - \omega_{Ba}- z_{{\parallel} a} s_{{\parallel} a}/ \sigma_a}, \end{gather}

respectively. Here, the local magnetic drift frequency is $\omega _{Ba} = \alpha _a ( x_a + 2 s_{\parallel a}^2 )$ (with $\alpha _a = \tau _a k_\perp /q_a$) and $z_{\parallel a} = \sqrt {2 \tau _a} k_\parallel / \sigma _a$.

The electromagnetic GK dispersion relation is obtained by inserting (A1) into the GK quasineutrality condition and making use of the GK Ampere's law, given by (2.3) and (2.4), respectively. This yields the GK dispersion relation

(A4)\begin{align} D(\omega; k_\perp, k_\parallel, R_N, R_{Ta}, \beta_e) & = \left(\sum_{a} \frac{q_a^2}{\tau_a}\left( 1 - \varGamma_0(a_a)\right) - \sum_a q_a \sum_{j=1}^3 \delta n_{a\phi}^{(j)}\right) \nonumber\\ & \quad \times \left( \frac{k^2}{\beta_e} + \sum_a\frac{q_a^2}{\sigma_a^2} \varGamma_0(a_a) - \sum_a q_a\sum_{j=1}^3 \delta u_{a\psi}^{(j)}\right)\nonumber\\ & \quad - \left(\sum_a q_a\sum_{j=0}^3 \delta n_{a\psi}^{(j)} \right) \left(\sum_{a'} q_{a'} \sum_{j'=1}^3 \delta u_{a'\phi}^{(j')} \right) =0, \end{align}

where the zeroth and first-order velocity moments of ${g}_a$ are defined by $\delta n_{a \phi }^{(j)} = \int \,{\rm d} \boldsymbol{v} {\rm J}_0(b_a\sqrt {x_a}) {g}_{a \phi }^{(j)}$, $\delta n_{a \psi }^{(j)} = \int \,{\rm d} \boldsymbol{v} {\rm J}_0(b_a\sqrt {x_a}) {g}_{a \psi }^{(j)}$, $\delta u_{a \phi }^{(j)} = \int \,{\rm d} \boldsymbol{v} {\rm J}_0(b_a\sqrt {x_a}) s_{\parallel a}{g}_{a \phi }^{(j)}$ and $\delta u_{a \psi }^{(j)} = \int \,{\rm d} \boldsymbol{v} {\rm J}_0(b_a\sqrt {x_a}) s_{\parallel a}{g}_{a \psi }^{(j)}$. In order to solve $D(\omega ) = 0$ for the mode complex frequency $\omega$, we consider the following transformation of the velocity resonant term for the unstable modes when $\text {Im}(\omega ) > 0$ (Frei et al. Reference Frei, Hoffmann and Ricci2022b)

(A5)\begin{equation} \frac{1}{ \omega - \omega_{Ba} - z_{{\parallel} a} s_{{\parallel} a}/ \sigma_a} ={-} {\rm i} \int_0^\infty \,{\rm d} \tau \exp\left({ {\rm i} \tau (\omega - \omega_{Ba} - z_{{\parallel} a} s_{{\parallel} a})}\right). \end{equation}

Equation (A5) allows us to perform analytically the velocity integrals, i.e. the zeroth and first velocity moments of ${g}_a$ (e.g. in $\delta n_{a \phi }^{(j)}$ and $\delta n_{a \psi }^{(j)}$). Using (A5), we derive the analytical expressions of the zeroth and first-order velocity moments of ${g}_a$

(A6a)\begin{align} \delta n_{a \phi}^{(1)} & ={-} \frac{q_a}{\tau_a} \varGamma_0(a_a), \end{align}
(A6b)\begin{align} \delta n_{a \phi}^{(2)} & ={-} \frac{ {\rm i}q_a}{\tau_a} \omega \int_0^\infty \,{\rm d} \tau {\rm e}^{{\rm i} \tau \omega} I_\perp(\tau) I_\parallel(\tau), \end{align}
(A6c)\begin{align} \delta n_{a \phi}^{(3)} & = {\rm i} k_\perp \int_0^\infty \,{\rm d} \tau {\rm e}^{{\rm i} \tau \omega} \left[ R_N I_\parallel (\tau) I_\perp (\tau) \right. \nonumber\\ & \quad \left.+ R_{Ta} \left( I_{{\parallel} }^{(2)} (\tau) I_\perp (\tau) + {\rm I}_{{\parallel} } (\tau) I_{{\perp} }^{(1)} (\tau) - \frac{3}{2} I_\parallel (\tau) {\rm I}_\perp (\tau) \right) \right], \end{align}
(A6d)\begin{align} \delta n_{a \psi}^{(2)} & = {\rm i} \frac{\sqrt{2 }}{\sigma_a} \frac{q_a}{\sqrt{\tau_a}} \int_0^\infty \,{\rm d} \tau \omega {\rm e}^{{\rm i} \tau \omega} I_\perp(\tau) I_{{\parallel} }^{(1)}(\tau), \end{align}
(A6e)\begin{align} \delta n_{a \psi}^{(3)} & ={-}{\rm i} k _\perp\frac{\sqrt{2 \tau_a}}{\sigma_a} \int_0^\infty \,{\rm d} \tau {\rm e}^{{\rm i} \tau \omega} \left[ R_N I_\perp(\tau) I_{{\parallel} }^{(1)}(\tau) \right. \nonumber\\ & \quad \left. + R_{Ta} \left( I_{{\perp} }^{(1)}(\tau) I_{{\parallel} }^{(1)}(\tau ) + {\rm I}_\perp(\tau) I_{{\parallel} }^{(3)}(\tau) - \frac{3}{2} I_\perp(\tau) {\rm I}_{{\parallel} }^{(1)}(\tau) \right)\right], \end{align}

and

(A7a)\begin{align} \delta u_{a \psi}^{(1)} & = \frac{q_a}{\sigma_a^2 } \varGamma_0(a_a), \end{align}
(A7b)\begin{align} \delta u_{a \psi}^{(2)} & = {\rm i} \frac{2 q_a}{\sigma_a^2} \int_0^\infty \,{\rm d} \tau {\rm e}^{{\rm i} \tau \omega} \omega I_\perp(\tau) I_{{\parallel} }^{(2)}(\tau), \end{align}
(A7c)\begin{align} \delta u_{a \psi}^{(3)} & ={-} {\rm i} k_\perp \frac{2 \tau_a}{\sigma_a^2} \int_0^\infty \,{\rm d} \tau {\rm e}^{{\rm i} \tau \omega} \left[ R_N I_\perp(\tau) I_{{\parallel} }^{(2)}(\tau) \right. \nonumber\\ & \quad \left. + R_{Ta} \left( I_{{\perp} }^{(1)}(\tau) I_{{\parallel} }^{(2)}(\tau ) + I_\perp(\tau) I_{{\parallel} }^{(4)}(\tau) - \frac{3}{2} I_\perp(\tau) I_{{\parallel} }^{(2)}(\tau) \right)\right], \end{align}
(A7d)\begin{align} \delta u_{a \phi}^{(1)} & = 0, \end{align}
(A7e)\begin{align} \delta u_{a \phi}^{(2)} & ={-} \frac{{\rm i}q_a \sqrt{2 }}{\sigma_a \sqrt{\tau_a}} \int_0^\infty \,{\rm d} \tau \omega {\rm e}^{{\rm i} \tau \omega} I_\perp (\tau) I_{{\parallel} }^{(1)}(\tau), \end{align}
(A7f)\begin{align} \delta u_{a \phi}^{(3)} & = {\rm i} k_\perp \frac{\sqrt{2 \tau_a}}{\sigma_a} \int_0^\infty \,{\rm d} \tau {\rm e}^{{\rm i} \tau \omega } \left[ R_N I_\perp(\tau) I_{{\parallel} }^{(1)}(\tau) \right. \nonumber\\ & \quad \left. + R_{Ta} \left( I_{{\perp} }^{(1)}(\tau) I_{{\parallel} }^{(1)}(\tau ) + {\rm I}_\perp(\tau) I_{{\parallel} }^{(3)}(\tau) - \frac{3}{2} I_\perp(\tau) {\rm I}_{{\parallel} }^{(1)}(\tau) \right)\right]. \end{align}

The $\tau$-dependent complex functions appearing in (A6) and (A7), which arise from the $s_\parallel$ integration, are given by

(A8a)\begin{gather} I_{{\parallel}}(\tau) = \frac{1}{\sqrt{1 + 2{\rm i} \alpha_a \tau}} \exp\left({ - z_\parallel^2 \tau^2 /4 /(1 + 2 {\rm i} \alpha_a \tau)}\right), \end{gather}
(A8b)\begin{gather}I_{{\parallel}}^{(1)} (\tau) ={-} \frac{{\rm i} \tau z_\parallel}{2 (1 + 2 {\rm i} \tau \alpha_a )^{3/2}} \exp\left({- \tau^2 z_\parallel^2 /4 / (1 + 2 {\rm i} \tau \alpha_a) }\right) , \end{gather}
(A8c)\begin{gather}I_{{\parallel}}^{(2)} (\tau) = \frac{(2(1 + 2 {\rm i} \tau \alpha_a) - \tau^2 z_\parallel^2) }{4 (1 + 2 {\rm i} \tau \alpha_a)^{5/2}} \exp({- z_\parallel^2 \tau^2 / ( 4 (1 + 2 {\rm i} \tau \alpha_a))}). \end{gather}
(A8d)\begin{gather}I_{{\parallel}}^{(3)} (\tau) ={-} \frac{{\rm i} z_\parallel \tau ( 6(1 + 2{\rm i} \alpha_a \tau)- \tau^2 z_\parallel^2)}{8 (1 + 2{\rm i} \alpha_a \tau)^{7/2}} \exp\left({ - \tau^2 z_\parallel^2 / 4 (1 + 2{\rm i} \alpha_a \tau)}\right), \end{gather}
(A8e)\begin{gather}I_{{\parallel} }^{(4)}(\tau) = \frac{(12 (1 + 2 {\rm i} \tau \alpha_a)^2 - 12 (1 + 2 {\rm i} \tau \alpha_a) \tau^2 z_\parallel^2 + z_\parallel^4 \tau^4 )}{16 (1 + 2 {\rm i} \tau \alpha_a)^{9/2}} \exp\left({- \tau^2 z_\parallel^2 / 4 / (1 + 2 {\rm i} \tau \alpha_a)}\right), \end{gather}

while the functions associated with the $x_a$ integration are

(A9a)\begin{align} I_{{\perp}}(\tau) & = \frac{1}{1 + {\rm i}\alpha_a \tau} {\rm I}_0\left(\frac{a_a }{1 + {\rm i}\alpha_a \tau} \right) e ^{{-}a_a /(1 + {\rm i} \alpha_s\tau)}, \end{align}
(A9b)\begin{align} I_{{\perp}}^{(1)} (\tau) & =\frac{\exp({- a_a / (1 + {\rm i} \alpha_a \tau)})}{2 (1 + {\rm i} \alpha_s \tau)^3} \nonumber\\ & \quad \times \left[ (2 (1 + {\rm i} \alpha_s \tau) - 2 a_s) {\rm I}_0 \left( \frac{a_a}{ (1 + {\rm i} \alpha_a \tau)} \right)+ 2 a_a {\rm I}_1 \left( \frac{a_a}{ (1 + {\rm i} \alpha_s \tau)} \right) \right]. \end{align}

The GK dispersion relation given in (A4), with the definitions in (A6) and (A7), constitutes the generalization of the ITG dispersion relation derived in Frei et al. (Reference Frei, Hoffmann and Ricci2022b) to the case of kinetic electrons and electromagnetic effects. We remark that, while the ${\rm I}_0$ and ${\rm I}_1$ functions can be expanded in the case of the electrons using the fact that $a_e \ll a_i \sim 1$, the electron FLR effects are kept here at arbitrary order in $a_e$.

The transformation performed in (A5) restricts the validity of the GK dispersion relation, (A4), to the case of unstable modes, while generalized plasma dispersion functions (Gürcan Reference Gürcan2014; Xie et al. Reference Xie, Li, Lu, Ou and Li2017a; Gültekin & Gürcan Reference Gültekin and Gürcan2018) can be used to include stable modes located in the negative quadrant of the complex plane where $\gamma < 0$.

A.2. Local limit of ITG and KBM

We now solve numerically the local dispersion relation, given in (A4), focusing on the case of electrostatic ITG and KBM. We compare the solution of the GK dispersion relation with the results obtained by solving the GM hierarchy equation, given in (2.25), in the same limit as a function of the number of GMs $(P,J)$.

We first focus on the ITG mode with kinetic electrons in the electrostatic limit. We consider the same values of the density and temperature gradients as in figure 6, and fix the parallel wavenumber at $k_\parallel = 0.1$. We scan over the perpendicular wavenumber $k_\perp = k_y$ and show the results in the top panels of figure 26. It is observed that, while the ITG mode converges with $(P,J) \simeq (16,8)$ for long perpendicular wavelengths, the GM approach requires larger values of $(P,J)$ to resolve FLR effects and magnetic gradient drift effect at smaller perpendicular scales (Frei et al. Reference Frei, Hoffmann and Ricci2022b). An excellent agreement with the local dispersion relation is found for $(P,J) \gtrsim (32,16)$. Additionally, we remark that the case of adiabatic electrons is in good agreement with the local GK dispersion relation with fewer GMs (i.e. $(P,J) = (16,8)$) than the case of non-adiabatic electrons with the same parameters. A scan over the parallel wavenumber at fixed $k_y = 0.4$, displayed in the bottom panels of figure 26, shows that a larger number of GMs is necessary to resolve localized modes in the parallel direction due to Landau damping.

Figure 26. The ITG growth rate $\gamma$ (a,c) and mode frequency $\omega _r$ (b,d) as a function of the binormal wavenumber $k_y$ at $k_\parallel = 0.1$ (a,b) and of the parallel wavenumber $k_\parallel$ at $k_y = 0.4$ (c,d) in the local limit for different numbers of GMs $(P,J)$ (coloured lines). The solution of the collisionless GK dispersion relation, (A4), is plotted (dashed lines). The case of adiabatic electrons (ae) is also shown for comparison. Here, the gradients are the same as in figure 6.

We now consider the case of KBM mode in the local limit by solving (A4) at finite electron plasma pressure, $\beta _e$. The same values of the temperature and density gradients as in figure 11 are used. The top panels of figure 27 shows the KBM growth rate $\gamma$ and mode frequency $\omega _r$ as a function of $\beta _e$ for different number of GMs at $k_y = 0.25$. The solution from the local GK dispersion relation is correctly retrieved by the GM approach and, consistently with the observations made in § 4.3, a fewer number of GMs $(P,J)$ is required than in the ITG case (see figure 26) to achieve convergence. The KBM mode growth rate and frequency are well approached with $(P,J) = (8,4)$. The same can be observed at smaller perpendicular wavelengths by varying the binormal wavenumber $k_y$ at fixed $\beta _e$, as shown in the results plotted in the bottom panels of figure 27. Finally, we remark that the ITG stabilization and KBM onset occurs at an electron plasma pressure (i.e. $\beta _e^c \simeq 0.002$, see figure 27), which is well below the MHD critical value $\beta _e^{{\rm MHD}}$ critical value observed in figure 11 with the same parameters (i.e. $\beta _e^{{\rm MHD}} \simeq 0.013$). This difference in the KBM onset is due to the absence of trapped electrons in the local dispersion relation, which destabilize the ITG mode to values of $\beta _e$ close to the MHD critical value (Weiland & Hirose Reference Weiland and Hirose1992).

Figure 27. The KBM growth rate $\gamma$ (a,c) and real mode frequency $\omega _r$ (b,d) as a function of $\beta _e$ at $k_y = 0.25$ (a,b) and of $k_y$ at $\beta _e = 0.008$ (c,d) obtained from the GM hierarchy (coloured lines) for different $(P,J)$. The analytical results from the collisionless GK dispersion relation, (A4), are shown by the dashed blacked lines. Here, $k_\parallel =0.1$ and the gradients are the same as figure 11.

References

Ajay, C.J., Brunner, S. & Ball, J. 2021 Effect of collisions on non-adiabatic electron dynamics in ITG-driven microturbulence. Phys. Plasmas 28 (9), 092303.Google Scholar
Aleynikova, K. & Zocco, A. 2017 Quantitative study of kinetic ballooning mode theory in simple geometry. Phys. Plasmas 24 (9), 092106.CrossRefGoogle Scholar
Applegate, D.J., Roach, C.M., Connor, J.W., Cowley, S.C., Dorland, W., Hastie, R.J. & Joiner, N. 2007 Micro-tearing modes in the mega ampere spherical tokamak. Plasma Phys. Control. Fusion 49 (8), 1113.CrossRefGoogle Scholar
Ball, J. & Brunner, S. 2021 A non-twisting flux tube for local gyrokinetic simulations. Plasma Phys. Control. Fusion 63 (6), 064008.CrossRefGoogle Scholar
Beer, M.A., Cowley, S.C. & Hammett, G.W. 1995 Field-aligned coordinates for nonlinear simulations of tokamak turbulence. Phys. Plasmas 2 (7), 2687.CrossRefGoogle Scholar
Beer, M.A. & Hammett, G.W. 1996 Toroidal gyrofluid equations for simulations of tokamak turbulence. Phys. Plasmas 3 (11), 4046.CrossRefGoogle Scholar
Belli, E.A. & Candy, J. 2010 Fully electromagnetic gyrokinetic eigenmode analysis of high-beta shaped plasmas. Phys. Plasmas 17 (11), 112314.CrossRefGoogle Scholar
Belli, E.A. & Candy, J. 2011 Full linearized Fokker–Planck collisions in neoclassical transport simulations. Plasma Phys. Control. Fusion 54 (1), 015015.CrossRefGoogle Scholar
Bufferand, H., Bucalossi, J., Ciraolo, G., Falchetto, G., Gallo, A., Ghendrih, P., Rivals, N., Tamain, P., Yang, H., Giorgiani, G., et al. 2021 Progress in edge plasma turbulence modelling–hierarchy of models from 2d transport application to 3D fluid simulations in realistic tokamak geometry. Nucl. Fusion 61 (11), 116052.CrossRefGoogle Scholar
Catto, P.J. & Rosenbluth, M.N. 1981 Trapped electron modifications to tearing modes in the low collision frequency limit. Phys. Fluids 24 (2), 243.CrossRefGoogle Scholar
Chapman, S. & Cowling, T.G. 1941 The velocity of diffusion in a mixed gas; the second approximation. Proc. R. Soc. Lond. A 179 (977), 159.Google Scholar
Churchill, R.M., Chang, C.S., Ku, S. & Dominski, J. 2017 Pedestal and edge electrostatic turbulence characteristics from an XGC1 gyrokinetic simulation. Plasma Phys. Control. Fusion 59 (10), 105014.CrossRefGoogle Scholar
Citrin, J., Garcia, J., Görler, T., Jenko, F., Mantica, P., Told, D., Bourdelle, C., Hatch, D.R., Hogeweij, G.M.D., Johnson, T., et al. 2014 Electromagnetic stabilization of tokamak microturbulence in a high-$\beta$ regime. Plasma Phys. Control. Fusion 57 (1), 014032.CrossRefGoogle Scholar
Connor, J.W., Hastie, R.J. & Helander, P. 2006 Stability of the trapped electron mode in steep density and temperature gradients. Plasma Phys. Control. Fusion 48 (6), 885.CrossRefGoogle Scholar
Connor, J.W., Hastie, R.J. & Taylor, J.B. 1978 Shear, periodicity, and plasma ballooning modes. Phys. Rev. Lett. 40 (6), 396.CrossRefGoogle Scholar
De Oliviera, D.S., Body, T.A., Galassi, D., Theiler, C., Laribi, E., Tamain, P., Stegmeir, A., Giacomin, M., Zholobenko, W. & Ricci, P. 2022 Validation of edge turbulence codes against the TCV-X21 diverted L-mode reference case. Nucl. Fusion 62 (9), 096001.CrossRefGoogle Scholar
Diamond, P.H., Itoh, S.I., Itoh, K. & Hahm, T.S. 2005 Zonal flows in plasma—a review. Plasma Phys. Control. Fusion 47 (5), R35.CrossRefGoogle Scholar
Dickinson, D., Roach, C.M., Saarelma, S., Scannell, R., Kirk, A. & Wilson, H.R. 2012 Kinetic instabilities that limit $\beta$ in the edge of a tokamak plasma: a picture of an h-mode pedestal. Phys. Rev. Lett. 108 (13), 135002.CrossRefGoogle ScholarPubMed
Dickinson, D., Roach, C.M., Saarelma, S., Scannell, R., Kirk, A. & Wilson, H.R. 2013 Microtearing modes at the top of the pedestal. Plasma Phys. Control. Fusion 55 (7), 074006.CrossRefGoogle Scholar
Dimits, A.M., Bateman, G., Beer, M.A., Cohen, B.I., Dorland, W., Hammett, G.W., Kim, C., Kinsey, J.E., Kotschenreuther, M., Kritz, A.H., et al. 2000 Comparisons and physics basis of tokamak transport models and turbulence simulations. Phys. Plasmas 7 (3), 969.CrossRefGoogle Scholar
Doerk, H., Jenko, F., Görler, T., Told, D., Pueschel, M.J. & Hatch, D.R. 2012 Gyrokinetic prediction of microtearing turbulence in standard tokamaks. Phys. Plasmas 19 (5), 055907.CrossRefGoogle Scholar
Dominski, J., Brunner, S., Görler, T., Jenko, F., Told, D. & Villard, L. 2015 How non-adiabatic passing electron layers of linear microinstabilities affect turbulent transport. Phys. Plasmas 22 (6), 062303.CrossRefGoogle Scholar
Dorland, W. 2000 Electron temperature gradient turbulence. Phys. Rev. Lett. 85 (26), 5579.CrossRefGoogle ScholarPubMed
Dorland, W. & Hammett, G.W. 1993 Gyrofluid turbulence models with kinetic effects. Phys. Fluids B 5 (3), 812.CrossRefGoogle Scholar
Dougherty, J.P. 1964 Model Fokker–Planck equation for a plasma and its solution. Phys. Fluids 7 (11), 1788.CrossRefGoogle Scholar
Ernst, D.R., Lang, J., Nevins, W.M., Hoffman, M., Chen, Y., Dorland, W. & Parker, S. 2009 Role of zonal flows in trapped electron mode turbulence through nonlinear gyrokinetic particle and continuum simulation. Phys. Plasmas 16 (5), 055906.CrossRefGoogle Scholar
Frei, B.J., Ball, J., Hoffmann, A.C.D., Jorge, R., Ricci, P. & Stenger, L. 2021 Development of advanced linearized gyrokinetic collision operators using a moment approach. J. Plasma Phys. 87, 905870501.CrossRefGoogle Scholar
Frei, B.J., Ernst, S. & Ricci, P. 2022 a Numerical implementation of the improved sugama collision operator using a moment approach. Phys. Plasmas 29 (9), 093902.CrossRefGoogle Scholar
Frei, B.J., Hoffmann, A.C.D. & Ricci, P. 2022 b Local gyrokinetic collisional theory of the ion-temperature gradient mode. J. Plasma Phys. 88 (3), 905880304.CrossRefGoogle Scholar
Frei, B.J., Jorge, R. & Ricci, P. 2020 A gyrokinetic model for the plasma periphery of tokamak devices. J. Plasma Phys. 86 (2), 905860205.CrossRefGoogle Scholar
Fulton, D.P., Lin, Z., Holod, I. & Xiao, Y. 2014 Microturbulence in DIII-D tokamak pedestal. I. Electrostatic instabilities. Phys. Plasmas 21 (4), 042110.CrossRefGoogle Scholar
Galassi, D., Theiler, C., Body, T., Manke, F., Micheletti, P., Omotani, J., Wiesenberger, M., Baquero-Ruiz, M., Furno, I., Giacomin, M., et al. 2022 Validation of edge turbulence codes in a magnetic X-point scenario in TORPEX. Phys. Plasmas 29 (1), 012501.CrossRefGoogle Scholar
Gao, Z. 2013 Collisional damping of the geodesic acoustic mode. Phys. Plasmas 20 (3), 032501.CrossRefGoogle Scholar
Garcia, J., De La Luna, E., Sertoli, M., Casson, F.J., Mazzi, S., Štancar, Ž., Szepesi, G., Frigione, D., Garzotti, L., Rimini, F., et al. 2022 New H-mode regimes with small ELMs and high thermal confinement in the Joint European Torus. Phys. Plasmas 29 (3), 032505.CrossRefGoogle Scholar
Giacomin, M., Stenger, L.N. & Ricci, P. 2020 Turbulence and flows in the plasma boundary of snowflake magnetic configurations. Nucl. Fusion 60 (2), 024001.CrossRefGoogle Scholar
Görler, T., Lapillonne, X., Brunner, S., Dannert, T., Jenko, F., Aghdam, S.K., Marcus, P., McMillan, B.F., Merz, F., Sauter, O., et al. 2011 Flux-and gradient-driven global gyrokinetic simulation of tokamak turbulence. Phys. Plasmas 18 (5), 056103.CrossRefGoogle Scholar
Gradshteyn, I.S. & Ryzhik, I.M. 2014 Table of Integrals, Series, and Products. Academic Press.Google Scholar
Grant, F.C. & Feix, M.C. 1967 Fourier–Hermite solutions of the Vlasov equations in the linearized limit. Phys. Fluids 10 (4), 696.CrossRefGoogle Scholar
Gültekin, Ö. & Gürcan, Ö.D. 2018 Stable and unstable roots of ion temperature gradient driven mode using curvature modified plasma dispersion functions. Plasma Phys. Control. Fusion 60 (2), 025021.CrossRefGoogle Scholar
Gürcan, Ö.D. 2014 Numerical computation of the modified plasma dispersion function with curvature. J. Comput. Phys. 269, 156.CrossRefGoogle Scholar
Hallatschek, K. & Dorland, W. 2005 Giant electron tails and passing electron pinch effects in tokamak-core turbulence. Phys. Rev. Lett. 95 (5), 055002.CrossRefGoogle ScholarPubMed
Hammett, G.W., Beer, M.A., Dorland, W., Cowley, S.C. & Smith, S.A. 1993 Developments in the gyrofluid approach to tokamak turbulence simulations. Plasma Phys. Control. Fusion 35 (8), 973.CrossRefGoogle Scholar
Han, M.K., Wang, Z.-X., Dong, J.Q. & Du, H. 2017 Multiple ion temperature gradient driven modes in transport barriers. Nucl. Fusion 57 (4), 046019.CrossRefGoogle Scholar
Hatch, D.R., Kotschenreuther, M., Mahajan, S., Valanju, P., Jenko, F., Told, D., Görler, T. & Saarelma, S. 2016 Microtearing turbulence limiting the JET-ILW pedestal. Nucl. Fusion 56 (10), 104003.CrossRefGoogle Scholar
Hazeltine, R.D. & Meiss, J.D. 2003 Plasma Confinement. Courier Corporation.Google Scholar
Hazeltine, R.D. & Strauss, H.R. 1976 Tokamak heat transport due to tearing modes. Phys. Rev. Lett. 37 (2), 102.CrossRefGoogle Scholar
Helander, P. & Sigmar, D.J. 2002 Collisional Transport in Magnetized Plasmas. Cambridge University Press.Google Scholar
Hinton, F.L. & Rosenbluth, M.N. 1999 Dynamics of axisymmetric ($E \times B$) and poloidal flows in tokamaks. Plasma Phys. Control. Fusion 41 (3A).CrossRefGoogle Scholar
Hirshman, S.P. & Sigmar, D.J. 1976 Approximate Fokker–Planck collision operator for transport theory applications. Phys. Fluids 19 (10), 1532.CrossRefGoogle Scholar
Hoffmann, A.C.D., Frei, B. J. & Ricci, P. 2023 Gyrokinetic simulations of plasma turbulence in a Z-pinch configuration using a moment approach and advanced collision operators. J. Plasma Phys. 89 (2), 905890214.CrossRefGoogle Scholar
Holland, C., Schmitz, L., Rhodes, T.L., Peebles, W.A., Hillesheim, J.C., Wang, G., Zeng, L., Doyle, E.J., Smith, S.P., Prater, R., et al. 2011 Advances in validating gyrokinetic turbulence models against l-and h-mode plasmas. Phys. Plasmas 18 (5), 056113.CrossRefGoogle Scholar
Idomura, Y., Ida, M., Kano, T., Aiba, N. & Tokuda, S. 2008 Conservative global gyrokinetic toroidal full-f five-dimensional Vlasov simulation. Comput. Phys. Commun. 179 (6), 391.CrossRefGoogle Scholar
Jenko, F., Dorland, W., Kotschenreuther, M. & Rogers, B.N. 2000 Electron temperature gradient driven turbulence. Phys. Plasmas 7 (5), 1904.CrossRefGoogle Scholar
Jorge, R., Frei, B.J. & Ricci, P. 2019 Nonlinear gyrokinetic Coulomb collision operator. J. Plasma Phys. 85 (6), 905850604.CrossRefGoogle Scholar
Jorge, R., Ricci, P. & Loureiro, N.F. 2017 A drift-kinetic analytical model for scrape-off layer plasma dynamics at arbitrary collisionality. J. Plasma Phys. 83 (6), 905830606.CrossRefGoogle Scholar
Jorge, R., Ricci, P. & Loureiro, N.F. 2018 Theory of the drift-wave instability at arbitrary collisionality. Phys. Rev. Lett. 121 (16), 165001.CrossRefGoogle ScholarPubMed
Kotschenreuther, M., Hatch, D.R., Mahajan, S., Valanju, P., Zheng, L. & Liu, X. 2017 Pedestal transport in h-mode plasmas for fusion gain. Nucl. Fusion 57 (6), 064001.CrossRefGoogle Scholar
Landau, L. 1936 Die kinetische Gleichung für den Fall Coulombscher Wechselwirkung. Phys. Z. Sowjetunion 10, 154.Google Scholar
Lapillonne, X., Brunner, S., Dannert, T., Jolliet, S., Marinoni, A., Villard, L., Görler, T., Jenko, F. & Merz, F. 2009 Clarifications to the limitations of the s-$\alpha$ equilibrium model for gyrokinetic computations of turbulence. Phys. Plasmas 16 (3), 032308.CrossRefGoogle Scholar
Lebedev, V.B., Yushmanov, P.N., Diamond, P.H., Novakovskii, S.V. & Smolyakov, A.I. 1996 Plateau regime dynamics of the relaxation of poloidal rotation in tokamak plasmas. Phys. Plasmas 3 (8), 3023.CrossRefGoogle Scholar
Lee, W.W. 1987 Gyrokinetic particle simulation model. J. Comput. Phys. 72 (1), 243.CrossRefGoogle Scholar
Li, B. & Ernst, D.R. 2011 Gyrokinetic Fokker–Planck collision operator. Phys. Rev. Lett. 106 (19), 195002.CrossRefGoogle ScholarPubMed
Li, Y. & Gao, Z. 2015 Comparison of collision operators for the geodesic acoustic mode. Nucl. Fusion 55 (4), 043001.CrossRefGoogle Scholar
Lin, Z., Nishimura, Y., Xiao, Y., Holod, I., Zhang, W.L. & Chen, L. 2007 Global gyrokinetic particle simulations with kinetic electrons. Plasma Phys. Control. Fusion 49 (12B), B163.CrossRefGoogle Scholar
Loureiro, N.F., Schekochihin, A.A. & Zocco, A. 2013 Fast collisionless reconnection and electron heating in strongly magnetized plasmas. Phys. Rev. Lett. 111 (2), 025002.CrossRefGoogle ScholarPubMed
Madsen, J. 2013 Full-f gyrofluid model. Phys. Plasmas 20 (7), 072301.CrossRefGoogle Scholar
Mandell, N., Dorland, W., Abel, I., Gaur, R., Kim, P., Martin, M. & Qian, T. 2022 GX: a GPU-native gyrokinetic turbulence code for tokamaks and stellarators. arXiv:2209.06731.Google Scholar
Mandell, N.R., Dorland, W. & Landreman, M. 2018 Laguerre–Hermite pseudo-spectral velocity formulation of gyrokinetics. J. Plasma Phys. 84 (1), 905840108.CrossRefGoogle Scholar
Mandell, N.R., Hakim, A., Hammett, G.W. & Francisquez, M. 2020 Electromagnetic full-gyrokinetics in the tokamak edge with discontinuous Galerkin methods. J. Plasma Phys. 86 (1), 905860109.CrossRefGoogle Scholar
Merlo, G., Sauter, O., Brunner, S., Burckel, A., Camenen, Y., Casson, F.J., Dorland, W., Fable, E., Görler, T., Jenko, F., et al. 2016 Linear multispecies gyrokinetic flux tube benchmarks in shaped tokamak plasmas. Phys. Plasmas 23 (3), 032104.CrossRefGoogle Scholar
Michels, D., Stegmeir, A., Ulbl, P., Jarema, D. & Jenko, F. 2021 Gene-x: a full-f gyrokinetic turbulence code based on the flux-coordinate independent approach. Comput. Phys. Commun. 264, 107986.CrossRefGoogle Scholar
Navarro, A.B., Happel, T., Görler, T., Jenko, F., Abiteboul, J., Bustos, A., Doerk, H., Told, D. & ASDEX Upgrade Team 2015 Gyrokinetic studies of core turbulence features in ASDEX Upgrade H-mode plasmas. Phys. Plasmas 22 (4), 042513.CrossRefGoogle Scholar
Neiser, T.F., Jenko, F., Carter, T.A., Schmitz, L., Told, D., Merlo, G., Bañón Navarro, A., Crandall, P.C., McKee, G.R. & Yan, Z. 2019 Gyrokinetic GENE simulations of DIII-D near-edge L-mode plasmas. Phys. Plasmas 26 (9), 092510.CrossRefGoogle Scholar
Novakovskii, S.V., Liu, C.S., Sagdeev, R.Z. & Rosenbluth, M.N. 1997 The radial electric field dynamics in the neoclassical plasmas. Phys. Plasmas 4 (12), 4272.CrossRefGoogle Scholar
Pan, Q. & Ernst, D.R. 2019 Gyrokinetic landau collision operator in conservative form. Phys. Rev. E 99 (2), 023201.CrossRefGoogle ScholarPubMed
Pan, Q., Ernst, D.R. & Crandall, P. 2020 First implementation of gyrokinetic exact linearized Landau collision operator and comparison with models. Phys. Plasmas 27 (4), 042307.CrossRefGoogle Scholar
Pan, Q., Ernst, D.R. & Hatch, D.R. 2021 Importance of gyrokinetic exact Fokker–Planck collisions in fusion plasma turbulence. Phys. Rev. E 103 (5), L051202.CrossRefGoogle ScholarPubMed
Peeters, A.G., Camenen, Y., Casson, F.J., Hornsby, W.A., Snodin, A.P., Strintzi, D. & Szepesi, G. 2009 The nonlinear gyro-kinetic flux tube code GKW. Comput. Phys. Commun. 180 (12), 2650.CrossRefGoogle Scholar
Pueschel, M.J., Hatch, D.R., Ernst, D.R., Guttenfelder, W., Terry, P.W., Citrin, J. & Connor, J.W. 2019 On microinstabilities and turbulence in steep-gradient regions of fusion devices. Plasma Phys. Control. Fusion 61 (3), 034002.CrossRefGoogle Scholar
Rafiq, T., Pankin, A.Y., Bateman, G., Kritz, A.H. & Halpern, F.D. 2009 Simulation of electron thermal transport in H-mode discharges. Phys. Plasmas 16 (3), 032505.CrossRefGoogle Scholar
Rosenbluth, M.N. & Hinton, F.L. 1998 Poloidal flow driven by ion-temperature-gradient turbulence in tokamaks. Phys. Rev. Lett. 80 (4), 724.CrossRefGoogle Scholar
Schekochihin, A.A., Parker, J.T., Highcock, E.G. & Dellar, P.J. 2016 Phase mixing versus nonlinear advection in drift-kinetic plasma turbulence. J. Plasma Phys. 82 (2), 905820212.CrossRefGoogle Scholar
Schmitz, L., Holland, C., Rhodes, T.L., Wang, G., Zeng, L., White, A.E., Hillesheim, J.C., Peebles, W.A., Smith, S.P., Prater, R., et al. 2012 Reduced electron thermal transport in low collisionality H-mode plasmas in DIII-D and the importance of TEM/ETG-scale turbulence. Nucl. Fusion 52 (2), 023003.CrossRefGoogle Scholar
Snyder, P.B., Groebner, R.J., Leonard, A.W., Osborne, T.H. & Wilson, H.R. 2009 Development and validation of a predictive model for the pedestal height. Phys. Plasmas 16 (5), 056118.CrossRefGoogle Scholar
Staebler, G.M., Kinsey, J.E. & Waltz, R.E. 2005 Gyro–Landau fluid equations for trapped and passing particles. Phys. Plasmas 12 (10), 102508.CrossRefGoogle Scholar
Stegmeir, A., Ross, A., Body, T., Francisquez, M., Zholobenko, W., Coster, D., Maj, O., Manz, P., Jenko, F., Rogers, B.N., et al. 2019 Global turbulence simulations of the tokamak edge region with grillix. Phys. Plasmas 26 (5), 052517.CrossRefGoogle Scholar
Sugama, H., Matsuoka, S., Satake, S., Nunami, M. & Watanabe, T.-H. 2019 Improved linearized model collision operator for the highly collisional regime. Phys. Plasmas 26 (10), 102108.CrossRefGoogle Scholar
Sugama, H., Watanabe, T.-H. 2006 Collisionless damping of geodesic acoustic modes. J. Plasma Phys. 72 (6), 825.CrossRefGoogle Scholar
Sugama, H., Watanabe, T.-H. & Nunami, M. 2009 Linearized model collision operators for multiple ion species plasmas and gyrokinetic entropy balance equations. Phys. Plasmas 16 (11), 112503.CrossRefGoogle Scholar
Tang, W.M., Connor, J.W. & Hastie, R.J. 1980 Kinetic-ballooning-mode theory in general geometry. Nucl. Fusion 20 (11), 1439.CrossRefGoogle Scholar
Told, D., Jenko, F., Xanthopoulos, P., Horton, L.D., Wolfrum, E. & ASDEX Upgrade Team 2008 Gyrokinetic microinstabilities in ASDEX Upgrade edge plasmas. Phys. Plasmas 15 (10), 102306.CrossRefGoogle Scholar
Tronko, N., Bottino, A., Görler, T., Sonnendrücker, E., Told, D. & Villard, L. 2017 Verification of gyrokinetic codes: theoretical background and applications. Phys. Plasmas 24 (5), 056115.CrossRefGoogle Scholar
Wan, W., Parker, S.E., Chen, Y., Yan, Z., Groebner, R.J. & Snyder, P.B. 2012 Global gyrokinetic simulation of tokamak edge pedestal instabilities. Phys. Rev. Lett. 109 (18), 185004.CrossRefGoogle ScholarPubMed
Weiland, J. & Hirose, A. 1992 Electromagnetic and kinetic effects on the ion temperature gradient mode. Nucl. Fusion 32 (1), 151.CrossRefGoogle Scholar
Winsor, N., Johnson, J.L. & Dawson, J.M. 1968 Geodesic acoustic waves in hydromagnetic systems. Phys. Fluids 11 (11), 2448.CrossRefGoogle Scholar
Xanthopoulos, P. & Jenko, F. 2006 Clebsch-type coordinates for nonlinear gyrokinetics in generic toroidal configurations. Phys. Plasmas 13 (9), 092301.CrossRefGoogle Scholar
Xiao, Y. & Catto, P.J. 2006 Short wavelength effects on the collisionless neoclassical polarization and residual zonal flow level. Phys. Plasmas 13 (10), 102311.CrossRefGoogle Scholar
Xiao, Y., Catto, P.J. & Molvig, K. 2007 Collisional damping for ion temperature gradient mode driven zonal flow. Phys. Plasmas 14 (3), 032302.CrossRefGoogle Scholar
Xie, H.-S. & Li, B. 2016 Global theory to understand toroidal drift waves in steep gradient. Phys. Plasmas 23 (8), 082513.CrossRefGoogle Scholar
Xie, H.S., Li, Y.Y., Lu, Z.X., Ou, W.K. & Li, B. 2017 a Comparisons and applications of four independent numerical approaches for linear gyrokinetic drift modes. Phys. Plasmas 24 (7), 072106.CrossRefGoogle Scholar
Xie, H.-S., Lu, Z.-X. & Li, B. 2018 Kinetic ballooning mode under steep gradient: high order eigenstates and mode structure parity transition. Phys. Plasmas 25 (7), 072106.CrossRefGoogle Scholar
Xie, H.-S. & Xiao, Y. 2015 Unconventional ballooning structures for toroidal drift waves. Phys. Plasmas 22 (9), 090703.CrossRefGoogle Scholar
Xie, H.-S., Xiao, Y. & Lin, Z. 2017 b New paradigm for turbulent transport across a steep gradient in toroidal plasmas. Phys. Rev. Lett. 118 (9), 095001.CrossRefGoogle ScholarPubMed
Zeiler, A., Drake, J.F. & Rogers, B. 1997 Nonlinear reduced Braginskii equations with ion thermal dynamics in toroidal plasma. Phys. Plasmas 4 (6), 2134.CrossRefGoogle Scholar
Figure 0

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

Figure 1

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.

Figure 2

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$.

Figure 3

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

Figure 4

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.

Figure 5

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$.

Figure 6

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 ac). 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}$.

Figure 7

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}$ (ac) and of the TEM at $k_y = 1.3$ (df), 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 (cf) for different values of $(P,J)$ and GENE. The same parameters as in figure 7 are used.

Figure 8

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.

Figure 9

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.

Figure 10

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.

Figure 11

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).

Figure 12

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.

Figure 13

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 )$.

Figure 14

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).

Figure 15

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

Figure 16

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.

Figure 17

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$.

Figure 18

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.

Figure 19

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$.

Figure 20

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$.

Figure 21

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$.

Figure 22

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

Figure 23

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

Figure 24

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

Figure 25

Figure 26. The ITG growth rate $\gamma$ (a,c) and mode frequency $\omega _r$ (b,d) as a function of the binormal wavenumber $k_y$ at $k_\parallel = 0.1$ (a,b) and of the parallel wavenumber $k_\parallel$ at $k_y = 0.4$ (c,d) in the local limit for different numbers of GMs $(P,J)$ (coloured lines). The solution of the collisionless GK dispersion relation, (A4), is plotted (dashed lines). The case of adiabatic electrons (ae) is also shown for comparison. Here, the gradients are the same as in figure 6.

Figure 26

Figure 27. The KBM growth rate $\gamma$ (a,c) and real mode frequency $\omega _r$ (b,d) as a function of $\beta _e$ at $k_y = 0.25$ (a,b) and of $k_y$ at $\beta _e = 0.008$ (c,d) obtained from the GM hierarchy (coloured lines) for different $(P,J)$. The analytical results from the collisionless GK dispersion relation, (A4), are shown by the dashed blacked lines. Here, $k_\parallel =0.1$ and the gradients are the same as figure 11.