# Local Gyrokinetic Collisional Theory of the Ion-Temperature Gradient Mode

B. J. Frei1<EMAIL_ADDRESS>A. C. R. Hoffmann 1  P. Ricci1  1 École
Polytechnique Fédérale de Lausanne (EPFL), Swiss Plasma Center (SPC), CH-1015
Lausanne, Switzerland

###### Abstract

We present a study of the linear properties of ion temperature gradient (ITG)
modes with collisions modelled by the linearized gyrokinetic (GK) Coulomb
collision operator (Frei et al., 2021) 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 gyro-
moments. We explore analytically the collisionless and high-collisional limits
of the gyro-moment hierarchy. Parameter scans revealing the dependence of the
ITG growth rate on the collisionality are reported, showing strong damping at
small scales as the collisionality increases. These properties are compared
with the predictions based on the Sugama, the momentum-conserving pitch-angle
scattering, the Hirshman-Sigmar-Clarke, and the Daugherty collision operators.
The importance of finite Larmor radius (FLR) terms in the collision operators
is pointed out by the appearance of a short wavelength (SW) ITG branch when
collisional FLR terms are neglected, this branch being completely suppressed
by collisional FLR effects. We demonstrate that energy diffusion is important
at high collisionality and small scale lengths and that, among the collision
operators considered in this work, the GK Sugama collision operator yields, in
general, the smallest deviation on the ITG growth rate compared to the GK
Coulomb collision operator. Convergence studies of the gyro-moment method are
reported.

##  1 Introduction

Belonging to the class of instabilities that develop on the ion-gyroscale, the
ion temperature gradient (ITG) modes are widely recognised as the main
candidate to explain the experimental observations of anomalous ion heat
turbulent transport in the tokamak core (Garbet et al., 2004). The ITG modes
can be driven unstable by the parallel plasma compression (slab ITG) and by
the presence of magnetic drifts (toroidal ITG). In the core region, since the
plasma collision frequency is considerably smaller than the typical ITG mode
frequency, it is usually assumed that collisions do not affect significantly
the ITG dynamics. In fact, the GK collisionless theory of the ITG modes is
well documented in the literature (see, e.g., Romanelli, 1989; Hahm & Tang,
1989; Romanelli & Briguglio, 1990; Chen et al., 1991). Because of the
temperature drop, the effects of collisions are expected to be important in
the tokamak boundary, i.e. the region that encompasses the edge and the
scrape-off layer (SOL). Collisions are known to modify the linear and
nonlinear properties of micro-instabilities. For instance, collisions
contribute to smear out small scale structures in velocity space (Watanabe &
Sugama, 2004; Ajay et al., 2021), to affect the linear properties of the ITG
mode and the transition to the trapped electron (TEM) mode (Manas et al.,
2015; Belli & Candy, 2017), to suppress short wavelength structures (Barnes et
al., 2009), and can ultimately modify the damping of the zonal flow and,
hence, the saturation mechanism of turbulent fluxes (Pan et al., 2020; Frei et
al., 2021). Also, the collisionality dependence of the ITG critical gradient
is important in the prediction of the transition from resistive ballooning
mode and ITG in the edge region (Zeiler et al., 1998). In general, despite the
fact that the boundary region is particularly important since it sets the heat
exhaust and the thermal load on the vessel walls, determines the overall
confinement capabilities of the device, controls the level of impurities, and
regulates the removal of helium ash, the role of the ITG mode is less clear in
this region than in the core.

The effects of collisions 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 models (Romanelli & Briguglio,
1990; Romanelli et al., 1991; Pusztai et al., 2009). Among the documented
studies of discrepancies produced by the collision operator models (Catto &
Ernst, 2009; Belli & Candy, 2017; Pan & Ernst, 2019; Pan et al., 2020) with
respect to the linearized Fokker-Planck collision operator (Rosenbluth et al.,
1957; Hazeltine & Meiss, 2003) which we refer to as Coulomb collision operator
in the present study, only a few have reported parameter scans and direct
operator comparisons when applied to ITG modes. In general, the limited number
of investigations of the impact of collisions on the ITG mode, and the related
limitations of the understanding of the role of ITG mode on the boundary
region stems from the frequent use of drift-reduced Braginskii fluid models
(Zeiler et al., 1997) for the simulation of this region (Dudson et al., 2009;
Tamain et al., 2016; Halpern et al., 2016; Paruta et al., 2018; Zhu et al.,
2018; Stegmeir et al., 2019; Giacomin et al., 2020). These fluid models are
based on the assumptions that $k_{\parallel}\lambda_{mfp}\ll 1$, i.e. that the
particle mean-free path $\lambda_{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 ion temperature gradient 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_{mfp}\ll 1$ assumption is well established at the
separatrix and in the SOL (with temperature $T\lesssim 10$ eV and density
$n\sim 10^{18}$ m-3 yielding typically $k_{\parallel}\lambda_{mfp}\lesssim
10^{-2}$ in, e.g., DIII-D and ITER), the mean free path is larger inside the
separatrix where the temperature increases. Hence, for typical H-mode pedestal
values ($T\sim 10^{2}-10^{3}$ eV and $n\sim 10^{20}$ m-3), an estimate of
$k_{\parallel}\lambda_{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.

To retain collisionless kinetic and finite Larmor radius (FLR) effects,
gyrofluid models have been developed by taking explicitly a finite number of
velocity space moments of the GK equation (Hammett et al., 1992; Dorland &
Hammett, 1993; Beer & Hammett, 1996; Snyder & Hammett, 2001). Collisionless
kinetic effects are then introduced in these models by designing ad hoc
closures for the high-order moments that mimic the kinetic linear response in
the absence of collisions. While usually neglected, collisional effects are
introduced in, e.g., 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 electrons
collisions. However, BGK-like operators and the absence of energy diffusion in
the pitch-angle scattering operator can potentially lead to significant
deviations when compared to the Coulomb collision operator.

In this work, we aim to study 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 Coulomb
collision operator. This is done by leveraging previous works that led to a
formulation of the GK model based on a Hermite-Laguerre expansion of the
full-F gyrocenter distribution function (Frei et al., 2020), and the
development of a full-F GK Coulomb collision operator derived using the same
technique expansion, as presented in Jorge et al. (2019a). More precisely, the
ITG model we use in this work is based on the linearized and local limit of
this formulation of the GK model, coupled with a linearized GK Coulomb
collision operator, 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 gyrocenter distribution function, that we
refer to as gyro-moments. The Coulomb collision operator we use is valid for
arbitrary mass and temperature ratios of the colliding species and arbitrary
values of the perpendicular wavenumber. By using the gyro-moment approach, we
first retrieve the collisionless and the high-collisional ITG limits. Using
the linearized GK Coulomb collision operator, wide parameters scans of the ITG
growth rate are then performed. We also compare 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). The gyro-moment approach is systematically applied to all the
collision operators considered in this work, allowing their numerical
implementation. We show the importance of energy diffusion and FLR terms in
the collisional damping of small-scale modes. 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 demonstrate the importance of retaining FLR terms in the
collision operator models and, in particular, that a short wavelength ITG
(SWITG) mode at high collisionality and steep ion temperature gradients can be
destabilized when the DK limit of the collision operator 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; Ulbl et al., 2021;
Francisquez et al., 2021). Finally, we investigate the convergence properties
of the gyro-moment approach as a function of relevant physical parameters in
different collisionality regimes. Additionally, we analyse the Hermite-
Laguerre spectrum of the ITG distribution function and show that the magnetic
drifts broaden significantly the gyro-moment spectrum in the absence of
collisions.

The structure of the paper is the following. In Sec. 2, we introduce the gyro-
moment hierarchy equation. Then, in Sec. 3, we derive the collisionless limit
of the ITG dispersion relation by leveraging the gyro-moment hierarchy
equation, and Sec. 4 focuses on the high-collisional limits. In Sec. 5, we
evaluate the ITG growth rate over a wide range of parameters with collisions
modeled by the GK Coulomb collision operator. Sec. 6 is focused on the
comparison of collision operator models with the GK Coulomb collision
operator. Finally, we analyze the Hermite-Laguerre spectrum and report
convergence studies in Sec. 7 and Sec. 8, respectively. The conclusion follow
in Sec. 9. The Hermite-Laguerre expansions of the collision operators used in
this work are reported in Appendices A, B and C, while analytical results are
detailed in Appendix D.

##  2 Gyrokinetic Model

Focusing on a shearless slab configuration and adopting a local GK approach
where quantities do not vary along the magnetic field line (Kadomtsev &
Shafranov, 2012), we consider the linearized electrostatic GK Boltzmann
equation expressed in gyrocenter phase-space coordinates
$(\bm{R},\mu,v_{\parallel},\theta)$ and assume that the electrons are
adiabatic. Here, $\bm{R}=\bm{r}-\bm{\rho}$ is the ion gyrocenter position
(with $\bm{r}$ being the particle position and
$\bm{\rho}(\bm{R},\mu,\theta)=\bm{b}\times\bm{v}/\Omega$ the ion gyroradius
with $\Omega=qB/m$ the ion gyrofrequency, $q$ the ion charge and
$\bm{b}=\bm{B}/B$), $\mu=mv^{2}/[2B(\bm{R})]$ is the magnetic moment,
$v_{\parallel}=\bm{b}\cdot\bm{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$, temperature $T$ and magnetic field strength $B$ by
$L_{N}=\left|\gradient_{\perp}\ln N\right|^{-1}$,
$L_{T}=\left|\gradient_{\perp}\ln T\right|^{-1}$ and
$L_{B}=\left|\gradient_{\perp}\ln B\right|^{-1}$, respectively. We assume that
the gyroaveraged ion gyrocenter distribution function,
$F=F(\bm{R},\mu,v_{\parallel})$, is a perturbed Maxwellian, that is
$F=F_{M}+\textsl{g}$, being g the perturbed part ($\textsl{g}/F_{M}\ll 1$) and
$F_{M}$ a local Maxwellian distribution function defined by
$F_{M}=N(\bm{R})e^{-s_{\parallel}^{2}-x}/(\pi^{3/2}v_{T}(\bm{R})^{3})$ with
$s_{\parallel}=v_{\parallel}/v_{T}(\bm{R})$, $x=\mu B(\bm{R})/T(\bm{R})$, and
$v_{T}(\bm{R})^{2}=2T(\bm{R})/m$ the ion thermal velocity. The linearized
electrostatic GK Boltzmann equation describing the time evolution of a single
Fourier harmonic $\textsl{g}(\bm{k},t)=\int
d\bm{R}\textsl{g}(\bm{R},t)e^{-i\bm{R}\cdot\bm{k}}$ is then given by

$\displaystyle\frac{\partial}{\partial t}\textsl{g}$
$\displaystyle+\left(i\omega_{B}+ik_{\parallel}v_{\parallel}\right)\textsl{g}-i\omega_{*}^{T}J_{0}(b\sqrt{x})\frac{e\phi}{T_{e}}F_{M}+i\omega_{\gradient
B}J_{0}(b\sqrt{x})\frac{q\phi}{T}F_{M}$
$\displaystyle+\frac{q}{T}ik_{\parallel}v_{\parallel}J_{0}(b\sqrt{x})\phi
F_{M}=\mathcal{C}.$ (1)

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

$\displaystyle\mathcal{C}=\left<C\right>_{\bm{R}}=\frac{1}{2\pi}\int_{0}^{2\pi}d\theta
C,$ (2)

where $C$ is the linearized ion collision operator acting on the ion perturbed
particle distribution function, and
$\left<\dots\right>_{\bm{R}}=\int_{0}^{2\pi}d\theta\dots/(2\pi)$ denotes the
gyro-average operator evaluated at constant $\bm{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 Eq. 2. The linearized GK
Boltzmann equation is closed by the GK quasi-neutrality condition that
determines the Fourier component of the electrostatic potential
$\phi(\bm{k},t)=\int d\bm{r}\phi(\bm{r},t)e^{i\bm{k}\cdot\bm{r}}$, that is

$\displaystyle\left[1+\frac{q^{2}n_{0}}{T}\left(1-\Gamma_{0}(a)\right)\right]\phi=q\left<\delta
n\right>,$ (3)

where $\left<\delta n\right>=2\pi\int dv_{\parallel}d\mu
B/mJ_{0}(b)\textsl{g}$ is the gyro-averaged ion gyrocenter perturbed density,
$a=b^{2}/2$ and $\Gamma_{0}(a)=I_{0}(a)e^{-a}$, with $I_{0}(a)$ the modified
Bessel function. The first term on the left hand-side of Eq. 3 is the
adiabatic electron response, while the second term constitutes the ion
polarization term.

We now normalize sections 2 and 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}/\Omega$ the ion sound Larmor radius. The ion-ion collision
frequency, $\nu=4Nq^{4}\ln\Lambda/[3\sqrt{\pi}\sqrt{m}T^{3/2}]$ that enters in
the ion-ion collision operator $\mathcal{C}$, is normalized to
$c_{s}/(NL_{N})$ and the gyrocenter density to $N$. The same dimensionless
units are used in the rest of the present paper. The normalized GK Boltzmann
equation, Sec. 2, then becomes

$\displaystyle\frac{\partial}{\partial t}\textsl{g}+i\left(\omega_{\gradient
B}+k_{\parallel}\sqrt{2\tau}s_{\parallel}\right)\textsl{g}-i\omega_{*}^{T}J_{0}(b\sqrt{x})\phi
F_{M}+\frac{q}{\tau}i\omega_{\gradient B}J_{0}(b\sqrt{x})\phi F_{M}$
$\displaystyle+\frac{\sqrt{2}}{\sqrt{\tau}}ik_{\parallel}s_{\parallel}J_{0}(b\sqrt{x})\phi
F_{M}=\mathcal{C},$ (4)

where $\omega_{\gradient B}=\tau\omega_{B}(x+2s_{\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,

$\left[1+\frac{q^{2}}{\tau}\left(1-\Gamma_{0}(a)\right)\right]\phi=q\left<\delta
n\right>.$ (5)

Sections 2 and 5 constitute a closed system that evolves self-consistently g
and $\phi$, in presence of the temperature and magnetic equilibrium gradients,
once the collision operator is specified (see Sec. 2.2).

###  2.1  Gyro-Moment Hierarchy Equation

To approach the solution of the GK model, given by sections 2 and 5, we use an
Hermite-Laguerre moment expansion of the ion distribution function. This
allows us to reduce the dimensionality of the linearized GK equation, Sec. 2,
by projecting it onto a velocity-space Hermite-Laguerre basis polynomials.
More precisely, we decompose the perturbed gyrocenter distribution function as
(Jorge et al., 2017, 2019a; Frei et al., 2020),

$\displaystyle\textsl{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},$
(6)

where we define the Hermite-Laguerre velocity moments of g,

$N^{pj}=\frac{1}{N}\int d\mu
dv_{\parallel}d\theta\frac{B}{m}\textsl{g}\frac{H_{p}(s_{\parallel})L_{j}(x)}{\sqrt{2^{p}p!}},$
(7)

with $N=\int d\mu dv_{\parallel}d\theta B\textsl{g}/m$ the gyrocenter density.
The Hermite-Laguerre coefficients $N^{pj}$ in Eq. 6 are referred to as the
gyro-moments. In Eq. 6, we introduce the Hermite and Laguerre polynomials,
$H_{p}$ and $L_{j}$, via their Rodrigues’ formulas (Gradshteyn & Ryzhik,
2014), that is

$\displaystyle H_{p}(x)$
$\displaystyle=(-1)^{p}e^{x^{2}}\frac{d^{p}}{dx^{p}}\left(e^{-x^{2}}\right),$
(8)

and

$\displaystyle L_{j}(x)$
$\displaystyle=\frac{e^{x}}{j!}\frac{d^{j}}{dx^{j}}\left(e^{-x}x^{j}\right),$
(9)

and their orthogonality relations, that is

$\displaystyle\int_{-\infty}^{\infty}dxH_{p}(x)H_{p^{\prime}}(x)e^{-x^{2}}$
$\displaystyle=2^{p}p!\sqrt{\pi}\delta_{p}^{p^{\prime}},$ (10)

and

$\displaystyle\int_{0}^{\infty}dxL_{j}(x)L_{j^{\prime}}(x)e^{-x}$
$\displaystyle=\delta_{j}^{j^{\prime}},$ (11)

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

$\displaystyle
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),$
(12)

where we introduce the $nth$-order kernel function

$\displaystyle\mathcal{K}_{n}=\frac{1}{n!}\left(\frac{b}{2}\right)^{2n}e^{-b^{2}/4},$
(13)

with argument $b=\sqrt{2\tau}k_{\perp}$ (Frei et al., 2020). The projections
of sections 2 and 5 onto the Hermite-Laguerre basis yields the linearized
gyro-moment hierarchy equation that states the time evolution of the gyro-
moments, $N^{pj}$, that is

$\displaystyle\frac{\partial}{\partial t}N^{pj}+$ $\displaystyle
ik_{\parallel}\sqrt{\tau}\sum_{p^{\prime}}\mathcal{H}^{\parallel
p}_{p^{\prime}}N^{p^{\prime}j}+i\frac{\tau}{q}\omega_{B}\sum_{p^{\prime}j^{\prime}}\left(\mathcal{H}_{Bp^{\prime}j^{\prime}}^{\parallel
pj}+\mathcal{H}_{Bp^{\prime}j^{\prime}}^{\perp
pj}\right)N^{p^{\prime}j^{\prime}}$
$\displaystyle=\left(\mathcal{F}^{pj}-\mathcal{F}_{B}^{\parallel
pj}\right)\phi+\mathcal{C}^{pj},$ (14)

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

$\displaystyle\mathcal{H}^{\parallel p}_{p^{\prime}}$
$\displaystyle=\sqrt{p+1}\delta_{p^{\prime}}^{p+1}+\sqrt{p}\delta_{p^{\prime}}^{p-1}$
(15a) $\displaystyle\mathcal{H}_{Bp^{\prime}j^{\prime}}^{\parallel pj}$
$\displaystyle=\delta_{j^{\prime}}^{j}\left(\sqrt{(p+1)(p+2)}\delta_{p^{\prime}}^{p+2}+\sqrt{2p+1}\delta_{p^{\prime}}^{p}+\sqrt{p(p-1)}\delta_{p^{\prime}}^{p-2}\right)$
(15b) $\displaystyle\mathcal{H}_{Bp^{\prime}j^{\prime}}^{\perp pj}$
$\displaystyle=\delta_{p^{\prime}}^{p}\left((2j+1)\delta_{j^{\prime}}^{j}-j\delta_{j^{\prime}}^{j-1}-(j+1)\delta_{j^{\prime}}^{j+1}\right),$
(15c)

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

$\displaystyle\mathcal{F}^{pj}$
$\displaystyle=i\omega_{*}\left[\mathcal{K}_{j}\delta_{p}^{0}+\eta\left(\mathcal{K}_{j}\frac{\delta_{p}^{2}}{\sqrt{2}}+\delta_{p}^{0}\left(2j\mathcal{K}_{j}-j\mathcal{K}_{j-1}-(j+1)\mathcal{K}_{j+1}\right)\right)\right],$
(16a) $\displaystyle\mathcal{F}_{B}^{\parallel pj}$
$\displaystyle=\frac{\delta_{p}^{1}}{\sqrt{\tau}}ik_{\parallel}\mathcal{K}_{j}+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]+i\omega_{B}\mathcal{K}_{j}\sqrt{2}\delta_{p}^{2}.$
(16b)

In Sec. 2.1, the Hermite-Laguerre projection of the collision operator,
$\mathcal{C}$, is defined by

$\displaystyle\mathcal{C}^{pj}=\int d\theta d\mu
dv_{\parallel}\frac{B}{m}\frac{H_{p}(s_{\parallel})L_{j}(x)}{\sqrt{2^{p}p!}}\mathcal{C}.$
(17)

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

$\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}.$
(18)

