Local gyrokinetic collisional theory of the ion-temperature gradient mode

Abstract

We present a study of the linear properties of the ion-temperature gradient (ITG) modes with collisions modelled for the first time by the linearized gyrokinetic (GK) Coulomb collision operator (Frei et al., J. Plasma Phys., vol. 87, issue 5, 2021, 905870501) in the local limit. The study is based on a Hermite–Laguerre polynomial expansion of the perturbed ion distribution function applied to the linearized GK Boltzmann equation, yielding a hierarchy of coupled equations for the expansion coefficients, referred to as gyromoments. We explore the collisionless and high-collisional limits of the gyromoment hierarchy analytically. Parameter scans revealing the dependence of the ITG growth rate on the collisionality modelled using the GK Coulomb operator are reported, showing strong damping at small scales as the collisionality increases and, therefore, the need for a steeper gradient for the ITG onset at high collisionality to overcome the finite Larmor radius (FLR) collisional stabilization. The predictions on the ITG growth rate by the GK Coulomb operator are compared with other collision operator models, such as the Sugama, the Dougherty, as well as the momentum-conserving pitch-angle scattering and the Hirshman–Sigmar–Clarke collision operators derived for the first time in terms of gyromoments. The importance of FLR terms in the collision operators is pointed out by the appearance of a short wavelength ITG branch when collisional FLR terms are neglected, this branch being completely suppressed by FLR collisional effects. Energy diffusion is shown to be important at high collisionality and at small scale lengths. Among the GK collision operators considered in this work, the GK Sugama collision operator yields, in general, the smallest deviation compared with the GK Coulomb collision operator, while the largest deviations are found with the GK Dougherty operator. Convergence studies of the gyromoment method are reported and show that the drifts associated with the gradient and curvature of the magnetic field increase the required number of gyromoments at low collisionality. Nevertheless, the low number of gyromoments necessary for convergence at high collisionality constitutes an attractive numerical and analytical feature of the gyromoment approach to study the plasma dynamics in the boundary of fusion devices.

1. Introduction

Belonging to the class of instabilities that develop on the ion gyroscale, the ion-temperature gradient (ITG) modes, driven unstable by the parallel plasma compression (slab ITG (sITG)) or by the presence of magnetic drifts (toroidal ITG), are widely recognized as the main candidate to explain the experimental observations of anomalous ion heat turbulent transport in the tokamak core (Garbet et al. 2004). Because of its important role and since the plasma collision frequency is considerably smaller than the typical ITG mode frequency in the core region such that collisions can be neglected, the gyrokinetic (GK) collisionless theory of the ITG modes is well documented in the literature (see, e.g. Hahm & Tang 1989; Romanelli 1989; Romanelli & Briguglio 1990; Chen, Briguglio & Romanelli 1991). On the other hand and because of the temperature drop, the effects of collisions on the ITG modes are expected to be important in the tokamak boundary, the region that encompasses the edge and the scrape-off layer (SOL). In fact, and more generally, collisions are known to modify the linear and nonlinear properties of microinstabilities. For instance, collisions contribute to smearing out small scale structures in velocity space (Watanabe & Sugama 2004; Ajay, Brunner & Ball 2021), to affect their linear properties (Manas et al. 2015; Belli & Candy 2017), to suppress short-wavelength structures (Barnes et al. 2009) and can, ultimately, affect their saturation mechanism via the damping of the zonal flow (Pan, Ernst & Crandall 2020; Frei et al. 2021). Since the boundary region sets the heat exhaust and the thermal load on the vessel walls and determines the overall confinement capabilities of the device, it is of primary importance to assess the role of the ITG modes in this key region that is characterized by a wide range of plasma collisionality. For instance, the collisionality dependence of the ITG mode is important in the prediction of, e.g. the transition from resistive ballooning to ITG driven turbulence regime in the boundary region (Zeiler et al. 1998; Dorf & Dorr 2022).

In the past, the effects of finite collisionality on the ITG linear properties have been introduced in kinetic models using simplified collision operators such as, e.g. an energy-conserving Krook and a pitch-angle operator model (Romanelli & Briguglio 1990; Romanelli, Chen & Briguglio 1991; Pusztai et al. 2009). More generally, the effect of collisions on microinstabilities (including the trapped electron and ITG modes) have been documented in previous studies using GK codes (Catto & Ernst 2009; Belli & Candy 2017; Pan & Ernst 2019; Crandall et al. 2020; Pan et al. 2020) showing discrepancies produced by the choice of collision operator models when compared with the linearized Fokker–Planck collision operator (Rosenbluth, MacDonald & Judd 1957; Hazeltine & Meiss 2003), that we refer to as the Coulomb collision operator in the present study. However, studies based on the GK Coulomb collision operator remain very limited (Pan et al. 2020), due to the very recent availability of such an operator (Frei et al. 2021). Additionally, we lack a comparison between this operator and simplified collision operator models for the ITG case.

The limited number of collisional GK investigations stems from the frequent use of drift-reduced Braginskii fluid models (Zeiler, Drake & Rogers 1997) for the simulation of the boundary region (Dudson et al. 2009; Halpern et al. 2016; Tamain et al. 2016; Paruta et al. 2018; Zhu, Francisquez & Rogers 2018; Stegmeir et al. 2019; Giacomin, Stenger & Ricci 2020). These fluid models are based on the assumptions that $k_\parallel \lambda _{{\rm mfp}} \ll 1$, i.e. that the particle mean free path $\lambda _{{\rm mpf}}$ is much shorter than the typical parallel length scales, described by the parallel wavenumber $k_\parallel$, and that $k_\perp \rho _i \ll 1$, i.e. that the modes develop a perpendicular wavenumber $k_\perp$ associated with a scale length much larger than the ion gyroradius $\rho _i$. Therefore, while the drift-reduced Braginskii models retain finite ITG effects (Hallatschek & Zeiler 2000; Mosetto et al. 2015), their applicability to the study of ITG modes is questionable since, for instance, these modes develop on a $k_\perp \rho _i \sim 1$ scale. In addition, while the $k_\parallel \lambda _{{\rm mfp}} \ll 1$ assumption is well established at the separatrix and in the SOL (with temperature $T \lesssim 50$ eV and particle density $n \sim 10^{18}$ m$^{-3}$ yielding typically $k_\parallel \lambda _{{\rm mfp}} \lesssim 10^{-2}$ in, e.g. JET, ITER and smaller tokamaks such as TCV (Frassinetti et al. 2019)), the mean free path is larger inside the separatrix, where the temperature increases. Hence, for typical H-mode pedestal values ($T \sim 10^{2}\text{--}10^{3}$ eV and $n \sim 10^{20}$ m$^{-3}$ (Leyland et al. 2014; Schneider et al. 2021)), an estimate of $k_\parallel \lambda _{{\rm mfp}} \sim 1$ can be obtained. Thus, one expects a level of collisionality sufficient to justify the use of the drift-reduced Braginskii fluid models only in the SOL and, possibly, in the outermost region of the edge in typical tokamak conditions.

Gyrofluid models have been developed by considering a finite number of velocity space moments of the GK equation to extend the validity of Braginskii-like fluid models (Hammett, Dorland & Perkins 1992; Dorland & Hammett 1993; Beer & Hammett 1996; Snyder & Hammett 2001). Collisionless kinetic effects are introduced in these models by designing ad hoc closures for the high-order moments that mimic the linear kinetic response in the absence of collisions. While usually neglected, collisional effects are introduced in, for example, the gyrofluid model proposed by Snyder & Hammett (2001) by considering a particle, momentum and energy-conserving BGK-like operator for ion collisions and a pitch-angle scattering operator for the collisions of the electrons. However, BGK-like operators and the absence of energy diffusion in the pitch-angle scattering operator can potentially lead to significant deviations compared with the Coulomb collision operator.

In this work, we aim to present a study of the linear properties of the ITG mode within a GK framework able to capture accurately the collisionless kinetic physics and, at the same time, the collisional effects described through the GK Coulomb collision operator for the first time. This is done by leveraging previous work that led to a formulation of the GK model for the boundary region based on a Hermite–Laguerre expansion of the gyrocentre distribution function (Frei, Jorge & Ricci 2020), and the development of a GK Coulomb collision operator derived using the same technique expansion, as presented in Jorge, Frei & Ricci (2019a). The present study represents the first physical application of this model, which has the potential of providing an efficient description of plasma turbulence in the boundary region. In fact, the Hermite–Laguerre expansion technique shows its ability to treat the ITG mode, numerically and analytically, for a wide range of plasma collisionality and with sophisticated collision operators (Jorge, Ricci & Loureiro 2018; Jorge et al. 2019b). More precisely, the ITG model we use in this work is based on the linearized (also known as $\delta f$) and local limits of the gyromoment formulation of the GK model coupled with a linearized GK Coulomb collision operator, recently derived and numerically implemented in Frei et al. (2021). The Hermite–Laguerre expansion allows us to express the self-consistent GK model in terms of moments of the gyrocentre distribution function that we refer to as gyromoments.

By using the gyromoment approach, we first retrieve the collisionless and the high-collisional limits of the ITG mode. In particular, the ability of the gyromoment method to be reduced to simpler fluid models at high collisionality is analytically demonstrated and numerically tested. Then, using the linearized GK Coulomb collision operator, wide parameter scans of the ITG growth rate are performed. It is demonstrated that a steeper temperature gradient is required for the ITG onset at high collisionality (although higher than typical SOL experimental values) to overcome the finite Lramor radius (FLR) collisional stabilization effects. We also compare, for the first time, the GK Coulomb collision operator with the drift-kinetic (DK) Coulomb collision operator, a GK/DK momentum-conserving pitch-angle scattering operator (Helander & Sigmar 2002), the zeroth-order DK Hirshman–Sigmar–Clarke (HSC) collision operator (Hirshman & Sigmar 1976a) and the DK/GK Dougherty collision operator (Dougherty 1964), in addition to the DK/GK Sugama collision operator implemented in Frei et al. (2021). Among the considered collision operator models, the GK Sugama operator yields the smallest deviation with respect to the GK Coulomb collision, while the largest deviations are found when the GK Dougherty collision operator is considered. We show the importance of energy diffusion and FLR terms in the collisional damping of ITG at small scales. Also, the importance of retaining FLR terms in the collision operator models is demonstrated by revealing that a short wavelength ITG (SWITG) mode at high collisionality and steep ITGs can be destabilized when the DK limit of the collision operators is used. We remark that the assessment of the potential deviations between collision operator models is of primary importance, since simplified and DK collision operators are often used in turbulent GK codes (Bernard et al. 2019; Francisquez et al. 2021; Ulbl, Michels & Jenko 2021). Finally, we investigate the convergence properties of the gyromoment approach as a function of relevant physical parameters in different collisionality regimes, analysing the Hermite–Laguerre spectrum of the ITG distribution function and showing that magnetic drifts broaden the gyromoment spectrum significantly in the absence of collisions. The capability of reducing the number of gyromoments with increasing collisionality is illustrated by the present study, a key feature of the present approach to simulate the boundary region.

The structure of the paper is the following. In § 2 we introduce the gyromoment hierarchy equation. Then, in § 3, we derive the collisionless limit of the ITG dispersion relation by using the gyromoment hierarchy equation, and § 4 focuses on the high-collisional limits. In § 5, we evaluate the ITG growth rate over a wide range of parameters with collisions modelled by the GK Coulomb collision operator. Section 6 reports on the comparison of collision operator models with the GK Coulomb collision operator. Finally, we analyse the Hermite–Laguerre spectrum and report convergence studies in §§ 7 and 8, respectively. The conclusions follow in § 9. The Hermite–Laguerre expansions of the collision operators used in this work are reported in appendices A–C, while analytical results are detailed in Appendix D.

2. Gyrokinetic model