Subsections 2.1 and 18 constitute an infinite set of linear coupled fluid
equations for the gyro-moments $N^{pj}$, and depend only on the Fourier mode
wavevector $\bm{k}$ and time $t$. In order to solve numerically Sec. 2.1
together with Eq. 18, we apply a simple closure by truncation, i.e. we solve
subsections 2.1 and 18 up to $(p,j)\leq(P,J)$ and set $N^{pj}=0$ for all
$(p,j)>(P,J)$, with $0<P,J<\infty$. We remark that a high-collisional closure
of the gyro-moment hierarchy is derived in Sec. 4.1.

###  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 et al., 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 operator
(Sugama et al., 2009), the GK and DK momentum-conserving pitch-angle
scattering operator (Helander & Sigmar, 2002), the zeroth-order DK Hirshman-
Sigmar-Clarke collision operator (Hirshman & Sigmar, 1976a), and the GK and DK
Dougherty collision operator (Dougherty, 1964). The gyro-moment expansion of
the linearized GK/DK Coulomb and the GK/DK Sugama collision operators are
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 other
collision operators are reported in Appendices A, B and C, respectively. We
remark that the gyro-moment approach offers an ideal framework to study the
impact of collisions and collision operator models on ITG modes thanks to its
capability of describing the distribution function at different levels of
fidelity, as needed, and as a function of the collisionality.

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. A model of 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 $\bm{r}$ (rather than at the gyrocenter position $\bm{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
$(\bm{r},\bm{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

$\displaystyle C=C^{T}+C^{F},$ (19)

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}(\bm{r},\bm{v})$ the particle Maxwellian distribution
function and $f=f(\bm{r},\bm{v})$ a small perturbation. Conservation
properties of the collision operator constraint the test and field components.
In fact, while particle conservation is satisfied by both components, such
that

$\displaystyle\int d\bm{v}C^{T}=\int d\bm{v}C^{F}=0,$ (20)

the momentum and energy conservations require that

$\displaystyle\int d\bm{v}\bm{v}C^{T}=-\int d\bm{v}\bm{v}C^{F},$ (21)
$\displaystyle\int d\bm{v}v^{2}C^{T}=-\int d\bm{v}v^{2}C^{F},$ (22)

respectively. From the conservation properties of the linearized collision
operator $C$, Eqs. 20, 21 and 22, constraints on the coefficients
$\mathcal{C}^{pj}$ can be derived in the zero gyroradius limit, when the
gyrocenter and particle position correspond. In fact, using Eq. 17 combined
with Eqs. 20, 21 and 22 yields the relations

$\displaystyle\mathcal{C}^{00}=0,\quad\mathcal{C}^{10}=0,\quad\mathcal{C}^{20}=\sqrt{2}\mathcal{C}^{01},$
(23)

which express the conservation of gyrocenters, momentum and energy. The
constraints in Eq. 23 are satisfied by all linearized collision operator
models considered in the present work in the long wavelength limit. We have
checked that, indeed, Eq. 23 are fulfilled up to machine precision when
numerically implementing the gyro-moment expansions of the collision operators
considered in this work. Finally, we note that the relations given in Eq. 23
are used to derive the high-collisional limit of the gyro-moment hierarchy in
Sec. 4.1

##  3 Collisionless Limit

In this section, we derive the collisionless limit of the gyro-moment
hierarchy equation given in Sec. 2.1, showing that the ITG collisionless
results (see, e.g., Romanelli, 1989; Hahm & Tang, 1989; Romanelli & Briguglio,
1990; Chen et al., 1991) can be retrieved as an asymptotic limit of the gyro-
moment hierarchy. We first derive the collisionless ITG dispersion relation in
Sec. 3.1 by solving directly the GK equation, Sec. 2. Then, we obtain closed
semi-analytical expressions of the collisionless linear gyro-moment response
in Sec. 3.2. The expressions of the gyro-moments we obtain are valid at
arbitrary order in the normalized magnetic frequency, $\omega_{\gradient
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 Sec. 3.3 and thanks to the expressions of the gyro-moments
developed in Sec. 3.2, we show the equivalence between the gyro-moment and GK
dispersion relations by considering the infinite gyro-moment limit focusing on
the slab branch. Ultimately, the results of the present section allow us to
determine the scaling of the collisionless gyro-moments spectrum that will be
used for the convergence studies in Sections 7 and 8.

###  3.1 Collisionless GK ITG Dispersion Relation

We first derive the ITG dispersion relation by solving directly the GK
equation, Sec. 2. Fourier transforming Sec. 2 in time with $\mathcal{C}=0$,
such that $\partial_{t}\to-i\omega$, and introducing
$\omega=\omega_{r}+i\gamma$, with $\omega_{r}$ and $\gamma$ the real mode
frequency and the growth rate respectively, the normalized perturbed ion
gyrocenter distribution function, $\textsl{g}/\phi$, can be expressed as

$\displaystyle\frac{\textsl{g}}{\phi}=\sum_{l=1}^{3}\hat{\textsl{g}}_{l}$ (24)

with

$\displaystyle\hat{\textsl{g}}_{1}$
$\displaystyle=-\frac{q}{\tau}J_{0}(b\sqrt{x})F_{M},$ (25a)
$\displaystyle\hat{\textsl{g}}_{2}$
$\displaystyle=\frac{q}{\tau}\frac{\omega}{\omega-\omega_{\gradient
B}-z_{\parallel}s_{\parallel}}J_{0}(b\sqrt{x})F_{M},$ (25b)
$\displaystyle\hat{\textsl{g}}_{3}$
$\displaystyle=-\frac{\omega_{T}^{*}}{\omega-\omega_{\gradient
B}-z_{\parallel}s_{\parallel}}J_{0}(b\sqrt{x})J_{0}(b\sqrt{x})F_{M},$ (25c)

where we introduce $z_{\parallel}=\sqrt{2\tau}k_{\parallel}$. The $l=1$
contribution to the sum in Eq. 24 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 GK ITG dispersion relation is deduced from the linearized GK
quasineutrality condition given in Eq. 5,

$\displaystyle 1+\frac{q^{2}}{\tau}\left(1-\Gamma_{0}(a)\right)-q\hat{n}=0,$
(26)

where $\hat{n}=\sum_{l=1}^{3}\hat{n}_{l}$ is the density perturbation with
$\hat{n}_{l}=\int d\bm{v}J_{0}\hat{\textsl{g}}_{l}$ and $\hat{\textsl{g}}_{l}$
given in Eq. 25. 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_{\gradient
B}-z_{\parallel}s_{\parallel})$, in integral form (Beer & Hammett, 1996)

$\displaystyle\frac{1}{\omega-\omega_{\gradient
B}-z_{\parallel}s_{\parallel}}=-i\int_{0}^{\infty}d\tau
e^{i\tau(\omega-\omega_{\gradient B}-z_{\parallel}s_{\parallel})}.$ (27)

We remark that Eq. 27 allows us to retain the effects of magnetic drift
resonance in $\hat{N}_{l}^{pj}$ at arbitrary order in $\omega_{B}/\omega$,
since no expansion in the $\omega_{B}/\omega$ parameter is carried out to
evaluate the velocity integrals, in contrast to previous toroidal
collisionless ITG theories (Jarmén et al., 1987; Kim et al., 1993). Using Eq.
27, one derives

$\displaystyle\hat{n}_{1}$ $\displaystyle=-\frac{q}{\tau}\Gamma_{0}(a),$ (28a)
$\displaystyle\hat{n}_{2}$
$\displaystyle=-\frac{iq}{\tau}\omega\int_{0}^{\infty}d\tau
e^{i\tau\omega}I_{\perp}(\tau)I_{\parallel}(\tau),$ (28b)
$\displaystyle\hat{n}_{3}$ $\displaystyle=ik_{\perp}\int_{0}^{\infty}d\tau
e^{i\tau\omega}\left[I_{\parallel}(\tau)I_{\perp}(\tau)\right.$
$\displaystyle\left.+\eta\left(I_{\parallel
2}(\tau)I_{\perp}(\tau)+I_{\parallel}(\tau)I_{\perp
1}(\tau)-\frac{3}{2}I_{\parallel}(\tau)I_{\perp}(\tau)\right)\right],$ (28c)

where we introduce

$\displaystyle I_{\perp}(\tau)$
$\displaystyle=\frac{1}{1+i\alpha\tau}I_{0}\left(\frac{a}{1+i\alpha\tau}\right)e^{-a/(1+i\alpha\tau)},$
(29a) $\displaystyle I_{\parallel}(\tau)$
$\displaystyle=\frac{1}{\sqrt{1+2i\alpha\tau}}e^{-z_{\parallel}^{2}\tau^{2}/[4(1+2i\alpha\tau)]},$
(29b) $\displaystyle I_{\perp 1}(\tau)$
$\displaystyle=\frac{e^{-a/(1+i\alpha\tau)}}{2(1+i\alpha\tau)^{3}}$
$\displaystyle\times\left[(2(1+i\alpha\tau)-2a)I_{0}\left(\frac{a}{(1+i\alpha\tau)}\right)+2aI_{1}\left(\frac{a}{(1+i\alpha\tau)}\right)\right],$
$\displaystyle I_{\parallel 2}(\tau)$
$\displaystyle=\frac{(2(1+2i\tau\alpha)-\tau^{2}z_{\parallel}^{2})}{4(1+2i\tau\alpha)^{5/2}}e^{-z_{\parallel}^{2}\tau^{2}/(4(1+2i\tau\alpha))},$
(29c)

with $\alpha=q\omega_{B}/\tau$. Equation 26, with the definitions in Eq. 28
and the velocity integrals in Eq. 29 constitute the GK ITG collisionless
dispersion relation for the unstable modes $\gamma>0$. We note that the
exponential factors $e^{i\tau\omega}$, appearing in Eq. 28, ensure the
convergence of the integrals at $\tau\to\infty$ for $\gamma>0$. On the other
hand, the inclusion of stable modes ($\gamma<0$) in Eq. 26 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 Eq. 27 is that it yields one-dimensional integrals that can be
easily evaluated numerically.

We finally note that the GK ITG dispersion relation can be simplified in the
slab ITG (sITG) case. Neglecting $\omega_{B}$ in Eq. 25, the velocity
integrals appearing in $\hat{n}$ and given in Eq. 28, can be expressed in
terms of the plasma dispersion function. Then, the GK ITG dispersion relation,
Eq. 26, becomes

$\displaystyle 1+\frac{1}{\tau}+\frac{1}{\tau}\xi
Z(\xi)\Gamma_{0}(a)-\frac{k_{\perp}}{k_{\parallel}}\frac{1}{\sqrt{2\tau}}\left[\Gamma_{0}(a)Z(\xi)\right.$
$\displaystyle\left.+\eta\left(\Gamma_{0}(a)\xi\left(1+\xi
Z(\xi)\right)-\frac{\Gamma_{0}(a)}{2}Z(\xi)+aZ(\xi)\left(\Gamma_{1}(a)-\Gamma_{0}(a)\right)\right)\right]=0,$
(30)

where $\xi=\omega\sqrt{2\tau}/k_{\parallel}$ is the normalized mode phase
velocity,
$Z(\xi)=\int_{-\infty}^{\infty}ds_{\parallel}e^{-s_{\parallel}^{2}}/(s_{\parallel}-\xi)/\sqrt{\pi}$
is the plasma dispersion function, and $\Gamma_{1}(a)=I_{1}(a)e^{-a}$.

###  3.2 Collisionless Linear Gyro-Moment Response

We now obtain the collisionless expressions of the gyro-moments by projecting
Eq. 24 onto the Hermite-Laguerre basis. This yields

$\displaystyle\frac{N^{pj}}{\phi}=\sum_{l=1}^{3}\hat{N}_{l}^{pj},$ (31)

where we introduce

$\displaystyle\hat{N}_{l}^{pj}=2\pi\int
dv_{\parallel}d\mu\frac{B}{m}\frac{H_{p}(s_{\parallel})L_{j}(x)\hat{g}_{l}}{\sqrt{2^{p}p!}}.$
(32)

with $\hat{g}_{l}$ defined in Eq. 25. Performing the velocity integral in Eq.
32 using Eq. 27 for unstable mode with $\gamma>0$ and identities involving
hypergeometric functions and Hermite-Laguerre polynomials (described in
Appendix D), we obtain

$\displaystyle\hat{N}_{1}^{pj}$
$\displaystyle=-\frac{q}{\tau}\mathcal{K}_{j}\delta_{p}^{0},$ (33a)
$\displaystyle\hat{N}_{2}^{pj}$
$\displaystyle=-\frac{iq\omega}{\tau}\sum_{n=0}^{\infty}\frac{\mathcal{K}_{n}}{\sqrt{2^{p}p!}}\mathcal{I}_{n}^{pj},$
(33b) $\displaystyle\hat{N}_{3}^{pj}$
$\displaystyle=\frac{ik_{\perp}}{\sqrt{2^{p}p!}}\sum_{n=0}^{\infty}\mathcal{K}_{n}\left\\{\mathcal{I}_{n}^{pj}+\eta\left[(2j+1)\mathcal{I}_{n}^{pj}-j\mathcal{I}_{n}^{pj-1}-(j+1)\mathcal{I}_{n}^{pj+1}\right.\right.$
$\displaystyle\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\\}.$
(33c)

where $\mathcal{I}_{n}^{pj}$ is given in Appendix D as a finite sum of terms
that involve integrals of velocity-independent hypergeometric functions with
complex argument. Generalized plasma dispersion relations can also be used to
express the collisionless gyro-moments $\hat{N}_{l}^{pj}$, given in Eq. 33,
extending the formulas in Appendix D to the lower half complex plane of
$\omega$.

In the case of the slab ITG branch, the velocity integrals containing the
resonant term, $1/(\omega-z_{\parallel}s_{\parallel})$, being $\omega_{B}=0$,
can be performed without using Eq. 27. In fact, in this case, the
collisionless gyro-moments becomes

$\displaystyle\frac{N^{pj}}{\phi}$
$\displaystyle=-\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.$
$\displaystyle\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.$
$\displaystyle\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),$
(34)

where $Z^{(p)}(\xi)=d^{p}Z(\xi)/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}ds_{\parallel}H_{p}(s_{\parallel})e^{-s_{\parallel}^{2}}/(s_{\parallel}-\xi)/\sqrt{\pi}$.

As an aside, we also notie that Eq. 31 can be related to the kinetic response
of the monomial velocity-space moments, $M_{j,k}=v_{T}^{j+2k}\int
d\bm{v}\textsl{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

$\displaystyle 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!}}.$
(35)

We remark that Eq. 35 provide a direct link between the gyro-moment hierarchy
equation described here and the moment equations derived in Beer & Hammett
(1996) and associated closures expressions.

###  3.3 Gyro-Moment Dispersion Relation and Infinite Gyro-Moment Limit

The ITG dispersion relation expressed in terms of $N^{pj}$ is obtained from
the GK quasineutrality condition Eq. 18, yielding

$\displaystyle
1+\frac{q^{2}}{\tau}\left(1-\sum_{j=0}^{\infty}\mathcal{K}_{j}^{2}\right)-\sum_{j=0}^{\infty}\mathcal{K}_{j}\sum_{l=1}^{3}\hat{N}^{0j}_{l}=0,$
(36)

where $\hat{N}_{l}^{0j}$ is given by Eq. 33 with $p=0$.

Focusing on the sITG branch, we now demonstrate that, in the $P\to\infty$ and
$J\to\infty$ limit, the ITG gyro-moment dispersion relation, Eq. 36, is
equivalent to the GK ITG dispersion relation given in Sec. 3.1. In order to
solve the gyro-moment dispersion relation Eq. 36, we use the expressions for
the collisionless gyro-moments, $N^{pj}/\phi$ given in Sec. 3.2 in the sITG
case. Setting $p=0$ in Sec. 3.2 and plugging the resulting expression for
$N^{p0}/\phi$ into Eq. 36 yields

$\displaystyle 1+\frac{1}{\tau}+\frac{1}{\tau}\xi
Z(\xi)\hat{\Gamma_{0}}(a)-\frac{k_{\perp}}{k_{\parallel}}\frac{1}{\sqrt{2\tau}}\left[\hat{\Gamma_{0}}(a,J)Z(\xi)+\eta\left(\hat{\Gamma_{0}}(a,J)\xi\left(1+\xi
Z(\xi)\right)\right.\right.$
$\displaystyle\left.\left.-\frac{\hat{\Gamma_{0}}(a,J)}{2}Z(\xi)+Z(\xi)\hat{\Pi}_{1}(a,J)\right)\right]=0,$
(37)

where we introduce the quantities

$\displaystyle\hat{\Gamma_{0}}(a,J)$
$\displaystyle=\sum_{j=0}^{J}\mathcal{K}_{j}^{2},$ (38a)
$\displaystyle\hat{\Pi}_{1}(a,J)$
$\displaystyle=\sum_{j=0}^{J}\mathcal{K}_{j}\left(2j\mathcal{K}_{j}-j\mathcal{K}_{j-1}-(j+1)\mathcal{K}_{j+1}\right).$
(38b)

which are associated with FLR effects. By comparing Sec. 3.3 with Sec. 3.1,
one observes that they are equivalent in the infinite moment limits provided
that $\hat{\Gamma_{0}}(a,\infty)=\Gamma_{0}$ and
$\hat{\Pi}_{1}(a,\infty)=a(\Gamma_{1}-\Gamma_{0})$. This can be verified
explicitly by considering the $J\to\infty$ limit in Eq. 38 and using the
definition of the kernel function $\mathcal{K}_{j}$ in Eq. 13. In fact, we
derive

$\displaystyle\hat{\Gamma_{0}}(a,\infty)$
$\displaystyle=\sum_{j=0}^{\infty}\frac{1}{(j!)^{2}}\left(\frac{b}{2}\right)^{4j}e^{-b^{2}/2}$
$\displaystyle=e^{-a}I_{0}\left(a\right)=\Gamma_{0}(a),$ (39)

and

$\displaystyle\hat{\Pi}_{1}(a,\infty)$
$\displaystyle=2e^{-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]$
$\displaystyle=ae^{-a}\left(I_{1}(a)-I_{0}(a)\right)=a\left(\Gamma_{1}(a)-\Gamma_{0}(a)\right).$
(40)

Figure 1 displays the quantities, $\hat{\Gamma_{0}}(a,J)$ and
$\hat{\Pi}_{1}(a,J)$ defined in Eq. 38, for different values of $J$, as well
as their $J\to\infty$ limits. As observed, $\hat{\Gamma_{0}}(a,J)$ and
$\hat{\Pi}_{1}(a,J)$ converges to $\Gamma_{0}$ and $a(\Gamma_{1}-\Gamma_{0})$
when $J\to\infty$. We notice that convergence is faster for low values of
$a=b^{2}/2$, i.e. for long wavelength mode. This is in agreement with the fact
that the truncation of the sums in Sec. 3.3 and Sec. 3.3 to include only
$j\leq J$ terms yields an error $O(b^{4J})$.

Figure 1: $\hat{\Gamma_{0}}(a,J)$ and $\hat{\Pi}_{1}(a,J)$ as a function of
$a=b^{2}/2$ for increasing values of $J$. The solid black lines are their
asymptotic limits, $\Gamma_{0}$ and $-a(\Gamma_{1}-\Gamma_{0})$, respectively.

The impact on the ITG mode due to approximating $\hat{\Gamma_{0}}(a,\infty)$
and $\hat{\Pi}_{1}(a,\infty)$ with $\hat{\Gamma_{0}}(a,J)$ and
$\hat{\Pi}_{1}(a,J)$ can be illustrated by solving the ITG dispersion
relation, Sec. 3.3, for the mode complex frequency
$\omega=\omega_{r}+i\gamma$, as a function of $J$, and compare its prediction
with the one from the GK ITG dispersion relation, given in Sec. 3.1. Figure 2
shows (left panel) the maximum growth rate, $\gamma$, and (right panel) the
corresponding real frequency, $\omega_{r}$, solution of the gyro-moment sITG
dispersion relation, given in Sec. 3.3, plotted as a function of the
perpendicular wavenumber $k_{\perp}$, when $\eta=3$ and $k_{\parallel}=0.1$,
for increasing values of $J$. For comparison, the collisionless GK solution,
obtained from Sec. 3.1, is shown. Both $\gamma$ and $\omega_{r}$, obtained by
the dispersion relation expressed in terms of gyro-moments, Sec. 3.3, approach
the one predicted by Sec. 3.1, as $J\to\infty$. Convergence is achieved with
small $J$ at long wavelength (e.g., the $k_{\perp}\simeq 0.6$ case is well
reproduced with $J\simeq 5$). In general and as a rule of thumb, $J\gtrsim
k_{\perp}^{2}$ must be retained in order to retrieve accurately the GK
collisionless prediction. We remark that, if $J$ is not sufficiently large,
the truncation of the quantities, Eq. 38, 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 Sec. 8.

Figure 2: Collisionless sITG growth rate, $\gamma$, and frequency,
$\omega_{r}$, as a function of the perpendicular wavenumber $k_{\perp}$ for
different values of $J$. The color lines are the solution of the gyro-moment
hierarchy dispersion relation, given in Sec. 3.3, while the black lines are
solutions of Sec. 3.1, 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 gyro-moment hierarchy,
Sec. 2.1, focusing on two reduced models that retain only a finite number of
gyro-moments. In the high-collisionality limit, high-order gyro-moments are
damped by the presence of collisions and the ITG mode can be accurately
described by considering only the evolution of the lowest order gyro-moments.
Indeed, the perturbed gyro-center distribution function is expected not to
significantly deviate from a perturbed Maxwellian at high collisionalities. We
consider, first, a reduced model based on $6$ gyro-moments ($6$GM), obtained
by using a closure by truncation of the gyro-moment hierarchy. Second, we
consider a reduced model based on $4$ gyro-moments ($4$GM), rigorously derived
as an asymptotic limit of the $6$GM using the Chapman-Enskog closure to
express the highest-order gyro-moments as a function of the lower order ones
in the limit of $\epsilon\sim k_{\parallel}\lambda_{mfp}\ll 1$. While the
Chapman-Enskog procedure allows us to reduce the number of gyro-moments,
closures for the FLR terms associated with the $\mathcal{K}_{n}$ kernel
functions and inherent to gyro-fluid models (Brizard, 1992; Waltz et al.,
1997; Dorland & Hammett, 1993; Beer & Hammett, 1996; Madsen, 2013a; Held et
al., 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 Eq. 13. Appendix E is dedicated to explore possible FLR closures by
using a Padé approximation technique applied to the $4$GM 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 to avoid 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 Sec. 4.1, we derive the $6$GM
and $4$GM reduced models. The high-collisional limits are also verified
against the GK code GENE. In Sec. 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 branch of the ITG mode.

###  4.1 $6$GM and $4$GM Reduced Models

We start by deriving the 6GM model. We introduce the lowest order normalized
fluid gyro-moments, i.e. the gyrocenter 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 gyro-moments are obtained
by setting $(p,j)=(0,0)$, $(1,0)$, $(2,0)$, $(0,1)$, $(3,0)$, and $(1,1)$ in
Sec. 2.1, respectively, and by neglecting all higher order gyro-moments. This
yields

$\displaystyle\frac{\partial}{\partial t}N$
$\displaystyle+ik_{\parallel}\sqrt{2\tau}V+i\tau\omega_{B}\left(T_{\parallel}-T_{\perp}\right)+i\Gamma_{1}\phi=0$
(41a) $\displaystyle\sqrt{2}\frac{\partial}{\partial t}V$
$\displaystyle+\frac{ik_{\parallel}}{\sqrt{\tau}}\mathcal{K}_{0}\phi+ik_{\parallel}\sqrt{\tau}T_{\parallel}+i\tau\omega_{B}\left(4\sqrt{2}V-Q_{\perp}+\sqrt{6}Q_{\parallel}\right)=0,$
(41b) $\displaystyle\frac{\partial}{\partial t}T_{\parallel}$
$\displaystyle+ik_{\parallel}\sqrt{6\tau}Q_{\parallel}+i\tau\omega_{B}\left(5T_{\parallel}-4N+T_{\perp}\right)$
$\displaystyle+i(\Gamma_{2}-\Gamma_{1})\phi=-\nu\alpha\left(T_{\parallel}-T_{\perp}\right),$
(41c) $\displaystyle\frac{\partial}{\partial t}T_{\perp}$
$\displaystyle+ik_{\parallel}\sqrt{\tau}(\sqrt{2}V+Q_{\perp})-i\tau\omega_{B}\left(3N-T_{\parallel}-3T_{\perp}\right)$
$\displaystyle-i(\Gamma_{3}-\Gamma_{1})\phi=\nu\frac{\alpha}{2}\left(T_{\parallel}-T_{\perp}\right),$
(41d) $\displaystyle\frac{\partial}{\partial t}Q_{\parallel}$
$\displaystyle+ik_{\parallel}\sqrt{\frac{3\tau}{2}}(T_{\parallel}-N)+i\tau\omega_{B}\left(2\sqrt{3}V+8Q_{\parallel}\right)=-\nu\left(\frac{8}{5}\sqrt{\frac{2}{\pi}}Q_{\parallel}+\frac{8}{5\sqrt{3\pi}}Q_{\perp}\right),$
(41e) $\displaystyle\frac{\partial}{\partial t}Q_{\perp}$
$\displaystyle+\frac{ik_{\parallel}}{\sqrt{\tau}}\mathcal{K}_{1}\phi+ik_{\parallel}\sqrt{\tau}(N-T_{\perp})$
(41f)
$\displaystyle-i\tau\omega_{B}\left(\sqrt{2}V-6Q_{\perp}\right)=-\nu\left(\frac{28}{15}\sqrt{\frac{2}{\pi}}Q_{\perp}+\frac{8}{5\sqrt{3\pi}}Q_{\parallel}\right).$
(41g)

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

$\displaystyle\Gamma_{1}$
$\displaystyle=\mathcal{K}_{0}\left(2\omega_{B}-\omega_{*}\right)+\mathcal{K}_{1}\left(\omega_{*}\eta-\omega_{B}\right),$
(42a) $\displaystyle\Gamma_{2}$
$\displaystyle=\left(2\omega_{B}-\omega_{*}\eta\right)\mathcal{K}_{0},$ (42b)
$\displaystyle\Gamma_{3}$
$\displaystyle=\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),$
(42c)

Additionally, the electrostatic potential, $\phi$, is determined by the self-
consistent quasi-neutrality condition obtained from Eq. 18, i.e.

$\displaystyle\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},$
(43)