Focusing on a shearless configuration and adopting a local GK approach where quantities do not vary along the magnetic field lines (Kadomtsev & Shafranov 2012), we consider the linearized electrostatic GK Boltzmann equation expressed in gyrocentre phase-space coordinates $(\boldsymbol {R}, \mu, v_\parallel,\theta )$ and assume that the electrons are adiabatic. Here, $\boldsymbol {R} = \boldsymbol {r} - \boldsymbol \rho$ is the ion gyrocentre position (with $\boldsymbol {r}$ being the particle position and $\boldsymbol \rho (\boldsymbol {R}, \mu, \theta ) = \boldsymbol {b} \times \boldsymbol v/\varOmega$ the ion gyroradius, with $\varOmega = q B / m$ the ion gyrofrequency, $q = Z e$ the ion charge with $Z =1$ and $\boldsymbol {b} = \boldsymbol {B} / B$), $\mu = m v^{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 magnetic field and $\theta$ is the gyroangle. Since only the evolution of the ion distribution is considered, we drop the ion species label for simplicity. We denote the equilibrium gradient scale lengths of ion density $N$, ion temperature $T$ and magnetic field strength $B$, by $L_N = | \boldsymbol {\nabla }_\perp \ln N |^{-1}$, $L_{T} = | \boldsymbol {\nabla }_\perp \ln T |^{-1}$ and $L_B = | \boldsymbol {\nabla }_\perp \ln B |^{-1}$, respectively. We assume that the gyroaveraged ion gyrocentre distribution function, $F = F (\boldsymbol {R}, \mu, v_\parallel )$, is a perturbed Maxwellian, that is $F = F_{M} + {g}$, where ${g}$ is the perturbed part (${g}/ F_{M} \ll 1$) and $F_{M}$ a local Maxwellian distribution function defined by $F_{M} = N(\boldsymbol {R}) {\rm e}^{- s_\parallel ^{2} - x} / ({\rm \pi} ^{3/2} v_{T}(\boldsymbol {R})^{3})$ with $s_\parallel = v_\parallel / v_{T}(\boldsymbol {R})$, $x= \mu B(\boldsymbol {R}) / T(\boldsymbol {R})$ and $v_T(\boldsymbol {R})^{2} = 2 T(\boldsymbol {R})/m$ the ion thermal velocity. The linearized electrostatic GK Boltzmann equation describing the time evolution of a single Fourier harmonic ${g}(\boldsymbol k,t) = \int \,{\rm d} \boldsymbol {R} {g}(\boldsymbol {R} ,t) {\rm e}^{- {\rm i} \boldsymbol R \boldsymbol {\cdot } \boldsymbol k}$ is then given by

(2.1)\begin{gather} \frac{\partial}{\partial t } {g} + \left({\rm i} \omega_{B} + {\rm i} k_\parallel v_\parallel \right) {g} - {\rm i} \omega_{*}^{T} {\rm J}_0(b\sqrt{x}) \frac{e \phi}{T_e} F_{M} + {\rm i} \omega_{\boldsymbol{\nabla} B} {\rm J}_0(b\sqrt{x}) \frac{q \phi }{T} F_{M}\nonumber\\ + \frac{q}{T} {\rm i} k_\parallel v_\parallel {\rm J}_0(b\sqrt{x}) \phi F_{M} = \mathcal{C}, \end{gather}

where $e$ is the electron charge and $\tau = T / T_e$. In (2.1), we express the wavevector as $\boldsymbol k = k_\parallel \boldsymbol {b} + \boldsymbol k_\perp$, and introduce the velocity-dependent resonant frequencies $\omega _{\boldsymbol {\nabla } B} = v_{T}^{2} ( \boldsymbol b \times \boldsymbol {\nabla } B /B) \boldsymbol {\cdot } \boldsymbol k /(2 \varOmega )(x + 2 s_\parallel ^{2})$ and $\omega _{*}^{T} = \omega _*[1 + \eta ( x + s_\parallel ^{2} - 3/2)]$, with the ion diamagnetic frequency $\omega _* = T_e (\boldsymbol b \times \boldsymbol {\nabla } \ln N ) \boldsymbol {\cdot } \boldsymbol k /(eB) > 0$. We also define the normalized temperature gradient $\eta = L_N / L_{T}$, a measure of the relative strength between the density and temperature gradients. The FLR effects are taken into account through the zeroth-order Bessel function, ${\rm J}_0$, with argument $k_\perp v_\perp / \varOmega = b \sqrt {x}$ being $b = k_\perp v_{T} /\varOmega$. Collisional effects are introduced in (2.1) through the GK linearized ion collision operator defined by

(2.2)\begin{equation} \mathcal{C} = \left\langle {C} \right\rangle_{\boldsymbol{R}} = \frac{1}{2 {\rm \pi}}\int_0^{2 {\rm \pi}} \,{\rm d} \theta C, \end{equation}

where $C$ is the linearized ion collision operator acting on the ion perturbed particle distribution function, and $\left \langle {\dots } \right \rangle _{\boldsymbol {R}} =\int _0^{2 {\rm \pi}} \,{\rm d} \theta \dots /(2 {\rm \pi})$ denotes the gyroaverage operator evaluated at constant $\boldsymbol {R}$. As the effects of collisions between ions and electrons develop over a time scale that is longer by a factor $O(\sqrt {m_i/ m_e})$ than the ones due to ion–ion collisions, we neglect the former, i.e. we assume $C = C_{ii}$ in (2.2). The linearized GK Boltzmann equation is closed by the GK quasineutrality condition that determines the Fourier component of the electrostatic potential $\phi (\boldsymbol k, t) = \int \,{\rm d} \boldsymbol {r} \phi (\boldsymbol {r},t) {\rm e}^{{\rm i} \boldsymbol k \boldsymbol {\cdot } \boldsymbol {r}}$, that is

(2.3)\begin{equation} \left[ \frac{ {\rm e}^{2}}{T_e} + \frac{ q^{2} }{T} \left( 1 - \varGamma_0(a) \right) \right] \phi = q \frac{\left\langle {\delta n} \right\rangle}{N}, \end{equation}

where $\left \langle {\delta n} \right \rangle = 2 {\rm \pi}\int \,{\rm d} v_\parallel \,{\rm d} \mu B / m {\rm J}_0(b) {g}$ is the perturbed gyroaveraged ion gyrocentre density, $a= b^{2} / 2$ and $\varGamma _0(a) = {\rm I}_0(a) {\rm e}^{- a}$, with ${\rm I}_0(a)$ the modified Bessel function. The first term on the left-hand side of (2.3) is the adiabatic electron response assumed to be given by $n_e \simeq e \phi / T_e$, while the second term constitutes the ion polarization term. We remark that the contribution from the flux-surface-average of $\phi$ in the electron adiabatic response is neglected in (2.3) due to the linear and local limit adopted in the present study (Dorland & Hammett 1993).

We now normalize (2.1) and (2.3). The physical quantities $t$, $v_\parallel$, $k_\perp$, $k_\parallel$, $\phi$ and $q$, are normalized, respectively, to $L_N/c_{s}$, $c_s$, $1/ \rho _s$, $1/ L_N$, $T_{e}/e$ and $e$, with $c_{s} = \sqrt {T_{e}/m}$ the ion sound speed and $\rho _s = c_s / \varOmega$ the ion sound Larmor radius. The ion–ion collision frequency, $\nu = 4 N q^{4} \ln \varLambda / [ 3 \sqrt {{\rm \pi} } \sqrt {m} T^{3/2}]$ that enters in the ion–ion collision operator $\mathcal {C}$, is normalized to $c_s / (N L_N)$ and the gyrocentre density to $N$. The same dimensionless units are used in the rest of the present paper. The normalized GK Boltzmann equation, (2.1), then becomes

(2.4)\begin{gather} \frac{\partial }{\partial t} {g} + {\rm i} ( \omega_{\boldsymbol{\nabla} B} + k_\parallel \sqrt{2 \tau} s_{{\parallel} }) {g} - {\rm i} \omega_{*}^{T} {\rm J}_0(b \sqrt{x}) \phi F_{M} + \frac{q}{\tau}{\rm i} \omega_{\boldsymbol{\nabla} B} {\rm J}_0(b \sqrt{x}) \phi F_{M} \nonumber\\ + \frac{\sqrt{2}}{\sqrt{\tau}} {\rm i} k_\parallel s_\parallel {\rm J}_0(b\sqrt{x}) \phi F_{M} = \mathcal{C}, \end{gather}

where $\omega _{\boldsymbol {\nabla } B} = \tau \omega _{B} ( x + 2 s_\parallel ^{2}) /q$ (with $\omega _{ B} = k_\perp R_B$, $\tau = T_i / T_e$ and $R_B = L_N / L_B$ the normalized magnetic gradient), $\omega _*^{T} = \omega _{*} [1 + \eta (x + s_\parallel ^{2} - 3/2)]$ (with $\omega _* = k_\perp$ the normalized ion diamagnetic frequency). The argument of the Bessel function, $b \sqrt {x}$, is normalized to $b = k_\perp \sqrt {2 \tau }$. The self-consistent electrostatic potential, $\phi$, is obtained by the normalized linearized GK quasineutrality,

(2.5)\begin{equation} \left[ 1 + \frac{q^{2}}{\tau}\left( 1 - \varGamma_0(a) \right) \right] \phi = q \left\langle {\delta n} \right\rangle. \end{equation}

Equations (2.4) and (2.5) constitute a closed system that evolves self-consistently ${g}$ and $\phi$, in the presence of the temperature and magnetic equilibrium gradients, once the collision operator is specified (see § 2.2).

2.1. gyromoment hierarchy equation

To approach the solution of the GK model, given by (2.4) and (2.5), we use a Hermite–Laguerre moment expansion of the ion distribution function. This allows us to reduce the dimensionality of the linearized GK equation, (2.1), by projecting it onto velocity-space Hermite–Laguerre basis polynomials. More precisely, we decompose the perturbed gyrocentre distribution function as (Jorge, Ricci & Loureiro 2017; Mandell, Dorland & Landreman 2018; Jorge et al. 2019a; Frei et al. 2020)

(2.6)\begin{equation} {g} = \sum_{p = 0}^{\infty} \sum_{j = 0}^{\infty} N^{pj} \frac{H_p(s_{{\parallel}}) L_j(x)}{\sqrt{2^{p} p!}} F_{M}, \end{equation}

where we define the Hermite–Laguerre velocity moments of ${g}$, as follows:

(2.7)\begin{equation} N^{pj} = \frac{1}{N} \int \,{\rm d} \mu \,{\rm d} {v_{{\parallel}}} \,{\rm d} \theta\frac{B}{m} {g} \frac{H_p(s_\parallel) L_j(x)}{\sqrt{2^{p} p!}}, \end{equation}

with $N = \int \,{\rm d} \mu \,{\rm d} {v_{\parallel }} \,{\rm d} \theta B F_M / m$. The Hermite–Laguerre coefficients $N^{pj}$ in (2.6) are referred to as the gyromoments. In (2.6), we introduce the Hermite and Laguerre polynomials, $H_p$ and $L_j$, via their Rodrigues’ formulae (Gradshteyn & Ryzhik 2014), that is

(2.8)\begin{equation} H_p(y) = ({-}1)^{p} {\rm e}^{y^{2}} \frac{{\rm d}^{p}}{{{\rm d} y}^{p}} ( {\rm e}^{- y^{2}} ) \end{equation}

and

(2.9)\begin{equation} L_j(y) = \frac{{\rm e}^{y}}{j!} \frac{{\rm d}^{j} }{{{\rm d} y}^{j}} ({\rm e}^{- y } y^{j}) \end{equation}

and their orthogonality relations, that is

(2.10)\begin{equation} \int_{- \infty}^{\infty} {{\rm d} y} H_p(y) H_{p'}(y) {\rm e}^{- y^{2}} = 2^{p} p! \sqrt{\rm \pi} \delta_{p}^{p'} \end{equation}

and

(2.11)\begin{equation} \int_0^{\infty} {{\rm d} y} L_j(y) L_{j'}(y) {\rm e}^{{-}y} = \delta_j^{j'}, \end{equation}

respectively. We now project (2.4) and (2.5) onto the Hermite–Laguerre polynomial basis. For this purpose, we note that the Bessel function ${\rm J}_m$ appearing in both (2.4) and (2.5) and in the GK collision operators (Frei et al. 2021) can be conveniently projected by expanding the Bessel functions, ${\rm J}_m$, onto associated Laguerre polynomials, $L^{m}_n(x) = (-1)^{m} \,{\rm d}^{m} L_{n + m}(x) / {{\rm d} x}^{m}$,

(2.12)\begin{equation} {\rm J}_m(b \sqrt{x}) = \left(\frac{b \sqrt{x}}{2}\right)^{m}\sum_{n=0}^{\infty} \frac{n!\mathcal{K}_{n}}{(n + m)!} L^{m}_n(x), \end{equation}

where we introduce the $n$th-order kernel function

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

with argument $b = \sqrt {2 \tau } k_\perp$ (Frei et al. 2020). The projections of (2.4) and (2.5) onto the Hermite–Laguerre basis yields the linearized gyromoment hierarchy equation that states the time evolution of the gyromoments, $N^{pj}$, that is

(2.14)\begin{gather} \frac{\partial}{\partial t} N^{pj} + {\rm i} k_\parallel \sqrt{\tau} \sum_{p'} \mathcal{H}^{{\parallel} p}_{p'} N^{p'j} + {\rm i} \frac{\tau}{q}\omega_{ B } \sum_{p'j'} ( \mathcal{H}_{B p'j'}^{{\parallel} pj} + \mathcal{H}_{B p'j'}^{{\perp} pj} ) N^{p'j'} \nonumber\\ = (\mathcal{F}^{pj} - \mathcal{F}_{B}^{{\parallel} pj}) \phi + \mathcal{C}^{pj} , \end{gather}

where we introduce the phase-mixing operators, associated with the parallel and perpendicular drifts,

(2.15a)\begin{gather} \mathcal{H}^{{\parallel} p}_{p'} = \sqrt{p +1 } \delta_{p'}^{p+1}+ \sqrt{p} \delta_{p'}^{p-1}, \end{gather}

(2.15b)\begin{gather}\mathcal{H}_{B p'j'}^{{\parallel} pj} = \delta_{j'}^{j} (\sqrt{(p+1)(p+2)} \delta_{p'}^{p+2} + \sqrt{2p+1} \delta_{p'}^{p} + \sqrt{p(p-1)} \delta_{p'}^{p-2}), \end{gather}

(2.15c)\begin{gather}\mathcal{H}_{B p'j'}^{{\perp} pj} = \delta_{p'}^{p} ( (2j+1)\delta_{j'}^{j} - j \delta_{j'}^{j-1} - (j+1) \delta_{j'}^{j+1}), \end{gather}

and the operators that provide the instability drive and are proportional to the equilibrium gradients

(2.16a)\begin{gather} \mathcal{F}^{pj} = {\rm i} \omega_* \left[ \mathcal{K}_{j} \delta_p^{0} + \eta \left( \mathcal{K}_{j} \frac{\delta_p^{2}}{\sqrt{2}} + \delta_p^{0} \left(2 j\mathcal{K}_{j} - j \mathcal{K}_{j-1} - (j+1)\mathcal{K}_{j+1} \right) \right) \right] , \end{gather}

(2.16b)\begin{gather}\mathcal{F}_{B}^{{\parallel} pj} =\frac{\delta_p^{1}}{\sqrt{\tau} } {\rm i} k_\parallel \mathcal{K}_{j} + {\rm i} \omega_{ B} \delta_p^{0} \left[ 2 ( j +1) \mathcal{K}_{j} - j \mathcal{K}_{j-1} - (j+1)\mathcal{K}_{j+1} \right]+ {\rm i} \omega_B \mathcal{K}_{j} \sqrt{2} \delta_p^{2} . \end{gather}

In (2.14), the Hermite–Laguerre projection of the collision operator, $\mathcal {C}$, is defined by

(2.17)\begin{equation} \mathcal{C}^{pj} = \int \,{\rm d} \theta \,{\rm d} \mu \,{\rm d} v_\parallel \frac{B}{m} \frac{H_p(s_\parallel) L_j(x)}{ \sqrt{2^{p} p!}} \mathcal{C}. \end{equation}

In this work, different collision operator models for $\mathcal {C}$ are considered and detailed in § 2.2. The linearized gyromoment hierarchy equation is coupled with the self-consistent GK quasineutrality condition that, projected onto the Hermite–Laguerre basis, is given by

(2.18)\begin{equation} \left[ 1 + \frac{q^{2}}{\tau}\left( 1 - \sum_{n=0}^{\infty} \mathcal{K}_{n}^{2} \right)\right] \phi = q \sum_{n=0}^{\infty} \mathcal{K}_{n} N^{0n}. \end{equation}

Equations (2.14) and (2.18) constitute an infinite set of linear coupled fluid equations for the gyromoments $N^{pj}$, and depend only on the Fourier mode wavevector $\boldsymbol k$ and time $t$. In order to solve numerically (2.14) together with (2.18), we apply a simple closure by truncation, i.e. we solve (2.14) and (2.18) up to $(p,j) \leq (P,J)$ and set $N^{pj} = 0$ for all $(p,j) > (P,J)$ given $0 < P, J < \infty$. When applying the closure by truncation to the gyromoment hierarchy equation, the infinite sums appearing in the GK quasineutrality, (2.18), are consistently truncated up to $J$, that is

(2.19)\begin{equation} \left[ 1 + \frac{q^{2}}{\tau}\left( 1 - \sum_{n=0}^{J} \mathcal{K}_{n}^{2} \right)\right] \phi = q \sum_{n=0}^{J} \mathcal{K}_{n} N^{0n}. \end{equation}

Finally, we remark that a high-collisional closure of the gyromoment hierarchy is derived in § 4.1.

A gyromoment hierarchy equation, similar to the one given in (2.14), is obtained by Mandell et al. (2018). In their formulation, the probabilist's Hermite polynomials, defined by $\hat {H}_p(y) = 2^{- p/2} H_p( y / \sqrt {2})$ and satisfying the orthogonality relation given in (2.10) but weighted by ${\rm e}^{-y^{2} / 2}$, are used, while collisional effects are modelled by the GK Dougherty collision operator model because of the simplicity of its expression on the Hermite–Laguerre basis (see Appendix C). Since the collisional effects are the main focus of the present study, we consider here more realistic collision operator models.

2.2. Gyrokinetic linearized collision operator models

While numerous collision operator models are available in the literature (see, e.g. Hirshman & Sigmar 1976a; Abel et al. 2008; Sugama, Watanabe & Nunami 2009; Li & Ernst 2011; Madsen 2013b; Sugama et al. 2019), the main focus of this work is a linearized GK Coulomb collision operator, and we leverage its formulation derived and implemented in Frei et al. (2021). This advanced linearized GK collision operator is obtained from the nonlinear full-F GK Coulomb operator derived in Jorge et al. (2019a). To highlight the importance of FLR terms in the Coulomb collision operator, the DK limit of the GK Coulomb collision operator is also considered. In addition, collision operators found in the literature and widely used in numerical codes are also considered for comparison. In particular, we consider the GK and DK Sugama collision operators (Sugama et al. 2009), the GK and DK momentum-conserving pitch-angle scattering operators (Helander & Sigmar 2002), the zeroth-order DK HSC collision operator (Hirshman & Sigmar 1976a) and the GK and DK Dougherty collision operators (Dougherty 1964). The gyromoment expansion of the linearized GK/DK Coulomb and the GK/DK Sugama collision operators were reported in Frei et al. (2021) where the linearized DK Coulomb and GK/DK Sugama operators are successfully benchmarked against the GK continuum code GENE (Jenko et al. 2000). The Hermite–Laguerre expansion of the pitch-angle, HSC and Dougherty collision operators are reported in appendices A–C, respectively.

Before turning to the investigation of the impact of collisions on the ITG mode, here we recall some important conservation properties of linearized like-species collision operators. These properties are useful to derive the high-collisional limit of the gyromoment hierarchy (see § 4.1) and to verify the numerical implementation of the collision operators. A model of a linearized collision operator for ion–ion collisions, $C$, is obtained by linearizing a nonlinear collision operator, that we denote by $C^{NL}$. Since collisions act at the particle position $\boldsymbol {r}$ (rather than at the gyrocentre position $\boldsymbol {R}$), the nonlinear collision operator is most often derived as acting on the full particle distribution function and is expressed in the particle phase space $(\boldsymbol {r}, \boldsymbol v)$. While the functional form of $C^{NL}$ depends on the choice of collision model, the linearized collision operator, $C$, can always be expressed as

(2.20)\begin{equation} C = C^{T} + C^{F}, \end{equation}

where we introduce the test component of the linearized collision operator, $C^{T} = C^{NL}(f,f_{M})$, and the field component, $C^{F} = C^{NL}(f_{M},f)$, with $f_M = f_M(\boldsymbol r, \boldsymbol v)$ the particle Maxwellian distribution function and $f = f (\boldsymbol {r},\boldsymbol v)$ a small perturbation. Conservation properties of the collision operator constrain the test and field components. In fact, while particle conservation is satisfied by both components, such that

(2.21)\begin{equation} \int \,{\rm d} \boldsymbol{v} C^{T} = \int \,{\rm d} \boldsymbol{v} C^{F} = 0, \end{equation}

the momentum and energy conservations require that

(2.22)\begin{gather} \int \,{\rm d} \boldsymbol{v} \boldsymbol v C^{T} ={-} \int \,{\rm d} \boldsymbol{v}\boldsymbol v C^{F} , \end{gather}

(2.23)\begin{gather}\int \,{\rm d} \boldsymbol{v} v^{2} C^{T} ={-} \int \,{\rm d} \boldsymbol{v} v^{2} C^{F} , \end{gather}

respectively. From the conservation properties of the linearized collision operator $C$, (2.21)–(2.23), constraints on the coefficients $\mathcal {C}^{pj}$ can be derived in the zero gyroradius limit, when the gyrocentre and particle position correspond. In fact, using (2.17) combined with (2.21)–(2.23) yields the relations

(2.24a–c)\begin{equation} \mathcal{C}^{00} = 0, \quad \mathcal{C}^{10}= 0, \quad \mathcal{C}^{20} = \sqrt{2} \mathcal{C}^{01}, \end{equation}

which express the conservation of gyrocentres, momentum and energy. The constraints in (2.24a–c) are satisfied by all linearized collision operator models considered in the present work in the long wavelength limit and are fulfilled numerically to machine precision.

3. Collisionless limit

In this section, we derive the collisionless limit of the gyromoment hierarchy equation given in (2.14), showing that the ITG collisionless results (see, e.g. Hahm & Tang 1989; Romanelli 1989; Romanelli & Briguglio 1990; Chen et al. 1991) can be retrieved as an asymptotic limit of the gyromoment hierarchy. We first derive the collisionless ITG dispersion relation in § 3.1 by solving directly the GK equation, (2.4). Then, we obtain closed semianalytical expressions of the collisionless linear gyromoment response in § 3.2. The expressions of the gyromoments we obtain are valid at arbitrary order in the normalized magnetic frequency, $\omega _{\boldsymbol {\nabla } B}/\omega$, and are written in terms of linear combinations of definite integrals of velocity-independent hypergeometric functions, which reduce to linear combinations of derivatives of plasma dispersion functions in the slab limit. Finally, in § 3.3 and thanks to the expressions of the gyromoments developed in § 3.2, we show the equivalence between the gyromoment and GK dispersion relations by considering the infinite gyromoment limit, focusing on the slab branch. Ultimately, the results of the present section (in particular the expressions of the collisionless gyromoment response) allow us to investigate the collisionless gyromoments spectrum and to progress in the analysis of the convergence properties of the gyromoment approach, discussed in §§ 7 and 8.

3.1. Collisionless GK dispersion relation

We first derive the ITG dispersion relation by solving directly the GK equation, (2.4). Fourier transforming (2.1) in time with $\mathcal {C} =0$, such that $\partial _t \to - {\rm i} \omega$, and introducing $\omega = \omega _r + {\rm i} \gamma$, with $\omega _r$ and $\gamma$ the real mode frequency and the growth rate, respectively, the normalized perturbed ion gyrocentre distribution function, ${g}/ \phi$, can be expressed as

(3.1)\begin{equation} \frac{{g}}{\phi}= \sum_{l=1}^{3} \hat{{g}}_{l} \end{equation}

with

(3.2a)\begin{gather} \hat{{g}}_{1 } ={-} \frac{q}{\tau} {\rm J}_0(b \sqrt{x}) F_{M}, \end{gather}

(3.2b)\begin{gather}\hat{{g}}_{2 } = \frac{q}{\tau} \frac{ \omega }{ \omega - \omega_{\boldsymbol{\nabla} B}- z_\parallel s_\parallel } {\rm J}_0(b \sqrt{x}) F_{M}, \end{gather}

(3.2c)\begin{gather}\hat{{g}}_{3} ={-}\frac{\omega_{T}^{*}}{\omega - \omega_{\boldsymbol{\nabla} B}- z_\parallel s_\parallel }{\rm J}_0(b \sqrt{x}) F_{M}, \end{gather}

where we introduce $z_\parallel = \sqrt {2 \tau } k_\parallel$. The $l=1$ contribution to the sum in (3.1) is the adiabatic response, while the $l=2$ and $l=3$ terms are associated with the Landau damping and equilibrium gradient drive. We remark the presence of the velocity-dependent magnetic frequency $\omega _B$ in the resonant denominator of the $l=2$ and $l=3$ contributions.

The collisionless GK dispersion relation is deduced from the linearized GK quasineutrality condition given in (2.5),

(3.3)\begin{equation} 1 + \frac{q^{2}}{\tau} \left( 1 - \varGamma_0(a) \right) - q \hat{n} =0, \end{equation}

where $\hat {n} = \sum _{l=1}^{3} \hat {n}_l$ is the density perturbation with $\hat {n}_l = \int \,{\rm d} \boldsymbol {v} {\rm J}_0 \hat {{g}}_l$ and $\hat {{g}}_l$ given in (3.2). By considering an unstable mode, i.e. $\gamma > 0$, the velocity integrals appearing in $\hat {n}_l$ can be performed explicitly by writing the resonant term, $1/(\omega - \omega _{\boldsymbol {\nabla } B} - z_\parallel s_\parallel )$, in integral form (Beer & Hammett 1996)

(3.4)\begin{equation} \frac{1}{ \omega - \omega_{\boldsymbol{\nabla} B} - z_\parallel s_\parallel} ={-} {\rm i} \int_0^{\infty} \,{\rm d} \sigma \exp({ {\rm i} \sigma (\omega - \omega_{\boldsymbol{\nabla} B} - z_\parallel s_\parallel)}). \end{equation}

We remark that (3.4) allows us to evaluate the velocity integrals, appearing in $\hat {n}_l$, exactly and to retain the effects of magnetic drifts in $\hat {N}_{l}^{pj}$ at arbitrary order in $\omega _B / \omega$. Using (3.4), one derives

(3.5a)\begin{align} \hat{n}_{1} & ={-} \frac{q}{\tau} \varGamma_0(a), \end{align}

(3.5b)\begin{align} \hat{n}_{2} & ={-} \frac{ {\rm i}q}{\tau} \omega \int_0^{\infty} \,{\rm d} \sigma {\rm e}^{{\rm i} \sigma \omega} {\rm I}_\perp(\sigma) {\rm I}_\parallel(\sigma), \end{align}

(3.5d)\begin{align} \hat{n}_{3} & = {\rm i} k_\perp \int_0^{\infty} \,{\rm d} \sigma {\rm e}^{{\rm i} \sigma \omega} \left[\vphantom{\frac{3}{2}} {\rm I}_\parallel (\sigma) {\rm I}_\perp (\sigma) \right. \nonumber\\ & \quad \left.+ \eta \left( {\rm I}_{{\parallel} 2} (\sigma) {\rm I}_\perp (\sigma) + {\rm I}_{{\parallel} } (\sigma) {\rm I}_{{\perp} 1} (\sigma) - \frac{3}{2} {\rm I}_\parallel (\sigma) {\rm I}_\perp (\sigma) \right) \right], \end{align}

where we introduce

(3.6a)\begin{align} {\rm I}_\perp(\sigma)& = \frac{1}{1 + {\rm i}\alpha \sigma} {\rm I}_0\left(\frac{a }{1 + {\rm i}\alpha \sigma} \right) \exp({-a /(1 + {\rm i} \alpha\sigma)}), \end{align}

(3.6b)\begin{align} {\rm I}_\parallel(\sigma) & = \frac{1}{\sqrt{1 + 2{\rm i} \alpha \sigma}} \exp({ - z_\parallel^{2} \sigma^{2} /[4 (1 + 2 {\rm i} \alpha \sigma)]}), \end{align}

\begin{align} {\rm I}_{{\perp} 1}(\sigma) & = \frac{\exp({- a / (1 + {\rm i} \alpha \sigma)})}{2 (1 + {\rm i} \alpha \sigma)^{3}} \nonumber\\ & \quad \times \left[ (2 (1 + {\rm i} \alpha \sigma) - 2 a) {\rm I}_0 \left( \frac{a}{ (1 + {\rm i} \alpha \sigma)} \right)+ 2 a {\rm I}_1 \left( \frac{a}{ (1 + {\rm i} \alpha \sigma)} \right) \right], \end{align}

(3.6c)\begin{align} {\rm I}_{{\parallel} 2}(\sigma) & = \frac{(2(1 + 2 {\rm i} \sigma \alpha ) - \sigma^{2} z_\parallel^{2}) }{4 (1 + 2 {\rm i} \sigma \alpha)^{5/2}} \exp({- z_\parallel^{2} \sigma^{2} / ( 4 (1 + 2 {\rm i} \sigma \alpha))}), \end{align}

with $\alpha = q \omega _B / \tau$. Equation (3.3), with the definitions in (3.5) and the velocity integrals in (3.6) constitute the collisionless GK dispersion relation for the unstable modes $\gamma > 0$. We note that the exponential factors ${\rm e}^{{\rm i} \sigma \omega }$, appearing in (3.5), ensure the convergence of the integrals at $\sigma \to \infty$ for $\gamma > 0$. On the other hand, the inclusion of stable modes ($\gamma < 0$) in (3.3) can be obtained by expressing the velocity integrals in $\hat {n}$ using generalized plasma dispersion functions (Gültekin & Gürcan 2019) and efficient numerical algorithms can be used to solve the dispersion relation in the entire complex plane (Xie et al. 2017; Gültekin & Gürcan 2018). However, here we focus on the upper half of the complex plane since one of the advantages of the transformation (3.4) is that it yields one-dimensional integrals that can be easily evaluated numerically.

We finally note that the collisionless GK dispersion relation can be simplified in the sITG case. Neglecting $\omega _B$ in (3.2), the velocity integrals appearing in $\hat {n}$ and given in (3.5), can be expressed in terms of the plasma dispersion function. Then, the GK dispersion relation, (3.3), becomes

(3.7)\begin{align} & 1 + \frac{1}{\tau} + \frac{1}{\tau} \xi Z(\xi) \varGamma_0(a) - \frac{k_\perp}{k_\parallel} \frac{1}{\sqrt{2 \tau}} \left[ \varGamma_0(a) Z(\xi) \right. \nonumber\\ & \qquad \left. + \eta \left( \varGamma_0(a) \xi \left( 1 + \xi Z(\xi)\right) - \frac{\varGamma_0(a)}{2}Z(\xi) + a Z(\xi)\left(\varGamma_1(a) - \varGamma_0(a) \right) \right)\right] =0, \end{align}

where $\xi = \omega \sqrt {2 \tau }/ k_\parallel$ is the normalized mode phase velocity, $Z(\xi ) \!=\! \int _{-\infty }^{\infty } \,{\rm d} s_\parallel {\rm e}^{- s_\parallel ^{2}}/(s_\parallel \!-\! \xi ) /\sqrt {{\rm \pi} }$ is the plasma dispersion function and $\varGamma _1(a) = {\rm I}_1(a) {\rm e}^{- a}$.

3.2. Collisionless linear gyromoment response

We now obtain the collisionless expressions of the gyromoments by projecting (3.1) onto the Hermite–Laguerre basis. This yields

(3.8)\begin{equation} \frac{N^{pj}}{\phi} = \sum_{l=1}^{3} \hat{N}_{l}^{pj}, \end{equation}

where we introduce

(3.9)\begin{equation} \hat{N}_{l}^{pj} = 2 {\rm \pi}\int \,{\rm d} v_\parallel \,{\rm d} \mu \frac{B}{m} \frac{H_p(s_\parallel ) L_j(x) \hat{{g}}_l }{ \sqrt{2^{p} p !}}, \end{equation}

with $\hat {{g}}_l$ defined in (3.2). Performing the velocity integral in (3.9) using (3.4) for an unstable mode with $\gamma >0$ and identities involving hypergeometric functions and Hermite–Laguerre polynomials (described in Appendix D), we obtain

(3.10a)\begin{align} \hat{N}_1^{pj} & ={-} \frac{q}{\tau} \mathcal{K}_{j} \delta_p^{0}, \end{align}

(3.10b)\begin{align} \hat{N}_{2}^{pj}& ={-} \frac{{\rm i} q \omega}{\tau} \sum_{n=0}^{\infty} \frac{ \mathcal{K}_{n} }{\sqrt{2^{p} p!} }\mathcal{I}_{n}^{pj} , \end{align}

(3.10c)\begin{align} \hat{N}_{3}^{pj}& = \frac{{\rm i} k_\perp}{\sqrt{2^{p} p!}}\sum_{n=0}^{\infty}\mathcal{K}_{n} \left\{\vphantom{\frac{3}{2}} \mathcal{I}_n^{pj} + \eta \left[\vphantom{\frac{3}{2}} (2j+1) \mathcal{I}_n^{pj} - j \mathcal{I}_n^{pj-1}- (j+1) \mathcal{I}_n^{pj+1} \right.\right. \nonumber\\ & \quad \left. \left. + \frac{\mathcal{I}_{n}^{p+2j}}{4} + \left(p + \frac{1}{2} \right) \mathcal{I}_n^{pj} + p(p-1) \mathcal{I}_n^{p-2j}- \frac{3}{2}\mathcal{I}_n^{pj}\right] \right\}, \end{align}

where $\mathcal {I}_{n}^{pj}$ is given in (D8) as a finite sum of terms that involve velocity-independent definite integrals of hypergeometric functions with complex argument. The expressions in (3.10) are semianalytical because they still contain a real definite integrals that need to be computed numerically. We remark that generalized plasma dispersion relations can also be used to express the collisionless gyromoments $\hat {N}_{l}^{pj}$, given in (3.10), extending the formulae in Appendix D to the lower half of the complex plane of $\omega$. We use the collisionless expressions of the gyromoments in § 7 to investigate their magnitude in the presence of kinetic effects, such as Landau damping and magnetic gradient drift resonance.

In the case of the sITG branch, the velocity integrals containing the resonant term, $1 /(\omega - z_\parallel s_\parallel )$ (being $\omega _B = 0$), can be performed without using (3.4). In fact, in this case, the collisionless expression of the gyromoments becomes

(3.11)\begin{align} \frac{N^{pj}}{\phi}& ={-}\frac{q}{\tau} \mathcal{K}_{j} \delta_{p}^{0} + \frac{1}{\sqrt{2 \tau}}\frac{k_\perp}{k_\parallel} \frac{1 }{\sqrt{2^{p} p!}} \left( \mathcal{K}_{j} ({-}1)^{p} Z^{(p)} (\xi) \right. \nonumber\\ & \quad \left. + \eta \left[\!\mathcal{K}_{j} \left( \frac{({-}1)^{p+2}}{4} Z ^{(p+2)}(\xi) + (p + 2j) ({-}1)^{p} Z ^{(p)}(\xi) + p(p-1) ({-}1)^{p-2} Z ^{(p-2)}(\xi)\!\right) \right. \right.\nonumber\\ & \quad \left. \left. - \left[ j \mathcal{K}_{j-1} + (j+1) \mathcal{K}_{j+1} \right] ({-}1)^{p} Z ^{(p)}(\xi) \right) \right] - \frac{q}{\tau} \frac{ \xi \mathcal{K}_{j}}{\sqrt{2^{p} p!}} ({-}1)^{p} Z^{(p)}(\xi), \end{align}

where $Z^{(p)} (\xi )= {\rm d}^{p} Z(\xi )/ {\rm d} \xi ^{p}$ is the $p$th-order derivative of the plasma dispersion relation, which can be defined in terms of the Hermite polynomials, $H_p$, by $Z^{(p)} (\xi ) = (-1)^{p} \int _{-\infty }^{\infty } \,{\rm d} s_\parallel H_p(s_\parallel ) {\rm e}^{- s_\parallel ^{2}}/(s_\parallel - \xi )/\sqrt {{\rm \pi} }$.

As an aside, we also note that (3.8) can be related to the kinetic response of the monomial velocity-space moments, $M_{j,k} = v_T^{j + 2k}\int \,{\rm d} \boldsymbol {v} {g}s_\parallel ^{j} x^{k} / 2^{k}$, used in Beer & Hammett (1996) to derive general closure expressions for toroidal gyrofluid equations, in the $k_\perp = 0$ limit, such as

(3.12)\begin{equation} N^{pj} = \sum_{p_1=0}^{\lfloor p/2 \rfloor} \sum_{j_1=0}^{j} \frac{({-}1)^{p_1} 2^{p-2p_1}}{p_1!(p-2p_1)!} \frac{({-}1)^{j_1} j !}{(j - j_1)!(j_1!)^{2}} \frac{2^{j_1}}{v_T^{2j_1 + p - 2p_1}} \frac{M_{p-2p_1,j_1}}{\sqrt{2^{p} p!}}. \end{equation}

Equation (3.12) provides a direct link between the gyromoment hierarchy equation described here and the moment equations derived in Beer & Hammett (1996) and associated closure expressions.

3.3. gyromoment dispersion relation and infinite gyromoment limit

The ITG dispersion relation expressed in terms of the gyromoments $N^{pj}$ (with the closure up to $(p,j) \leq (P,J)$) is obtained from the truncated GK quasineutrality condition (2.19), yielding

(3.13)\begin{equation} 1 + \frac{q^{2}}{\tau}\left( 1 - \sum_{j=0}^{J} \mathcal{K}_{j}^{2} \right) - \sum_{j=0}^{J} \mathcal{K}_{j} \sum_{l=1}^{3 }\hat{N}^{0j}_l = 0, \end{equation}

where $\hat {N}_l^{0j}$ is given by (3.10) with $p=0$.

Focusing on the sITG branch, we now demonstrate that, in the $P \to \infty$ and $J \to \infty$ limits, the ITG gyromoment dispersion relation, (3.13), is equivalent to the GK dispersion relation given in (3.7). In order to solve the gyromoment dispersion relation (3.13), we use the expressions for the collisionless gyromoments, $N^{pj} / \phi$ given in (3.11) in the sITG case. Setting $p =0$ in (3.11) and plugging the resulting expression for $N^{p0} / \phi$ into (3.13) yields

(3.14)\begin{align} & 1 + \frac{1}{\tau} + \frac{1}{\tau} \xi Z(\xi) \hat{\varGamma_0}(a) - \frac{k_\perp}{k_\parallel} \frac{1}{\sqrt{2 \tau}} \left[\vphantom{\frac{\hat{\varGamma_0}(a,J)}{2}} \hat{\varGamma_0}(a,J) Z(\xi) + \eta\left(\vphantom{\frac{\hat{\varGamma_0}(a,J)}{2}}\hat{\varGamma_0}(a,J) \xi \left( 1 + \xi Z(\xi)\right) \right. \right. \nonumber\\ & \quad \left. \left. - \frac{\hat{\varGamma_0}(a,J)}{2}Z(\xi) +Z(\xi) \hat{\varPi}_1(a,J) \right)\right] =0, \end{align}

where we introduce the quantities

(3.15a)\begin{gather} \hat{\varGamma_0}(a,J)= \sum_{j=0}^{J} \mathcal{K}_{j}^{2}, \end{gather}

(3.15b)\begin{gather}\hat{\varPi}_1(a,J) = \sum_{j=0}^{J} \mathcal{K}_{j} \left(2j \mathcal{K}_{j} - j \mathcal{K}_{j-1} -(j+1)\mathcal{K}_{j+1} \right), \end{gather}

which are associated with FLR effects. We remark that the sITG dispersion relation, given in (3.14) obtained using (3.11) and (3.13), agrees with the one found in Mandell et al. (2018). By comparing (3.14) with (3.7), one observes that they are equivalent in the infinite moment limits provided that $\hat {\varGamma _0}(a,\infty ) = \varGamma _0$ and $\hat {\varPi }_1(a,\infty ) = a (\varGamma _1 - \varGamma _0)$. This can be verified explicitly by considering the $J \to \infty$ limit in (3.15) and using the definition of the kernel function $\mathcal {K}_{j}$ in (2.13). In fact, we derive

(3.16)\begin{align} \hat{\varGamma_0}(a,\infty) & = \sum_{j=0}^{\infty} \frac{1}{(j!)^{2}} \left( \frac{b}{2}\right)^{4j} {\rm e}^{- b^{2} /2} \nonumber\\ & = {\rm e}^{- a} {\rm I}_0\left( a\right) = \varGamma_0(a) \end{align}

and

(3.17)\begin{align} \hat{\varPi}_1(a,\infty) & = 2 {\rm e}^{ - b^{2} /2}\sum_{j=0}^{\infty} \left[ \frac{1 }{j!(j+1)!} \left( \frac{b}{2} \right)^{4(j+1)} - \frac{1}{(j!)^{2}} \left( \frac{b}{2}\right)^{4j + 2}\right] \nonumber\\ & = a {\rm e}^{- a} \left( {\rm I}_1(a) - {\rm I}_0(a)\right) = a \left( \varGamma_1(a)- \varGamma_0(a)\right). \end{align}

Figure 1 displays the quantities, $\hat {\varGamma _0}(a,J)$ and $\hat {\varPi }_1(a,J)$ defined in (3.15), for different values of $J$, as well as their $J \to \infty$ limits. It can be observed that $\hat {\varGamma _0}(a,J)$ and $\hat {\varPi }_1(a,J)$ converge to $\varGamma _0$ and $a(\varGamma _1 - \varGamma _0)$ when $J\to \infty$. We notice that convergence is faster for low values of $a = b^{2}/2$, i.e. at long wavelength. This is in agreement with the fact that the truncation of the sums in (3.16) and (3.17) to include only $j \leq J$ terms yields an error $O( b^{4J})$ (Mandell et al. 2018).

Figure 1. Plots of $\hat {\varGamma _0}(a,J)$ (a) and $\hat {\varPi }_1(a,J)$ (b) as a function of $a = b^{2}/2$ for increasing values of $J$. The solid black lines represent their asymptotic limits, $\varGamma _0$ and $- a (\varGamma _1 - \varGamma _0)$, respectively.

The impact on the ITG mode due to approximating $\hat {\varGamma _0}(a,\infty )$ and $\hat {\varPi }_1(a,\infty )$ with $\hat {\varGamma _0}(a,J)$ and $\hat {\varPi }_1(a,J)$ can be illustrated by solving the ITG dispersion relation, (3.14), for the mode complex frequency $\omega = \omega _r + {\rm i} \gamma$, as a function of $J$, and comparing its prediction with the one from the GK dispersion relation, given in (3.7). Figure 2 shows the ITG growth rate $\gamma$ and the corresponding real frequency $\omega _r$ solution of the gyromoment sITG dispersion relation given in (3.14), plotted as a function of the perpendicular wavenumber $k_\perp$ for increasing values of $J$. For comparison, the collisionless GK solution, obtained from (3.7), is shown. Both $\gamma$ and $\omega _r$, obtained by the dispersion relation expressed in terms of gyromoments, (3.14), approach the one predicted by (3.7) as $J \to \infty$. Convergence is achieved with small $J$ at long wavelengths (e.g. the growth rate peak at $k_\perp \simeq 0.6$ is well reproduced with $J \sim 5$). On the other hand, for small perpendicular wavelengths and as a rule of thumb, $J \gtrsim k_\perp ^{2}$ must be retained in order to retrieve accurately the GK collisionless prediction (e.g. the $k_\perp \simeq 5$ case is well reproduced with $J \gtrsim 25$). We remark that, if $J$ is not sufficiently large, the truncation of the quantities, used in (3.15), can yield significant deviations in the growth rate $\gamma$ and frequency $\omega _r$. The presence of collisions is expected to affect the convergence estimate reducing $J$ even for short wavelength modes, as numerically demonstrated in § 8.

Figure 2. Collisionless sITG growth rate, $\gamma$ (a), and frequency, $\omega _r$ (b), as a function of the perpendicular wavenumber $k_\perp$ for different values of $J$. The coloured lines are the solution of the gyromoment hierarchy dispersion relation, given in (3.14), while the black lines are solutions of (3.7), i.e. the $J = \infty$ solution. Here, $\eta =3$ and $k_\parallel = 0.1$.

4. High-collisional limit

We now investigate the high-collisional limit of the gyromoment hierarchy, (2.14), focusing on two reduced models that retain only a finite number of gyromoments. In the high-collisionality limit, high-order gyromoments are damped by the presence of collisions and the ITG mode can be accurately described by considering only the evolution of the lowest-order gyromoments. Indeed, the perturbed gyrocentre distribution function is expected not to significantly deviate from a perturbed Maxwellian at high collisionalities. We consider, first, a reduced model based on six gyromoments (6GM), obtained by using a closure by truncation of the gyromoment hierarchy. Second, we consider a reduced model based on four gyromoments (4GM), rigorously derived as an asymptotic limit of the 6GM using the Chapman–Enskog closure to express the highest-order gyromoments as a function of the lower-order ones in the limit of $\epsilon \sim k_\parallel \lambda _{{\rm mfp}} \ll 1$. While the Chapman–Enskog procedure allows us to reduce the number of gyromoments, closures for the FLR terms associated with the $\mathcal {K}_{n}$ kernel functions and inherent to gyrofluid models (Brizard 1992; Dorland & Hammett 1993; Beer & Hammett 1996; Waltz et al. 1997; Madsen 2013a; Held, Wiesenberger & Kendl 2020), are still an open issue in the literature. The closure issue appears explicitly in the Fourier exponential form of the kernel function defined in (2.13). Appendix E is dedicated to exploring possible FLR closures by using a Padé approximation technique applied to the 4GM model. To make analytical progress, we focus here on long wavelength modes. Consistently, we model collisions using the DK Coulomb operator (Frei et al. 2021). This has the advantage of avoiding infinite sums with FLR terms that need to be truncated depending on the value of the perpendicular wavenumber, $k_\perp$.

The present section is structured as follows. In § 4.1, we derive the 6GM and 4GM reduced models. The predictions are compared with the ones of the full gyromoment hierarchy and are also verified against the GK code GENE. In § 4.2, we explore the high-collisional limit by deriving a general algebraic collisional dispersion relation from which we discuss the limits of the slab and toroidal branches of the ITG mode.

4.1. 6GM and 4GM reduced models

We start by deriving the 6GM model. We introduce the lowest-order normalized fluid gyromoments, i.e. the gyrocentre density $N = N^{00}$, the parallel velocity $V = N^{10}/\sqrt {2}$, the parallel and perpendicular temperatures, $T_\parallel = \sqrt {2} N^{20} + N$ and $T_\perp = N - N^{01}$ and, finally, the parallel and perpendicular heat fluxes, associated with the deviations of ${g}$ from the Maxwellian distribution function, $Q_\parallel = N^{30}$ and $Q_\perp = N^{11}$. The evolution equations of these gyromoments are obtained by setting $(p,j) = (0,0)$, $(1,0)$, $(2,0)$, $(0,1)$, $(3,0)$ and $(1,1)$ in (2.14), respectively, and by neglecting all higher-order gyromoments. This yields

(4.1a)\begin{align} & \frac{\partial}{\partial t } N + {\rm i} k_\parallel \sqrt{2 \tau} V + {\rm i} \tau \omega_{B} \left( T_\parallel{-} T_\perp \right) + {\rm i} \varGamma_1 \phi = 0 \end{align}

(4.1b)\begin{align} & \sqrt{2} \frac{\partial}{\partial t }V + \frac{{\rm i} k_\parallel}{\sqrt{\tau}} \mathcal{K}_{0} \phi + {\rm i} k_\parallel \sqrt{\tau} T_\parallel{+} {\rm i} \tau \omega_B \left( 4 \sqrt{2} V- Q_\perp{+} \sqrt{6} Q_\parallel \right)= 0,\end{align}

(4.1c)\begin{align} & \frac{\partial}{\partial t } T_\parallel{+} {\rm i} k_\parallel \sqrt{6\tau} Q_\parallel{+} {\rm i} \tau \omega_B \left( 5 T_\parallel{-} 4 N + T_\perp \right) \nonumber\\ & \qquad + {\rm i} (\varGamma_2 - \varGamma_1) \phi ={-} \nu \alpha \left( T_\parallel{-} T_\perp \right),\end{align}

(4.1d)\begin{align} & \frac{\partial}{\partial t }T_\perp{+} {\rm i} k_\parallel \sqrt{\tau} ( \sqrt{2 } V+Q_\perp ) - {\rm i} \tau \omega_B \left( 3 N - T_\parallel{-}3 T_\perp \right) \nonumber\\ & \qquad - {\rm i} (\varGamma_{3}- \varGamma_1) \phi = \nu \frac{\alpha}{2}\left( T_\parallel{-} T_\perp\right), \end{align}

(4.1e)\begin{align} & \frac{\partial}{\partial t } Q_\parallel{+} {\rm i} k_\parallel \sqrt{\frac{3 \tau}{2}} (T_\parallel{-} N) + {\rm i} \tau \omega_B \left( 2 \sqrt{3} V + 8Q_\parallel\right)={-} \nu \left( \frac{8}{5}\sqrt{\frac{2}{\rm \pi}} Q_\parallel{+} \frac{8}{5 \sqrt{3 {\rm \pi}}} Q_\perp \right), \end{align}

(4.1f)\begin{align} & \frac{\partial}{\partial t }Q_\perp{+} \frac{{\rm i} k_\parallel }{\sqrt{\tau}} \mathcal{K}_{1} \phi + {\rm i} k_\parallel \sqrt{\tau} ( N - T_\perp)\nonumber\\ & \qquad - {\rm i} \tau \omega_B \left( \sqrt{2} V - 6 Q_\perp \right)={-} \nu \left( \frac{28}{15}\sqrt{\frac{2}{\rm \pi}} Q_\perp{+} \frac{8}{5 \sqrt{3 {\rm \pi}}} Q_\parallel \right), \end{align}

where we introduce the numerical coefficients $\alpha = 16/15 \sqrt {2 / {\rm \pi}}$. We remark that the collisional terms in (4.1) are obtained by using the lowest-order gyromoments of the DK Coulomb collision operator (Frei et al. 2021) and the energy constraints expressed in (2.24a–c). The quantities appearing in the evolution equations, (4.1), are defined by

(4.2a)\begin{gather} \varGamma_1 = \mathcal{K}_{0} \left( 2 \omega_B - \omega_* \right) + \mathcal{K}_{1} \left( \omega_* \eta - \omega_B\right) , \end{gather}

(4.2b)\begin{gather}\varGamma_2 = \left( 2 \omega_B- \omega_* \eta \right) \mathcal{K}_{0}, \end{gather}

(4.2c)\begin{gather}\varGamma_3 = \mathcal{K}_{0} \left( \omega_* \eta - \omega_B \right) + \mathcal{K}_{1} \left( 4 \omega_B- \omega_* ( 1 + 2 \eta)\right) + 2 \mathcal{K}_{2} \left( \eta \omega_* - \omega_B \right). \end{gather}

Additionally, the electrostatic potential, $\phi$, is determined by the self-consistent quasineutrality condition obtained from (2.18), i.e.

(4.3)\begin{equation} \left[ 1 + \frac{q^{2}}{\tau} \left( 1 - \sum_{j=0}^{1} \mathcal{K}_{j}^{2} \right) \right]\phi = q(\mathcal{K}_{0} + \mathcal{K}_{1} ) N - \mathcal{K}_{1}T_\perp, \end{equation}

where the infinite sum appearing in the polarization term in (2.18) is truncated to $j = 1$ consistently with the fact that the FLR terms are represented up to $\mathcal {K}_{1}$ in the right-hand side of the same equation. The system in (4.1) with (4.3) constitutes a closed set of fluid-like equations that we refer to as the 6GM model.

To derive the 4GM model, further reducing the number of gyromoments appearing in the 6GM model, we apply the Chapman–Enskog asymptotic closure scheme (Chapman & Cowling 1941; Jorge et al. 2017) to the 6GM model. We introduce the dimensionless small parameter $\epsilon \sim \lambda _{{\rm mfp}} k_\parallel \sim \sqrt {\tau } k_\parallel / \nu \ll 1$, which quantifies the importance of the non-Maxwellian part of ${g}$. Thus, we expect that the Chapman–Enskog closure and the reduced models derived in this section to be valid when typically $\nu \gg \sqrt {\tau } k_\parallel$. The small parameter $\epsilon$ allows us to express the heat fluxes $Q_\parallel$ and $Q_\perp$ as a function of the lowest-order gyromoments, by applying the ordering $Q_\parallel \sim Q_\perp \sim \epsilon N$ (with $T_\parallel \sim T_\perp \sim V \sim N$) and $\partial _t \sim \gamma \sim \epsilon \nu$. We neglect the terms proportional to $\omega _B$, which are smaller than the terms proportional to $k_\parallel T_\parallel$ and $k_\parallel T_\perp$ since $\omega _B / k_\parallel \lesssim 1$ in the boundary region. With these assumptions, (4.1c) and (4.1d) can be solved at the leading order in $\epsilon$, yielding

(4.4a)\begin{gather} Q_{{\parallel} } = {\rm i} k_\parallel \frac{\chi_{{\perp}}^{{\parallel}} }{\tau} \mathcal{K}_{1} \phi + {\rm i} k_\parallel \chi_{{\perp}}^{{\parallel}} ( N - T_\perp) - {\rm i} k_\parallel \chi_{{\parallel}}^{{\parallel}} (T_\parallel{-} N ), \end{gather}

(4.4b)\begin{gather}Q_\perp{=}{-} {\rm i} k_\parallel \frac{\chi_\perp^{{\perp}} }{\tau} \mathcal{K}_{1} \phi - {\rm i} k_\parallel \chi_\perp^{{\perp}} ( N - T_\perp) + {\rm i} k_\parallel \chi^{{\perp}}_{{\parallel}} (T_\parallel{-} N ), \end{gather}

where we introduce the thermal conductivities,

(4.5a–d)\begin{align} \chi_{{\perp}}^{{\parallel}} = \frac{5}{16}\sqrt{\frac{\rm \pi}{3}} \frac{\sqrt{\tau}}{\nu}, \quad \chi_{{\parallel}}^{{\parallel}} = \frac{35}{32} \sqrt{\frac{2 {\rm \pi}}{3}}\frac{\sqrt{\tau}}{\nu}, \quad \chi_\perp^{{\perp}} = \frac{5\sqrt{2{\rm \pi}}} {16}\frac{\sqrt{\tau}}{\nu}, \quad \chi_\parallel^{{\perp}} = \frac{5\sqrt{\rm \pi}}{16} \frac{\sqrt{\tau}}{\nu}. \end{align}

We remark that $\chi _{\parallel }^{\parallel } > \chi _{\perp }^{\parallel }$ and that $\chi _\perp ^{\perp } > \chi _\parallel ^{\perp }$. Equation (4.1), with the closure relations expressed in (4.4), yields the 4GM for the lowest-order gyromoments fluid quantities $N$, $V$, $T_\parallel$ and $T_\perp$ and provides the rigorous $\nu \gg \sqrt {\tau } k_\parallel$ asymptotic limit of the gyromoment hierarchy. We remark that an analytical procedure to close the gyromoment hierarchy similar to the one considered here for the ITG mode, can be applied in the more general cases that include variation of equilibrium quantities along the magnetic field line and electron dynamics.

To study the limits of validity of the 6GM and 4GM models, we compare in figure 3 the estimate of the ITG growth rate $\gamma$ and real frequency $\omega _r$ with the results of the gyromoment hierarchy using the DK Coulomb collision operator (Frei et al. 2021) and the GK collisionless dispersion relation given in (3.3). To compute $\gamma$ and $\omega _r$ from the gyromoment hierarchy equation (see (2.14)) and from the 6GM and 4GM models (see (4.1) and (4.4)), an initial-value calculation as well as a direct eigenvalue solver are used, and their results are compared for consistency. The same numerical methods are used to evaluate all the results reported in the present paper. For verification purposes, the results of GENE simulations using the same operator (Jenko et al. 2000) are also shown. Here, we consider $(P,J) = (18,6)$ with a closure by truncation (convergence studies are reported in § 8) and $k_\perp = 0.5$, which corresponds approximately to the fastest growing ITG mode. The results are shown as a function of the collisionality $\nu$, for two different normalized temperature gradients, $\eta =3$ and $\eta =5$. First, an excellent agreement is observed between the gyromoment approach and the GENE simulations at arbitrary collisionality $\nu$, using the same DK Coulomb collision operator, confirming the preliminary study reported in Frei et al. (2021). Second, both GENE and the gyromoment ITG growth rates and real frequencies agree well with the collisionless limit (i.e. $\nu \ll 1$) and with the high-collisional predictions (i.e. $\nu \gg \sqrt {\tau } k_\parallel \sim 0.1$), obtained from the 6GM and 4GM models. Overall, the 6GM provides a better approximation to the low-collisionality solution than the 4GM, while they both agree in the high-collisional limits. It is remarkable that the 6GM retrieves surprisingly well the ITG growth rate $\gamma$ when $\eta =5$, but fails to predict correctly the mode frequency $\omega _r$. Finally, we remark that the 4GM yields an unphysical $\omega _r < 0$, when the collisionality is below $\nu \lesssim 10^{-2}$ with $\eta =3$. This is mainly due to the fact that, in the low-collisionality limit, the $1/\nu$ dependence of the thermal conductivities given in (4.5a–d) produces arbitrarily large heat fluxes. The good agreement observed in figure 3 illustrates the multi-fidelity nature of the gyromoment approach that allows the derivation of reduced models able to address both the collisionless and high-collisional limits of the GK equation, an attractive feature in the assessment of the plasma dynamics in the boundary region.

Figure 3. The sITG growth rate $\gamma$ (a) and real frequency $\omega _r$ (b) as a function of collisionality $\nu$ obtained using the gyromoment approach with the DK Coulomb collision operator (solid blue line), GENE with the same operator (red markers), the collisionless GK dispersion relation (dotted black) and the 4GM and 6GM high-collisional limits (red dashed–dotted and red dotted lines, respectively). The solution of the dispersion relation given by (4.6), derived in the case of $Q_{\parallel, \perp } = 0$, is also plotted for comparison by the solid thin black lines. Here, we consider $k_\parallel = 0.1$, $k_\perp = 0.5$ and $\tau = 1$.

4.2. High-collisional dispersion relation

We now derive an algebraic dispersion relation using the 4GM equations and discuss the limits of the slab and toroidal branch of the ITG mode. As collisional effects enter through the parallel and perpendicular thermal conductivities such that $\chi ^{\parallel }_{\perp } \sim \chi ^{\parallel }_{\parallel } \sim \chi ^{\perp }_{\parallel } \sim \chi ^{\perp }_{\perp } \sim 1/\nu$ (see (4.5a–d)), we assume that the heat fluxes become negligible when $\nu \gg k_\parallel$, i.e. $Q_{\parallel } \simeq 0$ and $Q_\perp \simeq 0$ (see (4.4)). Hence, a dispersion relation can be derived from the 4GM equations, (4.1), yielding

(4.6)\begin{equation} a \omega^{4} + b \omega^{3} + c \omega^{2} + d \omega + e = 0, \end{equation}

where

(4.7a)\begin{align} a & = \varGamma_\phi , \end{align}

(4.7b)\begin{align} b & = \frac{3 }{2} {\rm i} \alpha \nu \varGamma_\phi - \varGamma_1 \mathcal{K}_{0} - \varGamma_3 \mathcal{K}_{1}- 16 \varGamma_\phi \tau \omega_{B}, \end{align}

(4.7c)\begin{align} c & = \frac{ {\rm i} \alpha \nu}{2} \left( \mathcal{K}_{1}\left( \varGamma_2 - 2 \varGamma_3\right) - 32 \tau \varGamma_\phi \omega_B \varGamma_1 \mathcal{K}_{0} \right) \nonumber\\ & \quad - 3 k_\parallel^{2} \tau \varGamma_\phi - k_\parallel^{2} \mathcal{K}_{0}^{2} + \tau \omega_B \left( \mathcal{K}_{0} \left( 14 \varGamma_1 - \varGamma_2 + \varGamma_3 \right) + \mathcal{K}_{1} \left(\varGamma_1 + 12 \varGamma_3 \right) \right), \end{align}

(4.7d)\begin{align} d& = {\rm i} \alpha \nu \left[ k_\parallel^{2} \left( \mathcal{K}_{0} \mathcal{K}_{1} - \frac{5}{2} \tau \varGamma_\phi - \frac{3}{2} \mathcal{K}_{0}^{2} \right)\right.\nonumber\\ & \quad \left.+ \tau \omega_B \left( \varGamma_1 \left(13 \mathcal{K}_{0} + 2 \mathcal{K}_{1} \right) - (\varGamma_2 - 2 \varGamma_3)(\mathcal{K}_{0} + 3 \mathcal{K}_{1}) \right)\vphantom{\left( \mathcal{K}_{0} \mathcal{K}_{1} - \frac{5}{2} \tau \varGamma_\phi - \frac{3}{2} \mathcal{K}_{0}^{2} \right)} \right] \nonumber\\ & \quad + \tau k_\parallel^{2} \left[ \omega _{B } \left(18 \tau \varGamma _{\phi }+\mathcal{K}_{0} \left(8 \mathcal{K}_{0} +\mathcal{K}_{1} \right)\right)+ \left(2 \varGamma _1-\varGamma _2\right) \mathcal{K}_{0} +3 \varGamma _3 \mathcal{K}_{1} \right], \end{align}

(4.7e)\begin{align} e & = {\rm i} \alpha \nu \left[ 5 k_\parallel^{2} \tau \omega_B \left( \tau \varGamma _{\phi } + \mathcal{K}_{0}^{2} \right) +k_\parallel^{2} \tau \left( \varGamma _1- \frac{\varGamma _2}{2}+ \varGamma _3\right) \left(\mathcal{K}_{0}+ \mathcal{K}_{1}\right] \right) \nonumber \\ &+ k_\parallel^{2} \tau^{2} \omega_B \left(\mathcal{K}_{1} \left(\varGamma _2-6 \varGamma _3 -2 \varGamma _1 \right) -2 \mathcal{K}_{0} \left(4 \varGamma _1-2 \varGamma _2+\varGamma _3\right) \right), \end{align}

being terms proportional to $\omega _B^{2}$ neglected ($\omega _B \ll 1$) and having defined $\varGamma _\phi = 1 + q^{2} (1 - \mathcal {K}_{0}^{2} - \mathcal {K}_{1}^{2}) / \tau$. In (4.6), collisional effects are represented through the terms proportional to ${\rm i}\alpha \nu$. The solution of (4.6) is plotted in figure 3 and approaches the 6GM and 4GM models when $\nu \gg k_\parallel$, such that $Q_\parallel = Q_\perp \simeq 0$.

We now consider the slab and toroidal limits of (4.6) separately in detail. An estimate of the sITG peak can be derived from (4.6) as a function of $k_\parallel$, $\nu$ and $\eta$. Besides neglecting the terms proportional to $\omega _B$ and FLR effects (i.e. imposing $\mathcal {K}_{0} \simeq 1$, $\mathcal {K}_{1} \simeq 0$, $\mathcal {K}_{2} \simeq 0$, such that $\varGamma _\phi \simeq 1$), we assume that the sITG mode is far from the marginal stability limit $\eta \sim \nu \gg 1$. Then, the dispersion relation in (4.6) reduces to

(4.8)\begin{align} \hat{\omega}_*^{4} + \frac{3}{2} {\rm i} \alpha \hat{\nu}_* \hat{\omega}_*^{3} - \hat{k}_{{\parallel} *}^{2} \left( 1 + 3 \tau\right) \hat{\omega}_*^{2} + \hat{k}_{{\parallel} *}^{2} \left[ \tau - \frac{3}{2} {\rm i} \alpha \hat{\nu}_* \left( 1 + \frac{5}{2} \tau\right)\right] \hat{\omega}_* + \frac{3}{2} {\rm i} \alpha \hat{k}_{{\parallel} *} \hat{\nu}_* \tau =0, \end{align}

where we introduce $\hat {\omega }_* = \omega / (\eta \omega _*)$, $\hat {\nu }_* = \nu / (\eta \omega _*)$ and $\hat {k}_{\parallel *} = \hat {k}_\parallel / (\eta \omega _*)$. Given the normalization of (4.8), the maximum growth rate is expected to be proportional to $\gamma \simeq g(\tau, \hat {\nu }_*) \eta \omega _*$ and occurs at $k_\parallel \simeq f(\tau, \hat {\nu }_*) \eta \omega _*$. In figure 4, we show the maximum ITG growth rate, $\hat {\gamma }_* = \gamma / (\eta \omega _*)$ (i.e. the imaginary part of $\hat {\omega }_*$) obtained by solving (4.8) as a function of $\hat {k}_{\parallel *}$ and $\hat {\nu }_*$. The effect of increasing $\hat {\nu }_*$ is primarily to shift the growth rate peak to larger values of $\hat {k}_\parallel ^{*}$ for all values of $\tau$. Additionally, while $g(\tau,\nu )$ increases with $\tau$ for all $\nu$, it is found that $f(\tau, \hat {\nu }_*)$ is an increasing function of $\tau$ at large values of $\hat {\nu }_*$.

Figure 4. Estimates of the normalized sITG peak, $\hat {\gamma }_*$, obtained from (4.8), for increasing values of $\hat {\nu }_* = 0.015$ (from (a) to (b)), as a function of $\hat {k}_{\parallel *}$ and the temperature ratio $\tau$. The local maxima of the growth rate is indicated by the dashed black lines. The growth rates $\hat {\gamma }_*$ are normalized to their maximal values.

We remark that the collisionality dependence of the functions $f$ and $g$ is not present in previous fluid ITG models based on the drift-reduced fluid equations (Mosetto et al. 2015). This dependence arises here because the 4GM model (and more generally the gyromoment hierarchy equation in (2.14)) allows for a different evolution of the perpendicular and parallel temperature fluctuations assuming, in general, $T_\parallel \neq T_\perp$. On the other hand, the drift-reduced Braginskii fluid model in, e.g. Zeiler et al. (1997), assumes that $T_\parallel = T_\perp$. In fact, the $\nu$ dependence in (4.8) enters in the dispersion relation, given in (4.8), through terms proportional to $\nu \alpha (T_\parallel - T_\perp )$ which arise from the gyromoment $\mathcal {C}^{20}$ and $\mathcal {C}^{01}$ components of the DK Coulomb collision operator.

To investigate the effect of anistotropic temperature fluctuations, we derive the evolution equation of the temperature, $T = (T_\parallel + 2 T_\perp ) /3 - N$. Using (4.1c) and (4.1d) and the continuity equation, (4.1a), this can be expressed as

(4.9)\begin{equation} \frac{\partial }{\partial t} T - {\rm i} k_\parallel \frac{ \sqrt{2 \tau}}{3}V + \frac{{\rm i}}{3} \left( \varGamma_2 - 2 \varGamma_3 - 2 \varGamma_1\right) \phi =0. \end{equation}

From (4.9) with (4.1a) and (4.1b), and considering the limits used to derive (4.8), we obtain the estimate of the maximum sITG growth rate,

(4.10)\begin{equation} \hat{\omega}_*^{3} - \hat{k}_{{\parallel} *}^{2} \left( 1 + \frac{2}{3}\tau \right) \hat{\omega}_* + \hat{k}_{{\parallel}*}^{2} \tau = 0. \end{equation}

From (4.10), one finds that the peak of the sITG growth rate is proportional to $\gamma \simeq g(\tau ) \eta \omega _*$ and occurs at $k_{\parallel } \simeq f(\tau ) \eta \omega _*$. It is found that $g(\tau )$ is an increasing function of $\tau$, while $f(\tau )$ decreases with $\tau$. We remark that (4.10) agrees with Mosetto et al. (2015) by replacing the numerical factor $2/3$ by $5/3$, a difference due to the additional degree of freedom assumed in the 4GM model.

We now turn to the toroidal ITG mode, driven unstable by the magnetic drifts and by the presence of a finite temperature gradient. We derive a simple dispersion relation for the toroidal ITG from (4.6) indicative of its properties. Since, in contrast to the sITG, the toroidal ITG subsists at $k_\parallel = 0$, we set $k_\parallel =0$ thus removing the coupling with sound waves, and we adopt the long perpendicular wavelength limit, keeping terms up to $O(k_\perp ^{2})$, such that $\mathcal {K}_{0} \simeq 1 - \tau k_\perp ^{2}/2$, $\mathcal {K}_{1} \simeq \tau k_\perp ^{2}/2$, and $\mathcal {K}_{2} \simeq 0$. Since $\omega _B \ll \omega _*$ in the boundary region, we obtain

(4.11)\begin{equation} \omega^{3} + b_t \omega^{2} + c_t \omega + d_t = 0, \end{equation}

where

(4.12a)\begin{align} b_t & = \frac{3}{2} {\rm i} \alpha \nu + \omega_* - 2 \omega_B \left( 1 + 8 \tau \right) + k_\perp^{2} \left[ 2 \omega_B \left( 1 + \frac{3}{2} \tau \right) - \omega_* \left( 1 + \tau + \eta \tau \right) \right], \end{align}

(4.12b)\begin{align} c_t & = \frac{3}{2 } {\rm i} \alpha \nu \left[ \omega_* - 2 \omega_B - \frac{32}{3} \tau \omega_B + k_\perp^{2} \left( 2 \omega_B \left( 1 + \frac{19}{12} \tau \right) -\omega_* (1 + \tau + \tau \eta) \right) \right] \nonumber\\ & \quad + \omega_B \omega_* \left[ 2 \tau ( \eta -7) + k_\perp^{2} \left(10 \eta \tau \left(\tau - \frac{1}{5}\right) +\tau\left(13 \tau +14\right) \right) \right], \end{align}

(4.12c)\begin{align} d_t & ={\rm i} \alpha \nu \omega _{B } \omega _* \tau \left[ \left(3 \eta -13 \right) + k_\perp^{2} \left(3 \eta \left(2 \tau-1 \right) +11 \tau +13 \right)\right]. \end{align}

An estimate of the value of $k_\perp$ that yields the largest growing mode can be derived from (4.11) by neglecting terms proportional to $\nu$, for simplicity. Then, the dispersion relation in (4.11) becomes a second-order equation for $\omega$ presenting the largest solution $\omega$ when the coefficient $b_t$ is minimized. Hence, by minimizing $b_t$ and considering $\omega _B \ll \omega _*$, the growth rate is maximized when $k_\perp ^{2} \simeq 1 / [ 1 + \tau (1 + \eta )]$. Figure 5 shows the toroidal ITG growth rate obtained by solving (4.11), as a function of the perpendicular wavenumber $k_\perp$ and temperature ratio $\tau$. We observe that the decrease of the perpendicular wavenumber of the fastest growing mode with $\tau$ is present also at finite collisionality, $\nu$. We also remark that collisions tend to destabilize long wavelength modes.

Figure 5. Toroidal ITG growth rate as a function of the perpendicular wavenumber, $k_\perp$ and the temperature ratio $\tau$ obtained from (4.11), for a low (a) and a high (b) value of collisionality. The local maxima of the growth rate are indicated by the dashed black lines. The growth rates are normalized to the maximum value. Here, $\eta = 7$ and $R_B = 0.1$.

5. ITG mode with GK coulomb collision operator

We now perform parameter scans to explore the dependence of the ITG linear growth rate on the perpendicular and parallel wavenumbers, on the level of collisionality and on the temperature and magnetic gradient strengths. Collisional effects are modelled with the GK Coulomb collision operator (Frei et al. 2021). As a reference value, we consider $k_\parallel \sim 0.1$. This is justified by the strong ballooning character of the ITG mode in the boundary region where $k_\parallel$ is of the order of $L_N / q_0 R_0$ (with $q_0$ the safety factor and $R_0$ the major radius of the tokamak devices) and typical H-mode edge values of $q_0 \sim 4$ and $L_N / R_0 \lesssim 1/50$. In addition, we note that $1 \lesssim \tau \lesssim 4$ in the boundary. Therefore, we consider a temperature ratio of $\tau = 1$ for simplicity, and note that the effect of increasing $\tau$ can be inferred from figures 4 and 5 for the slab and toroidal cases, respectively. We remark, that given the value of $k_\parallel$ and $\tau$ used in this section, the high-collisional regime is reached when $\nu \gg \sqrt {\tau } k_\parallel =0.1$ in agreement with § 4.1. For the present numerical calculations, we fix $(P,J) = (18,6)$. As we show in § 8, this choice ensures good convergence of the numerical results in the range of parameters explored in this section.

Figure 6 displays the ITG growth rate in the $(k_\parallel,k_\perp )$ parameter space when $\eta =3$ for increasing collisionality $\nu$ (from figure 6a to 6c and from figure 6d to 6f) and magnetic gradient $R_B$ (slab in figure 6a–c and toroidal in figure 6a–c). For comparison purposes, we show the stability boundaries obtained numerically from the gyromoment hierarchy (dotted white) and from the 6GM and 4GM reduced models (red dotted and dashed–dotted). First, we observe that, as collisionality increases, FLR collisional stabilization effects become important and eventually suppress the ITG modes when $k_\perp \gtrsim 1$ and $\nu \gtrsim 0.5$. In fact, we notice that, at intermediate and high collisionalities, the gyromoment hierarchy stability boundary follows closely the ones predicted by the 6GM and 4GM models for long parallel ($k_\parallel \lesssim 0.15$) and perpendicular ($k_\perp \lesssim 0.2$) wavelengths where the FLR collisional effects are negligible and kinetic effects are small. However, while the ITG mode is completely suppressed for $k_\perp \gtrsim 0.5$ when $\nu > 0.5$ according to the full gyromoment calculation, the 6GM and 4GM stability boundaries extend to $k_\perp \gtrsim 1$. This highlights the importance of FLR terms in the Coulomb collision operator model. In fact, despite the lack of collisionless FLR terms in the 6GM and 4GM models, these models agree with the gyromoment hierarchy at high collisionality when collisional effects are modelled by the DK Coulomb even when $k_\perp \sim 2$ (see, e.g. figure 12). We remark that the complete stabilization of the ITG mode by FLR collisional damping requires a level of collisionality (i.e. $\nu \gtrsim 0.5$) that is even higher than the typical SOL experimental values, i.e. $T \gtrsim 10$ eV, $N \sim 10^{19}$ m$^{-3}$ and $L_N \sim 5$ cm yields a simple estimate of $\nu \lesssim 0.3$. Second, at a given collisionality, the toroidal ITG mode extends to shorter perpendicular wavelength with respect to the slab case. At the same time, the peak of the growth rate of the toroidal ITG (which occurs near $k_\perp \sim 0.5$) is shifted at longer wavelengths in the parallel direction compared with the slab branch, consistently with the fact that the toroidal mode subsists in the $k_\parallel = 0$ limit. In all cases, the stabilization of the growth rate at large $k_\parallel$ is caused by Landau damping. Third, the peak growth rate of the toroidal ITG is approximately twice as large as the slab one at all collisionalities.

Figure 6. The ITG linear growth rate, $\gamma$, in the $(k_\perp, k_\parallel )$ parameter space, for increasing collisionality $\nu$ (from (a) to (c) and from ( d) to ( f)) and magnetic gradient strength $R_B$ (from (a–c) to (d–f)) at a temperature gradient strength of $\eta =3$. The stability boundary of the gyromoment hierarchy (dotted white line) and 6GM and 4GM (dotted and dot–dashed lines, respectively) are shown for comparison. Here, we fix $(P,J) = (18,6)$.

We now assess the effects of collisions on the critical value of the temperature gradient above which the ITG develops. For this purpose, we compute the ITG growth rate, maximized over $k_\parallel$ and $k_\perp$, as a function of $\nu$ and $\eta$. This allows us to evaluate the stability boundary shown in figure 7(a) for the slab and in figure 7(b) for the toroidal ITG modes and compared them with the predictions of the 6GM and 4GM models. The slab and toroidal ITG modes are destabilized by increasing $\eta$ at all collisionality and are damped (and/or suppressed) by FLR collisional effects. We remark that the gyromoment stability boundary agrees with the collisionless one as obtained from (3.3). The sITG stability boundary is highly sensitive to the level of collisionality, in particular for $\nu \gtrsim 0.1$. In fact, the critical temperature gradient for the sITG mode onset is increased from $\eta \sim 1$, when $\nu \lesssim 1$, to $\eta \sim 3$, when $\nu \gtrsim 1$, in order to overcome the strong FLR collisional stabilization effects. A similar observation can be made when using other GK collision operators, but not shown (see § 6). However, we remark that the increase of $\eta$ for the mode onset due to FLR collisional damping is expected to be much weaker in experiments because of lower values of $\nu$. On the other hand, the toroidal ITG mode stability boundary is not significantly affected by collisions and remains near the collisionless threshold. The reason of this difference can be inferred from figure 6 which shows a peak of the growth rate near $k_\perp \simeq 0.3$ when $\nu = 5$ in the toroidal case (figure 6d–f), while the mode is completely suppressed in the slab case (figure 6a–c). Finally, we notice that the 6GM and 4GM stability boundaries converge to a constant value near $\eta \simeq 0.5$ in the high-collisional limit. At low collisionality, $\nu \lesssim 0.1$, the ITG mode is completely suppressed according to the 4GM model because of kinetic effects associated with the parallel streaming.

Figure 7. The ITG growth rate maximized over $k_\parallel$ and $k_\perp$, denoted by $\max _{k_\parallel, k_\perp }(\gamma )$, as a function of the collisionality $\nu$ and normalized temperature gradient $\eta$ for the slab (a) and toroidal (b) ITG modes. The stability boundaries (which separates the unstable and stable modes) predicted by the gyromoment approach using the GK Coulomb collision operator and the 6GM and 4GM models are plotted by the white dotted, red dotted and red dotted–dashed lines, respectively. The values of the growth rates in the low collisionality limit are in agreement with (3.3).

6. Comparisons between collision operator models

In this section, we investigate the differences between the collision operator models introduced in § 2.2. In particular, we compare them as a function of collisionality and perpendicular wavenumber and measure their relative deviation with respect to the GK Coulomb collision operator in § 6.1. We consider different collisionality regimes and temperature gradient strengths. We demonstrate that energy diffusion and FLR effects are important in the collisional damping of ITG modes, and find that the deviations, when compared with the GK Coulomb collision operator, decrease with temperature gradient strength at all collisionalities. One of the main results of this comparison is that the smallest relative deviation (${\lesssim }15\,\%$) is produced by the GK Sugama operator, while the largest (${\lesssim }50\,\%$) is obtained from the GK Dougherty operator. The importance of FLR effects is highlighted in § 6.2 by investigating a SWITG typically peaking near $k_\perp \sim 1.5$ that can be destabilized if collisional FLR terms are neglected. The present results are particularly important since DK collision operators are often used in the continuum GK codes to simulate the plasma boundary. In this section, we truncate the gyromoment hierarchy, (2.14), at $(P,J) = (18,6)$ for all cases. We have checked that convergence is achieved in all cases, in agreement with the convergence studies performed in § 8.

6.1. Comparisons between the GK coulomb and collision operators models

We focus on the differences between the GK Coulomb collision operator, the GK/DK Sugama collision operator (Sugama et al. 2009), the GK/DK momentum-conserving pitch-angle scattering (pitch) operator (Helander & Sigmar 2002), the zeroth-order DK HSC collision operator (Hirshman & Sigmar 1976a) and the GK/DK Dougherty collision operator (Dougherty 1964). The gyromoment expansion of the Sugama collision operator is detailed in Frei et al. (2021), while the gyromoment expansions of the other collision operators are reported in appendices A–C.

Figure 8 shows the slab (figure 8a–c) and the toroidal (figure 8d–f) ITG linear growth rate at a relative temperature gradient strength of $\eta =3$ from low (figure 8a,d) to high (figure 8c, f) collisionality as a function of the perpendicular wavenumber $k_\perp$ using the GK collision operators. The collisionless and high-collisional limits are plotted for comparison. It is remarkable that the ITG growth rate is sensitive to both the collision operator models and to the presence of FLR related terms. In fact, the GK Dougherty collision operator produces the results that mostly differ with respect to the other GK operator models, yielding the strongest ITG damping. The effect of energy diffusion can be observed by comparing the GK pitch-angle and GK Sugama collision operators. While the former systematically overestimates the ITG growth rate (compared with the GK Coulomb operator) due to the absence of energy diffusion, the latter damps stronger the ITG mode with growth rates smaller (but closer) than the GK Coulomb. Finally, we remark that the GK Coulomb operator predictions approach the ones of the 6GM and 4GM models in the limit of long perpendicular wavelengths for all collisionalities, since the FLR collisional damping becomes negligible when $k_\perp \ll 1$.

Figure 8. Comparisons between GK (solid lines) collision operator models in the case of the slab (a–c) and toroidal (d–f) ITG growth rate as a function of the perpendicular wavenumber $k_\perp$, from the low (a,d) to high (c, f) collisionality regime. The collisionless (dashed black) and 6GM and 4GM (red dotted and dashed–dotted lines, respectively) growth rates are plotted for comparisons. Here, $\eta = 3$, $k_\parallel = 0.1$.

By neglecting the FLR terms, the DK collision operators lead to a larger ITG growth rate than the GK operators as shown in figure 9. Contrary to the latter, the DK collision operators can eventually destabilize a SWITG peaking near $k_\perp \simeq 1.5$. The presence of this SWITG can be inferred from figure 9 from the increase of the ITG growth rate occurring for $k_\perp \gtrsim 1$ when $\nu \gtrsim 0.1$ (see § 6.2). This behaviour is particularly pronounced in the pitch angle and HSC collision operators at short wavelengths, and it is related to the absence of energy diffusion in these operators (Ricci, Rogers & Dorland 2006; Barnes et al. 2009), which is retained in the Coulomb, Sugama and Dougherty collision operators. We also notice that the absence of energy diffusion in the pitch angle and HSC collision operators affects the high-collisional limit. While the deviations between the DK Coulomb, Sugama and Dougherty collision operators are reduced as the collisionality increases and, ultimately, becomes negligible when $\nu \gg 1$, converging to the 6GM and 4GM solutions, a difference subsists in the pitch-angle and HSC collision operators. Finally, among the DK collision operators considered, the DK Sugama approaches better the DK Coulomb collision operator.

Figure 9. Same as figure 8 using the DK collision operator models.

The deviations of the collision operator models with respect to the GK Coulomb operator can be quantified by computing their signed relative difference of the growth rate, $\sigma _C = (\gamma _C - \gamma ) / \gamma _C$ where $\gamma _C$ is the growth rate using the GK Coulomb and $\gamma$ is the growth rate obtained using the other collision operator models considered in figures 8 and 9. The results are shown in figure 10 where we plot $\sigma _C$ as a function of collisionality, $\nu$, and of temperature gradient strength, $\eta$, for the toroidal branch ($R_B = 0.1$). The GK Sugama collision operator provides the results that best approach the GK Coulomb collision operator. Among the other GK operators considered, the largest deviation is produced by the GK Dougherty where the relative difference is larger than $30\,\%$ when $\nu \gtrsim 0.5$. Surprisingly, while the GK pitch-angle operator systematically overestimates the ITG growth rate (as observed from $\sigma _C< 0$), it yields an error smaller than (or of the order of) $30\,\%$ in almost all the parameter space considered. In particular, the relative signed difference of the GK pitch-angle collision operator is similar in absolute value to the one of the GK Sugama at high $\eta$ and all collisionalities, but with opposite sign. We remark also that the negative sign of $\sigma _C$ for the DK collision operators is the result of ignoring the FLR collisional damping. For all the considered operators, the relative deviation with GK Coulomb operator decreases at steeper temperature gradients while it increases with collisionality.

Figure 10. Signed relative difference, $\sigma _C$ (defined in the main text), of the growth rate obtained using the collision operator models, used in figures 8 and 9, with respect to the GK Coulomb operator as a function of collisionality $\nu$ and temperature gradient strength $\eta$. The colourbars are saturated at a maximal relative deviation of $\sigma _C = \pm 0.4$. Here, $k_\parallel = 0.1$, $k_\perp = 0.5$ with $R_B = 0.1$ are considered, corresponding approximately to the peak ITG growth rate.

The gyromoment approach allows for the investigation of the ITG eigenvalue spectrum at finite collisionality. An example is shown in figure 11 for the toroidal case when $\nu = 0.1$ using the DK (figure 11a) and GK (figure 11b) pitch angle, Dougherty, Sugama and Coulomb collision operators near to the peak growth rate occurring at $k_\perp = 0.6$. First, focusing on the spectrum of the DK collision operators, we note that, while the DK Sugama reproduces qualitatively well the eigenvalues of the DK Coulomb, the DK Dougherty collision operator displays an eigenvalue spectrum with strongly damped subdominant modes characterized by $\gamma < 0$. Similar differences between the DK Dougherty and the DK Coulomb collision operators in the eigenvalue spectrum is reported in the study of electron plasma waves at arbitrary collisionality by Jorge et al. (2019b). In addition, despite the absence of energy diffusion, the DK pitch-angle operator features similar eigenvalues than the DK Coulomb. Focusing now on the spectrum of the GK collision operators, we remark that a large difference exists between the GK Dougherty operator and the other GK collision models. It is noticeable that the GK Coulomb collision operator produces subdominant modes that are less damped than the ones predicted by the GK Sugama and pitch-angle operators but, nevertheless, yields a similar eigenvalue spectrum. We remark that the eigenvalues have all a positive real frequency, $\omega _r > 0$, corresponding to the ion diamagnetic direction.

Figure 11. Toroidal ITG eigenvalue spectrum when $\nu = 0.1$ and $k_\perp = 0.6$ using the DK (a) and GK (b) collision operator models. Here, $R_B = 0.1$ and $\eta = 3$.

6.2. Destabilization of the SWITG mode

An initial analysis of the effects of FLR terms in the collision operators can be carried out by observing, first, the discrepancies between the GK Coulomb and the 6GM and 4GM stability boundaries in the $\nu =1$ case (see figure 6) and, second, from the comparison between GK and DK collision operators (see figures 8 and 9). A clear deviation between GK and DK Coulomb operators occurs at sufficiently steep temperature gradients, $\eta \gtrsim 3$, and intermediate collisionality, i.e. when $\nu \gtrsim 0.1$. For a more detailed analysis, we plot in figure 12 the ITG growth rate evaluated with the DK Coulomb in the slab case for two different values of temperature gradient strength.

Figure 12. Slab SWITG growth rate (unstable when the DK Coulomb model is used) when $\eta =3$ (a–c) and $\eta =5$ (d–f). The stability boundaries of the gyromoment hierarchy (dotted white line), 6GM and 4GM (red dotted and dot–dashed lines, respectively) are plotted for comparison. The main ITG branch is identified near $k_\perp \simeq 0.5$ and the SWITG near $k_\perp \simeq 2$, appearing as the collisionality increases when $\nu \gtrsim 0.5$. The colourbar is saturated at the maximum of $\gamma$ when $\nu = 0.05$.

The comparison between figures 6 and 12 reveals that a short wavelength branch of the ITG, referred to as the SWITG mode and typically peaking at $k_\perp \sim 2$, is present when $\nu \gtrsim 0.5$ if the DK Coulomb operator is used. We remark that the SWITG is also captured by the 6GM and 4GM models. This mode has been identified in the collisionless limit as a continuous extension of the main branch of the ITG mode in the $k_\perp > 1$ region. In fact, the collisionless ITG growth rate reveals a typical ‘double-humped’ behaviour when plotted as a function of $k_\perp$ (Smolyakov, Yagi & Kishimoto 2002; Gao et al. 2003, 2005). Additionally, figure 12 shows the presence of a stability region that increases with collisionality, isolating the conventional ITG and the SWITG modes. We note that, in the toroidal case, the SWITG is stabilized (or even suppressed) by increasing magnetic gradients, a similar feature observed also in the collisionless limit (Gao et al. 2005).

The SWITG mode is also predicted by all the other DK collision operators considered in this work, but it is suppressed by FLR collision damping in the GK collision operators. This questions the applicably of DK collision operator models to correctly predict the turbulent transport level in ITG driven turbulence in the boundary region, since the SWITG can produce an additional (yet spurious) contribution in addition to the main ITG mode peaking at larger $k_\perp$ (Gao et al. 2005). However, further investigations are required to go beyond the local approximation used in the present work.

7. Gyromoment spectrum

We now study the ITG gyromoment spectrum. In particular, we numerically evaluate the collisionless gyromoment spectrum using the semianalytical expressions obtained in § 3.2. This reveals that magnetic drift resonance effects broaden the gyromoment spectrum compared with the slab case. The effect of collisions is also considered by deriving asymptotic expressions of the spectrum in the limit of a large number of gyromoments. We first focus on the slab limit and, then, consider the purely toroidal ($k_\parallel = 0$) case. The results presented in this section are used to analyse the convergence properties in § 8.

Considering the sITG mode, we focus on the parameters $k_\parallel = 0.1$, $\eta =5$ and $k_\perp = 0.5$ (i.e. near the growth rate peak). We first note that the solution of the GK dispersion relation, given in (3.3), yields the solution $\omega = 0.1754 + 0.1504 {\rm i}$. Given $\omega$, the normalized collisionless gyromoments spectrum, $N^{pj} / \phi$, can be evaluated using (3.11) (or equivalently (3.3) with $R_B = 0$). This is shown in figure 13 as a function of $(p,j)$ on a logarithmic scale. As observed, the collisionless gyromoment spectrum decreases rapidly with $j$, while the decrease of the magnitude of the gyromoments is slower along the $p$ direction, as highlighted in figure 13(b) that shows $N^{p0} / \phi$ and $N^{0j} / \phi$. The finite amplitude of the gyromoments $N^{pj} / \phi$ at $j \neq 0$ (and small $p$) is due to the presence of FLR effects, yielding terms proportional to $\mathcal {K}_{j}$ in the collisionless gyromoment response (see (3.11)) and, to a smaller extent, to the presence of temperature gradients, yielding terms proportional to $\eta \mathcal {K}_{j \pm 1}$ in (3.11). As a consequence, the gyromoment amplitude rapidly decreases with $j$ because of the behaviour of the FLR kernel $\mathcal {K}_{j}$ for large $j$, such that $\mathcal {K}_{j} \sim (b/2)^{2j}/j!$. We remark that the decrease of modulus of $\hat {N}^{pj}$ with $j$ is slower at large $b$ and follows, in general, $1/j!$. On the other hand, the gyromoments with $p \neq 0$ arise because of the parallel resonance terms proportional to the $p$th derivative of the plasma dispersion functions, $Z^{(p)}(\xi ) / \sqrt {2^{p} p!}$, as illustrated by figure 13. The extend of the spectrum in figure 13 demonstrates that the sITG mode requires a larger number of Hermite than Laguerre gyromoments, i.e. $P > J$, to be resolved.

Figure 13. Modulus of the collisionless gyromoments spectrum $N^{pj} / \phi$ (normalized to the maximum value) in the case of sITG, as a function of $p$ and $j$ plotted on a logarithmic scale and artificially saturated at $10^{-4}$. The modulus of the spectrum $N^{p0} / \phi$ and $N^{0j} / \phi$ are shown in panel (b) (solid blue and red lines, respectively). Here, $k_\parallel = 0.1$, $\eta =5$ and $k_\perp = 0.5$.

We now consider the collisionless gyromoment spectrum in the presence of magnetic gradients. As in the slab case, we focus on $k_\perp = 0.5$ and $\eta = 5$, but we choose $R_B = 0.5$ and $k_\parallel = 0$. We remark that setting $k_\parallel =0$ allows us to investigate the kinetic effects associated with the ${\rm i} \omega _{\boldsymbol {\nabla } B}$ term in (2.4) independently of Landau damping. The solution of the GK dispersion relation (3.3), for the complex mode frequency, yields $\omega = 0.7972 + 0.4364 {\rm i}$. The expressions given in (3.10) are used to evaluate the gyromoment spectrum. More precisely, the quantities $\mathcal {I}_n^{pj}$ are evaluated numerically using (D10) and the infinite sums in (3.10) are truncated when machine precision is reached. Similarly to figure 13, we plot the toroidal collisionless gyromoment spectrum, $N^{pj} / \phi$, on a logarithmic scale in figure 14. Compared with the slab case, the gyromoment spectrum is significantly broader. This is attributed to the velocity-dependence of the magnetic drift resonant coefficients, proportional to $1/[\omega - \omega _{\boldsymbol {\nabla } B}]$, that drives high-order gyromoments. In fact, while the dependence on $j$ of the gyromoments stems mostly from terms proportional to $\mathcal {K}_{j}$ in the slab case, the definite integrals of Gauss hypergeometric functions, appearing in $\mathcal {I}_{n}^{pj}$ (see (D10)), extend the gyromoment spectrum, producing a considerably slower decay than $1/j!$. We also notice that, since $k_\parallel =0$, the gyromoment spectrum vanishes for odd values of $p$, where $\mathcal {I}_n^{pj} = 0$ (see (D10)). Finite $k_\parallel$ yields non-vanishing odd gyromoments in $p$ due to the term proportional to $(1 + (-1)^{p-2p_1 +1})$ in (D8).

Figure 14. Modulus of the collisionless gyromoments spectrum $N^{pj} / \phi$ (normalized to the maximum value) in the case of purely toroidal ($k_\parallel = 0$) ITG with $R_B = 0.5$, as a function of $p$ and $j$ plotted on a logarithmic scale and artificially saturated at $10^{-4}$. The moduli of the spectrum of $N^{p0} / \phi$ and $N^{0j} / \phi$ are shown in panel (b) (solid blue and red lines, respectively). Here, $k_\perp = 0.5$, $\eta =5$.

Despite the fact that the gyromoment spectrum is significantly broader with respect to the slab limit when the magnetic drifts are important (large $R_B$), we numerically demonstrate the convergence of the gyromoment representation of the perturbed distribution function, ${g} /\phi$. Given the complex frequency $\omega$, a truncated distribution function ${g}^{PJ}$ can be defined from the Hermite–Laguerre expansion in (2.6) according to

(7.1)\begin{equation} {g}^{PJ} = \sum_{p=0}^{P} \sum_{j=0}^{J} \frac{\hat{N}^{pj}}{\phi} \frac{H_p(s_{{\parallel}}) L_j(x)}{\sqrt{2^{p} p!}} F_{M}. \end{equation}

From ${g}^{PJ}$ and ${g} /\phi$ obtained from (3.1), the error norm $e(P,J)$ can introduced by

(7.2)\begin{equation} e(P,J) = \int_{- \infty}^{\infty} \,{\rm d} s_\parallel \int_0^{\infty} x \left| {g}^{PJ} - \frac{{g}}{\phi} \right|^{2}, \end{equation}

and can be computed as a function of $(P,J)$. The results are shown in figure 15 for the slab and toroidal cases. Convergence is observed in both cases. However, while the magnitude of $e(P,J)$ becomes negligibly small when $(P,J) \gtrsim (20,10)$ (corresponding to $e(20,10) \lesssim 10^{-3}$) in the slab case, a larger number of gyromoments is necessary to achieve a similar error, i.e. $(P,J ) \gtrsim (50,20)$.

Figure 15. Error norm $e(P,J)$, (7.2), on a logarithmic scale for the slab (a) and toroidal (b) ITG cases. The same parameters as in figures 13 and 14 are used, respectively. The colourbars are saturated artificially at $\log _{10} e(P,J) = -4$.

We now turn to the collisional effects, which provide a natural cutoff of the gyromoment hierarchy. We aim to derive an ad hoc proxy of the scaling of $(P,J)$ at which the gyromoment hierarchy can be truncated as a function of the collisionality. In particular, we consider an equation for the generalized gyromoment amplitude $G ^{pj} = \sqrt {p} (g^{pj})^{2} /2$ (where $g^{pj} = (-1)^{j} i^{ p +j } \text {sgn}(k_\parallel )^{p} N^{pj}$ is introduced to remove the oscillatory behaviour of $N^{pj}$) using the gyromoment hierarchy equation, (2.14) in the limit $p,j \gg 1$ (Zocco & Schekochihin 2011; Loureiro, Schekochihin & Zocco 2013; Jorge et al. 2018). In this limit, the amplitude $g^{pj}$ can be considered as a continuous complex function of $p$ and $j$. This allows us to approximate $g^{p \pm 1j} \simeq g^{pj} \pm \partial _p g^{pj}$ (and $g^{p j \pm 1} \simeq g^{pj} \pm \partial _j g^{pj}$). Assuming that an unstable mode is considered, $\partial _t G^{pj} = 2 \gamma G^{pj}$, we derive the equation for $G^{pj}$, up to first order in $1/p$ and $1/j$, i.e.

(7.3)\begin{align} |k_\parallel |\sqrt{\tau} \partial_{s} G^{pj} + \frac{ \tau}{q} \omega_{B} ( 2 j +1) \partial_j G^{pj}+ \frac{2 \tau}{q} \omega_B \left( 1 + {\rm i} \left( 2j +1\right)\right) G^{pj} = 2 \nu (f_{p,j} + b^{\alpha} h_{p,j} - \gamma ) G^{pj}, \end{align}

with $s = \sqrt {p}$. In deriving (7.3), collisional effects are modelled using the ansatz $\mathcal {C}^{pi} \simeq \nu [ f_{p,j} + b^{\alpha } h_{p,j}] N^{pj}$, where the functions $f_{p,j} = f(p,j)$ and $h_{p,j} = h(p,j)$, and the constant $\alpha$ depends on the choice of the collision operator. We remark that the coupling along the Hermite direction associated with the curvature drifts vanishes at the leading order in $1/p$ in (7.3). We solve (7.3) for the slab and toroidal $k_\parallel = 0$ ITG cases using the GK Dougherty and GK Coulomb collision operators. While the functions $f_{p,j}$ and $h_{p,j}$ are evaluated numerically for the GK Coulomb collision operator, the sparse gyromoment expansion of the GK Dougherty collision operator, given in Appendix C, allows us to solve (7.3) analytically.

Using the gyromoment expansion of the GK Dougherty collision operator, we solve first (7.3) in the case of the sITG mode. With $f_{p,j} = - p - 2j$, $h_{p,j} = -1/2$ and $\alpha = 2$ (see (C3)), the solution of (7.3) is

(7.4)\begin{equation} \left| N^{pj}\right| \sim \frac{1}{p^{1/4}} \exp\left\{ - \left( \frac{p}{p_\gamma}\right)^{1/2} - \left( \frac{p}{p_{j \perp}}\right)^{1/2}- \left(\frac{p}{p_{D\nu}} \right)^{3/2} \right\}, \end{equation}

where we introduce $p_\gamma = \tau / \gamma _\parallel ^{2}$ (with $\gamma _\parallel = \gamma / k_\parallel$), $p_{j\perp } = \tau /[ \nu _\parallel (2 j + b^{2}/2)]^{2}$ ($\nu _\parallel = \nu / k_\parallel$) and $p_{D \nu } = (9 \tau / \nu _\parallel ^{2})^{1/3}$. The quantities $p_\gamma$, $p_{j \perp }$ and $p_{D \nu }$ introduce different decay scale lengths in the gyromoment spectrum. They are associated with the growth of mode ($p_\gamma$), with the $j$ dependence of $\mathcal {C}^{pj}$ ($p_{j \perp }$) and with the $p$ dependence of $\mathcal {C}^{pj}$ ($p_{D \nu }$). The $1 / \nu _\parallel$ dependence of $p_{D \nu }$ indicates that the exponential decay of (7.4) is slower at high $k_\parallel$ and faster at high $\nu$. At high $k_\parallel$, higher-order $p$ gyromoments are excited due to parallel resonance effects. At high $\nu$, the Chapman–Enskog closure procedure ensures that $\nu _\parallel \sim \sqrt {\tau } / \epsilon \gg 1$, yielding $p_{D\nu } \sim \nu _\parallel ^{-2/3} \sim \epsilon ^{2/3} / \tau ^{1/3}$ such that that the ITG dynamics can be described by the lowest-order gyromoments at high collisionality when $\epsilon \ll 1$ (see § 4.1). Moreover, when $\nu \gg \gamma$ (such that $p_{D\nu }$ precedes $p_\gamma$) and in the DK limit (such that the $b^{2}$ term in $p_{j \perp }$ is neglected), the collisional damping provides a truncation at $P \gtrsim p_{D \nu }$. We remark that FLR terms in the collision operator have a small effect on the scaling of (7.4) at finite collisionality. In fact, near the ITG peak growth rate (i.e. $k_\perp \sim 0.5$), one typically has $j \gtrsim b^{2}/4$. Finally, we remark that the $j$ dependence in (7.3) becomes important when the second exponent becomes comparable in magnitude to the last one. Despite that $p_{j \perp } / p_{ D\nu } \sim \epsilon ^{4/3} / j^{2} \ll 1$ at large collisionality and at large $j$, the different exponents $1/2$ and $3/2$ in (7.4) imply that this happens when

(7.5)\begin{equation} \frac{p_{D \nu}^{3/2}}{p (p_{j \perp})^{1/2}} \sim 1, \end{equation}

which is satisfied when $j \sim p / 10$. This shows that the decay of the gyromoments becomes faster with $j$ because of the increase of the second exponent in magnitude.

We now derive the solution of (7.3) for the toroidal $k_\parallel = 0$ ITG case, still considering the GK Dougherty collision operator. It yields

(7.6)\begin{equation} \left| N^{pj}\right| \sim \frac{1}{p^{1/4}} \exp\left\{ - \nu_B j- \frac{1}{2}\ln\left( 1 + 2j \right) \left[ \gamma_B + 1 + \nu_B \left( p - 1+ b^{2}/2 \right) \right]\right\}, \end{equation}

where we introduce $\gamma _B = q \gamma /( \tau \omega _B)$ and $\nu _B = q \nu / (\tau \omega _B)$. The $\nu _B$ dependence highlights the broadening of the gyromoment spectrum associated with the presence of magnetic gradient drifts (see figure 14), an effect reduced at high collisionality. Similar to the slab case, FLR terms do not affect the scaling of (7.4), since the linear term dominates over the second term at large $j$ values, and $p-1 \gg b^{2} /2\sim k_\perp ^{2}$ for ITG modes.

Because the Coulomb collision operator has a non-trivial block matrix structure when projected onto the Hermite–Laguerre basis (Frei et al. 2021), an approximate expression of $\mathcal {C}^{pj}$ needs to be determined numerically. Based on the results derived from the GK Dougherty collision operator, we assume that $\mathcal {C}^{pj} \simeq \nu f_{p,j} N^{pj}$, thus neglecting the FLR effects in the spectrum scaling. Fitting numerically the values of $\mathcal {C}^{pj}$ using the GK Coulomb collision operator yields $f_{p,j} \simeq - A j - B \sqrt {p}$ ($A \simeq 0.3$ and $B \simeq 0.8$). Similar values (within a $10\,\%$ accuracy) are typically found in the case of the DK Coulomb and show that FLR terms have a small effect on the ITG spectrum scaling when $k_\perp \sim 1$. We remark that the $\sqrt {p}$ dependence of $f_{p,j}$ is consistent with previous DK moment hierarchy using the DK Coulomb collision operator (Jorge et al. 2018). Given $f_{p,j}$, we again solve (7.3) for the case of the slab and toroidal $k_\parallel =0$ ITG modes yielding, respectively,

(7.7)\begin{equation} \left| N^{pj} \right| \sim \frac{1}{p^{1/4}} \exp\left\{ - \left( \frac{p}{p_\gamma}\right)^{1/2} - \left( \frac{p}{p_j} \right)^{1/2} - \frac{p}{p_\nu} \right\}, \end{equation}

with $p_j = \tau / \nu _\parallel A j^{2}$ and $p_\nu = 2 \sqrt {\tau } / (B \nu _\parallel )$ and

(7.8)\begin{equation} \left| N^{pj} \right| \sim \frac{1}{p^{1/4}} \exp\left\{ - \frac{ A}{2} \nu_B j - \frac{1}{2}\ln\left( 1 + 2j \right) \left[ \gamma_B + 1 + \nu_B \left( B \sqrt{p} - A /2 \right) \right] \right\}. \end{equation}

The main differences between the results of the GK Coulomb collision operator, given in (7.7) and (7.8) and the results of the GK Dougherty operator, in (7.4) and (7.6), are that the scaling in $p$ predicted by the GK Coulomb collision operator is weaker in the slab case (compared with a $3/2$ exponent in (7.4)), and that the numerical factor in the linear term in $j$ in the toroidal case, $A / 2 \simeq 0.15$, is smaller.

We compare the scaling of the gyromoment spectra obtained using the GK Dougherty and GK Coulomb collision operators in figure 16 for the slab and toroidal ITG modes for different values of $\nu _\parallel$ (at $k_\parallel = 0.1$) and $\nu _B$ (at $R_B = 0.1$). The differences between the $p$ exponent and the numerical factors between GK Dougherty and GK Coulomb, do not significantly impact the gyromoment spectra that feature the same qualitative behaviour. The differences between the GK Dougherty and GK Coulomb observed when $\nu _\parallel = 1$ and $\nu _B =2$ stem from the difference in the exponents and numerical factors in the spectrum. Finally, we remark that the asymptotic dependence of the expansion coefficients, i.e. $N^{pj}$, on their indices $p$ and $j$ (see, e.g. (7.7) and (7.8)) shows that the gyromoment method features a spectral convergence rate that increases with collisionality (Boyd 2001).

Figure 16. Normalized slab (a) and toroidal (b) gyromoment spectra as a function of $p$ and $j$, respectively, obtained using the GK Coulomb collision operator from (7.7) and (7.8) (solid lines) and using the GK Dougherty collision operator from (7.4) and (7.6) (dashed lines). The slab and toroidal growth rates are estimated at $\gamma \simeq 0.1$ and $\gamma \simeq 0.2$ (see figure 8), respectively, near the ITG peaks ($k_\perp \simeq 0.5$).

8. Convergence studies

Following the analysis of the gyromoment spectrum presented in § 7, we perform convergence studies of the ITG growth rate as a function of collisionality in § 8.1, perpendicular wavenumber in § 8.2 and magnetic gradient in § 8.3 for various numbers of gyromoments. Collisions are modelled using the GK and DK Coulomb collision operators (Frei et al. 2021), since similar convergence properties are observed with the other collision operator models considered in this work. We show that the convergence of the gyromoment approach improves with collisionality.

8.1. Collisionality convergence

We first consider the convergence of the ITG linear growth mode as a function of the collisionality, $\nu$, using the GK and DK Coulomb collision operators. The collisionless and high-collisional limits are obtained by solving the collisionless GK dispersion relation, (3.3), and the 6GM and 4GM models (see § 4.1). The results are shown in figure 17(a) for the slab and figure 17(b) for the toroidal ITG modes. Different numbers of gyromoments, $(P,J)$, are used. The collisionless GK solution is retrieved for $(P,J) \simeq (18,6)$ with both the GK (dashed) and DK (solid) Coulomb collision operators. In the slab case, convergence is achieved with $P \lesssim 8$, when $\nu \simeq 0.1$, and with $P \gtrsim 18$, when $\nu \lesssim 10^{-3}$. This is in qualitative agreement with the estimate of $P$ that can be obtained from (7.7). In fact, solving the gyromoment hierarchy up to $P \gtrsim p_\nu$ yields $p_\nu \sim 3$ for $\nu = 0.1$ and $p_\nu \sim 20$ for $\nu = 10^{-3}$. Also consistent with the observations made in § 7 is the fact that only a few numbers of the Laguerre gyromoments are necessary to resolve the sITG, as illustrated by the overlap between the $(10,4)$ and $(10,6)$ lines in figure 17. This is because the sITG mode is mainly driven by the parallel dynamics, resolved by the Hermite part of the gyromoment spectrum (see figure 13), and the relatively small value of $k_\perp$ considered in this case (see figure 18). This remains also true for the toroidal ITG case considered in figure 17 because of the relatively small value of $R_B$ used. In fact, larger values of $R_B$ deteriorates the convergence as discussed below. While the GK Coulomb yields a strong FLR collisional damping as $\nu$ increases, the gyromoment hierarchy with DK Coulomb agrees with the 6GM and 4GM when $\nu \gg 1$. Nevertheless, we notice that, because of finite magnetic drifts (see § 3), agreements with the 6GM and the 4GM models occur at higher collisionality in the toroidal case. In fact, in the slab case, the high-collisional limits are retrieved when $\nu \gtrsim 1$, while for the toroidal case ($R_B = 0.1$) when $\nu \gtrsim 10^{2}$.

Figure 17. Convergence of the ITG linear growth rate, $\gamma$, as a function of the collisionality $\nu$ for different values of $(P,J)$, using the GK (dashed lines) and DK (solid lines) Coulomb operators, in the case of slab (a) and toroidal ((b) with $R_B = 0.1$) ITG branches. The collisionless limit (dotted black lines) and the high-collisional limits, 6GM and 4GM models (red dotted and dash–dot lines, respectively) are shown for comparison. Similar plots are obtained for the other GK and DK operator models. The parameters are $k_\perp = 0.5$, $k_\parallel = 0.1$, $\eta = 3$.

Figure 18. Slab ITG growth rate $\gamma$ as a function of $k_\perp$ for different values of $J$, in the low (a), intermediate (b) and high (c) collisionality regimes, with the GK (dashed lines) and DK (solid lines) Coulomb collision operators. Similar plots are obtained for the other GK and DK operator models.

8.2. Perpendicular wavenumber convergence

We now investigate the convergence associated with the FLR effects due to the gyroaverage of the electrostatic potential, i.e. with the number of Laguerre gyromoments. The coupling between the Laguerre gyromoments driven by magnetic gradients are neglected by focusing on the sITG. We scan the growth rate $\gamma$ as a function of the perpendicular wavenumber, $k_\perp$, for different values of $J$ and collisionality $\nu$. The number of the Hermite gyromoment is fixed at $P = 18$. The results are shown in figure 18. We observe that convergence occurs at lower $J$ at high collisionalities. Nevertheless, the value of $J$ necessary for convergence increases with $k_\perp$ at all collisionalities. This is in agreement with the unfavourable scaling of the FLR kernel with $b$, i.e. $\mathcal {K}_{j} \sim (b/2)^{2j}/j!$ (see § 7). We also notice that the ITG mode is strongly damped by the GK Coulomb collision operator that departs rapidly from the DK Coulomb as $k_\perp$ and $\nu$ increases, in agreement when comparing figure 8 with figure 9. The same qualitative observations hold in the toroidal case.

8.3. Magnetic gradient convergence

Finally, we consider the convergence of the gyromoment hierarchy as a function of the normalized magnetic gradient, $R_B$. As discussed in § 7, the presence of finite magnetic drifts broadens significantly the collisionless gyromoment spectrum. However, even at low collisionality, the gyromoment approach converges correctly to the collisionless solution $\gamma _{\text {GK}}$ obtained from (3.3). This is shown in figure 19 where the ratio of the ITG growth rate obtained using the gyromoment hierarchy to $\gamma _{\text {GK}}$ is plotted as a function of $P$ and $J$. It is remarkable that, consistently with § 7, a larger number of gyromoments is needed as $R_B$ increases.

Figure 19. Ratio between the ITG growth rate $\gamma$ obtained using the gyromoment hierarchy and the collisionless solution $\gamma _{\textrm {GK}}$ as a function of $P$ with $J = 10$ (a) and as a function of $J$ with $P=32$ (b) for two different values of $R_B$. Here, $k_\perp =0.5$, $k_\parallel = 0.1$ and $\nu = 0.001$.

Collisions can compete against the strong kinetic drive of the magnetic drifts by damping higher-order gyromoments (see figure 16). To investigate the convergence of the gyromoment approach at finite $R_B$ in the presence of collisions, we focus on the toroidal ITG with $k_\parallel =0$ and, hence, neglecting the coupling between the Hermite gyromoments associated with the parallel streaming. We scan the toroidal ITG linear growth rate $\gamma$ as a function of $R_B$ and increase the collisionality. The results are reported in figure 20. We observe that, at all collisionalities, the convergence of the gyromoment approach deteriorates as $R_B$ increases. Nevertheless and consistently with § 7, the number of gyromoment to achieve convergence is reduced as the collisionality increases, demonstrating the competition between the kinetic effects driven by the magnetic gradient drifts and the collisional effects. We also remark that the 4GM limit is retrieved at higher collisionality when $R_B$ increases.

Figure 20. Toroidal ($k_\parallel = 0$) ITG growth rate $\gamma$ as a function of the normalized magnetic gradient, $R_B$, for different $(P,J)$ and increasing collisionality (from (a) to (c)). Here, the GK (dashed lines) and DK (solid lines) Coulomb collision operators are shown, with the collisionless (black dotted) and high-collisional 4GM limit (red dotted–dashed lines).

9. Conclusion

In this paper, the properties of the ITG mode are studied using the GK Coulomb collision operator for the first time. The investigation is based on a Hermite–Laguerre moment expansion of the perturbed gyrocentre distribution function, which we refer to as the gyromoment approach. By projecting the GK Boltzmann equation onto the Hermite–Laguerre basis, a gyromoment hierarchy equation is deduced, which retains the effects of like-species collisions between ions thanks to the numerical implementation of the GK Coulomb collision operator described in Frei et al. (2021). Using the gyromoment hierarchy equation, the collisionless and high-collisional limits are derived analytically. In the collisionless limit, we find that the magnetic drift resonance effects significantly broaden the gyromoment spectrum compared with the slab case. In the high-collisional limit, the gyromoment hierarchy is reduced to a fluid model retaining only a finite number of gyromoments where the Chapman–Enskog asymptotic closure scheme is used. Using these analytical results, we show that the gyromoment hierarchy can retrieve the collisionless limit, where kinetic features are essential, and, at the same time, the high-collisional limit using a reduced number of gyromoments.

Numerical experiments are performed to study the ITG linear growth rate using the GK Coulomb collision operator and investigate the importance of the FLR collisional terms. We compare the GK Coulomb collision operator with other collision operator models such as the Sugama (Sugama et al. 2009), the momentum-conserving pitch-angle scattering operator (Helander & Sigmar 2002), the zeroth-order DK HSC collision operator (Hirshman & Sigmar 1976a) and the Dougherty collision operator (Dougherty 1964). The gyromoment expansions of the pitch angle and HSC collision operators are also derived for the first time here. We find that the ITG mode is strongly damped (or even suppressed) as the collisionality increases and that a steeper temperature gradient for the mode onset is necessary to overcome the FLR collisional stabilization when using the GK Coulomb collision operator. We reveal the importance of FLR terms in the Coulomb collision operator by demonstrating that neglecting these FLR terms destabilizes a short wavelength branch of the ITG mode, peaking near $k_\perp \rho _s \sim 1.5$. These observations are also found when using the other collision operator models considered in this work. In addition to these findings, the main outcome of the systematic comparison between the collision operators considered in this work is that the GK Sugama slightly underestimates the ITG growth rate compared to the GK Coulomb operator and that the GK pitch-angle scattering operator systematically yields a larger ITG linear growth because of the absence of energy diffusion in the latter. The largest deviations with respect to the GK Coulomb collision operator are observed when the GK Dougherty operator is used. In general, we observe that energy diffusion in the collision operators is important at high collisionality and in the damping of ITG at small scales. Finally, the numerical efficiency of the gyromoment approach is highlighted by performing convergence studies and show that the number of the gyromoment can be reduced as the collisionality increases. In particular, the ability of the presented method to naturally reduce its dimensionality (i.e. the number of gyromoments) without compromising the accuracy of kinetic effects and the details of advanced collision operators makes it an ideal framework to develop high-fidelity simulations of the boundary region, a possible advantage over present continuum GK codes.

Complementary to previous works (see, e.g. Jorge et al. 2017, 2018; Frei et al. 2020, 2021), the present study of the ITG mode at arbitrary collisionality shows the ability of the gyromoment approach to capture the relevant physics in the kinetic and collisional regimes, with a fidelity cost that depends on the collisionality. The present study can be extended to more realistic conditions, including variation along the parallel direction, electron dynamics, electromagnetic effects and nonlinear physics. In particular, nonlinear simulations using the gyromoment approach are required to further investigate and to extrapolate the effects of collision operator models, FLR terms and convergence reported in the present work to, for example, ITG driven turbulence. For instance, the role of FLR collisional damping by the GK collision operators of the ITG (see, e.g. figure 11) is a mechanism that is expected to affect the level of the saturated turbulent fluxes. Additionally, the choice of the collision operator models (e.g. Coulomb, Dougherty and pitch-angle operators) used in the simulations can play an important role in the nonlinear saturated state, especially in the high-collisional regime where their analytical details can no longer be ignored. Nonlinear simulations are being carried out and will be subject to a future publication. Finally, the $\delta f$ formulation adopted in this work might not be valid for the study of the turbulent dynamics in the boundary region, where the large amplitude fluctuations require a full-F gyromoment hierarchy equation, a subject of future work.