where the infinite sum appearing in the polarization term in Eq. 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. Equation 41 constitutes a closed set of fluid-like equations that we
refer to as the $6$GM model.

To derive the $4$GM model, further reducing the number of gyro-moments
appearing in the $6$GM model, we apply the Chapman-Enskog asymptotic closure
scheme (Chapman & Cowling, 1941; Jorge et al., 2017) to the $4$GM model. We
introduce the dimensionless small parameter
$\epsilon\sim\lambda_{mfp}k_{\parallel}\sim\sqrt{\tau}k_{\parallel}/\nu\ll 1$,
which quantifies the importance of the non-Mawellian part of g. This allows us
to express the heat fluxes $Q_{\parallel}$ and $Q_{\perp}$ as a function of
the lowest order gyro-moments, 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 neglecting 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, Eqs. 41e and 41g can be solved at the leading order in $\epsilon$
yielding

$\displaystyle
Q_{\parallel}=ik_{\parallel}\frac{\chi_{\perp}^{\parallel}}{\tau}\mathcal{K}_{1}\phi+ik_{\parallel}\chi_{\perp}^{\parallel}(N-T_{\perp})-ik_{\parallel}\chi_{\parallel}^{\parallel}(T_{\parallel}-N),$
(44a) $\displaystyle
Q_{\perp}=-ik_{\parallel}\frac{\chi_{\perp}^{\perp}}{\tau}\mathcal{K}_{1}\phi-
ik_{\parallel}\chi_{\perp}^{\perp}(N-T_{\perp})+ik_{\parallel}\chi^{\perp}_{\parallel}(T_{\parallel}-N).$
(44b)

where we introduce the thermal conductivities,

$\displaystyle\chi_{\perp}^{\parallel}=\frac{5}{16}\sqrt{\frac{\pi}{3}}\frac{\sqrt{\tau}}{\nu},\quad\chi_{\parallel}^{\parallel}=\frac{35}{32}\sqrt{\frac{2\pi}{3}}\frac{\sqrt{\tau}}{\nu},\quad\chi_{\perp}^{\perp}=\frac{5\sqrt{2\pi}}{16}\frac{\sqrt{\tau}}{\nu},\quad\chi_{\parallel}^{\perp}=\frac{5\sqrt{\pi}}{16}\frac{\sqrt{\tau}}{\nu}.$
(45)

We remark that $\chi_{\parallel}^{\parallel}>\chi_{\perp}^{\parallel}$ and
that $\chi_{\perp}^{\perp}>\chi_{\parallel}^{\perp}$. Equation 41, with the
closure relations expressed in Eq. 44, yield the $4$GM for the lowest-order
gyro-moments fluid quantities $N$, $V$, $T_{\parallel}$ and $T_{\perp}$,
providing the rigorous $\nu\gg 1$ asymptotic limit of the gyro-moment
hierarchy. We remark that an analytical procedure to close the gyro-moment
hierarchy similar to the one considered here for the ITG mode, can be applied
in more the general cases that include variation of equilibrium quantities
along the magnetic field line and electron dynamics.

To study the limits of validity of the $6$GM and $4$GM models, we compare in
Fig. 3 the estimate of the ITG growth rate $\gamma$ and real frequency
$\omega_{r}$ with the results of the gyro-moment hierarchy using the DK
Coulomb collision operator (Frei et al., 2021) and the GK ITG collisionless
dispersion relation, given in Eq. 26. 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)$ (convergence studies are reported in
Sec. 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 gyro-moment
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 gyro-moment ITG growth rates
and real frequencies agree well with the high-collisional ($\nu\gg 1$)
predictions obtained from the $6$GM and $4$GM models, in addition to being in
agreement between them in the collisionless limit, i.e. $\nu\ll 1$. Overall,
the $6$GM provide a better approximation to the low-collisionality solution
than the $4$GM, while they both agree in the high-collisional limits. It is
remarkable that the $6$GM 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 $4$GM 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 Eq. 45 produces
arbitrarily large heat fluxes. The good agreement observed in Fig. 3
illustrates the multi-fidelity nature of the gyro-moment approach that allows
the derivation of reduced models able to address both the collisionless and
high-collisional limits of the GK equation.

Figure 3: (Left) slab ITG growth rate $\gamma$ and (right) real frequency
$\omega_{r}$ as a function of collisionality $\nu$ obtained using the gyro-
moment 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 $4$GM and $6$GM high collisional limits (red dashed-dotted
and red dotted lines, respectively). The solution of the dispersion relation
in the case of $Q_{\parallel,\perp}=0$, given by Eq. 46, is also plotted for
comparison by the solid thin black lines. Here, we consider
$k_{\parallel}=0.1$ and $k_{\perp}=0.5$.

###  4.2 High-Collisional Dispersion Relation

We now derive an algebraic dispersion relation using the $4$GM 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 Eq. 45), we assume that the heat fluxes become negligible when
$\nu\gg 1$, i.e. $Q_{\parallel}\simeq 0$ and $Q_{\perp}\simeq 0$ from Eq. 44.
Hence, a dispersion relation can be derived from the $4$GM equations, Eq. 41,
yielding

$\displaystyle a\omega^{4}+b\omega^{3}+c\omega^{2}+d\omega+e=0,$ (46)

where

$\displaystyle a$ $\displaystyle=\Gamma_{\phi},$ (47a) $\displaystyle b$
$\displaystyle=\frac{3}{2}i\alpha\nu\Gamma_{\phi}-\Gamma_{1}\mathcal{K}_{0}-\Gamma_{3}\mathcal{K}_{1}-16\Gamma_{\phi}\tau\omega_{B},$
(47b) $\displaystyle c$
$\displaystyle=\frac{i\alpha\nu}{2}\left(\mathcal{K}_{1}\left(\Gamma_{2}-2\Gamma_{3}\right)-32\tau\Gamma_{\phi}\omega_{B}\Gamma_{1}\mathcal{K}_{0}\right)$
$\displaystyle-3k_{\parallel}^{2}\tau\Gamma_{\phi}-k_{\parallel}^{2}\mathcal{K}_{0}^{2}+\tau\omega_{B}\left(\mathcal{K}_{0}\left(14\Gamma_{1}-\Gamma_{2}+\Gamma_{3}\right)+\mathcal{K}_{1}\left(\Gamma_{1}+12\Gamma_{3}\right)\right),$
(47c) $\displaystyle d$
$\displaystyle=i\alpha\nu\left[k_{\parallel}^{2}\left(\mathcal{K}_{0}\mathcal{K}_{1}-\frac{5}{2}\tau\Gamma_{\phi}-\frac{3}{2}\mathcal{K}_{0}^{2}\right)+\tau\omega_{B}\left(\Gamma_{1}\left(13\mathcal{K}_{0}+2\mathcal{K}_{1}\right)-(\Gamma_{2}-2\Gamma_{3})(\mathcal{K}_{0}+3\mathcal{K}_{1})\right)\right]$
$\displaystyle+\tau
k_{\parallel}^{2}\left[\omega_{B}\left(18\tau\Gamma_{\phi}+\mathcal{K}_{0}\left(8\mathcal{K}_{0}+\mathcal{K}_{1}\right)\right)+\left(2\Gamma_{1}-\Gamma_{2}\right)\mathcal{K}_{0}+3\Gamma_{3}\mathcal{K}_{1}\right],$
(47d) $\displaystyle e$
$\displaystyle=i\alpha\nu\left[5k_{\parallel}^{2}\tau\omega_{B}\left(\tau\Gamma_{\phi}+\mathcal{K}_{0}^{2}\right)+k_{\parallel}^{2}\tau\left(\Gamma_{1}-\frac{\Gamma_{2}}{2}+\Gamma_{3}\right)\left(\mathcal{K}_{0}+\mathcal{K}_{1}\right]\right)$
$\displaystyle+k_{\parallel}^{2}\tau^{2}\omega_{B}\left(\mathcal{K}_{1}\left(\Gamma_{2}-6\Gamma_{3}-2\Gamma_{1}\right)-2\mathcal{K}_{0}\left(4\Gamma_{1}-2\Gamma_{2}+\Gamma_{3}\right)\right),$
(47e)

and where terms proportional to $\omega_{B}^{2}$ are neglected (being
$\omega_{B}\ll 1$). We also define
$\Gamma_{\phi}=1+q^{2}(1-\mathcal{K}_{0}^{2}-\mathcal{K}_{1}^{2})/\tau$. In
Eq. 46, collisional effects are represented through the terms proportional to
$i\alpha\nu$. The solution of Eq. 46 is plotted in Fig. 3 and compared with
the $6$GM and $4$GM models, showing a good agreement when $\nu\gg 1$, such
that $Q_{\parallel}=Q_{\perp}\simeq 0$.

We now consider the slab and toroidal limits of Eq. 46 separately in details.
An estimate the slab ITG peak can be derived from Eq. 46 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
$\Gamma_{\phi}\simeq 1$), we assume that the slab ITG mode is far from
marginal stability limit $\eta\sim\nu\gg 1$. Then, the dispersion relation in
Eq. 46 reduces to

$\displaystyle\hat{\omega}_{*}^{4}+\frac{3}{2}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}i\alpha\hat{\nu}_{*}\left(1+\frac{5}{2}\tau\right)\right]\hat{\omega}_{*}+\frac{3}{2}i\alpha\hat{k}_{\parallel*}\hat{\nu}_{*}\tau=0.$
(48)

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 Eq. 48, the maximum growth rate is expected to be
proportional to $\gamma\simeq g(\tau,\hat{\nu}_{*})\eta\omega_{*}$ and occur
at $k_{\parallel}\simeq f(\tau,\hat{\nu}_{*})\eta\omega_{*}$. In Fig. 4, we
show the maximum ITG growth rate, $\hat{\gamma}_{*}=\gamma/(\eta\omega_{*})$
(i.e., the imaginary part of $\hat{\omega}_{*}$) obtained by solving Eq. 48,
as a function of $\hat{k}_{\parallel*}$ and $\hat{\nu}_{*}$. The effects 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}_{*}$.

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
$4$GM model (and more generally the gyro-moment hierarchy equation in Sec.
2.1) 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 Eq. 48 enters in the dispersion relation, given in Eq. 48,
through terms proportional to $\nu\alpha(T_{\parallel}-T_{\perp})$ which arise
from the gyro-moment $\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}+2T_{\perp})/3-N$.
Using Eqs. 41c and 41d and the continuity equation Eq. 41a, this can be
expressed as

$\displaystyle\frac{\partial}{\partial
t}T-ik_{\parallel}\frac{\sqrt{2\tau}}{3}V+\frac{i}{3}\left(\Gamma_{2}-2\Gamma_{3}-2\Gamma_{1}\right)\phi=0.$
(49)

From Eq. 49 with Eq. 41a and Eq. 41b, and considering the limits used to
derive Eq. 48, we obtain the estimate of the maximum sITG growth rate,

$\displaystyle\hat{\omega}_{*}^{3}-\hat{k}_{\parallel*}^{2}\left(1+\frac{2}{3}\tau\right)\hat{\omega}_{*}+\hat{k}_{\parallel*}^{2}\tau=0.$
(50)

From Eq. 50, 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 Eq. 50 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 $4$GM model.

Figure 4: Estimates of the normalized slab ITG peak, $\hat{\gamma}_{*}$,
obtained from Eq. 48, for increasing values of $\hat{\nu}_{*}=0.015$ (from
left to right), as a function of $\hat{k}_{\parallel*}$ and temperature ratio
$\tau$. The trajectory of the maximum growth rate is indicated by the dashed
black line. The growth rates $\hat{\gamma}_{*}$ are normalized to their
maximal values.

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 Eq. 46 indicative of its
properties. Since, in contrast to the sITG, the toroidal ITG subsists at
$k_{\parallel}=0$, we set $k_{\parallel}=0$ removing the coupling with sound
waves, and 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

$\displaystyle\omega^{3}+b_{t}\omega^{2}+c_{t}\omega+d_{t}=0,$ (51)

where

$\displaystyle b_{t}$
$\displaystyle=\frac{3}{2}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],$
(52a) $\displaystyle c_{t}$
$\displaystyle=\frac{3}{2}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]$
$\displaystyle+\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],$
(52b) $\displaystyle d_{t}$
$\displaystyle=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].$
(52c)

An estimate of the value of $k_{\perp}$ that yields the largest growing mode
can be derived from Eq. 51 by neglecting terms proportional to $\nu$, for
simplicity. Then, the dispersion relation in LABEL:{eq:cubictITG} 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 Eq. 51, 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 temperature ratio $\tau$ obtained from Eq. 51, for
a low (left) and a high (right) value of collisionality. The maximum growth
rate is indicated by the dashed black line. 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

Figure 6: ITG linear growth rate, $\gamma$, in the $(k_{\perp},k_{\parallel})$
parameter space, for increasing collisionality, $\nu$ (from left to right) and
magnetic gradient strength $R_{B}$ (from top to bottom) at a temperature
gradient strength of $\eta=3$. The stability boundary of the gyro-moment
hierarchy (dotted white line) and $6$GM and $4$GM (dotted and dot-dashed
lines, respectively) are shown for comparison. Here, we fix $(P,J)=(18,6)$.

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 reference value, we consider $k_{\parallel}\sim 0.1$. This
is because of the strong ballooning character of the ITG mode in boundary
region such that $k_{\parallel}$ is of the order of $L_{N}/qR_{0}$, and
typical H-mode edge values are $q\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 Figs. 4 and 5 for the slab
and toroidal cases, respectively. For the present numerical calculations, we
fix $(P,J)=(18,6)$. As we show in Sec. 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 left
to right) and magnetic gradient $R_{B}$ (slab at the top and toroidal at
bottom). For comparison purposes, we show the stability boundaries obtained
numerically from the gyro-moment hierarchy (dotted white) and from the $6$GM
and $4$GM 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 gyro-moment hierarchy stability boundary follows closely
the ones predicted by the $6$GM and $4$GM 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\gtrsim 0.5$ according to the full gyro-moment
calculation, the $6$GM and $4$GM stability boundaries extend to
$k_{\perp}\gtrsim 1$. This highlights the importance of FLR terms in the
Coulomb collision operator model since the $6$GM and $4$GM express the $\nu\gg
1$ limits of the gyro-moment hierarchy when collisions are modelled with the
DK Coulomb (see Sec. 6.2), i.e. neglecting FLR terms in the collision
operator. 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 reported at longer wavelength in the parallel
direction, compared to 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 about twice as
large as the slab one at all collisionalities.

Figure 7: ITG growth rate maximized over $k_{\parallel}$ and $k_{\perp}$,
denoted by $\max_{k_{\parallel},k_{\perp}}(\gamma)$, a a function of the
collisionality $\nu$ and normalized temperature gradient $\eta$ for the slab
(left) and toroidal (right) ITG modes. The stability boundaries predicted by
the gyro-moment approach using the GK Coulomb collision operator (which
separates the unstable and stable modes) and the $6$GM and $4$GM models are
plotted by the white dotted, red dotted and red dotted-dashed lines,
respectively. The collisionless temperature gradient critical value is
retrieved when $\nu\ll 1$.

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 Fig. 7 for the slab (left) and toroidal (right) ITG modes
and compared with the predictions of the $6$GM and $4$GM models. The slab and
toroidal ITG mode are destabilized by increasing $\eta$ at all collisionality
and are damped (and/or suppressed) by FLR collisional effects. We remark that
the gyro-moment stability boundary matches the collisionless one when $\nu\ll
1$. 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\simeq 1$
when $\nu\lesssim 1$ to $\eta\gtrsim 3$ when $\nu\gtrsim 1$ with an abrupt
transition) in order to overcome the strong FLR collisional stabilization
effects that suppress all unstable modes. 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 Fig. 6 which shows a peak near $k_{\perp}\simeq 0.3$ when
$\nu=5$ in the toroidal case (bottom row), while the mode is completely
suppressed in the slab case (top row). Finally, we notice that the $6$GM and
$4$GM 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 $4$GM model because of
kinetic effects associated with the parallel streaming.

##  6 Comparisons Between Collision Operator Models

In this section, we investigate the differences between the collision operator
models introduced in Sec. 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 Sec. 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
to the GK Coulomb collision operator, decrease with temperature gradient
strength at all collisionalities. One of the main result 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 Sec. 6.2
by investigating a short wavelength ITG branch (SWITG) typically peaking near
$k_{\perp}\sim 1.5$ that can be destabilized if collisional FLR terms are
neglected. The present results are particular important since DK collision
operators are often used in the continuum GK codes to simulate the plasma
boundary. In this section, we truncate the gyro-moment hierarchy, Sec. 2.1, 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 Sec. 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 Hirshman-Sigmar-Clarke (HSC) collision operator (Hirshman
& Sigmar, 1976a) and the GK/DK Dougherty collision operator (Dougherty, 1964).
The gyro-moment expansion of the Sugama collsion operator are detailed in Frei
et al. (2021), while the gyro-moment expansion of the other collision
operators are reported in Appendices A, B and C.

Figure 8 shows the slab (top) and the toroidal (bottom) ITG linear growth rate
at a relative temperature gradient strength of $\eta=3$ from low (left) to
high (right) collisionality, as a function of the perpendicular wavenumber
$k_{\perp}$. The DK collision operators are displayed by the colored solid
lines, while their GK operators by the colored square markers. The
collisionless and high-collisional limits are plotted for comparison. Is is
remarkable that the ITG growth rate is sensitive to both the collision
operator models and to the presence of FLR related terms. The GK Dougherty
collision operator produces the results that most differ with respect to the
other GK operator models, and the GK pitch-angle operator yields a systematic
overestimation of the ITG growth rate, an effect of the absence of energy
diffusion in the latter operator. On the other hand, the DK collision
operators lead to a larger ITG growth rate than the GK operators, and can
eventually destabilize a SWITG peaking near $k_{\perp}\simeq 1.5$. The
presence of this SWITG can be inferred in Fig. 8 from the increase of the ITG
growth rate occurring for $k_{\perp}\gtrsim 1$ when $\nu\gtrsim 0.1$ (see Sec.
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 (Barnes et al., 2009; Ricci et al., 2006),
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 $6$GM and $4$GM 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 8: Comparisons between DK (solid lines) and GK (square markers)
collision operator models in the case of the slab (top) and toroidal (bottom)
ITG growth rate as a function of the perpendicular wavenumber $k_{\perp}$,
from the low (left) to high (right) collisionality regime. The collisionless
(dashed black) and $6$GM and $4$GM (red dotted and dashed-dotted lines,
respectively) growth rates are plotted for comparisons. Here, $\eta=3$,
$k_{\parallel}=0.1$. Figure 9: Toroidal ITG eigenvalue spectrum when $\nu=0.1$
and $k_{\perp}=0.6$ using the DK (left) and GK (right) collision operator
models. Here, $R_{B}=0.1$ and $\eta=3$.

The deviations of the collision operator models can be quantified by computing
their relative difference, $|\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. The results are shown in
Fig. 10 where we plot $|\gamma_{C}-\gamma|/\gamma_{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, it yields an error smaller than $30\%$ in all the
parameter space considered. In particular, the relative difference between the
GK Coulomb and the GK pitch-angle collision operators is similar to the one of
GK Sugama at low $\eta$ and high collisionality. For all operators, this
deviation decreases at steeper temperature gradients while it increases with
collisionality.

Figure 10: Relative deviations in the growth rate obtained using the collision
operator models, used in Fig. 8, with respect to the GK Coulomb results, as a
function of collisionality $\nu$ and temperature gradient strength $\eta$. The
colorbars are saturated at a maximal relative deviation of $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 gyro-moment approach allows the investigation of ITG eigenvalue spectrum
at finite collisionality. An example is shown in Fig. 9 for the toroidal case
when $\nu=0.1$ using the DK (left) and GK (right) 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 has been reported in the study
of electron plasma waves at arbitrary collisionality (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, yield a similar eigenvalue spectrum. The
differences observed in the ITG spectrum shown in Fig. 9 can potentially yield
large deviations in nonlinear ITG simulations where subdominant mode can
interact nonlinearly. We remark that the eigenvalues have all a positive real
frequency, $\omega_{r}>0$, corresponding to the ion diamagnetic direction.

###  6.2 Destabilization of the Short Wavelength ITG Mode

Figure 11: Slab SWITG branch driven unstable by the DK Coulomb model when
$\eta=3$ (top) and $\eta=5$ (bottom). The stability boundaries of the gyro-
moment hierarchy (dotted white line), $6$GM and $4$GM (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 colobar
is saturated at the maximum of $\gamma$ when $\nu=0.05$.

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 $6$GM and $4$GM stability boundaries in the $\nu=1$ case (see Fig. 6)
and, second, from the comparison between GK and DK collision operators (see
Fig. 8). A clear deviation between GK and DK Coulomb operators occurs at
sufficiently steep temperature gradient, $\eta\gtrsim 3$, and large
collisionality, $\nu\gtrsim 0.1$. For a more detailed analysis, we present in
Fig. 11 the ITG growth rate evaluated with the DK Coulomb in the slab case for
two different values of temperature gradient strength.

The comparison between Fig. 6 and Fig. 11 reveal 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 $6$GM and $4$GM
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 et al.,
2002; Gao et al., 2003, 2005). Additionally, Fig. 11 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 the other DK collision operators
considered in this work, while it is suppressed by all the other GK collision
operators. This questions the applicably of DK collision operator models for
the prediction of ITG driven turbulence level in edge conditions, since the
SWITG can produce significant, yet spurious, level of transport 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 Gyro-Moment Spectrum

Using the collisionless results obtained in Sec. 3, we now study the structure
of the ITG gyro-moment spectrum. In particular, we compute numerically the
collisionless gyro-moment spectrum using the semi-analytical expressions
obtained in Sec. 3.2 for the slab and toroidal branches. We find that magnetic
drift resonance effects broaden the collisionless spectrum compared to the
slab case. The semi-analytical expression of the collisionless gyro-moment,
given in Eq. 31, allows us to evaluate numerically the collisionless gyro-
moment spectrum of the normalized ion perturbed distribution function,
$\textsl{g}/\phi$. We first focus on the slab limit case and, then, on the
purely toroidal ($k_{\parallel}=0$) case. Finally, we consider the effect of
collisions by deriving an approximate expression of the cut-off of the gyro-
moment spectrum.

Considering the slab ITG 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 Eq. 26, yields the
mode complex frequency $\omega=0.1754+0.1504i$. Given $\omega$, the normalized
collisionless gyro-moments spectrum, $N^{pj}/\phi$, can be evaluated using
Sec. 3.2 (or equivalently Eq. 26 with $R_{B}=0$). This is shown in Fig. 12 as
a function of $(p,j)$ on a logarithmic scale. The collisionless gyro-moment
spectrum decreases rapidly with $j$, while the decrease of the magnitude of
the gyro-moments is slower along the $p$ direction, as highlighted on the
right panel of Fig. 12 that shows $N^{p0}/\phi$ and $N^{0j}/\phi$. The finite
amplitude of the gyro-moments $N^{pj}/\phi$ at $j\neq 0$ is due to the
presence of FLR effects, yielding terms proportional to $\mathcal{K}_{j}$ in
the collisionless gyro-moment and, to a smaller extent, to the presence of
temperature gradients, yielding terms proportional to $\eta\mathcal{K}_{j\pm
1}$ (see Sec. 3.2). The rapid decay of the gyro-moment amplitude with $j$
stems from the behaviour of the FLR kernel $\mathcal{K}_{j}$ for large $j$,
such that $\mathcal{K}_{j}\sim(b/2)^{2j}/j!$. In fact, as Fig. 12 shows, the
decrease of the $\hat{N}^{pj}$ coefficients with $j$ is slower at large $b$
and it follows, in general, $1/j!$. On the other hand, the gyro-moments of
$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 the black solid and dashed
dotted lines in Fig. 12.

Figure 12: The modulus of the normalized collisionless gyro-moments spectrum,
$N^{pj}/\phi$, in the case of slab ITG, as a function of $p$ and $j$, is shown
on the left panel (plotted on a logarithmic scale and artificially saturated
at $10^{-6}$). The spectrum $N^{p0}/\phi$ and $N^{0j}/\phi$ are shown on the
right panel (solid blue and red lines, respectively). Here,
$k_{\parallel}=0.1$, $\eta=5$ and $k_{\perp}=0.5$. Dashed lines in the right
plot correspond to the case of $k_{\perp}=2$.

We now consider the collisionless gyro-moment spectrum in the presence of
magnetic gradients. As in the slab case, we focus on $k_{\perp}=0.5$,
$k_{\parallel}=0.1$ and $\eta=5$, but we choose $R_{B}=0.5$. The solution of
the GK dispersion relation Eq. 26, for the complex mode frequency, yields
$\omega=0.7972+0.4364i$. Equation 33 is used to evaluate the definite
integrals in $\mathcal{I}_{n}^{pj}$, given in Appendix D numerically, while
the infinite sums in Eq. 33 are truncated when machine precision is reached.
Similarly to Fig. 12, we plot the toroidal collisionless gyro-moment spectrum,
$N^{pj}/\phi$, on a logarithmic scale in Fig. 13. Compared with the slab case,
the gyro-moment spectrum is significantly broader. This is attributed to the
velocity-dependence of the magnetic drift resonant coefficients, proportional
to $1/[\omega-\omega_{\gradient B}]$, that drives high-order gyro-moments. In
fact, while the high $j$ dependence of the gyro-moments 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 Appendix D), extend the gyro-moment spectrum to
higher $j$, producing a considerably slower decay than $1/j!$. We also notice
that, since $k_{\parallel}=0$, the gyro-moment spectrum vanishes for odd
values of $p$, where $\mathcal{I}_{n}^{pj}=0$ (see Appendix D). Finite
$k_{\parallel}$ yields non vanishing odd gyro-moments in $p$ due to the term
proportional to $(1+(-1)^{p-2p_{1}+1})$ in Appendix D.

Figure 13: The modulus of the normalized collisionless gyro-moments spectrum,
$N^{pj}/\phi$, in the case of toroidal ITG ($R_{B}=0.5$), as a function of $p$
and $j$, is shown on the left panel (plotted on a logarithmic scale and
artificially saturated at $10^{-6}$). The spectrum $N^{p0}/\phi$ and
$N^{0j}/\phi$ are shown on the right panel (solid blue and red lines,
respectively). Here, $k_{\perp}=0.5$, $k_{\parallel}=0$, $\eta=5$.

Despite the fact that the gyro-moment 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 gyro-moment
representation of the perturbed distribution function, $\textsl{g}/\phi$.
Given the complex frequency $\omega$, a truncated distribution function
$g^{PJ}$ is defined from the Hermite-Laguerre expansion in Eq. 6 according to

$\displaystyle\textsl{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}.$
(53)

Then, from $g^{PJ}$ and $\textsl{g}/\phi$ obtained from Eq. 24, the error norm
$e(P,J)$ can be defined by

$\displaystyle
e(P,J)=\int_{-\infty}^{\infty}ds_{\parallel}\int_{0}^{\infty}x\left|\textsl{g}^{PJ}-\frac{\textsl{g}}{\phi}\right|^{2},$
(54)

and can be computed as a function of $(P,J)$. The results are shown in Fig. 14
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 gyro-moments is necessary to achieve a similar
error, i.e. $(P,J)\gtrsim(50,20)$.

Figure 14: Error norm $e(P,J)$, Eq. 54, on a logarithmic scale for the slab
(left) and toroidal (right) ITG cases. The same parameters as in Fig. 12 and
Fig. 13 are used respectively. The colorbars are saturated artificially at
$\log_{10}e(P,J)=-6$.

We now turn to the collisional effects, which provide a natural cutoff of the
gyro-moment hierarchy. We aim to derive an ad-hoc proxy of the scaling of
$(P,J)$ at which the gyro-moment hierarchy can be truncated as a function of
the collisionality. In particular, we consider an equation for the generalized
gyro-moment 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 gyro-moment hierarchy
equation, Sec. 2.1 in the limit $p,j\gg 1$ (Zocco & Schekochihin, 2011;
Loureiro et al., 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^{pj\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.

$\displaystyle|k_{\parallel}|\sqrt{\tau}\partial_{s}G^{pj}+\frac{\tau}{q}\omega_{B}(2j+1)\partial_{j}G^{pj}+\frac{2\tau}{q}\omega_{B}\left(1+i\left(2j+1\right)\right)G^{pj}=2\nu\left(f_{p,j}+b^{\alpha}h_{p,j}-\gamma\right)G^{pj},$
(55)

with $s=\sqrt{p}$. In deriving Eq. 55, collisional effects are modelled using
the ansatz
$\mathcal{C}^{pi}\simeq\nu\left[f_{p,j}+b^{\alpha}h_{p,j}\right]N^{pj}$, where
the functions $f_{p,j}=f(p,j)$ and $h_{p,j}=h(p,j)$, and the constant $\alpha$
depend 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 Eq. 55. We solve Eq. 55 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 gyro-moment
expansion of the GK Dougherty collision operator, given in Appendix C, allows
us to solve analytically Eq. 55.

Using the gyro-moment expansion of the GK Dougherty collision operator, we
solve first Eq. 55 in the case of the slab ITG mode. With $f_{p,j}=-p-2j$,
$h_{p,j}=-1/2$ and $\alpha=2$, the solution of Eq. 55 is

$\displaystyle\left|N^{pj}\right|\sim\frac{1}{p^{1/4}}\exp\left\\{-\left(\frac{p}{p_{\gamma}}\right)^{1/2}-\left(\frac{p}{p_{\perp}}\right)^{1/2}-\left(\frac{p}{p_{D\nu}}\right)^{3/2}\right\\},$
(56)

where we introduce $p_{\gamma}=\tau/\gamma_{\parallel}^{2}$ (with
$\gamma_{\parallel}=\gamma/k_{\parallel}$),
$p_{D\nu}=(9\tau/\nu_{\parallel}^{2})^{1/3}$
($\nu_{\parallel}=\nu/k_{\parallel}$) and
$p_{\perp}=\tau/[\nu_{\parallel}(2j+a)]^{2}$. The $1/\nu_{\parallel}$ of
$p_{D\nu}$ indicates that the exponential decay of Eq. 56 is slower at high
$k_{\parallel}$ and faster at high $\nu$. At high $k_{\parallel}$, higher-
order $p$ gyro-moments are excited due to parallel resonance effects. Indeed,
in the high-collisional limit, the Chapman-Enskog closure procedure ensures
that $\nu_{\parallel}\sim\sqrt{\tau}/\epsilon\gg 1$, such that
$p_{D\nu}\sim\nu_{\parallel}^{-2/3}\sim\epsilon^{2/3}/\tau^{1/3}$ and,
therefore, that the ITG dynamics can be described by the lowest order gyro-
moments at high collisionality (see Sec. 4.1). Moreover, when $\nu\gg\gamma$
(such that $p_{D\nu}$ is small that $p_{\gamma}$) and in the DK limit, the
collisional damping provide a truncation at $P\gtrsim p_{D\nu}$. We remark
that FLR terms have a small effect on the scaling of Eq. 56 since, for ITG
modes, $p_{\perp}\gtrsim p_{D\nu}$ and the term with the exponent $3/2$ is
dominant at large $p$. Since Eq. 55 does not couple the Laguerre gyro-moment,
the evolution at large $j$ is dominated by collisions, and, thus, are damped
in the high-collisional regime.

We now derive the solution of Eq. 55 for the toroidal $k_{\parallel}=0$ ITG
case, still considering the GK Dougherty collision operator. It yields

$\displaystyle\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+a\right)\right]\right\\}.$
(57)

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 gyro-moment spectrum associated with the presence of
magnetic gradient drifts (see Fig. 13), an effect reduced at high
collisionality. Similar to the slab case, FLR terms do not affect the scaling
of Eq. 56, since the linear term dominates over the second term at large $j$
values, and $p-1\gg a\sim k_{\perp}^{2}$ for ITG modes. We remark that, in the
case of toroidal ITG with $k_{\parallel}=0$, the Hermite gyro-moment evolution
is dominated by collisional damping.

Figure 15: Normalized slab (left) and toroidal (right) gyro-moment spectra as
a function of $p$ and $j$, respectively. The solid lines are the gyro-moment
spectrum obtained using the GK Coulomb collision operator, given in Eq. 58 and
Eq. 59, and the dotted lines the ones obtained from the GK Dougherty collision
operator, given in Eq. 56 and Eq. 57. The slab and toroidal growth rate are
estimated at $\gamma\simeq 0.1$ and $\gamma\simeq 0.4$, respectively near the
ITG peaks ($k_{\perp}\simeq 0.5$).

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}^{pi}\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-Aj-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 effects 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 solve again Eq. 55 for the
case of the slab and toroidal $k_{\parallel}=0$ ITG modes yielding,
respectively,

$\displaystyle\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\\},$
(58)

with $p_{j}=\tau/\nu_{\parallel}Aj^{2}$ and
$p_{\nu}=2\sqrt{\tau}/(B\nu_{\parallel})$, and

$\displaystyle\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\\},$
(59)

The main differences between the results of the GK Coulomb collision operator,
given in Eq. 58 and Eq. 59 and the results of the GK Dougherty operator, in
Eqs. 56 and 57, are that the scaling in $p$ predicted by the GK Coulomb
collision operator is weaker in the slab case (compared to a $3/2$ exponent in
Eq. 56), 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 gyro-moment spectra obtained using the GK
Dougherty and GK Coulomb collision operators in Fig. 15 for the slab (left)
and toroidal (right) ITG modes for different values of $\nu_{\parallel}$ (at
$k_{\parallel}=0.1$) and $\nu_{B}$ (at $\nu=0.1$). The differences between the
$p$ exponent and the numerical factors between GK Dougherty (dotted lines) and
GK Coulomb (solid lines), do not significantly impact the gyro-moment spectra
that feature the same qualitative behaviour. The larger deviation appear in
the toroidal case when $\nu_{B}$ is small (large $R_{B}$ and/or small $\nu$).
We compare the predictions of the values of $P$ and $J$, at which the gyro-
moment hierarchy can be truncated deduced from Fig. 15, with the numerical
convergences results obtained in Sec. 8.

##  8 Convergence Studies

Following the theoretical analysis of the gyro-moment spectrum presented in
Sec. 7, we perform convergence studies of the ITG growth rate as a function of
collisionality, perpendicular wavenumber and magnetic gradient. Collisions are
modelled using the GK and DK Coulomb collision operator (Frei et al., 2021).
Similar convergence properties are observed with other collision operator
models considered in this work. We show that the convergence of the gyro-
moment approach improves with collisionality.

Figure 16: 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
(left) and toroidal (right), with $R_{B}=0.1$ branches. The collisionless
limit (dotted black lines) and the high collisional limits, $6$GM and $4$GM
models (red dotted and dash-dot lines, respectively) are shown for
comparisons. 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$.

First, convergence on the ITG linear growth mode using the GK and DK Coulomb
collision operators are performed numerically as function of the
collisionality $\nu$. The collisionless and high-collisional limits are
considered for comparison by solving the collisionless dispersion relation,
Eq. 26, and the $6$GM and $4$GM models. The results are shown in Fig. 16 for
the slab (left) and for the toroidal (right) ITG modes. Different number of
gyro-moments, $(P,J)$, are used. The collisionless solution is retrieved for
$(P,J)\simeq(18,6)$ with the GK (dashed) and DK (solid) Coulomb collision
operators. In the slab case, for $\nu\simeq 0.1$ and $\nu\lesssim 10^{3}$,
convergence is achieved with $P\lesssim 10$ and $P\gtrsim 18$. This is in
qualitative agreement with the estimates deduced from Fig. 15 with the
corresponding values of $\nu_{\parallel}=\nu/k_{\parallel}\simeq 1$ and
$\nu_{\parallel}\lesssim 0.01$. Also consistent with the observations made in
Sec. 7, only a few number of Laguerre gyro-moments are necessary to resolve
the slab ITG, as illustrated by the overlapping between the $(10,4)$ and
$(10,6)$ lines in Fig. 16 as the slab ITG mode is mainly driven by the
parallel dynamics which is resolved by the Hermite gyro-moments. This remains
also true for the toroidal ITG case considered in Fig. 16. 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 gyro-moment
hierarchy, with DK Coulomb, agrees well with the $6$GM and $4$GM when $\nu\gg
1$. Nevertheless, we notice that, because of finite magnetic drifts (see Sec.
3), agreements with the $6$GM 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}$.

We investigate the convergence associated with the perpendicular dynamics in
velocity space, i.e. with the Laguerre gyro-moments and associated with the
FLR effects due to the gyro-average of the electrostatic potential. The
coupling between the Laguerre gyro-moments driven by magnetic gradients are
neglected by focusing on the slab ITG branch. We scan the growth rate $\gamma$
as a function of the perpendicular wavenumber, $k_{\perp}$, for different
values of $J$ and collisionality $\nu$, while we fix $P=18$. The results are
shown in Fig. 17. Convergence occurs at lower $J$ for high collisionalities,
but in all cases the values of $J$ for convergence increases with $k_{\perp}$.
This is in agreement with the unfavorable scaling of the FLR kernel with $b$,
i.e. $\mathcal{K}_{j}\sim(b/2)^{2j}/j!$ (see Sec. 7).

Figure 17: Slab ITG growth rate $\gamma$ as a function of $k_{\perp}$ for
different values of $J$, in the low (left), intermediate (center), and high
(right) collisionality regimes, with the GK (square markers) and DK (solid
lines) Coulomb collision operators. Similar plots are obtained for the other
GK and DK operator models. Figure 18: Ratio between the ITG growth rate
$\gamma$ and the collisionless solution $\gamma_{\textrm{collisionless}}$ as a
function of $P$ with $J=10$ (left) and as a function of $J$ with $P=32$
(right) for two different values of $R_{B}$. Here, $k_{\perp}=0.5$,
$k_{\parallel}=0.1$ and $\nu=0.001$.

Finally, we consider the convergence of the gyro-moment hierarchy as a
function of the normalized magnetic gradient, $R_{B}$. As illustrated in Sec.
7, the presence of finite magnetic drifts broadens significantly the
collisionless gyro-moment spectrum compared to the slab case. However, in all
cases, the gyro-moment approach converge correctly to the collisionless
solution of Eq. 26, as shown in Fig. 18, where the ratio of the ITG growth
rate to the collisionless value is plotted as a function of $P$ and $J$. It is
remarkable that a larger number of gyro-moment is needed as $R_{B}$ increases.
Collisions can compete against this strong kinetic drive by damping higher
order gyro-moments (see Fig. 15). To investigate the convergence of the gyro-
moment approach at finite $R_{B}$ in the presence of collisions, we focus on
the toroidal ITG with $k_{\parallel}=0$, i.e. neglecting the coupling between
the Hermite gyro-moments 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 Fig. 19. At all
collisionalities, the convergence deteriorates as $R_{B}$ increases.
Nevertheless, convergence is observed with the increase of the collisionality,
showing the competition between the magnetic gradient kinetic and collisional
effects. We also remark that the $4$GM limit is retrieved at higher
collisionality when $R_{B}$ increases.

Figure 19: 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 left to right). Here, the GK (dashed lines)
and DK (solid lines) Coulomb collision operators are shown, with the
collisionless (black dotted) and high collisional $4$GM 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 gyrocenter distribution
function, which we refer to as gyro-moment approach. By projecting the GK
Boltzmann equation onto the Hermite-Laguerre basis, a gyro-moment 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 gyro-moment
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 gyro-moment spectrum compared to
the slab case. In the high-collisional limit, the gyro-moment hierarchy can be
reduced to a fluid model retaining only a finite number of gyro-moments using
the Chapman-Enskog asymptotic closure scheme. Using these analytical results,
we show that the gyro-moment 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 gyro-moments.

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 Hirshman-Sigmar-Clarke collision operator (Hirshman &
Sigmar, 1976a), and the Dougherty collision operator (Dougherty, 1964). The
gyro-moment expansions of some of these 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 gyro-moment expansion developed in
this work allows us to perform a systematic comparison between collision
operator models and assess their accuracy. The main outcome of this comparison
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 for the GK Dougherty operator. In
general, we observe that energy diffusion is important at high collisionality
and in the damping of small-scale modes.

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 gyro-moment approach to capture the
relevant physics in the kinetic and collisional regimes, with a fidelity cost
that depends on the collisionality. Finally, we remark that the present study
can be extended to more realistic conditions including the parallel direction,
electron dynamics, and electromagnetic effects.

## Acknowledgement

The authors acknowledge helpful discussions with S. Brunner, P. Donnel, and M.
Held.

This research was supported in part by the Swiss National Science Foundation,
and has been carried out within the framework of the EUROfusion Consortium and
has received funding from the Euratom research and training programme 2014 -
2018 and 2019 - 2020 under grant agreement No 633053. The views and opinions
expressed herein do not necessarily reflect those of the European Commission.

##  Appendix A Pitch-Angle Scattering Collision Operator

The linearized pitch-angle scattering operator model is the sum of the pitch-
angle scattering operator for the test component coupled and an ad hoc field
component such that the momentum conservation is satisfied. This operator
neglects energy diffusion, contrary to the Sugama and Coulomb collision
operators, but conserves particle, momentum and energy. The test component of
the GK pitch-angle scattering operator is given by (Helander & Sigmar, 2002)

$\displaystyle\mathcal{C}^{T}=\nu^{D}(v)\left\\{-\mathcal{L}^{2}h-\frac{b^{2}}{2v_{T}^{2}}(2v_{\parallel}^{2}+v_{\perp}^{2})h\right\\},$
(60)

and the field component is defined by

$\displaystyle\mathcal{C}^{F}=2F_{M}\nu^{D}(v)\left(v_{\perp}J_{1}\overline{U}_{\perp}^{D}[h]+v_{\parallel}J_{0}\overline{U}_{\parallel}^{D}[h]\right),$
(61)

where we introduce

$\displaystyle\overline{U}_{\perp}^{D}[h]=\frac{3\int
d\bm{v}v_{\perp}J_{1}\nu^{D}(v)h}{\int
d\bm{v}v^{2}\nu^{D}(v)F_{M}},\quad\overline{U}_{\parallel}^{D}[h]=\frac{3\int
d\bm{v}v_{\parallel}J_{0}\nu^{D}(v)h}{\int d\bm{v}v^{2}\nu^{D}(v)F_{M}}.$ (62)

with $h=\textsl{g}+qJ_{0}\phi F_{M}/\tau$ the nonadiabatic part of the
perturbed gyrocenter distribution function g. In Eq. 60, $\mathcal{L}^{2}$ is
the spherical angular operator defined by
$v^{2}\gradient_{\bm{v}}^{2}f=\partial_{v}(v^{2}\partial_{v}f)-\mathcal{L}^{2}f$
and
$\nu^{D}(v)=\nu\left[\erf(s)-(\erf(s)-s\erf^{\prime}(s))/2s^{2}\right]/s^{3}$
(with $s=v/v_{T}$) is the pitch-angle scattering (or deflection) velocity-
dependent frequency.

Following Frei et al. (2021), we now project Eqs. 60 and 61 onto the Hermite-
Laguerre basis in terms of the nonadiabatic moments, $n^{pj}$, defined by

$\displaystyle n^{pj}=\frac{1}{N}\int d\mu
dv_{\parallel}d\theta\frac{B}{m}h\frac{H_{p}(s_{\parallel})L_{j}(x)}{\sqrt{2^{p}p!}},$
(63)

and expressed in terms of the gyro-moments, $N^{pj}$, by
$n^{pj}=N^{pj}+q\mathcal{K}_{j}\phi\delta_{p}^{0}/\tau$. The gyro-moment
expansion of the test and field components, Eq. 60 and Eq. 61, are given by

$\displaystyle\mathcal{C}^{Tlk}=\mathcal{C}^{Dlk}+\overline{R}^{Dlk},$ (64)

and

$\displaystyle\mathcal{C}^{Flk}=\overline{R}_{\parallel}^{Dlk}+\overline{R}_{\perp}^{Dlk},$
(65)

respectively. In Eqs. 64 and 65, we introduce

$\displaystyle\mathcal{C}^{Dlk}=-\sum_{j=0}^{\infty}\sum_{p=0}^{l+2k}\sum_{h=0}^{k+\lfloor
l/2\rfloor}\sum_{g=0}^{p+2j}\sum_{t=0}^{j+\lfloor
l/2\rfloor}\sum_{d=0}^{h}\frac{2^{p}(p!)^{2}\sqrt{2^{g}g!}}{(2p)!}\frac{\left(T^{-1}\right)^{ph}_{lk}}{\sigma_{p}^{j}\sqrt{2^{l}l!}}\frac{p(p+1)}{(2p+1)}T_{pj}^{gt}L_{hd}^{p}\overline{\nu}_{*}^{Dpjd}n^{gt},$
(66)

with

$\displaystyle\overline{\nu}_{*}^{Dpjd}=\frac{4}{\sqrt{\pi}}\sum_{j_{1}=0}^{j}L_{jj_{1}}^{p}\overline{\nu}^{Dp+d+j_{1}+1},$
(67)

and

$\displaystyle\overline{R}^{Dlk}$
$\displaystyle=-\frac{b^{2}}{\sqrt{\pi}}\sum_{p,j}\sum_{g=0}^{l+2k}\sum_{h=0}^{k+\lfloor
l/2\rfloor}\sum_{r=0}^{p+2j}\sum_{s=0}^{j+\lfloor
p/2\rfloor}\sum_{s_{1}=0}^{s}\sum_{h_{1}=0}^{h}\frac{\left(T^{-1}\right)^{gh}_{lk}}{\sqrt{2^{l}l!}}\frac{\left(T^{-1}\right)^{rs}_{pj}}{\sqrt{2^{p}p!}}L^{r}_{ss_{1}}L^{g}_{hh_{1}}\mathcal{\alpha}_{\perp}^{grs_{1}h_{1}}n^{pj},$
(68a) $\displaystyle\overline{R}_{\perp}^{lk}$
$\displaystyle=\overline{U}_{\perp}[h]\frac{4}{\sqrt{\pi}}\sum_{n=0}^{\infty}\sum_{s=0}^{n+k+1}\sum_{h=0}^{s+\lfloor
l/2\rfloor}\sum_{h_{1}=0}^{h}L_{hh_{1}}^{0}\frac{d_{nks}^{1}b\mathcal{K}_{n}}{(n+1)}\frac{\left(T^{-1}\right)_{ls}^{gh}}{\sqrt{2^{l}l!}}\overline{\nu}^{Dh_{1}+1},$
(68b) $\displaystyle\overline{R}_{\parallel}^{lk}$
$\displaystyle=\overline{U}_{\parallel}[h]\frac{8}{3\sqrt{\pi}}\sum_{n=0}^{\infty}\sum_{s=0}^{n+k}\sum_{g=0}^{l+2s}\sum_{h=0}^{s+\lfloor
l/2\rfloor}\sum_{h_{1}=0}^{h}d_{nks}^{0}\mathcal{K}_{n}\frac{L_{hh_{1}}^{1}\left(T^{-1}\right)_{ls}^{gh}}{\sqrt{2^{l}l!}}\overline{\nu}^{Dh_{1}+2}[l>0\cup
s>0],$ (68c)

where $[\cdot]$ denotes the Iverson bracket. In Eq. 68, we define

$\displaystyle\mathcal{\alpha}_{\perp}^{grs_{1}h_{1}}$
$\displaystyle=\frac{2\delta_{g}^{r}}{(2g+1)}\overline{\nu}^{Dg+2+s_{1}+h_{1}}+d_{g}^{r}\overline{\nu}^{D(g+r)/2+2+s_{1}+h_{1}},$
(69)

with the numerical coefficient

$\displaystyle
d_{g}^{r}=\begin{cases}\dfrac{2(r+1)(r+2)}{(2r+1)(2r+3)(2r+5)}&\mbox{for
}g=r+2\\\ \dfrac{2(2r^{2}+2r-1)}{(2r-1)(2r+1)(2r+3)}&\hbox{for }g=r\\\
\dfrac{2r(r-1)}{(2r-3)(2r-1)(2r+1)}&\hbox{for }g=r-2,\end{cases}$ (70)

and

$\displaystyle\overline{U}_{\parallel}[h]$
$\displaystyle=\frac{1}{\overline{\nu}^{D2}}\sum_{p,j}\sum_{n=0}^{\infty}\sum_{s=0}^{j+n}\sum_{g=0}^{p+2s}\sum_{h=0}^{s+\lfloor
p/2\rfloor}\sum_{h_{1}=0}^{h}d_{njs}^{0}\mathcal{K}_{n}\frac{L_{hh_{1}}^{1}\left(T^{-1}\right)_{ps}^{gh}}{\sqrt{2^{p}p!}}\overline{\nu}^{Dh_{1}+2}n^{pj}[p>0\cup
s>1],$ (71a) $\displaystyle\overline{U}_{\perp}[h]$
$\displaystyle=\frac{3}{2\overline{\nu}^{D2}}\sum_{p,j}\sum_{n=0}^{\infty}\sum_{s=0}^{n+j+1}\sum_{h=0}^{s+\lfloor
p/2\rfloor}\sum_{h_{1}=0}^{h}L_{hh_{1}}^{0}\frac{b\mathcal{K}_{n}}{(n+1)}\frac{d_{njs}^{1}\left(T^{-1}\right)_{ps}^{0h}}{\sqrt{2^{p}p!}}\overline{\nu}^{Dh_{1}+1}n^{pj},$
(71b)

In Eqs. 68 and 71, we introduce the speed integrated deflection frequency by
$\overline{\nu}^{Dk}=\int dse^{-s^{2}}s^{2k}\nu^{D}(v)$. Using the definitions
of $\nu^{D}(v)$, we derive

$\displaystyle\overline{\nu}^{Dk}$
$\displaystyle=\frac{1}{2}\nu\left(2E^{k-2}-E^{k-3}+e^{k-2}\right),$ (72)

where the definitions of $E^{k}$ and $e^{k}$ can be found in (Frei et al.,
2021).

The DK limit of the pitch-angle scattering operator is obtained from Eqs. 64
and 65 by taking the zero gyroradius limit, $b\to 0$. Noticing that
$n^{pj}\simeq N^{pj}$, we derive the DK test and the field components given by

$\displaystyle\mathcal{C}^{Tlk}$ $\displaystyle=\mathcal{C}^{Dlk}$ (73)
$\displaystyle\mathcal{C}^{Flk}$
$\displaystyle=U_{\parallel}[\textsl{g}]\frac{8}{3\sqrt{\pi}}\sum_{g=0}^{l+2k}\sum_{h=0}^{k+\lfloor
l/2\rfloor}\sum_{h_{1}=0}^{h}\frac{L_{hh_{1}}^{1}\left(T^{-1}\right)_{lk}^{gh}}{\sqrt{2^{l}l!}}\overline{\nu}^{Dh_{1}+2}[l>0\cup
k>0],$ (74)

respectively, where $C^{Dlk}$ is given in Eq. 66 and with

$\displaystyle U_{\parallel}^{D}[\textsl{g}]$
$\displaystyle=\frac{1}{\overline{\nu}^{D2}}\sum_{p,j}\sum_{g=0}^{p+2j}\sum_{h=0}^{j+\lfloor
p/2\rfloor}\sum_{h_{1}=0}^{h}N^{pj}\frac{L_{hh_{1}}^{1}\left(T^{-1}\right)_{pj}^{gh}}{\sqrt{2^{p}p!}}\overline{\nu}^{Dh_{1}+2}[p>0\cup
k>j].$ (75)

##  Appendix B Hirshman-Sigmar-Clarke Collision Operator

We consider the gyro-moment expansion of the Hirshman-Sigmar-Clarke (HSC)
collision operator (Hirshman et al., 1976), also referred to as the zeroth-
order Hirshman-Sigmar collision operator (Hirshman & Sigmar, 1976b). The HSC
collision operator model has been widely used for neoclassical transport
calculations (Hirshman et al., 1976; Belli & Candy, 2008). This operator is
only considered in the DK limit. Additionally, it conserves particle, momentum
and energy. While the test component of the HSC operator is the pitch-angle
scattering operator, i.e.

$\displaystyle\mathcal{C}^{T}=-\nu^{D}(v)\mathcal{L}^{2}\textsl{g},$ (76)

the field component of the HSC operator is defined by

$\displaystyle\mathcal{C}^{F}=v_{\parallel}\left[\nu^{S}(v)\frac{2r[\textsl{g}]}{v_{T}^{2}}+\frac{2u[\textsl{g}]}{v^{2}}\left(\nu^{D}(v)-\nu^{S}(v)\right)\right]F_{M},$
(77)

where

$\displaystyle r[\textsl{g}]$ $\displaystyle=\dfrac{3\int
d\bm{v}v_{\parallel}\nu^{S}(v)\textsl{g}}{\int d\bm{v}s^{2}\nu^{S}(v)F_{M}},$
(78) $\displaystyle u[\textsl{g}]$ $\displaystyle=\frac{3}{4\pi}\int d\Omega
v_{\parallel}\frac{\textsl{g}}{F_{M}},$ (79)

with $\int d\Omega=\int_{0}^{2\pi}d\theta\int_{-1}^{1}d\xi$ the integral over
the velocity-space solid angle, and
$\nu^{S}(v)=2\nu(\erf(s)-s\erf^{\prime}(s))/s^{3}$.

We project the HSC collision operator, given in Eq. 77, onto the Hermite-
Laguerre basis. As the Hermite-Laguerre projection of
$-\nu^{D}(v)\mathcal{L}^{2}g$ is given by Eq. 66, we focus on the projection
of Eq. 77. Performing the Hermite-Laguerre projection yields

$\displaystyle\mathcal{C}^{Flk}=R^{Slk}_{1}+R^{Dlk},$ (80)

where we introduce

$\displaystyle R_{1}^{Slk}$
$\displaystyle=\sum_{g=0}^{l+2k}\sum_{h=0}^{k+\lfloor
l/2\rfloor}\sum_{h_{1}=0}^{h}L_{hh_{1}}^{1}\frac{\left(T^{-1}\right)^{gh}_{lk}}{\sqrt{2^{l}l!}}\frac{8}{3\sqrt{\pi}}\frac{r[\textsl{g}]}{v_{T}}\overline{\nu}^{Sh_{1}+2}[l>0\cup
k>0],$ (81a) $\displaystyle R^{SDlk}$
$\displaystyle=\frac{8}{3\sqrt{\pi}}\sum_{p,j}\sum_{g=0}^{p+2j}\sum_{h=0}^{j+\lfloor
p/2\rfloor}\sum_{s=0}^{l+2k}\sum_{t=0}^{k+\rfloor
l/2\lfloor}\sum_{h_{1}=0}^{h}\sum_{t_{1}=0}^{t}\frac{\left(T^{-1}\right)^{st}_{lk}}{\sqrt{2^{l}l!}}\frac{\left(T^{-1}\right)^{gh}_{pj}}{\sqrt{2^{p}p!}}$
$\displaystyle\times
L_{tt_{1}}^{1}L_{hh_{1}}^{1}\left(\overline{\nu}^{D2+t_{1}+h_{1}}-\overline{\nu}^{S2+t_{1}+h_{1}}\right)N^{pj}[p>0\cup
j>0][l>0\cup k>0],$ (81b)

and

$\displaystyle\frac{r[\textsl{g}]}{v_{T}}$
$\displaystyle=\frac{1}{\overline{\nu}^{S2}}\sum_{p,j}\sum_{s=0}^{p+2j}\sum_{t=0}^{j+\lfloor
p/2\rfloor}\sum_{t_{1}=0}^{t}\left(T^{-1}\right)^{st}_{pj}L_{tt_{1}}^{1}\frac{N^{pj}}{\sqrt{2^{p}p!}}\overline{\nu}^{St_{1}+2}[p>0\cup
j>0].$ (82)

In, Eqs. 81a and 82, we introduce the speed integrated frequency by
$\overline{\nu}^{Sk}=\int dse^{-s^{2}}s^{2k}\nu^{S}(v)$. Using the definitions
of $\nu^{S}(v)$, we derive

$\displaystyle\overline{\nu}^{Sk}$
$\displaystyle=2\nu\left(E^{k-2}-e^{k-1}\right).$ (83)

##  Appendix C Dougherty Collision Operator

The linearized GK Dougherty collision operator is given by the test and field
components (Dougherty, 1964)

$\displaystyle\mathcal{C}^{T}$
$\displaystyle=\nu\left[\frac{T}{m}\frac{\partial^{2}}{\partial
v_{\parallel}^{2}}h+2\frac{T}{B}\frac{\partial}{\partial\mu}\left(\mu\frac{\partial}{\partial\mu}h\right)-\frac{1}{2}b^{2}h+3h+v_{\parallel}\frac{\partial}{\partial
v_{\parallel}}h+2\mu\frac{\partial}{\partial\mu}h\right],$ (84)

and

$\displaystyle\mathcal{C}^{F}$
$\displaystyle=\nu\left[2J_{0}T[h]\left(x+s_{\parallel}^{2}-\frac{3}{2}\right)+2u_{\parallel}[h]s_{\parallel}J_{0}(b)+J_{1}(b)2u_{\perp}[h]\sqrt{x}\right]F_{M},$
(85)

respectively. Here, the particle fluid quantities are defined by
$T[h]=(T_{\parallel}[h]+2T_{\perp}[h])/3-n[h]$ where $T_{\parallel}[h]=\int
d^{3}\bm{v}mv_{\parallel}^{2}J_{0}h/T$, $T_{\perp}[h]=\int d^{3}\bm{v}J_{0}\mu
Bh/T$, $n[h]=\int d^{3}\bm{v}J_{0}h/N$, and $u_{\parallel}[h]=\int
d^{3}\bm{v}J_{0}s_{\parallel}h$, $u_{\perp}[h]=\int
d^{3}\bm{v}\sqrt{x}J_{1}h$.

Projecting the GK Dougherty collision operator, Eqs. 84 and 85, onto the
Hermite-Laguerre basis yields,

$\displaystyle\mathcal{C}^{Tlk}=-\nu\left(2k+l+\frac{b^{2}}{2}\right)n^{lk}$
(86)

and

$\displaystyle\mathcal{C}^{Flk}=$
$\displaystyle\nu\left(T[h]\left[2\delta_{l}^{0}\left(2k\mathcal{K}_{k}-(k+1)\mathcal{K}_{k+1}-k\mathcal{K}_{k-1}\right)+\mathcal{K}_{k}\sqrt{2}\delta_{l}^{2}\right]\right.$
$\displaystyle\left.+u_{\parallel}[h]\mathcal{K}_{k}\delta_{l}^{1}+u_{\perp}[h]b\left(\mathcal{K}_{k}-\mathcal{K}_{k-1}\right)\delta_{l}^{0}\right),$
(87)

where the particle fluid quantities are expressed as

$\displaystyle n[h]$ $\displaystyle=\sum_{j=0}^{\infty}\mathcal{K}_{j}n^{0j},$
(88a) $\displaystyle T_{\parallel}[h]$
$\displaystyle=\sum_{j=0}^{\infty}\mathcal{K}_{j}\left[\sqrt{2}n^{2j}+n^{0j}\right],$
(88b) $\displaystyle T_{\perp}[h]$
$\displaystyle=\sum_{j=0}^{\infty}n^{0j}\left[(2j+1)\mathcal{K}_{j}-(j+1)\mathcal{K}_{j+1}-j\mathcal{K}_{j-1}\right],$
(88c) $\displaystyle u_{\parallel}[h]$
$\displaystyle=\sum_{j=0}^{\infty}\mathcal{K}_{j}n^{1j},$ (88d) $\displaystyle
u_{\perp}[h]$
$\displaystyle=b\sum_{j=0}^{\infty}\frac{1}{2}\left[\mathcal{K}_{j}-\mathcal{K}_{j-1}\right]n^{0j}.$
(88e)

We remark that the GK Dougherty collision operator, in Eq. 84 and Eq. 85, is
used in a previous Hermite-Laguerre velocity-space pseudo-spectral GK
formulation because of its sparse representation on this basis (Mandell et
al., 2018). In fact, $\mathcal{C}^{Tlk}=-\nu\left(2k+l+b^{2}/2\right)n^{lk}$
and $\mathcal{C}^{Flk}=0$ for $l>2$.

The DK Dougherty collision operator is obtained in the zero gyroradius limit
of Eqs. 86 and C, and is given by

$\displaystyle\mathcal{C}^{lk}=\nu\left[-(2k+l)N^{lk}-(\sqrt{2}N^{20}-2N^{01})\left(\frac{2}{3}\delta_{l}^{0}\delta_{k}^{1}-\frac{\sqrt{2}}{3}\delta_{l}^{2}\delta_{k}^{0}\right)+N^{10}\delta_{l}^{1}\delta_{k}^{0}\right].$
(89)

We remark that the DK Dougherty operator, given in Eq. 89, is equivalent to
the one used in Jorge et al. (2018).

##  Appendix D Collisionless Gyro-Moment Expressions

This appendix reports on the derivations of the collisionless gyro-moment
expressions, defined in Eq. 31. We start from Eq. 24 that we multiply by the
Hermite-Laguerre basis yielding

$\displaystyle\frac{N^{pj}}{\phi}=\sum_{l=1}^{3}\hat{N}_{l}^{pj},$ (90)

where we introduce the velocity integrals,

$\displaystyle\hat{N}_{1}^{pj}$
$\displaystyle=-\frac{q}{\tau}\int_{0}^{\infty}dx\int_{-\infty}^{\infty}\frac{ds_{\parallel}}{\sqrt{\pi}}\frac{H_{p}(s_{\parallel})L_{j}(x)}{\sqrt{2^{p}p!}}J_{0}(b\sqrt{x})e^{-s_{\parallel}^{2}-x},$
(91a) $\displaystyle\hat{N}_{2}^{pj}$
$\displaystyle=\frac{q}{\tau}\int_{0}^{\infty}dx\int_{-\infty}^{\infty}\frac{ds_{\parallel}}{\sqrt{\pi}}\frac{\omega
J_{0}(b\sqrt{x})e^{-x-s_{\parallel}^{2}}}{\omega-\omega_{\gradient
B}-\sqrt{2\tau}k_{\parallel}s_{\parallel}}\frac{H_{p}(s_{\parallel})L_{j}(x)}{\sqrt{2^{p}p!}},$
(91b) $\displaystyle\hat{N}_{3}^{pj}$
$\displaystyle=-\int_{0}^{\infty}dx\int_{-\infty}^{\infty}\frac{ds_{\parallel}}{\sqrt{\pi}}\omega_{T}^{*}\frac{J_{0}(b\sqrt{x})e^{-x-s_{\parallel}^{2}}}{\omega-\omega_{\gradient
B}-\sqrt{2\tau}k_{\parallel}s_{\parallel}}\frac{H_{p}(s_{\parallel})L_{j}(x)}{\sqrt{2^{p}p!}}.$
(91c)

We now perform analytically the velocity integrals that contain the resonant
term proportional to $1/(\omega-\omega_{\gradient
B}-\sqrt{2\tau}k_{\parallel}s_{\parallel})$ using the transformation Eq. 27
valid for unstable modes. The adiabatic part of the collisionless gyro-moment
response is deduced from the orthogonality relations, Eqs. 10 and 12, such
that

$\displaystyle\hat{N}_{1}^{pj}=-\frac{q}{\tau}\mathcal{K}_{j}\delta_{p}^{0}.$
(92)

Expanding the Bessel function in terms of Laguerre polynomials using Eq. 12,
$\hat{N}^{pj}_{2}$, defined in Eq. 91b, can be written as

$\displaystyle\hat{N}_{2}^{pj}$
$\displaystyle=-i\frac{q}{\tau}\sum_{n=0}^{\infty}\mathcal{K}_{n}\int_{0}^{\infty}d\tau\omega
e^{i\tau\omega}\int_{-\infty}^{\infty}\frac{ds_{\parallel}}{\sqrt{\pi}}e^{-(1+i\tau
2\alpha)s_{\parallel}^{2}-i\tau
z_{\parallel}s_{\parallel}}\frac{H_{p}(s_{\parallel})}{\sqrt{2^{p}p!}}$
$\displaystyle\times\int_{0}^{\infty}dxL_{n}(x)L_{j}(x)e^{-(1+i\tau\alpha)x},$
(93)

where $\alpha=\tau\omega_{B}/q$. The $x$-integration can be performed using
the identity (Gradshteyn & Ryzhik, 2014)

$\displaystyle\int_{0}^{\infty}dxe^{-\beta x}L_{n}(x)L_{m}(x)$
$\displaystyle=\frac{(m+n)!}{m!n!}\frac{(\beta-1)^{n+m}}{\beta^{m+n+1}}\mathcal{F}^{m}_{n}\left[\frac{\beta(\beta-2)}{(\beta-1)^{2}}\right],$
(94)

for $\text{Re }\beta>0$ and where $\mathcal{F}^{a}_{b}[z]=F[-a,-b;-a-b,z]$ is
the Gauss hypergeometric function. We remark that, in the slab limit
corresponding to $\alpha=0$ in Appendix D (and thus $\beta=1$ in Eq. 94)), Eq.
94 reduces to the orthogonality relation between Laguerre polynomials, Eq. 11.
The $s_{\parallel}$-integration in Appendix D is evaluated by expanding
$H_{p}(s_{\parallel})$ such that

$\displaystyle H_{p}(s_{\parallel})=\sum_{p_{1}=0}^{\lfloor
p/2\rfloor}\frac{(-1)^{p_{1}}2^{p-2p_{1}}p!}{p_{1}!(p-2p_{1})!}s_{\parallel}^{p-2p_{1}}.$
(95)

Using the above identities, we derive

$\displaystyle\hat{N}_{2}^{pj}$
$\displaystyle=-\frac{iq\omega}{\tau}\sum_{n=0}^{\infty}\frac{\mathcal{K}_{n}}{\sqrt{2^{p}p!}}\mathcal{I}_{n}^{pj},$
(96)

where we introduce the one-parameter integral

$\displaystyle\mathcal{I}_{n}^{pj}$ $\displaystyle=\sum_{p_{1}=0}^{\lfloor
p/2\rfloor}\frac{(-1)^{p_{1}}2^{p-2p_{1}-1}p!}{p_{1}!(p-2p_{1})!}\int_{0}^{\infty}\frac{d\tau}{\sqrt{\pi}}e^{i\tau\omega}(1+2i\tau\alpha)^{p_{1}-\frac{p}{2}-\frac{1}{2}}$
$\displaystyle\times\frac{(n+j)!}{n!j!}\frac{(i\alpha\tau)^{n+j}}{(1+i\alpha\tau)^{j+n+1}}\mathcal{F}^{j}_{n}\left[\frac{(1+i\alpha\tau)(i\tau\alpha-1)}{(i\alpha\tau)^{2}}\right]$
$\displaystyle\times\left\\{\left(1+(-1)^{p-2p_{1}}\right)\left(\frac{p}{2}-p_{1}-\frac{1}{2}\right)!\mathcal{M}^{{p/2-p_{1}+1/2}}_{1/2}\left[-\frac{z_{\parallel}^{2}\tau^{2}}{4(1+2\alpha
i\tau)}\right]\right.$
$\displaystyle\left.-2i\frac{z_{\parallel}\tau\left(1+(-1)^{p-2p_{1}+1}\right)}{2\sqrt{(1+2i\tau\alpha)}}\left(\frac{p}{2}-p_{1}\right)!\mathcal{M}^{\frac{p}{2}-p_{1}+1}_{3/2}\left[-\frac{z_{\parallel}^{2}\tau^{2}}{4(1+2\alpha
i\tau)}\right]\right\\},$ (97)

where $\mathcal{M}^{a}_{b}[z]=M(a;b;z)$ is the Kummer confluent hypergeometric
function, which stems from the $s_{\parallel}$-integration.

Finally, using the recursive properties of Hermite and Laguerre polynomials,
we deduce that

$\displaystyle\hat{N}_{3}^{pj}$
$\displaystyle=\frac{ik_{\perp}}{\sqrt{2^{p}p!}}\sum_{n=0}^{\infty}\mathcal{K}_{n}\left\\{\mathcal{I}_{n}^{pj}+\eta\left[(2j+1)\mathcal{I}_{n}^{pj}-j\mathcal{I}_{n}^{pj-1}-(j+1)\mathcal{I}_{n}^{pj+1}\right.\right.$
$\displaystyle\left.\left.+\frac{\mathcal{I}_{n}^{p+2j}}{4}+(p+\frac{1}{2})\mathcal{I}_{n}^{pj}+p(p-1)\mathcal{I}_{n}^{p-2j}-\frac{3}{2}\mathcal{I}_{n}^{pj}\right]\right\\}.$
(98)

With Eqs. 92, 96 and D, the closed semi-analytical expressions for the
collisionless gyro-moment, Eq. 31, can be evaluated. We remark that a simpler
expression for $\mathcal{I}_{n}^{pj}$ can be obtained in the purely toroidal
limit, i.e. $k_{\parallel}=0$, therefore, neglecting the resonance due to the
parallel streaming. Hence, from Appendix D, one obtains

$\displaystyle\mathcal{I}_{n}^{pj}$ $\displaystyle=\sum_{p_{1}=0}^{\lfloor
p/2\rfloor}\frac{(-1)^{p_{1}}2^{p-2p_{1}}p!}{p_{1}!(p-2p_{1})!}\frac{(n+j)!}{n!j!\sqrt{\pi}}\left(\frac{p}{2}-p_{1}-\frac{1}{2}\right)!$
$\displaystyle\times\int_{0}^{\infty}d\tau
e^{i\tau\omega}(1+2i\tau\alpha_{s})^{p_{1}-\frac{p}{2}-\frac{1}{2}}\frac{(i\alpha_{s}\tau)^{n+j}}{(1+i\alpha_{s}\tau)^{j+n+1}}\mathcal{F}^{j}_{n}\left[\frac{1+\alpha_{s}^{2}\tau^{2}}{\alpha_{s}^{2}\tau^{2}}\right],$
(99)

if $p$ is even and $0$ otherwise.

##  Appendix E FLR Closures

We explore possible FLR closures for the terms, proportional to
$\mathcal{K}_{0}$, $\mathcal{K}_{1}$ and $\mathcal{K}_{2}$ appearing in the
$6$GM and $4$GM using the Padé approximation technique. In order to deduce
Padé approximation of the kernel functions, we rewrite the exact form of
$\mathcal{K}_{n}$ defined in Eq. 13 as

$\displaystyle\mathcal{K}_{j}=\frac{(-1)^{j}}{j!}\frac{\partial^{j}}{\partial\beta^{j}}\Upsilon(a\beta)\bigg{\rvert}_{\beta=1},$
(100)

where we introduce the basic FLR operator
$\Upsilon(a\beta)=\sum_{l=0}^{\infty}(\beta a)^{l}/(l!2^{l})=e^{-\beta a/2}$.
From Eq. 100, approximation to $\mathcal{K}_{j}$ can be constructed by
specifying the basic FLR operator $\Upsilon$. We remark that in contrast to
the choice $\Upsilon(a\beta)=e^{-\beta a/2}$ introduced in (Brizard, 1992),
Dorland & Hammett (1993) proposed a modified basic FLR operator such as
$\Upsilon(a\beta)=\sqrt{\Gamma_{0}(a\beta)}$. The latter approximation is
motivated by the fact that basic FLR operator, $\Upsilon(a\beta)=e^{-\beta
a/2}$, yields a large damping at short wavelength with respect to Dorland &
Hammett (1993) when applied to, e.g., the collisionless ITG case. As
collisions damp short wavelength modes (see Fig. 8), we focus on the basic FLR
operator $\Upsilon(a\beta)=e^{-\beta a/2}$, since little difference with
respect to $\sqrt{\Gamma_{0}(a\beta)}$ is expected in the high collisional
limits where short wavelength ITG modes are damped.

Figure 20: (Left) lowest order kernel functions $\mathcal{K}_{j}$ (solid
line), and the correspinding Padé approximants,
${}_{q}^{p}\check{\mathcal{K}_{j}}$ for $j=0,1,2$ when $(p,q)=(1,2)$ (dashed
line) and $(p,q)=(1,4)$ (dotted line). (Right) ITG linear growth rate $\gamma$
as a function of the perpendicular wavenumber $k_{\perp}$, obtained by the
$4$GM with the exact kernel functions $\mathcal{K}_{j}$ (solid line) and with
the Padé approximants ${}_{q}^{p}\check{\mathcal{K}_{j}}$, with $(p,q)=(1,2)$
(dashed line) and $(p,q)=(1,4)$ (dotted line), respectively.

We now aim to construct a Padé approximant, of order $(p,q)$, of the basic FLR
operator, that we denote by
${}_{q}^{p}\check{\Upsilon}(x)=R^{(p)}(x)/S^{(q)}(x)$ where $R^{(p)}(x)$ and