2. Gyrokinetic model equation

We consider the electrostatic GK Boltzmann equation in the presence of magnetic field, density and temperature gradients. We use the gyrocentre phase-space coordinates $\boldsymbol {Z} = (\boldsymbol {R}, \mu , v_\parallel ,\theta )$, where $\boldsymbol {R} = \boldsymbol {r} - \boldsymbol {\rho }_a$ is the gyrocentre position (with $a$ the species index), with $\boldsymbol {r}$ the particle position and $\boldsymbol {\rho }_a(\boldsymbol {R}, \mu , \theta ) = \boldsymbol {b} \times \boldsymbol {v}/\varOmega _a$ the gyroradius of species $a$ ($\boldsymbol {b} = \boldsymbol {B}/ B$ and $\varOmega _a = q_a B / m_a$), $\mu = m_a v_\perp ^2/[2 B(\boldsymbol {R})]$ is the magnetic moment, $v_\parallel = \boldsymbol {b} \boldsymbol {\cdot } \boldsymbol {v}$ is the component of the velocity parallel to the magnetic field, and $\theta$ is the gyroangle. We linearize the GK Boltzmann equation by assuming that the gyrocentre distribution function of species $a$, $F_a = F_a (\boldsymbol {R}, \mu , v_\parallel )$, is a perturbed Maxwellian, that is $F_{a} = F_{Ma} + g_a$, with ${g}_a = {g}_a(\boldsymbol {R} ,\mu , v_\parallel , t)$, the perturbation with respect to the local Maxwellian distribution function $F_{Ma} = N \exp ({- s_{\parallel a}^2 - x_a}) / ({\rm \pi} ^{3/2} v_{Ta}^3)$ such that ${g}_{a}/ F_{Ma} \ll 1$, with $N = N_i(\boldsymbol {R}) = N_e(\boldsymbol {R})$ the background gyrocentre density, $s_{ \parallel a} = v_\parallel / v_{Ta}(\boldsymbol {R})$, $x_a= \mu B(\boldsymbol {R}) / T_{a}(\boldsymbol {R})$ and $v_{T_a}^2(\boldsymbol {R}) = 2 T_{a}(\boldsymbol {R})/m_a$ the thermal velocity based on the equilibrium temperature $T_a(\boldsymbol {R})$. The linearized electrostatic GK Boltzmann equation that describes the time evolution of ${g}_a(\boldsymbol {k},t) = \int \,\textrm {d} \boldsymbol {R} {g}_a(\boldsymbol {R} , t) \exp ({- \textrm {i} \boldsymbol {R} \boldsymbol {\cdot } \boldsymbol {k}})$ can be written as (Hazeltine & Meiss Reference Hazeltine and Meiss2003)

(2.1)\begin{equation} \frac{\partial}{\partial t } {g}_a + \textrm{i} \omega_{Ba} h_a + v_\parallel \boldsymbol{\nabla}_\parallel h_a - \frac{\mu}{m_s} (\boldsymbol{b} \boldsymbol{\cdot} \boldsymbol{\nabla} B)\frac{\partial }{\partial v_\parallel} h_a - \textrm{i} \omega_{Ta}^* \textrm{J}_0(b_a \sqrt{x_a})\phi F_{Ma} = \sum_{b} \mathcal{C}_{ab}, \end{equation}

where $\phi = \phi (\boldsymbol {k}, t)$ is the Fourier component of the electrostatic potential and we introduce the non-adiabatic part, $h_a = h_a (\boldsymbol {k}, \mu , v_\parallel ,t)$, of the perturbed gyrocentre distribution function ${g}_a$, that is

(2.2)\begin{equation} h_a = {g}_a + \frac{q_a}{T_{a}} F_{Ma} \textrm{J}_0(b_a \sqrt{x_a}) \phi. \end{equation}

In (2.1), we define $\boldsymbol {\nabla }_\parallel = \boldsymbol {b} \boldsymbol {\cdot } \boldsymbol {\nabla }$, $\omega _{Ba} = v_{Ta}^2 ( x_a+ 2 s_{\parallel a}^2) \omega _B/ (2 \varOmega _a)$, with $\omega _B = ( \boldsymbol {b} \times \boldsymbol {\nabla } \ln B ) \boldsymbol {\cdot } \boldsymbol {k}$, and $\omega _{Ta}^* = [\omega _N + \omega _{T_a}( x_a + s_{\parallel a}^2 -3/2) ]$, with $\omega _N = \boldsymbol {b} \times \boldsymbol {\nabla } \ln N \boldsymbol {\cdot } \boldsymbol {k}/ B$ and $\omega _{T_a} = \boldsymbol {b} \times \boldsymbol {\nabla } \ln T_a \boldsymbol {\cdot } \boldsymbol {k}/B$. Finite Larmor radius effects are taken into account through the zeroth-order Bessel function, $\textrm {J}_0(b_a \sqrt {x_a})$, where $b_a = k_\perp v_{T_a} /\varOmega _a$ is the normalized perpendicular wavevector, with $k_\perp = |\boldsymbol {k} - (\boldsymbol {b} \boldsymbol {\cdot } \boldsymbol {k}) \boldsymbol {b}|$. On the right-hand side of (2.1), we define the linearized GK collision operator between species $a$ and $b$, $\mathcal {C}_{ab} = \mathcal {C}_{ab}(\boldsymbol {R}, \mu , v_\parallel )$,

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

that describes the effects of small-angle Coulomb collisions between particle species $a$ and $b$. While the linearized collision operator $C_{ab}$ is the focus of § 3, we remark that $\left \langle \dots \right \rangle _{\boldsymbol {R}} =\int _0^{2 {\rm \pi}} \,\textrm {d} \theta \ldots /(2 {\rm \pi})$ denotes the gyroaverage operator evaluated at the gyrocentre position (the subscript $\boldsymbol {R}$ indicates that the integral over the gyroangle should be computed holding $\boldsymbol {R}$ constant). We remark that the linearized GK collision operator is defined as a function of the gyrocentre phase space, i.e. $\mathcal {C}_{ab} = \mathcal {C}_{ab}(\boldsymbol {R}, \mu , v_\parallel )$, and it is therefore gyrophase independent.

The linearized Boltzmann equation is closed by the GK quasineutrality condition that determines self-consistently the electrostatic potential,

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

where $a_a= b_a^2 / 2$, and $\varGamma _0(x) = \textrm {I}_0(x) \textrm {e}^{- x}$ with $\textrm {I}_0$ the modified Bessel function.

In order to approach the solution of the GK system, (2.1) and (2.4), we use a Hermite–Laguerre moment expansion of the gyrocentre perturbed distribution function ${g}_a$ (Mandell, Dorland & Landreman Reference Mandell, Dorland and Landreman2018; Frei et al. Reference Frei, Jorge and Ricci2020). Then, by projecting (2.1) onto the complete Hermite–Laguerre basis polynomials, we reduce its dimensionality removing the $(\mu , v_\parallel )$ dependence. More precisely, we decompose the perturbed gyrocentre distribution function, ${g}_a$, onto a set of Hermite–Laguerre polynomials (Jorge et al. Reference Jorge, Ricci and Loureiro2017, Reference Jorge, Frei and Ricci2019; Frei et al. Reference Frei, Jorge and Ricci2020),

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

where the Hermite–Laguerre velocity moments of ${g}_a$ are defined as

(2.6)\begin{equation} N_a^{pj} = \frac{1}{N_a} \int \,\textrm{d} \mu\, \textrm{d} v_{{\parallel}} \,\textrm{d} \theta\frac{B}{m_a} {g}_a \frac{H_p(s_{{\parallel} a}) L_j(x_a)}{\sqrt{2^p p!}}, \end{equation}

with $N_a = \int \,\textrm {d} \mu \,\textrm {d} v_{\parallel } \,\textrm {d} \theta B {g}_a / m_a$ the gyrocentre density. In (2.5), we introduce the Hermite and Laguerre polynomials, $H_p$ and $L_j$, via their Rodrigues’ formulae (Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik2014)

(2.7)\begin{gather} H_p(x) = ({-}1)^p \textrm{e}^{x^2} \frac{\textrm{d}^p}{{\textrm{d}\kern0.07pt x}^p} \left( \textrm{e}^{- x^2} \right), \end{gather}

(2.8)\begin{gather}L_j(x) = \frac{\textrm{e}^{x}}{j!} \frac{\textrm{d}^j }{{\textrm{d}\kern0.07pt x}^j} \left ( \textrm{e}^{- x } x^j\right) \end{gather}

and we note their orthogonality relations over the intervals, $] -\infty , \infty [$ weighted by $\textrm {e}^{-x^2}$, and $[0,+ \infty [$ weighted by $\textrm {e}^{-x}$, respectively,

(2.9)\begin{gather} \int_{- \infty}^\infty\,{\textrm{d}\kern0.07pt x} H_p(x) H_{p'}(x) \textrm{e}^{- x^2} = 2^p p! \sqrt{\rm \pi} \delta_{p}^{p'}, \end{gather}

(2.10)\begin{gather}\int_0^\infty \,{\textrm{d}\kern0.07pt x} L_j(x) L_{j'}(x) \textrm{e}^{{-}x} = \delta_j^{j'}. \end{gather}

We denote the Hermite–Laguerre coefficients, $N_a^{pj}$ in (2.5), as the gyromoments. We remark that any function, $f = f (\mu , v_\parallel )$, that satisfies

(2.11)\begin{equation} \int \,\textrm{d} \mu \,\textrm{d} v_\parallel \,\textrm{d} \theta \frac{B}{m_a} |\,f|^2 \textrm{e}^{{-}s_{{\parallel} a}^2 - x_a} <{+} \infty, \end{equation}

can be decomposed onto the orthogonal basis defined by the Hermite–Laguerre polynomials (Wong Reference Wong1998). This is, in particular, fulfilled by the perturbed gyrocentre distribution function $g_a$, and also by its non-adiabatic part, $h_a$.

We now project (2.1) and (2.4) onto the Hermite–Laguerre polynomial basis. For this purpose, we note that the Bessel function $\textrm {J}_0$ (appearing in both (2.1) and (2.4) and arising from FLR effects) and, more generally $\textrm {J}_m$, with $m > 0$, can be conveniently expanded onto associated Laguerre polynomials, $L^m_n(x) = (-1)^m d^m L_{n + m}(x) / {\textrm {d}\kern 0.05em x}^m$, such as (Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik2014)

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

where we introduce the velocity-independent expansion coefficients

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

We refer to $\mathcal {K}_{n}(b_a)$ as the $n$th-order kernel function, with velocity-independent argument $b_a$ (Frei et al. Reference Frei, Jorge and Ricci2020). In fact, (2.12) allows us to isolate the velocity dependence, $\sqrt {x_a}$, appearing in the argument of the Bessel function $\textrm {J}_m$ from fluid quantities such as $b_a$.

Multiplying the GK Boltzmann equation, (2.1), by the basis element $B H_p(s_{\parallel a}) L_j(x_a) / \sqrt {2^p p!}$, and integrating over the velocity space yields the gyromoment hierarchy equation,

(2.14)$$\begin{align} & \frac{\partial}{\partial t} N_a^{pj} + v_{Ta} \boldsymbol{\nabla}_\parallel \left( \sqrt{\frac{p+1}{2}} n_a^{p+1j}+ \sqrt{\frac{p}{2}} n_a^{p-1j} \right) \nonumber\\ &\qquad - v_{Ta}\boldsymbol{\nabla}_\parallel \ln B \left[ (\,j\!+\!1)\sqrt{p+1} n_a^{p+1j}+ (\,j+1)\sqrt{p} n_a^{p-1j} - j \sqrt{p+1} n_a^{p+1j-1}- j \sqrt{p} n_a^{p-1j-1 } \right] \nonumber\\ &\qquad + \frac{ T_a \textrm{i} \omega_B}{q_a B} \left[ \sqrt{(p+1)(p+2)} n_a^{p+2j} +(2p+1) n_a^{pj}+ \sqrt{p(p-1)} n_a^{p-2j} \right. \nonumber\\ & \qquad\left. + (2 j +1) n_a^{pj} - j n_a^{pj-1} - (\,j+1) n_a^{pj+1} \right] \nonumber\\ &\qquad + v_{Ta} \boldsymbol{\nabla}_\parallel \ln B \sqrt{ \frac{p}{2}} \left[ (2 j +1) n_a^{p-1j} - (\,j +1) n_a^{p-1j+1} - j n_a^{p-1 j-1} \right] \nonumber\\ &\qquad - \textrm{i} \phi \left\{ \mathcal{K}_{j} \delta_p^0 \omega_N + \omega_{T_a}\left[\mathcal{K}_{j} \left(\frac{\delta_p^2}{\sqrt{2}} - \delta_p^0\right) + \delta_p^0 \left((2 j+1) \mathcal{K}_{j} - j \mathcal{K}_{j-1 } - (\,j+1)\mathcal{K}_{j+1}\right) \right] \right\} \nonumber\\ &\quad = \sum_{b} \mathcal{C}_{ab}^{pj}, \end{align}$$

where we define $\mathcal {C}_{ab}^{pj} =\mathcal {C}_{ab}^{pj}(\boldsymbol {R},t)$ as the Hermite–Laguerre gyromoment expansion of the linearized GK collision operator $\mathcal {C}_{ab}$,

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

In (2.14), the non-adiabatic gyromoments $n_a^{pj}$, are defined by

(2.16)\begin{equation} n_a^{pj} = \frac{1}{N_a} \int \,\textrm{d} \mu\, \textrm{d} v_{{\parallel}}\, \textrm{d} \theta\frac{B}{m_a} h_a \frac{H_p(s_{{\parallel} a}) L_j(x_a)}{\sqrt{2^p p!}}, \end{equation}

and can be expressed in terms of the gyromoments $N_a^{pj}$ by projecting (2.2) onto the Hermite–Laguerre basis, yielding

(2.17)\begin{equation} n_a^{pj} = N^{pj}_a + \frac{q_a}{ T_a} \mathcal{K}_{j} \phi \delta_p^0. \end{equation}

Finally, the GK quasineutrality condition projected onto the Hermite–Laguerre basis is expressed as

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

The linearized gyromoment hierarchy, (2.14), with the GK quasineutrality condition (2.18), describes the evolution of the gyromoments as it results from the interplay of Landau damping, magnetic gradient drifts, magnetic trapping and, finally, the driving density and temperature gradients, $\omega _N$ and $\omega _{Ta}$. Collisional effects are expressed by the collisional gyromoment of the GK collision operator $\mathcal {C}_{ab}^{pj}$ defined in (2.15). Contrary to the GK Boltzmann equation, given in (2.1), that depends on the phase-space coordinates, the gyromoment hierarchy, (2.14), constitutes an infinite system of fluid-like equations for the variables $N_a^{pj}$, which are only functions of the Fourier mode wavevector $\boldsymbol {k}$ and time $t$.

We remark that (2.14) is equivalent to the previous Hermite–Laguerre moment hierarchy derived by Mandell et al. (Reference Mandell, Dorland and Landreman2018), where the probabilist's Hermite polynomials, $He_{p}(x)$, are used instead, defined by $He_{p}(x) = 2^{-p/2} H_p(x/2)$ (with $H_p(x)$ in (2.7)), albeit with collisional effects modelled by a linearized GK Dougherty operator (Dougherty Reference Dougherty1964), which is characterized by a sparse Hermite–Laguerre representation.

3. Linearized collision operator models

In general, a linearized collision operator between particles of species $a$ and $b$, $C_{ab}$, is obtained by linearizing a nonlinear collision operator, $C_{ab}^{NL}$. Nonlinear collision operators are usually defined on the full particle distribution function expressed in the particle phase space $(\boldsymbol {r}, \boldsymbol {v})$. This is because collisions occur at the particle position $\boldsymbol {r}$ (rather than at the gyrocentre position $\boldsymbol {R}$). Given $f_{Ma} = f_{Ma}(\boldsymbol {r}, \boldsymbol {v})$, the particle Maxwellian distribution function of species $a$, and $f_a = f _a(\boldsymbol {r},\boldsymbol {v})$ its small amplitude perturbation, $f_a / f_{Ma} \ll 1$, the linearized collision operator, $C_{ab} = C_{ab}(\boldsymbol {r}, \boldsymbol {v}) = C_{ab} (f_a,f_b)$, can be expressed as

(3.1)\begin{equation} C_{ab} = C_{ab}^T + C_{ab}^F, \end{equation}

where we introduce the test component of the linearized collision operator, $C_{ab}^T = C_{ab}^T(\boldsymbol {r}, \boldsymbol {\boldsymbol {v}}) = C_{ab}^T(f_a(\boldsymbol {r}, \boldsymbol {\boldsymbol {v}})) = C^{NL}(f_a,f_{Mb})$, and the field component, $C_{ab}^F = C_{ab}^F(\boldsymbol {r}, \boldsymbol {\boldsymbol {v}}) = C_{ab}^F (f_b(\boldsymbol {r}, \boldsymbol {\boldsymbol {v}})) = C^{NL}(f_{Ma},f_b)$. Local conservation properties of the collision operator constraint the test and field components. In fact, while particle conservation is satisfied by both test and field components (Xu & Rosenbluth Reference Xu and Rosenbluth1991), such that

(3.2)\begin{equation} \int \,\textrm{d} \boldsymbol{v} C_{ab}^T = \int \,\textrm{d} \boldsymbol{v} C_{ab}^F = 0, \end{equation}

the momentum and energy conservations require that

(3.3)\begin{gather} \int \,\textrm{d} \boldsymbol{v} m_a \boldsymbol{v} C_{ab}^T + \int \,\textrm{d} \boldsymbol{v} m_b \boldsymbol{v} C_{ba}^F = 0, \end{gather}

(3.4)\begin{gather}\int \,\textrm{d} \boldsymbol{v} m_a v^2 C_{ab}^T + \int \,\textrm{d} \boldsymbol{v} m_b v^2 C_{ba}^F = 0, \end{gather}

respectively. We remark that the velocity integrals in (3.2)–(3.4) are evaluated at constant $\boldsymbol {r}$, and that local conservation laws do not hold in the gyrocentre coordinates $(\boldsymbol {R}, \mu , v_\parallel , \theta )$, because of the difference between the gyrocentre and particle positions, being $\boldsymbol {r} - \boldsymbol {R} = \boldsymbol {\rho }_a(\boldsymbol {R},\mu ,\theta )$. As an example of the implication of the difference between $\boldsymbol {r}$ and $\boldsymbol {R}$, we consider the zeroth-order gyromoment of the GK linearized collision operator. In Fourier space and performing the velocity at constant $\boldsymbol {R}$, one has

(3.5)$$\begin{align} \int \,\textrm{d} \boldsymbol{v} \left\langle C^T_{ab}(\boldsymbol{r},\boldsymbol{v}) \right\rangle_{\boldsymbol{R}} & = 2 {\rm \pi}\int\, \textrm{d} \boldsymbol{k} \int \textrm{d} v_\parallel \,\textrm{d} \mu \frac{B}{m_a} \left\langle \exp({ \textrm{i} \boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{\rho}_a})C^T_{ab}(\boldsymbol{k},\mu, v_\parallel, \theta) \right\rangle_{\boldsymbol{R}} \textrm{e}^{\textrm{i} \boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{R}} \nonumber\\ & = 2 {\rm \pi}\int\, \textrm{d} \boldsymbol{k} \sum_{j=0}^\infty \int \textrm{d} v_\parallel \,\textrm{d} \mu \frac{B}{m_a} \left\langle \frac{(\textrm{i}\boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{\rho}_a)^j}{j!} C^T_{ab}(\boldsymbol{k},\mu, v_\parallel, \theta) \right\rangle_{\boldsymbol{R}} \textrm{e}^{\textrm{i} \boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{R}}. \end{align}$$

Since $C_{ab}^T$ conserves particles, i.e. $\int \,\textrm {d} \boldsymbol {v} C_{ab}^T = 0$, we expect only the $j = 0$ term in the $j$ sum to vanish when the gyroaverage is performed, and the $j \neq 0$ terms to give finite contributions as the particle gyroradius $\boldsymbol {\rho }_a =\boldsymbol {\rho }_a (\boldsymbol {R}, \mu ,\theta )$ depends on the gyroangle $\theta$. In fact, (3.5) does not vanish, and the linearized GK collision operator yields classical gyrodiffusion in gyrocentre phase space. On the other hand, in the DK limit, where finite $\boldsymbol {k} \boldsymbol {\cdot } \boldsymbol {\rho }_a$ effects in the collision operator are neglected, i.e. $\left \langle C_{ab}(\boldsymbol {r},\boldsymbol {v}) \right \rangle _{\boldsymbol {R}} \simeq \left \langle C_{ab}(\boldsymbol {R},\boldsymbol {v}) \right \rangle _{\boldsymbol {R}}$, the local conservation laws, in (3.2)–(3.4), for the gyrocentre positions are satisfied.

In the following, of the present section, we derive the Hermite–Laguerre expansion of two advanced linearized GK collision operators, that are the linearized GK Coulomb collision operator (Jorge et al. Reference Jorge, Frei and Ricci2019) and the linearized GK Sugama collision operator (Sugama et al. Reference Sugama, Watanabe and Nunami2009). First, following the procedure outlined in Jorge et al. (Reference Jorge, Frei and Ricci2019) in § 3.1, we derive the linearized GK Coulomb collision operator using a spherical harmonic moment expansion. This allows us to perform the gyroaverage and the expansion of the Coulomb collision operator in terms of gyromoments analytically at arbitrary order in the perpendicular wavenumber. We also obtain the DK limit by explicitly taking the zeroth gyroradius limit, i.e. $\boldsymbol {k} \boldsymbol {\cdot } \boldsymbol {\rho }_a \to 0$. Similarly, we project the GK Sugama collision operator (Sugama et al. Reference Sugama, Watanabe and Nunami2009) onto the Hermite–Laguerre basis in § 3.2, and derive its DK limit.

3.1. Linearized GK Coulomb collision operator

Starting from the nonlinear Coulomb collision operator, the test component of the linearized Coulomb collision operator around a local particle Maxwellian distribution function, $f_{Mb}$, is given in the particle pitch-angle coordinates $(\boldsymbol {r}, v,\xi , \theta )$ (with $v$ the modulus of $\boldsymbol {v}$ and $\xi = v_\parallel /v$ the pitch angle) by (Rosenbluth et al. Reference Rosenbluth, MacDonald and Judd1957; Helander & Sigmar Reference Helander and Sigmar2002; Hazeltine & Meiss Reference Hazeltine and Meiss2003)

(3.6)$$\begin{align} C^T_{ab}(f_a) & = \frac{m_a \nu_{ab} v_{Ta}}{n_b} \left[ \left( \frac{2 \partial_v G(f_{Mb})}{ v^2} + \left( 1 - \frac{m_a}{m_b} \right) \partial_v H(f_{Mb}) \right) \partial_v f_{a} \right. \nonumber\\ &\quad \left. - \frac{1}{v^3} \partial_v G(f_{Mb}) \mathcal{L}^2 f_a + \partial_{v}^2 G(f_{Mb}) \partial^2_v f_{a} + \frac{m_a}{m_b} 8 {\rm \pi}f_{Mb} f_a\right]. \end{align}$$

In (3.6), we introduce the Rosenbluth potentials $H(f) = 2 \int \,\textrm {d} \boldsymbol {v}' f(\boldsymbol {v}') /u$ and $G(f) = \int \,\textrm {d} \boldsymbol {v}' f(\boldsymbol {v}')u$ (with $u = |\boldsymbol {v} - \boldsymbol {v}'|$), the collision frequency between species $a$ and $b$, $\nu _{ab} = 4 \sqrt {{\rm \pi} } N_b q_a^2 q_b^2 \ln \varLambda /(3 m_a^{1/2} T_a^{3/2})]$ and the spherical angular operator $\mathcal {L}^2 f = \partial _{v}(v^2 \partial _v f) - v^2 \boldsymbol {\nabla }_{\boldsymbol {v}}^2 f$. Finally, we note that the thermal velocity is defined as $v_{Ta}^2 = 2 T_a / m_a$. When evaluated with a Maxwellian distribution function, analytical expressions of the Rosenbluth potentials can be obtained, and are given by $H(f_{Mb}) = 2 n_b \mathrm {erf}(s_b)/v$ and $G(f_{Mb}) = n_b v_{Tb} [ ( 1 + 2 s_b^2 )\mathrm {erf}(s_b)/s_b + \mathrm {erf}'(s_b)]/2$, with $s_b = v / v_{Tb}$, the error function, $\mathrm {erf}(x) = 2 \int _0^x \,\textrm {d} s \textrm {e}^{-s^2}/\sqrt {{\rm \pi} }$ and its derivative, $\mathrm {erf}'(x) = 2 \textrm {e}^{-x^2}/\sqrt {{\rm \pi} }$.

The field component of the linearized Coulomb collision operator, expressed in the particle coordinates $(\boldsymbol {r}, v,\xi , \theta )$, is

(3.7)\begin{equation} C^F_{ab}(f_b) = \frac{2 \nu_{ab} v_{Ta} f_{Ma}}{n_b} \left[ 2 s_a^2 \partial_v^2 G(f_b) - H(f_b) - \left( 1 - \frac{m_a}{m_b} \right) v \partial_v H(f_b) + \frac{m_a}{m_b} 4 {\rm \pi}v_{Ta}^2 f_b \right]. \end{equation}

We refer to (3.6) and (3.7), respectively, as the exact test and field components of the Coulomb collision operator between species $a$ and $b$ (Rosenbluth et al. Reference Rosenbluth, MacDonald and Judd1957). We remark that, while the test component, (3.6), contains velocity derivatives evaluated at constant $\boldsymbol {r}$, the field component, (3.7), has additional velocity-space integrals contained in the Rosenbluth potential to evaluate holding $\boldsymbol {r}$ constant.

Since the gyrocentre coordinate transformation mixes spatial and velocity coordinates, such as $\boldsymbol {r} = \boldsymbol {R} + \boldsymbol {\rho }_a(\boldsymbol {R},\mu , \theta )$, and challenges the analytical and numerical treatments of the exact field component, ad hoc field components have been developed to implement collisional effects in GK codes (Abel et al. Reference Abel, Barnes, Cowley, Dorland and Schekochihin2008; Sugama et al. Reference Sugama, Watanabe and Nunami2009, Reference Sugama, Matsuoka, Satake, Nunami and Watanabe2019). Following closely the method adopted in Jorge et al. (Reference Jorge, Frei and Ricci2019), we overcome the difficulty associated with the evaluation of the exact Coulomb collision operator by using a spherical harmonic moment expansion of the particle distribution function. Ultimately, this allows us to evaluate the exact test and field components of the linearized gyroaveraged Coulomb collision operator, $\mathcal {C}_{ab} = \left \langle C_{ab} \right \rangle _{\boldsymbol {R}}$, at arbitrary order in the perpendicular wavenumber.

3.1.1. Spherical harmonics moment expansion and gyroaverage

In order to derive an expression of the linearized Coulomb collision operator ready for an expansion in the Hermite–Laguerre basis, following Ji & Held (Reference Ji and Held2006) and Jorge et al. (Reference Jorge, Frei and Ricci2019), the perturbed particle distribution function of species $a$, $f_a = f_a(\boldsymbol {x}, \boldsymbol {v})$, is expanded in the spherical harmonics basis as

(3.8)\begin{equation} f_a = \sum_{p=0}^\infty \sum_{j=0}^\infty \frac{f_{Ma}}{\sigma_j^p} \boldsymbol{\mathsf{M}}_a^{pj} \boldsymbol{\cdot} \boldsymbol{\mathsf{Y}}^{p}(\boldsymbol{s}_a) L_j^{p+1/2}(s_a^2), \end{equation}

with $\sigma _j^p = p!(p + j +1/2)!/[2^p (p +1/2)! j!]$ and $\boldsymbol {s}_a = \boldsymbol {v} / v_{Ta}$. We remark that the quantities $\boldsymbol{\mathsf{M}}_a^{pj}$ are defined below in (3.14). In (3.8), the velocity-space basis, $\boldsymbol{\mathsf{Y}}^{p}(\boldsymbol {s}_a) L_j^{p+1/2}(s_a^2)$, is the product between associated Laguerre polynomials, $L_j^{p+1/2}(s_s^2)$, defined as (Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik2014)

(3.9)\begin{equation} L_j^{p+1/2}(x) = \sum_{l=0}^j L_{jl}^{p} x^l, \end{equation}

with

(3.10)\begin{equation} L_{jl}^{p} = \frac{({-}1)^l (p + j + 1/2)!}{(\,j - l)!(l +p + 1/2 )! l!} \end{equation}

and the $p$th-order traceless symmetric spherical harmonic tensors, $\boldsymbol{\mathsf{Y}}^p(\boldsymbol {s}_a)$. The spherical harmonic tensors have the property that $\boldsymbol {\nabla }_{\boldsymbol {v}}^2 \boldsymbol{\mathsf{Y}}^p =0$, such that $\mathcal {L}^2 \boldsymbol{\mathsf{Y}}^p = p(p+1)\boldsymbol{\mathsf{Y}}^p$ and $\boldsymbol{\mathsf{Y}}^p(\boldsymbol {v}) = v^p \boldsymbol{\mathsf{Y}}^p(\hat {\boldsymbol {v}})$ ($\hat {\boldsymbol {v}} = \boldsymbol {v} / v$) (Jorge et al. Reference Jorge, Frei and Ricci2019), and can be explicitly defined by introducing the spherical harmonic basis, $\boldsymbol {e}^{pm}$, that satisfies the orthogonality relation $\boldsymbol {e}^{pm} \boldsymbol {\cdot }\boldsymbol {e}^{pm'} = (-1)^m \delta _{-m}^{m'}$ (Snider Reference Snider2017). The expression of the spherical harmonic tensors, $\boldsymbol{\mathsf{Y}}^p(\boldsymbol {s}_a)$, on that basis is

(3.11)\begin{equation} \boldsymbol{\mathsf{Y}}^p (\boldsymbol{s}_a)= s_a^p \sqrt{\frac{2 {\rm \pi}^{3/2} p!}{2^p (p+1/2)!}}\sum_{m={-}p}^{p } Y_p^m(\xi,\theta) \boldsymbol{e}^{pm}, \end{equation}

where

(3.12)\begin{equation} Y_p^m(\xi,\theta) = \sqrt{\frac{(2p+1)(p-|m|)!}{4 {\rm \pi}(p+|m|)!}} P_p^m(\xi) \textrm{e}^{\textrm{i}m\theta}, \end{equation}

are scalar harmonic functions, with $P_p^m(\xi ) = (-1)^m (1 - \xi ^2)^{{m/2}} \textrm {d}^{m} P_p(\xi ) / \textrm {d} \xi ^m$ the associated Legendre polynomials, being with $P_p(\xi ) = \textrm {d}^p [(\xi ^2 -1)^p] / \textrm {d} \xi ^p / [2^p p!]$ the Legendre polynomial (Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik2014). We remark that (3.11) is particularly useful for performing analytically the gyroaverage of the Coulomb collision operator since the $\theta$ dependence of $\boldsymbol{\mathsf{Y}}^{p}$ is isolated in $Y_p^m(\xi , \theta )$. The spherical harmonic basis, used in the expansion, given in (3.8), satisfies the orthogonality relation (Snider Reference Snider2017)

(3.13)\begin{equation} \int\, \textrm{d} \boldsymbol{s} \textrm{e}^{- s^2} L_j^{p+1/2}(s^2) \boldsymbol{\mathsf{Y}}^{p'}(\boldsymbol{s}) L_{j'}^{p'+1/2}(s^2) \boldsymbol{\mathsf{Y}}^{p}(\boldsymbol{s}) L_j^{p+1/2}(s^2) \boldsymbol{\cdot} \boldsymbol{\mathsf{T}}^{p} = \delta_{pp'} \delta_{ll'} {\rm \pi}^{3/2} \sigma_p^j \boldsymbol{\mathsf{T}}^{p}, \end{equation}

with $\boldsymbol{\mathsf{T}}^p$ an arbitrary $p$th-order tensor. We note that the dot product, appearing in (3.8) and (3.13), between two $p$th-order tensors, $\boldsymbol{\mathsf{T}}^p$ and $\boldsymbol {R}^p$, i.e. $\boldsymbol{\mathsf{T}}^p \boldsymbol {\cdot } \boldsymbol {R}^p$, yields a scalar.

Projecting the spherical harmonic basis and using (3.8), one deduces that the expansion coefficient $\boldsymbol{\mathsf{M}}_a^{pj} =\boldsymbol{\mathsf{M}}_a^{pj}(\boldsymbol {r})$, defined as the spherical harmonic particle moments of species $a$, are given by

(3.14)\begin{equation} \boldsymbol{\mathsf{M}}_a^{pj} = \frac{1}{n_a }\int \,\textrm{d} \boldsymbol{v} f_a \boldsymbol{\mathsf{Y}}^{p}(\boldsymbol{s}_a) L_j^{p+1/2}(s_a^2), \end{equation}

with $n_a$ the particle density background of species $a$, with $n_a = N_a$.

Using the spherical harmonics particle moment expansion of $f_b$, one can obtain the spherical particle moment expansion of the Rosenbluth potentials, $H(f_b)$ and $G(f_b)$ (Jorge et al. Reference Jorge, Frei and Ricci2019),

(3.15a)$$\begin{align} H(f_b) & = \frac{4 n_b}{v_{Tb}} \sum_{p=0}^\infty \sum_{j=0}^\infty\sum_{m=0}^j \frac{L_{jm}^p }{ \sigma_p^j} \frac{ \boldsymbol{\mathsf{Y}}^{p}(\hat{\boldsymbol{s}}_b) \boldsymbol{\cdot} \boldsymbol{\mathsf{M}}_b^{pj}}{2p+1} \left[ \frac{\textrm{I}_{-}^{2p + 2m +2}}{s_b^{p+1}} + s_b^p \textrm{I}_{+}^{2m +1} \right],\end{align}$$

(3.15b)$$\begin{align} G(f_b) & = 2 n_b v_{Tb} \sum_{p=0}^\infty \sum_{j=0}^\infty\sum_{l=0}^j \frac{L_{jl}^p}{\sigma_j^p}\frac{ \boldsymbol{\mathsf{M}}_b^{pj} \boldsymbol{\cdot} \boldsymbol{\mathsf{Y}}^p(\hat{\boldsymbol{s_b}}) }{(2p+1)} \nonumber\\ &\quad \times \left[ \frac{1}{(2p +3)} \left( \frac{\textrm{I}_-^{2p+2l +4}}{s_b^{p+1}} + s_b^{p+2} \textrm{I}_+^{2l+1}\right) -\frac{1}{(2p-1)} \left( \frac{\textrm{I}_-^{2p+2l +2}}{s_b^{p-1}} + s_b^{p} \textrm{I}_+^{2l+3}\right) \right], \end{align}$$

with $\hat {\boldsymbol {s}_a} = \boldsymbol {s}_a / s_a$ and the upper and lower incomplete gamma functions defined by $\textrm {I}_-^k = \int _0^x \,\textrm {d} s s^{(k-1)/2} \textrm {e}^{-s} / \sqrt {{\rm \pi} }$ and $\textrm {I}_+^k = \int _x^\infty \, \textrm {d} s s^{(k-1)/2} \textrm {e}^{-s} / \sqrt {{\rm \pi} }$, respectively. In deriving (3.15), the associated Laguerre polynomials, $L_j^{p+1/2}(s_b^2)$, are expanded using (3.9).

The velocity derivatives of the Rosenbluth potentials appearing in the field component of the Coulomb collision operator, (3.7), can be analytically evaluated using (3.15), and are given by

(3.16a)$$\begin{align} \frac{\partial}{\partial s_b} H_b(f_b) &= \frac{4 n_b}{v_{Tb}} \sum_{p=0}^\infty \sum_{j=0}^\infty \sum_{m=0}^j \frac{L_{jm}^p }{ \sigma_p^j} \frac{ \boldsymbol{\mathsf{Y}}^{p}(\hat{\boldsymbol{s}}_b) \boldsymbol{\cdot} \boldsymbol{\mathsf{M}}_b^{pj}}{2p+1} \left[ -(p+1)\frac{\textrm{I}_{-}^{2p + 2m +2}}{s_b^{p+2}} + p s_b^{p-1} \textrm{I}_{+}^{2m +1} \right],\end{align}$$

(3.16b)$$\begin{align} \frac{\partial^2}{\partial s_b^2} G_b(f_b ) & = 2 n_b v_{Tb} \sum_{p=0}^\infty \sum_{j=0}^\infty \sum_{l=0}^j \frac{L_{jl}^p}{\sigma_j^p}\frac{ \boldsymbol{\mathsf{M}}_b^{pj} \boldsymbol{\cdot} \boldsymbol{\mathsf{Y}}^p(\hat{\boldsymbol{s_b}}) }{(2p+1)} \left[ \frac{(p+1)(p+2)}{(2p +3)} \left( \frac{\textrm{I}_-^{2p+2l +4}}{s_b^{p+3}} + s_b^{p} \textrm{I}_+^{2l+1}\right) \right. \nonumber\\ &\quad \left. -\frac{p(p-1)}{(2p-1)} \left( \frac{\textrm{I}_-^{2p+2l +2}}{s_b^{p+1}} + s_b^{p-2} \textrm{I}_+^{2l+3}\right) \right], \end{align}$$

respectively. The spherical harmonic expansion of $f_a$, (3.8), can also be inserted into the test component of the Coulomb operator, given in (3.6). This finally yields the spherical harmonic moment expansion of the test and field components of the linearized Coulomb collision operator $C_{ab}$,

(3.17a)\begin{gather} C^T_{ab} = \sum_{p=0}^\infty \sum_{j=0}^\infty \frac{f_{Ma}}{\sigma_j^p} \boldsymbol{\mathsf{M}}_a^{pj}\boldsymbol{\cdot} \boldsymbol{\mathsf{Y}}^p(\hat{\boldsymbol{v} })\nu_{ab}^{Tpj}, \end{gather}

(3.17b)\begin{gather}C^F_{ab} = \sum_{p=0}^\infty \sum_{j=0}^\infty \frac{f_{Mb}}{\sigma_j^p} \boldsymbol{\mathsf{M}}_b^{pj}\boldsymbol{\cdot} \boldsymbol{\mathsf{Y}}^p(\hat{\boldsymbol{v}})\nu_{ab}^{Fpj}, \end{gather}

where the test and field speed functions, $\nu _{ab}^{Tpj}$ and $\nu _{ab}^{Fpj}$, are defined in Appendix A. We remark that the $v$ dependence of $C_{ab}^T$ and $C_{ab}^F$ is entirely contained in the speed functions, while the angular $(\theta ,\xi )$ dependence is isolated in the spherical harmonic tensor $\boldsymbol{\mathsf{Y}}^{p}(\hat {\boldsymbol {v}})$.

We now focus on evaluating the gyroaverage of the test component of the linearized Coulomb operator in (3.17a). Since the gyroaverage of $C_{ab}$, which appears in the linearized GK equation (2.1), is performed at constant $\boldsymbol {R}$, the spatial dependence of $\boldsymbol{\mathsf{M}}^{pj}_a(\boldsymbol {r})$ is first expressed as a function of $\boldsymbol {R}$. Since $\boldsymbol {r} = \boldsymbol {R} + \boldsymbol {\rho }_a$, the gyroaverage can be carried out in Fourier space by observing that $\boldsymbol{\mathsf{M}}_a^{pj}(\boldsymbol {r}) = \int \,\textrm {d} \boldsymbol {k} \boldsymbol{\mathsf{M}}_a^{pj} \textrm {e}^{\textrm {i} \boldsymbol {k} \boldsymbol {\cdot } \boldsymbol {\rho }_a } \textrm {e}^{\textrm {i} \boldsymbol {k} \boldsymbol {\cdot } \boldsymbol {R}}$ and by using the Jacobi–Anger identity to express the phase factor difference $\exp ({ \textrm {i} \boldsymbol {k} \boldsymbol {\cdot } \boldsymbol {\rho }_a }) = \sum _{n} i^n \textrm {J}_n(b_a\sqrt {x_a}) \textrm {e}^{\textrm {i}n \theta }$ (Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik2014). Thus, in Fourier space, the linearized test component of the Coulomb collision operator at $\boldsymbol {R}$ can be written as

(3.18)\begin{equation} C_{ab}^T = \sum_{p=0}^\infty \sum_{p=j}^\infty \sum_{n={-} \infty}^\infty \frac{F_{Ma}}{\sigma_j^p} \boldsymbol{\mathsf{M}}_a^{pj} \boldsymbol{\cdot} \boldsymbol{\mathsf{Y}}^p(\hat{\boldsymbol{v} }) i^n \textrm{J}_n(b_a \sqrt{x_a}) \textrm{e}^{\textrm{i} n \theta}\nu_{ab}^{Tpj}. \end{equation}

By expanding the spherical tensors, $\boldsymbol{\mathsf{Y}}^{p}$ appearing in (3.18), as shown in (3.11), one obtains

(3.19)\begin{gather} \left\langle C_{ab}^T \right\rangle_{\boldsymbol{R}} = \sum_{p=0}^\infty \sum_{j=0}^\infty \sum_{m={-}p}^p\sum_{n ={-}\infty}^{\infty} \frac{1}{\sigma_p^j} \sqrt{\frac{2 {\rm \pi}^{3/2} p!}{2^p (p+1/2) !}} F_{Ma} \nu_{ab}^{Tpj} i^n \textrm{J}_n(b_a\sqrt{x_a}) \nonumber\\ \times \left\langle \textrm{e}^{\textrm{i} n \theta} Y_p^m(\xi,\theta) \right\rangle_{\boldsymbol{R}} \boldsymbol{\mathsf{M}}_a^{pj} \boldsymbol{\cdot} \boldsymbol{e}^{pm}. \end{gather}

The scalar harmonic functions $Y_p^m(\xi ,\theta )$ can then be expanded in terms of associated Legendre polynomials $P_p^{m}(\xi )$, as indicated by (3.12). Thus, (3.19) can be written as

(3.20)\begin{gather} \left\langle C_{ab}^T \right\rangle_{\boldsymbol{R}} = \sum_{p=0}^\infty \sum_{j=0}^\infty \sum_{m={-}p}^p \frac{1}{\sigma_p^j} F_{Ma} \nu_{ab}^{Tpj} i^m \textrm{J}_m(\sqrt{x_a} b_a) \sqrt{\frac{{\rm \pi}^{1/2} p!}{2^p (p-1/2) !}} \nonumber\\ \times\sqrt{\frac{(p-|m|)!}{ (p+|m|)!}} P_p^{m}(\xi) \boldsymbol{\mathsf{M}}_a^{pj} \boldsymbol{\cdot} \boldsymbol{e}^{pm}. \end{gather}

We remark that the same gyroaverage procedure can be applied to the field component. In fact, the gyroaverage of the field component is given by (3.20), having replaced $\nu _{ab}^{Tpj}(v)$ by $\nu _{ab}^{Fpj}(v)$ and $\boldsymbol{\mathsf{M}}^{pj}_a$ by $\boldsymbol{\mathsf{M}}^{pj}_b$, i.e.

(3.21)\begin{gather} \left\langle C_{ab}^F \right\rangle_{\boldsymbol{R}} = \sum_{p=0}^\infty \sum_{j=0}^\infty \sum_{m={-}p}^p \frac{1}{\sigma_p^j} F_{Ma} \nu_{ab}^{Fpj} i^m \textrm{J}_m(\sqrt{x_a} b_a) \sqrt{\frac{{\rm \pi}^{1/2} p!}{2^p (p-1/2) !}} \nonumber\\ \times\sqrt{\frac{(p-|m|)!}{ (p+|m|)!}} P_p^{m}(\xi) \boldsymbol{\mathsf{M}}_b^{pj} \boldsymbol{\cdot} \boldsymbol{e}^{pm}. \end{gather}

3.1.2. Gyrokinetic Hermite–Laguerre expansion

Equations (3.20) and (3.21) are in a suitable form for the projection onto the Hermite–Laguerre basis, $\mathcal {C}_{ab}^{lk} = \mathcal {C}_{ab}^{lk}(\boldsymbol {k},t)$, defined (2.15). We focus first on the expression for the test component in (3.20). This yields

(3.22)\begin{equation} \mathcal{C}_{ab}^{Tlk} = \mathcal{C}_{ab}^{Tlk}(\boldsymbol{k},t) = \sum_{p=0}^\infty \sum_{j=0}^\infty \sum_{m={-}p}^p \frac{1}{\sigma_p^j} i^m \sqrt{\frac{{\rm \pi}^{1/2} p!}{2^p (p-1/2) !}} \sqrt{\frac{(p-|m|)!}{ (p+|m|)!}} \boldsymbol{\mathsf{M}}_a^{pj} \boldsymbol{\cdot} \boldsymbol{e}^{pm} \frac{ \textrm{I}_{abm}^{Tlkpj}}{\sqrt{2^l l!}}, \end{equation}

where we introduce,

(3.23)\begin{equation} \textrm{I}_{abm}^{Tlkpj} = \int \,\textrm{d} \boldsymbol{v} F_{Ma} \nu_{ab}^{Tpj} \textrm{J}_m(\sqrt{x_a} b_a) P_p^{m}(\xi) H_l(s_{{\parallel} a}) L_k(x_{ a}). \end{equation}

A similar expression can be obtained for the field component, $\mathcal {C}_{ab}^{Flk}$, having replaced $\nu _{ab}^{Tpj}$ by $\nu _{ab}^{Fpj}$ in (3.23), and $\boldsymbol{\mathsf{M}}_a^{pj}$ by $\boldsymbol{\mathsf{M}}_b^{pj}$ in (3.22). In order to express the Hermite–Laguerre projection of the linearized GK Coulomb collision operator, given in (3.22), in a closed analytical form in terms of gyromoments, the spherical harmonic particle moments, $\boldsymbol{\mathsf{M}}_a^{pj}$, defined in (3.14), must be written as a function of the gyromoments $N_a^{pj}$, and the velocity integral contained in (3.23) must be evaluated.

As a first step, we relate the spherical harmonics particle moments, $\boldsymbol{\mathsf{M}}_a^{pj}$ in (3.14), to the gyromoments $N_a^{pj}$ defined in (2.6). To proceed, we transform the velocity-space integral in (3.14) to a phase-space integral, i.e.

(3.24)\begin{equation} \boldsymbol{\mathsf{M}}_a^{pj}(\boldsymbol{r}) = \frac{1}{n_a }\int \,\textrm{d} \boldsymbol{v} \int \,\textrm{d} \boldsymbol{r}' \delta (\boldsymbol{r} - \boldsymbol{r}') f_a (\boldsymbol{r}', \boldsymbol{v}) \boldsymbol{\mathsf{Y}}^{p}(\boldsymbol{s}_a) L_j^{p+1/2}(s_a^2). \end{equation}

The coordinate transformation $\boldsymbol {r}' = \boldsymbol {r}' (\boldsymbol {R}, \mu , \theta ) = \boldsymbol {R} + \boldsymbol {\rho }_a(\mu , \theta )$ in (3.25) allows us to write

(3.25)\begin{equation} \boldsymbol{\mathsf{M}}_a^{pj}(\boldsymbol{r}) = \frac{1}{n_a } \int \,\textrm{d} \theta \,\textrm{d} \mu \,\textrm{d} v_\parallel \,\textrm{d} \boldsymbol{R} \frac{B}{m_a} \delta (\boldsymbol{r} - \boldsymbol{R} - \boldsymbol{\rho}_a ) f_a (\boldsymbol{r}', \boldsymbol{v}) \boldsymbol{\mathsf{Y}}^{p}(\boldsymbol{s}_a) L_j^{p+1/2}(s_a^2), \end{equation}

where $f_a (\boldsymbol {r}', \boldsymbol {v}) = f_a(\boldsymbol {r}'(\boldsymbol {R}, \mu , \theta ) , \boldsymbol {v}(\mu , v_\parallel , \theta ))$. We now aim to express the perturbed particle distribution function, $f_a(\boldsymbol {r}',\boldsymbol {v})$, in terms of the perturbed gyrocentre distribution function, ${g}_a(\boldsymbol {R},\mu , v_\parallel )$. The scalar invariance of the total particle and gyrocentre distribution functions yields

(3.26)\begin{equation} f_a(\boldsymbol{r}', \boldsymbol{v}) ={g}_a(\boldsymbol{R}, \mu, v_\parallel) + F_{Ma}(\boldsymbol{R}, \mu, v_\parallel) - f_{Ma}(\boldsymbol{r}', \boldsymbol{v}) . \end{equation}

Since the gyrocentre coordinates $(\boldsymbol {R}, \mu , v_\parallel , \theta )$ are obtained through a perturbative coordinate transformation from the particle phase-space $(\boldsymbol {r}, \boldsymbol {v})$, within a perturbation approach known as Lie-transform perturbation method (Cary Reference Cary1981; Brizard & Mishchenko Reference Brizard and Mishchenko2009; Frei et al. Reference Frei, Jorge and Ricci2020), the functional form of any function $F$, defined on $(\boldsymbol {R}, \mu , v_\parallel , \theta )$, can be expressed in terms of the functional form of a function $f$, defined on $(\boldsymbol {r}, \boldsymbol {v})$, by $f = \mathcal {T} F$, where the operator $\mathcal {T}$ is known as the pull-back operator (Frei et al. Reference Frei, Jorge and Ricci2020). Taking the case of $f = f_{Ma}$ and $F = F_{Ma}$ such that $f_{Ma} = \mathcal {T} F_{Ma}$ in (3.26), and at the leading order in $\epsilon \sim e \phi / T_e \ll 1$, one derives that the perturbed particle and gyrocentre distribution functions are related by

(3.27)\begin{equation} f_a(\boldsymbol{r}',\boldsymbol{v}) = {g}_a(\boldsymbol{R},\mu,v_\parallel) - \frac{q_a}{T_a} \left( \phi(\boldsymbol{r}') -\left\langle \phi \right\rangle_{\boldsymbol{R}}(\boldsymbol{R})\right) F_{Ma} . \end{equation}

Since $\phi - \left \langle \phi \right \rangle _{\boldsymbol {R}} \sim (k_\perp \rho _s)^2 e \phi / T_e$, the last term in (3.27) can be omitted in DK collision operators but, in general, it cannot be neglected in GK collisional theories, that consider $k_\perp \rho _s \sim 1$. The difference between the particle and gyrocentre Maxwellian distribution functions (Madsen Reference Madsen2013) leads to polarization effects in the collision operator (Xu & Rosenbluth Reference Xu and Rosenbluth1991; Dimits & Cohen Reference Dimits and Cohen1994). On the other hand, we remark that polarization effects are neglected in the nonlinear full-F GK Coulomb collision operator developed in Jorge et al. (Reference Jorge, Frei and Ricci2019), where it is assumed that $\epsilon _\nu e \phi / T_e \ll \epsilon ^2$ (with $\epsilon _\nu \sim \nu _{ii} / \varOmega _i \sim \epsilon ^2$), an ordering that we do not consider here (the same ordering is considered in the nonlinear GK Dougherty collision operator derived in Frei et al. (Reference Frei, Jorge and Ricci2020)). We also note that polarization effects are also neglected in the linearized GK Coulomb operator developed in Li & Ernst (Reference Li and Ernst2011). Using (3.27) in the definition of the spherical harmonics particle moments given in (3.14), we can write

(3.28)\begin{equation} \boldsymbol{\mathsf{M}}_a^{pj}(\boldsymbol{r}) = \boldsymbol{\mathsf{m}}_a^{pj}(\boldsymbol{r}) + \sum_{i=1}^{2}\boldsymbol{\mathsf{m}}^{pj}_{a \phi i}(\boldsymbol{r}), \end{equation}

with

(3.29a)\begin{gather} \boldsymbol{\mathsf{m}}_a^{pj}(\boldsymbol{r}) =\frac{1}{N_a} \int \,\textrm{d} \theta\,\textrm{d} \mu \,\textrm{d} v_\parallel \frac{B}{m_a} \left\langle {g}_a \boldsymbol{\mathsf{Y}}^p(\boldsymbol{s}_a) \right\rangle_{\boldsymbol{r}}^{\dagger} L_j^{p+1/2}(s_a^2), \end{gather}

(3.29b)\begin{gather}\boldsymbol{\mathsf{m}}^{pj}_{ a \phi 1}(\boldsymbol{r}) ={-} \frac{q_a}{N_a T_a}\int \,\textrm{d} \theta\, \textrm{d} \mu\, \textrm{d} v_\parallel \frac{B}{m_a} \left\langle \phi(\boldsymbol{R} + \boldsymbol{\rho}_a) \boldsymbol{\mathsf{Y}}^p(\boldsymbol{s}_a) \right\rangle^{\dagger}_{\boldsymbol{r}} L_j^{p+1/2}(s_a^2) F_{Ma}, \end{gather}

(3.29c)\begin{gather}\boldsymbol{\mathsf{m}}^{pj}_{ a \phi 2}(\boldsymbol{r}) = \frac{q_a}{N_a T_a}\int\, \textrm{d} \theta\, \textrm{d} \mu\, \textrm{d} v_\parallel \frac{B}{m_a} \left\langle \left\langle \phi \right\rangle_{\boldsymbol{R}}(\boldsymbol{R}) \boldsymbol{\mathsf{Y}}^p(\boldsymbol{s}_a) \right\rangle^{\dagger}_{\boldsymbol{r}} L_j^{p+1/2}(s_a^2) F_{Ma}, \end{gather}

where we introduce the adjoint gyroaverage operator at the particle position $\boldsymbol {r}$,

(3.30)\begin{equation} \left\langle \xi \right\rangle_{\boldsymbol{r}}^{\dagger}{=} \int_0^{2 {\rm \pi}} \frac{\textrm{d} \theta}{2 {\rm \pi}} \int \,\textrm{d} \boldsymbol{R} \delta(\boldsymbol{R} + \boldsymbol{\rho}_a - \boldsymbol{r}) \xi(\boldsymbol{R}, v_\parallel,\mu,\theta), \end{equation}

defined over any gyrocentre phase-space function $\xi$. We remark that the operator in (3.30) is the adjoint of the gyroaverage operator at constant $\boldsymbol {R}$, such that $\int \,\textrm {d} \boldsymbol {R} \left \langle \chi \right \rangle _{\boldsymbol {R}} \xi = \int \,\textrm {d} \boldsymbol {x} \left \langle \xi \right \rangle ^{\dagger} _{\boldsymbol {r}}\chi$, and plays a central role in polarization and magnetization effects in the field equations of GK theories (Frei et al. Reference Frei, Jorge and Ricci2020).

We now focus on evaluating the velocity integrals in $\boldsymbol{\mathsf{m}}_a^{pj}$ defined in (3.29a). As a first step, we compute $\left \langle {g}_a \boldsymbol{\mathsf{Y}}^p(\boldsymbol {s}_a) \right \rangle _{\boldsymbol {r}}^{\dagger}$ in Fourier space, such that

(3.31)$$\begin{align} \boldsymbol{\mathsf{m}}_a^{pj}(\boldsymbol{r}) \!& =\! \frac{1}{N_a} \int\, \textrm{d} \boldsymbol{k} \int\, \textrm{d} \theta \,\textrm{d} \mu \textrm{d} v_\parallel \frac{B}{m_a} {g}_a(\boldsymbol{k}, \mu, v_\parallel) \int_0^{2 {\rm \pi}} \frac{\textrm{d} \theta'}{2 {\rm \pi}} \exp({- \boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{\rho}_a'})\boldsymbol{\mathsf{Y}}^p(\boldsymbol{s}_a') L_j^{p+1/2}(s_a^2) \textrm{e}^{\textrm{i} \boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{R}} \nonumber\\ & = \frac{1}{N_a} \sum_{n ={-} \infty}^\infty (-\textrm{i})^n \int\, \textrm{d} \boldsymbol{k} \int\, \textrm{d} \theta \,\textrm{d} \mu \textrm{d} v_\parallel \frac{B}{m_a} {g}_a(\boldsymbol{k}, \mu, v_\parallel) \textrm{J}_n(b_a \sqrt{x_a} ) \nonumber\\ &\quad \times \int_0^{2 {\rm \pi}} \frac{\textrm{d} \theta'}{2 {\rm \pi}} \textrm{e}^{ \textrm{i} n \theta'}\boldsymbol{\mathsf{Y}}^p(\boldsymbol{s}_a') L_j^{p+1/2}(s_a^2) \textrm{e}^{\textrm{i} \boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{R}}, \end{align}$$

where we use the Jacobi–Anger identity and where the $\boldsymbol {\rho }_a'$ and $\boldsymbol {s}_a'$ notation indicate that these quantities have to be evaluated at $\theta '$. The velocity integrals in (3.31) can be performed by expanding the Bessel functions $\textrm {J}_n$ into associated Laguerre polynomials, thanks to (2.12), and by using the definition of the gyromoment, given in (2.6). A basis transformation between associated Legendre and Laguerre polynomials to Hermite and Laguerre polynomials is necessary to express $\boldsymbol{\mathsf{m}}_a^{pj}$ as a function of the gyromoments $N_a^{pj}$. This basis transformation (and its inverse) can be expressed as (Jorge et al. Reference Jorge, Frei and Ricci2019)

(3.32a)\begin{gather} s_s^p P_p^{m}(\xi) L_j^{p+1/2}(s_s^2) = \sum_{l=0}^{p + m + 2j} \sum_{k=0}^{j + \lfloor (p +m)/2 \rfloor} T_{pjm}^{lk} H_{l}(s_{{\parallel} s}) L_k(x_s) x_s^{m/2}, \end{gather}

(3.32b)\begin{gather}H_{l}(s_{{\parallel} s}) L_k(x_s) x_s^{m/2} = \sum_{p=0}^{l + m +2k} \sum_{j=0}^{k + \lfloor (l+m)/2 \rfloor} (T^{{-}1})^{pjm}_{lk} s_s^p P_p^{m}(\xi) L_j^{p+1/2}(s_s^2), \end{gather}

where the closed analytical expression of the coefficients $T_{pjm}^{lk}$ can be found in Jorge et al. (Reference Jorge, Frei and Ricci2019), and the coefficients $(T^{-1})^{pjm}_{lk}$ of the inverse basis transformation, (3.32b), are given by

(3.33)\begin{equation} (T^{{-}1})^{pjm}_{lk} = \frac{\sqrt{\rm \pi} 2^l l j!(p+1/2)}{(\,j+p+ 1/2)!} \frac{(p-m)!}{(p+m)!} \times \sum_{k'=0}^{\min(\,j + \lfloor (p+m)/2 \rfloor,m +k)} T_{pjm}^{lk'} d^{m}_{0kk'} [ l \le p+ m +2j], \end{equation}

with $[{\cdot }]$ the Iverson bracket, defined by $[A] = 1$ if $A$ is true, and $0$ otherwise. The basis transformation given by (3.32a) yields the gyromoment expansion

(3.34)\begin{equation} \boldsymbol{\mathsf{m}}_a^{pj} =\boldsymbol{\mathsf{m}}_a^{pj}(\boldsymbol{k}) = \sqrt{\frac{ {\rm \pi}^{1/2} p!}{2^p (p-1/2)!}}\sum_{m={-}p}^{p } (-\textrm{i})^m \sqrt{\frac{(p-|m|)!}{ (p+|m|)!}} \mathcal{M}_a^{pjm} \boldsymbol{e}^{pm}, \end{equation}

where we introduce

(3.35)\begin{equation} \mathcal{M}_a^{pjm} = \sum_{g=0}^{p+|m|+2j} \sum_{h=0}^{j + \lfloor (p+|m|)/2\rfloor}\sum_{n=0}^\infty \sum_{s_1=0}^{n+|m|+h} d_{nhs_1}^{|m|} T_{pj|m|}^{gh} \frac{ n!\sqrt{2^g g!}}{(n+|m|)!} \left( \frac{b_a}{2}\right)^{ |m|} \mathcal{K}_{n} (b_a)N_a^{gs_1}. \end{equation}

In (3.35), the numerical coefficients $d_{nks_1}^{|m|}$, which arise from the product between Laguerre and associated Laguerre polynomials,

(3.36)\begin{equation} L_n^{m}(x) L_k(x) x^{m} =\sum_{s=0}^{n+m +k} d^{m}_{nks} L_s(x), \end{equation}

are given by the closed formula (Jorge et al. Reference Jorge, Frei and Ricci2019),

(3.37)\begin{equation} d_{nks}^{m} =\sum_{n_1=0}^{n} \sum_{k_1=0}^k \sum_{s_1=0}^s L_{kk_1}^{{-}1/2} L_{nn_1}^{m-1/2} L_{ss_1}^{{-}1/2} (n_1 + k_1 + s_1 + m)!, \end{equation}

with the coefficients $L_{jl}^{p}$ defined in (3.10).

We now turn to the polarization terms in (3.28), i.e. $\boldsymbol{\mathsf{m}}_{a \phi 1}^{pj}$ and $\boldsymbol{\mathsf{m}}_{a \phi 2}^{pj}$ defined in (3.29b) and (3.29c), respectively. Following the same steps considered in the evaluation of $\boldsymbol{\mathsf{m}}_s^{pj}$, we first compute the terms in $\left \langle \phi (\boldsymbol {R} + \boldsymbol {\rho }_a) \boldsymbol{\mathsf{Y}}^p(\boldsymbol {s}_a) \right \rangle _{\boldsymbol {r}}^{\dagger}$ and $\left \langle \left \langle \phi \right \rangle _{\boldsymbol {R}}(\boldsymbol {R}) \boldsymbol{\mathsf{Y}}^p(\boldsymbol {s}_a) \right \rangle _{\boldsymbol {r}}^{\dagger}$. Using the coordinate transformation $\boldsymbol {r} = \boldsymbol {R} + \boldsymbol {\rho }_a$, we obtain

(3.38)\begin{equation} \left\langle \phi(\boldsymbol{R} + \boldsymbol{\rho}_a) \boldsymbol{\mathsf{Y}}^p(\boldsymbol{s}_a) \right\rangle_{\boldsymbol{r}}^{\dagger}{=} \phi(\boldsymbol{r}) \left\langle \boldsymbol{\mathsf{Y}}^p(\boldsymbol{s}_a) \right\rangle_{\boldsymbol{r}}^{\dagger}, \end{equation}

while using the Jacobi–Anger identity

(3.39)\begin{equation} \left\langle \left\langle \phi \right\rangle_{\boldsymbol{R}}(\boldsymbol{R}) \boldsymbol{\mathsf{Y}}^p(\boldsymbol{s}_a) \right\rangle_{\boldsymbol{r}}^{\dagger}{=} \int\, \textrm{d} \boldsymbol{k} \sum_{n ={-} \infty}^\infty (-\textrm{i})^n \phi(\boldsymbol{k}) \textrm{J}_0(b_a \sqrt{x_a}) \textrm{J}_n( b_s \sqrt{x_s}) \left\langle \textrm{e}^{\textrm{i}n \theta} \boldsymbol{\mathsf{Y}}^p(\boldsymbol{s}_a) \right\rangle_{\boldsymbol{r}}^{\dagger} \textrm{e}^{\textrm{i} \boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{r}}. \end{equation}

Equations (3.38) and (3.39) allow us to perform analytically the velocity integrals contained in (3.29b) and (3.29c). We first use the orthogonality relation, given in (3.13), and expand $\boldsymbol{\mathsf{Y}}^p$ in associated Legendre polynomials with (3.11), while expressing the Bessel functions in associated Laguerre polynomials according to (2.12). Finally, we use the basis transformation given in (3.32a). This yields

(3.40a)\begin{gather} \boldsymbol{\mathsf{m}}_{a \phi 1}^{pj} ={-} \frac{q_a \phi}{T_a} \delta_{p}^0 \delta_j^0 \boldsymbol{e}^{p0}, \end{gather}

(3.40b)\begin{gather}\boldsymbol{\mathsf{m}}_{a \phi 2}^{pj} = \frac{q_a \phi}{T_a} \sum_{m'={-}p}^{p }(-\textrm{i})^{m'} \sqrt{\frac{ {\rm \pi}^{1/2} p!}{2^p (p-1/2)!}} \sqrt{\frac{(p-|m'|)!}{ (p+|m'|)!}} \mathcal{\varPi}^{pjm'}_a \boldsymbol{e}^{pm'}, \end{gather}

where we introduce the coefficients

(3.41)\begin{equation} \mathcal{\varPi}^{pjm'}_a = \sum_{h=0}^{ j + \lfloor(p+|m'| )/2 \rfloor} \sum_{s=0}^\infty \sum_{s'=0}^{s+|m'| +h} T^{0h}_{pj|m'|}\left(\frac{b_a}{2}\right)^{|m'|}\frac{s! \mathcal{K}_{s}(b_a)\mathcal{K}_{s'}(b_a)}{(s +|m'|)!} d^{|m'|}_{shs'}. \end{equation}

The electrostatic potential $\phi$, appearing in (3.40), can be self-consistently expressed as a function of $N_a^{pj}$ using the GK quasineutrality condition, (2.18).

As a last step, we compute the velocity integral in $\textrm {I}_{abm}^{Tlkpj}$ defined in (3.23). To proceed, we first expand the Bessel functions $\textrm {J}_m$ in terms of associated Laguerre polynomials, according to (2.12), and use the inverse basis transformation given in (3.32b), to obtain

(3.42)$$\begin{align} \textrm{I}_{abm}^{Tlkpj} &= \sum_{n=0}^\infty \sum_{f =0}^{n+k} \sum_{g=0}^{l + m +2f} \sum_{h=0}^{f + \lfloor (l+m)/2 \rfloor} (T^{{-}1})^{ghm}_{lf} \bar{d}^{m}_{nkf}\nonumber\\ &\quad \times \int \,\textrm{d} \boldsymbol{v} F_{Ma} s_a^g L_h^{g+1/2}(s_a^2) \nu_{ab}^{Tpj} P_p^{m}(\xi) P_g^{m}(\xi) , \end{align}$$

where the product of associated Laguerre and Laguerre polynomials is expressed as a single series of Laguerre polynomials thanks to (Jorge et al. Reference Jorge, Frei and Ricci2019),

(3.43)\begin{equation} L_n^{m}(x_a) L_k(x_a) = \sum_{f =0}^{n+k} \bar{d}^{m}_{nkf} L_f(x_a), \end{equation}

where

(3.44)\begin{equation} \bar{d}^{m}_{nkf} = \sum_{n_1 = 0}^n \sum_{k_1 =0}^k \sum_{f_1=0}^f L_{kk _1}^{{-}1/2} L_{nn_1}^{m-1/2} L_{f f_1}^{{-}1/2} (n_1 + k_1 + f_1)!. \end{equation}

Finally, using the orthogonality relations (Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik2014)

(3.45)\begin{equation} \int_{{-}1}^{1} \,\textrm{d} \xi P_p^{m}(\xi) P_g^{m}(\xi) = \frac{2 \delta_g^p}{(2p +1)} \frac{(p + m)!}{(p - m)!} \end{equation}

and $\boldsymbol {e}^{pm} \boldsymbol {\cdot }\boldsymbol {e}^{pm'} = (-1)^m \delta _{-m}^{m'}$ (Snider Reference Snider2017), the Hermite–Laguerre expansion of the linearized GK Coulomb collision operator, $\mathcal {C}_{ab}^{lk} = \mathcal {C}_{ab}^{lk}(\boldsymbol {k},t)$, is obtained,

(3.46)\begin{equation} \mathcal{C}_{ab}^{lk} = \mathcal{C}^{Tlk}_{ab0} +\mathcal{C}^{Flk}_{ab0} + \sum_{i=1}^{2} \left( \mathcal{C}^{Tlk}_{ab\phi i} + \mathcal{C}^{Flk}_{ab\phi i} \right), \end{equation}

with the Hermite–Laguerre expansion of the test and field components, $\mathcal {C}^{Tlk}_{ab0}$ and $\mathcal {C}^{Flk}_{ab0}$, given by

(3.47)$$\begin{align} \mathcal{C}_{ab0}^{Tlk} & =\sum_{p=0}^\infty \sum_{j=0}^\infty \sum_{n=0}^\infty \sum_{f=0}^{k+n} \sum_{m=0}^{p} \sum_{t=0}^{f + \lfloor (l +m)/2\rfloor} \ \frac{ a_m }{\sigma_p^j} \frac{2^{p} (p!)^2 }{(2 p)! } \frac{ \bar{d}_{nkf}^{m}}{(2p+1)}\nonumber\\ &\quad \times \frac{ \left( T ^{{-}1} \right)^{ptm}_{lf} }{\sqrt{2^l l!}} \frac{ n! \mathcal{K}_{n}(b_a) }{(n + m)! } \left( \frac{b_a}{2}\right)^{ m} \mathcal{M}_a^{pjm} \bar{\nu}_{ab}^{Tpjt} \left[ p \leq l + m + 2f\right] ,\end{align}$$

(3.48)$$\begin{align} \mathcal{C}_{ab0}^{Flk} & =\sum_{p=0}^\infty \sum_{j=0}^\infty \sum_{n=0}^\infty \sum_{f=0}^{k+n} \sum_{m=0}^{p} \sum_{t=0}^{f + \lfloor (l +m)/2\rfloor} \frac{ a_m }{\sigma_p^j} \frac{2^{p} (p!)^2 }{(2 p)! } \frac{ \bar{d}_{nkf}^{m}}{(2p+1)}\nonumber\\ &\quad \times \frac{ \left( T ^{{-}1} \right)^{ptm}_{lf} }{\sqrt{2^l l!}} \frac{ n! \mathcal{K}_{n}(b_a) }{(n + m)! } \left( \frac{b_a}{2}\right)^{ m} \mathcal{M}_b^{pjm} \bar{\nu}_{ab}^{Fpjt} \left[ p \leq l + m + 2f\right], \end{align}$$

respectively, and the associated polarization contributions,

(3.49a)$$\begin{align} \mathcal{C}^{Tlk}_{ab \phi 1 } & ={-} \frac{q_a}{T_a} \phi \sum_{n=0}^\infty \sum_{s=0}^{n+k} \sum_{h=0}^{s +\lfloor l/2 \rfloor}d_{nks}^0 \mathcal{K}_{n}(b_a) \frac{\left( T^{{-}1} \right)^{0h0}_{ls}}{\sqrt{2^l l!}} \bar{\nu}_{ab}^{T00h}, \end{align}$$

(3.49b)$$\begin{align} \mathcal{C}^{Flk}_{ab \phi 1 } & ={-} \frac{q_b}{T_b} \phi \sum_{n=0}^\infty \sum_{s=0}^{n+k} \sum_{h=0}^{s +\lfloor l/2 \rfloor}d_{nks}^0 \mathcal{K}_{n}(b_a) \frac{\left( T^{{-}1} \right)^{0h0}_{ls}}{\sqrt{2^l l!}} \bar{\nu}_{ab}^{F00h},\end{align}$$

(3.49c)$$\begin{align} \mathcal{C}_{ab \phi 2 }^{T lk} & = \frac{q_a}{T_a } \phi \sum_{p=0}^\infty \sum_{j=0}^\infty \sum_{m=0}^{p} \sum_{n=0}^\infty\sum_{f =0}^{n+k} \sum_{t=0}^{f + \lfloor (l+m)/2 \rfloor} \frac{a_m}{\sigma_p^j} \frac{2^{p} (p!)^2 }{(2 p)! } \frac{ 1}{( 2p+1)} \frac{\mathcal{\varPi}_a^{pjm}}{\sqrt{2^l l!}} \nonumber\\ &\quad \times \bar{d}^{m}_{nkf} \frac{n!\mathcal{K}_{n}(b_a)}{(n+m)!}\left( \frac{b_a}{2}\right)^{m} (T^{{-}1})^{ptm}_{lf} \bar{\nu}_{ab}^{Tpjt} \left[ p \leq l + m + 2f\right],\end{align}$$

(3.49d)$$\begin{align} \mathcal{C}_{ab \phi 2}^{Flk} & = \frac{q_b}{T_b } \phi \sum_{p=0}^\infty \sum_{j=0}^\infty \sum_{m=0}^{p} \sum_{n=0}^\infty\sum_{f =0}^{n+k} \sum_{t=0}^{f + \lfloor (l+m)/2 \rfloor} \frac{a_m}{\sigma_p^j} \frac{2^{p} (p!)^2 }{(2 p)! } \frac{ 1}{( 2p+1)} \frac{\mathcal{\varPi}_b^{pjm}}{\sqrt{2^l l!}} \nonumber\\ &\quad \times \bar{d}^{m}_{nkf} \frac{n!\mathcal{K}_{n}(b_a)}{(n+m)!}\left( \frac{b_a}{2}\right)^{m} (T^{{-}1})^{ptm}_{lf} \bar{\nu}_{ab}^{Fpjt} \left[ p \leq l + m + 2f\right] , \end{align}$$

where $a_m = 2$ for $m \geq 1$ while $a_0 =1$. In (3.47)–(3.49), we introduce the velocity- integrated speed functions

(3.50)\begin{gather} \bar{\nu}_{ab}^{Tpjt} = \sum_{d =0}^t L_{td}^p\bar{\nu}_{*ab}^{Tpjd}, \end{gather}

(3.51)\begin{gather}\bar{\nu}_{ab}^{Fpjt} = \sum_{d =0}^t L_{td}^p\bar{\nu}_{*ab}^{Fpjd}, \end{gather}

where $\bar {\nu }_{*ab}^{Tpjd} = \int \,\textrm {d} \boldsymbol {v} s_a^{p+2d} \nu _{ab}^{Tpj} F_{Ma}$ and $\bar {\nu }_{*ab}^{Fpjd} = \int \,\textrm {d} \boldsymbol {v} s_a^{p+2d} \nu _{ab}^{Fpj} F_{Ma}$. Closed analytical expressions for $\bar {\nu }_{ab*}^{Tpjd}$ and $\bar {\nu }_{ab*}^{Fpjd}$ are reported in Appendix A. We remark that the Hermite–Laguerre expansions of the test and field components, (3.47) and (3.48), respectively, can be further written in terms of the gyromoments using the definition of $\mathcal {M}_a^{pjm}$ and $\mathcal {M}_b^{pjm}$ given in (3.35). The remaining quantities in (3.47)–(3.49) are numerical coefficients.

Equation (3.46), with (3.47)–(3.49), define the gyromoment expansion of the linearized GK Coulomb collision operator. This operator is valid for arbitrary mass and temperature ratios. It retains FLR effects associated with the difference between the particle and gyrocentre positions, $\boldsymbol {r}$ and $\boldsymbol {R}$, at arbitrary perpendicular wavenumber. These enter explicitly in the linearized GK Coulomb collision operator through terms proportional to $b_s^m \mathcal {K}_{n}(b_s)$ in the test and field components. Coupled with the quasineutrality equation, (2.18), and the definitions of the different numerical factors appearing in (3.47) and (3.48), $\mathcal {C}_{ab}^{lk}$ can be numerically implemented, and be used as a collision operator model in the linearized gyromoment hierarchy, (2.14), to describe collisions between like and unlike species.

3.1.3. Drift-kinetic limit of Coulomb collision operator

The DK limit of the GK linearized Coulomb collision operator, given in (3.46), can be deduced by neglecting the FLR effects in the test and field components as well as the polarization terms. For instance, in the zero gyroradius limit, (3.27) reduces to $f_a = {g}_a$. In addition, the gyromoment expansion of the spherical harmonic particle moments, $\boldsymbol{\mathsf{M}}_a^{pj}$, given in (3.28), reduces to

(3.52)\begin{equation} \boldsymbol{\mathsf{M}}_a^{pj} = \sqrt{\frac{ {\rm \pi}^{1/2} p!}{2^p (p-1/2)!}} \mathcal{M}_a^{pj} \boldsymbol{e}^{p0}, \end{equation}

where

(3.53)\begin{equation} \mathcal{M}_a^{pj} = \sum_{g=0}^{p+2j} \sum_{h=0}^{j + \lfloor p/2\rfloor} T_{pj}^{gh} \sqrt{2^g g!} N_a^{gh}. \end{equation}

Neglecting also FLR effects in (3.20), such that $\left \langle C_{ab}(\boldsymbol {r},\boldsymbol {v}) \right \rangle _{\boldsymbol {R}} \simeq \left \langle C_{ab}(\boldsymbol {R},\boldsymbol {v}) \right \rangle _{\boldsymbol {R}}$, and projecting the resulting expression onto the Hermite–Laguerre basis, we obtain the linearized DK Coulomb collision operator with the test and field components given by

(3.54)\begin{gather} \mathcal{C}_{ab}^{Tlk} = \sum_{j=0}^{\infty} \sum_{p=0}^{l+2k} \sum_{t=0}^{k + \lfloor l/2\rfloor} \sum_{d=0}^t \frac{ 2^{p} (p!)^2 }{(2 p)! \sigma_p^j} \frac{ L_{td}^p}{(2p+1)} \frac{\left( T ^{{-}1} \right)^{pt}_{lk}}{\sqrt{2^l l!}} \mathcal{M}_a^{pj} \bar{\nu}_{*ab}^{Tpjd}, \end{gather}

(3.55)\begin{gather}\mathcal{C}_{ab}^{Flk} = \sum_{j=0}^{\infty} \sum_{p=0}^{l+2k} \sum_{t=0}^{k + \lfloor l/2\rfloor} \sum_{d=0}^t \frac{ 2^{p} (p!)^2 }{(2 p)! \sigma_p^j} \frac{ L_{td}^p}{(2p+1)} \frac{\left( T ^{{-}1} \right)^{pt}_{lk}}{\sqrt{2^l l!}} \mathcal{M}_b^{pj} \bar{\nu}_{*ab}^{Fpjd}, \end{gather}

respectively. Equations (3.54) can also be obtained by taking explicitly the zero gyroradius limit of (3.46). We remark that the DK linearized Coulomb collision operator obtained here is equivalent to the DK collision operator derived in Jorge, Ricci & Loureiro (Reference Jorge, Ricci and Loureiro2018).

3.2. Gyrokinetic Sugama collision operator

The Sugama collision operator model (Sugama et al. Reference Sugama, Watanabe and Nunami2009) is an advanced model collision operator that approximates the full linearized Coulomb collision operator such that it conserves particle, momentum and energy. It also satisfies the H-theorem and is subject to the self-adjoint relations for the test and field components, i.e.

(3.56a)\begin{gather} \int \,\textrm{d} \boldsymbol{v} \frac{f_a}{f_{Ma}} C_{ab}^T \left( g_a \right) = \int \,\textrm{d} \boldsymbol{v} \frac{g_a}{f_{Ma}} C_{ab}^T \left( f_a \right), \end{gather}

(3.56b)\begin{gather}T_a \int \,\textrm{d} \boldsymbol{v} \frac{f_a}{f_{Ma}} C_{ab}^F \left( f_b \right) = T_b \int \,\textrm{d} \boldsymbol{v} \frac{f_b}{f_{Mb}} C_{ba}^F \left( f_a \right). \end{gather}

The Sugama collision operator model extends the Abel et al. (Reference Abel, Barnes, Cowley, Dorland and Schekochihin2008) collision operator to the general case of collisions between different species, and reduces to it in the case of like species, i.e. $a = b$. We remark that, because the difference with the Sugama and Coulomb operator is expected to increase at high collisionality, the Sugama operator was recently improved and extended to the regime of high-collisionality in Sugama et al. (Reference Sugama, Matsuoka, Satake, Nunami and Watanabe2019). In this work, we focus on the original Sugama operator, presented in Sugama et al. (Reference Sugama, Watanabe and Nunami2009), but remark that the gyromoment expansion can be applied also to the improved Sugama by following the same analytical procedure detailed below.

In the particle coordinates $(\boldsymbol {r}, v,\xi , \theta )$, the test component of the Sugama collision operator, $C_{ab}^T(f_a) = C_{ab}^T(\boldsymbol {r},\boldsymbol {v})$, is defined by

(3.57)\begin{equation} C_{ab}^{T}(f_a) = C_{ab}^S(f_a)+ \sum_{n=1}^3 X^n_{ab}, \end{equation}

where

(3.58)\begin{equation} C_{ab}^S(f_a ) ={-} \nu_{ab}^D(v)\mathcal{L}^2 f_a + \frac{1}{v^2} \frac{\partial}{\partial v} \left( \nu_{ab}^{{\parallel}}(v) v^4 F_{aM} \frac{\partial}{\partial v} \left( \frac{f_a}{f_{aM}}\right)\right) \end{equation}

and

(3.59)\begin{gather} X_{ab}^1 = (\theta_{ab} -1 ) f_{Ma} \left[ \frac{m_a}{T_a} \frac{1}{n_a} \int \,\textrm{d} \boldsymbol{v}' \boldsymbol{v}' C_{ab}^S(f_a) \boldsymbol{\cdot} \boldsymbol{v} + \left( s_a^2 - \frac{3}{2}\right) \frac{1}{n_a} \int \,\textrm{d}\boldsymbol{v} \frac{2}{3} s_a^2 C_{ab}^S(f_a) \right], \end{gather}

(3.60)\begin{gather}X_{ab}^2 = (\theta_{ab} -1 ) \left[ \frac{m_a}{T_a} \boldsymbol{u}_a(f_a) \boldsymbol{\cdot} C_{ab}^S( f_{Ma} \boldsymbol{v}) + \frac{\delta T_a (f_a)}{T_a} C_{ab}^S(F_{aM}s_a^2) \right], \end{gather}

(3.61)\begin{gather}X_{ab}^3 ={-} \frac{4 \nu_{ab}}{3 \sqrt{\rm \pi}}f_{Ma} \frac{\chi(\theta_{ab} -1 )^2}{ (1 + \chi^2)^{1/2}} \left[ \frac{m_a}{T_a} \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{u}_a(f_a) + \frac{\delta T_a(f_a)}{T_a}\left( s_a^2 -\frac{3}{2}\right) \frac{2 }{ (1 + \chi^2)} \right], \end{gather}

with $\boldsymbol {u}_a(f_a) = \int \,\textrm {d} \boldsymbol {v} \boldsymbol {v} f_a / n_a$, $\delta T_a(f_a) = T_a \int \,\textrm {d} \boldsymbol {v} f_a ( 2s_a^2/3 - 1 ) /n_a$, $\theta _{ab} = \sqrt { \tau + \chi ^2/(1 + \chi ^2)}$, $\tau = T_a / T_b$, $\chi = v_{Ta} / v_{Tb} = \sqrt {\tau / \sigma }$ and $\sigma = m_a / m_b$. In (3.57), we introduce the pitch-angle scattering (also denoted as deflection) and energy diffusion frequencies defined by $\nu _{ab}^D(v ) = \nu _{ab} [ \mathrm {erf}(s_b) - \varPhi (s_b)]/s_a^3$ and $\nu _{ab}^{\parallel }(v) = 2 \nu _{ab} \varPhi (s_b) / s_a^3$, respectively, with $\varPhi (x) = [ \mathrm {erf}(x) - x \mathrm {erf}'(x)]/(2 x^2)$ the Chandrasekhar function. The closed analytical expressions for $C_{ab}^{S}(f_{Ma} \boldsymbol {v} )$ and $C_{ab}^{S}( s_a^2 f_{Ma})$ appearing in (3.60) are given by (Sugama et al. Reference Sugama, Watanabe and Nunami2009)

(3.62a)\begin{gather} C_{ab}^{S}(f_{Ma} \boldsymbol{v} ) ={-}\left( 1 + \chi^2\right) f_{Ma} \nu_{ab}^{{\parallel}}(v) s_a^2 \boldsymbol{v}, \end{gather}

(3.62b)\begin{gather}C_{ab}^{S}( s_a^2 f_{Ma}) ={-} \nu_{ab} f_{Ma} \frac{2}{\chi^2 s_a} \left[ \mathrm{erf}(s_b) - s_b \left( 1 + \chi^2 \right) \mathrm{erf}'(s_b)\right]. \end{gather}

We remark that $C_{ab}^{S}$, given in (3.58), is actually equivalent to the test component of the linearized Coulomb operator in (3.6), for the case of like-species collisions. More precisely, $C_{ab}^{S}$ in (3.58) is obtained by using the analytical equilibrium Rosenbluth potentials, $H(f_{Mb})$ and $G(f_{Mb})$, in (3.6) and neglecting the term proportional to $(1 - T_a/T_b)$.

The field component of the Sugama collision operator is constructed from the test component $C_{ab}^T$, in (3.57), such that the local conservation laws, (3.2)–(3.4), and the self-adjoint relations, (3.56), are satisfied. Sugama et al. (Reference Sugama, Watanabe and Nunami2009) shows that this is the case if the field component, $C_{ab}^F(f_b) = C_{ab}^F(\boldsymbol {r},\boldsymbol {v})$, is given by

(3.63)\begin{equation} C_{ab}^{F}(f_b) ={-} \frac{m_a \boldsymbol{V}_{ab}(f_b)}{T_a} \boldsymbol{\cdot} C_{ab}^{T}\left( f_{Ma} \boldsymbol{v}\right) - W_{ab}(f_b) C_{ab}^{T}( f_{Ma} s_a^2), \end{equation}

where $C_{ab}^{T}$ is test component of the Sugama collision operator, given in (3.57). From (3.57), one has

(3.64a)$$\begin{align} C_{ab}^{T}( f_{Ma} \boldsymbol{v})& ={-} \theta_{ab} \left( 1 + \chi^2\right) \frac{4 \nu_{ab}}{3 \sqrt{\rm \pi}} f_{Ma} \boldsymbol{v} \left[ \frac{3 \sqrt{\rm \pi}}{2} \frac{\varPhi(s_b)}{s_a} + \frac{\chi (\theta_{ab}-1)}{(1 + \chi^2)^{3/2}}\right], \end{align}$$

(3.64b)$$\begin{align} C_{ab}^{T}(f_{Ma} s_a^2) &={-} \frac{4 \nu_{ab}}{3 \sqrt{\rm \pi}} \theta_{ab} f_{Ma} \left[ \frac{3 \sqrt{\rm \pi}}{2} \frac{1}{\chi^2 s_a} \left\{ \mathrm{erf}(s_b) - s_b \mathrm{erf}'(s_b) (1 + \chi^2) \right\} \right. \nonumber\\ & \quad \left. + \frac{2 \chi (\theta_{ab} -1 )}{(1 + \chi^2)^{3/2}} \left( s_a^2 - \frac{3}{2}\right)\right]. \end{align}$$

The quantities, $\boldsymbol {V}_{ab}[f_b]$ and $W_{ab}[f_b]$, appearing in (3.63), are defined by

(3.65)\begin{gather} \boldsymbol{V}_{ab}[f_b] = \frac{m_b}{\gamma_{ab}} \int \,\textrm{d} \boldsymbol{v} \frac{f_b}{f_{Mb}} C_{ba}^{T}\left( f_{Mb} \boldsymbol{v} \right), \end{gather}

(3.66)\begin{gather}W_{ab}[f_b] = \frac{T_b}{\eta_{ab}} \int \,\textrm{d} \boldsymbol{v} \frac{f_b}{f_{Mb}} C_{ba}^{T}(f_{Mb} s_b^2), \end{gather}

where $\gamma _{ab}$ and $\eta _{ab}$ are given by

(3.67)\begin{gather} \gamma_{ab} ={-} \frac{4 \nu_{ab}}{3 \sqrt{\rm \pi}} \frac{n_a m_a \chi}{(1 + \chi^2)^{3/2}} \left( \frac{T_a}{T_b} + \chi^2\right), \end{gather}

(3.68)\begin{gather}\eta_{ab} ={-} \frac{4 \nu_{ab}}{3 \sqrt{\rm \pi}} \frac{3 \chi n_a T_a}{(1 + \chi^2)^{5/2}} \left( \frac{T_a}{T_b} + \chi^2\right). \end{gather}

We remark that $C_{ba}^{T}( f_{Mb} \boldsymbol {v})$ and $C_{ba}^{T}(f_{Mb} s_b^2)$ correspond to $C_{ab}^T(f_{Ma} \boldsymbol {v})$ and $C_{ab}^{T}(f_{Ma} s_a^2)$, given in (3.64), with $a$ replaced by $b$.

In order to express the linearized GK Sugama collision operator in the gyrocentre phase space coordinates $(\boldsymbol {R}, \mu , v_\parallel , \theta )$, Sugama et al. (Reference Sugama, Watanabe and Nunami2009) considers the non-adiabatic part $h_a$ of the perturbed particle distribution function $f_a$, which is defined by

(3.69)\begin{equation} f_a(\boldsymbol{r}, \boldsymbol{v}) = h_a(\boldsymbol{R},\mu,v_\parallel) - \frac{q_a }{T_a}\phi (\boldsymbol{r}) f_{Ma}. \end{equation}

Then, to evaluate the velocity derivatives contained (3.57) on $h_a$ and to perform the gyroaverage, the GK transformation (Xu & Rosenbluth Reference Xu and Rosenbluth1991; Abel et al. Reference Abel, Barnes, Cowley, Dorland and Schekochihin2008)

(3.70)\begin{equation} \mathcal{C}^T_{ab} = \left\langle \exp({ \textrm{i} \boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{\rho}_a}) C^T_{ab}( h_a \exp({- \textrm{i} \boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{\rho}_a})) \right\rangle_{\boldsymbol{R}} \end{equation}

is used. This yields the test component of the GK Sugama operator $\mathcal {C}_{ab}^T = \mathcal {C}_{ab}^T(\boldsymbol {k}, \mu , v_\parallel )$,

(3.71)\begin{equation} \mathcal{C}_{ab}^{T} = C_{ab }^{S}(h_a) + R^\perp_{ab} + \sum_{n=1}^3 \bar{X}_{ab}^n, \end{equation}

with $C_{ab }^{S}$ defined in (3.58) and

(3.72)\begin{gather} R_{ab}^\perp={-} \frac{k_\perp^2}{4 \varOmega_a^2} h_a \left[ \nu_{ab}^D(v) \left( 2 v_\parallel^2 + v_\perp^2 \right) + \nu_{ab}^{{\parallel}}(v) v_\perp^2 \right], \end{gather}

(3.73)\begin{gather}\bar{X}_{ab}^1 = \left( \theta_{ab}-1\right) F_{Ma} \left[ \textrm{J}_0 s_{{\parallel} a} \mathcal{R}_{ab}^{{\parallel} 1} + \textrm{J}_1 \sqrt{x_a} \mathcal{R}_{ab}^{{\parallel} 2} + \textrm{J}_0 \left(s_a^2 - \frac{3}{2} \right) \mathcal{R}_{ab}^{{\parallel} 3} \right] , \end{gather}

(3.74)\begin{gather}\bar{X}_{ab}^2 =\left( \theta_{ab}-1\right) \left[ v_{Ta} C_{ab}^S(F_{Ma} m_a v_\parallel{/}T_a ) \left(\textrm{J}_0 \bar{u}_{{\parallel} a} + \textrm{J}_1 \frac{v_\perp}{v_\parallel} \bar{u}_{{\perp} a } \right)+ \textrm{J}_0 C_{ab}^S(F_{Ma} s_a^2)\bar{T}_a \right] , \end{gather}

(3.75)\begin{gather}\bar{X}_{ab}^3 ={-} F_{Ma} \frac{8 \nu_{ab} (\theta_{ab} -1 )^2 \chi}{ 3 \sqrt{\rm \pi}(1 + \chi^2)^{1/2}} \left[ \textrm{J}_0 s_{{\parallel} a} \bar{u}_{{\parallel} a} + \textrm{J}_1 \sqrt{x_a} \bar{u}_{{\perp} a} + \frac{\textrm{J}_0}{1 + \chi^2} \left(s_a^2 - \frac{3}{2} \right) \bar{T}_a \right], \end{gather}

with

(3.76a)\begin{gather} \mathcal{R}_{ab}^{{\parallel} 1} = \frac{1}{n_a}\int \,\textrm{d} \boldsymbol{v} \textrm{J}_0 \frac{h_a}{F_{Ma}} C_{ab}^S(F_{Ma} m_a v_\parallel{/}T_a), \end{gather}

(3.76b)\begin{gather}\mathcal{R}_{ab}^{{\parallel} 2} =\frac{1}{n_a} \int \,\textrm{d} \boldsymbol{v} \textrm{J}_1 \frac{h_a}{F_{Ma}} \frac{v_\perp}{v_\parallel} C_{ab}^S(F_{Ma} m_a v_\parallel{/}T_a), \end{gather}

(3.76c)\begin{gather}\mathcal{R}_{ab}^{{\parallel} 3} = \frac{1}{n_a}\int \,\textrm{d} \boldsymbol{v} \textrm{J}_0 \frac{h_a}{F_{Ma}} \frac{2}{3} C_{ab}^S(F_{Ma} s_a^2) \end{gather}

and

(3.77a)\begin{gather} \bar{u}_{{\parallel} a} = \frac{1}{n_a}\int \,\textrm{d} \boldsymbol{v} \textrm{J}_0 h_a s_{{\parallel} a}, \end{gather}

(3.77b)\begin{gather}\bar{u}_{{\perp} a} = \frac{1}{n_a}\int \,\textrm{d} \boldsymbol{v} \textrm{J}_1 \sqrt{x_a} h_a, \end{gather}

(3.77c)\begin{gather}\bar{T}_{a} = \frac{1}{n_a} \int \,\textrm{d} \boldsymbol{v} \textrm{J}_0 h_a \frac{2}{3} \left( s_a^2 - \frac{3}{2} \right), \end{gather}

being $\textrm {J}_0 = \textrm {J}_0(b_a \sqrt {x_a})$ and $\textrm {J}_1 = \textrm {J}_1(b_a \sqrt {x_a})$.

The GK field component of the Sugama operator $\mathcal {C}_{ab}^T = \mathcal {C}_{ab}^T(\boldsymbol {k}, \mu , v_\parallel )$, expressed in terms of $h_b$, is also obtained by applying the transformation in (3.70) to (3.63). It yields,

(3.78)\begin{equation} \mathcal{C}_{ab}^{F} ={-} C_{ab}^{T}(F_{Ma} v_{{\parallel}}) \textrm{J}_0 \bar{V}_{{\parallel} ab} - C_{ab}^{T}( F_{Ma} v_{{\parallel}}) \textrm{J}_1 \frac{v_\perp}{v_\parallel}\bar{V}_{ab \perp} - C_{ab}^{T}( F_{Ma} s_a^2) \textrm{J}_0 \bar{W}_{ab}, \end{equation}

where

(3.79a)\begin{gather} \bar{V}_{{\parallel} ab} =\frac{m_a m_b}{\gamma_{ab} T_a} \int \,\textrm{d} \boldsymbol{v} \textrm{J}_0 \frac{h_b}{F_{Mb}} C_{ba}^{T} ( F_{Mb} v_\parallel), \end{gather}

(3.79b)\begin{gather}\bar{V}_{{\perp} ab} = \frac{ m_a m_b}{T_a \gamma_{ab}} \int \,\textrm{d} \boldsymbol{v} \textrm{J}_1 \frac{h_b}{F_{Mb}}\frac{v_\perp}{v_\parallel} C_{ba}^{T} \left( F_{bM} v_\parallel \right), \end{gather}

(3.79c)\begin{gather}\bar{W}_{ab} = \frac{T_b}{\eta_{ab}}\int \,\textrm{d} \boldsymbol{v} \textrm{J}_0 \frac{h_b}{F_{Mb}} C_{ba}^{T} \left( F_{Mb} s_b^2\right). \end{gather}

3.2.1. Gyrokinetic Hermite–Laguerre expansion

We now perform the Hermite–Laguerre expansion of the test and field components of the GK Sugama collision operator given in (3.71) and (3.78), respectively. Since the GK Sugama operator is formulated in terms of $h_a$, the Hermite–Laguerre expansion is written in terms of the non-adiabatic part of the gyromoments $n_a^{pj}$, defined in (2.16), which can be expressed as a function of the gyromoments $N_a^{pj}$, thanks to (2.17).

The Hermite–Laguerre expansion of the test part of the GK Sugama collision operator is obtained by projecting (3.71) onto the Hermite–Laguerre basis. The velocity integrals necessary for the projection are performed analytically in pitch-angle coordinates using the basis transformation from Hermite–Laguerre polynomials to Legendre and associated Laguerre polynomials given in (3.32a) and (3.32b). It is found that the gyromoment expansion of the test component of the GK Sugama operator is given by

(3.80)\begin{equation} \mathcal{C}_{ab}^{Tlk} = \mathcal{C}_{ab}^{Slk}+ R_{ab}^{ {\perp} lk} + \sum_{n=1}^3 \bar{X}^{nlk}_{ab}, \end{equation}

where

(3.81)\begin{equation} \mathcal{C}_{ab}^{Slk} = \sum_{j=0}^{\infty} \sum_{p=0}^{l+2k} \sum_{h=0}^{k + \lfloor l/2\rfloor} \sum_{g=0}^{p+2j} \sum_{r=0}^{j + \lfloor p/2\rfloor} \sum_{d=0}^h \frac{2^{p} (p!)^2}{(2p)!} \frac{L_{hd}^p \left(T^{{-}1} \right)^{ph}_{lk} }{\sigma_p^j \sqrt{2^l l!}} \frac{T_{pj}^{gr} \sqrt{2^g g!} }{(2p +1)} \bar{\nu}_{*ab}^{Spjd} n_{a}^{gr} \end{equation}

and

(3.82)\begin{equation} R_{ab }^{ {\perp} lk} ={-} \frac{b_a^2}{4} \sum_{p=0}^\infty \sum_{j=0}^\infty \sum_{g=0}^{l + 2 k} \sum_{h=0}^{k + \lfloor l/2\rfloor} \sum_{r=0}^{ p + 2j} \sum_{s=0}^{j + \lfloor p/2\rfloor} \frac{\left( T^{{-}1} \right)^{gh}_{lk}}{\sqrt{2^l l!}}\frac{\left( T^{{-}1} \right)^{rs}_{pj}}{\sqrt{2^p p!}} \mathcal{A}_{ab\perp}^{gr}n_a^{pj}. \end{equation}

The details of the calculations leading to the expressions of $\mathcal {C}_{ab}^{Slk}$ and $R_{ab}^{lk}$, given in (3.81) and (3.82), respectively, are reported in Appendix B. In (3.81), the Sugama velocity-integrated speed function is introduced, $\bar {\nu }_{*ab}^{Spjd} =\bar {\nu }_{*ab}^{\parallel pjd} - p (p+1) \bar {\nu }_{*ab}^{D pjd}$ (expressions for $\bar {\nu }_{*ab}^{D pjd}$ and $\bar {\nu }_{*ab}^{ \parallel pjd}$ can be found in (B10) and (B11), respectively). In (3.82), we introduce the coefficients

(3.83)\begin{equation} \mathcal{A}_{ab\perp}^{gr} = \sum_{s_1 =0}^s\sum_{h_1 =0}^h L^{r}_{ss_1} L^{g}_{hh_1} \mathcal{\alpha}_{ab\perp}^{grs_1 h_1} , \end{equation}

with

(3.84)$$\begin{align} \mathcal{\alpha}_{ab\perp}^{grs_1 h_1} &= \frac{\delta_g^r}{(2g+1)} \left( \bar{\nu}_{ab}^{{\parallel} g + 2 + s_1 + h_1} + \bar{\nu}_{ab}^{D g + 2 + s_1 + h_1 }\right) \nonumber\\ &\quad + \frac{d_{g}^r}{2} \left(\bar{\nu}_{ab}^{D(g+r)/2 + 2 + s_1 + h_1 } - \bar{\nu}_{ab}^{{\parallel} (g +r)/2 + 2 + s_1 + h_1} \right) , \end{align}$$

where

(3.85)\begin{equation} d_g^r = \begin{cases}\dfrac{2 (r + 1)(r+2)}{(2 r+1)(2r +3)(2r+5)}, & \text{for } g = r+2,\\ \dfrac{2(2r^2 + 2r -1)}{(2r-1)(2r+1)(2r+3)}, & \text{for } g=r, \\ \dfrac{2 r(r-1)}{(2r -3)(2r -1)(2r +1)}, & \text{for } g=r-2 \end{cases} \end{equation}

and the speed integrated energy and deflection functions, $\bar {\nu }_{ab}^{\parallel k} = 4 \int _0^\infty \,\textrm {d} s_a s_a^{2k} \nu _{ab}^{\parallel }(v) \textrm {e}^{- s_a^2} / \sqrt {{\rm \pi} }$ and $\bar {\nu }_{ab}^{Dk} =4 \int \,\textrm {d} s_a \textrm {e}^{- s^2} s_a^{2k} \nu _{ab}^D(v) / \sqrt {{\rm \pi} }$, with closed analytical expressions,

(3.86a)\begin{gather} \bar{\nu}_{ab}^{{\parallel} k} =\frac{4 \nu_{ab}}{\sqrt{\rm \pi}\chi^2} \left(E_{ab}^{k-3} - \chi e_{ab}^{k-2} \right), \end{gather}

(3.86b)\begin{gather}\bar{\nu}_{ab}^{Dk} = \frac{4 \nu_{ab} }{\sqrt{\rm \pi}} \left[ E_{ab}^{k-2} - \frac{1}{2 \chi^2} E_{ab}^{k-3} + \frac{1}{2 \chi} e_{ab}^{k-2}\right], \end{gather}

respectively. Finally, the $\bar {X}_{ab}^{nlk}$ terms in (3.80) are defined by

(3.87a)\begin{gather} \bar{X}^{1lk}_{ab} = \left( \theta_{ab} -1 \right) \left[ \mathcal{U}_{{\parallel}}^{lk}\mathcal{R}_{ab}^{{\parallel} 1} + \mathcal{U}_{{\perp}}^{lk} \mathcal{R}_{ab}^{{\parallel} 2} + \mathcal{T}_{{\perp}}^{lk} \mathcal{R}_{ab}^{{\parallel} 3} \right] , \end{gather}

(3.87b)\begin{gather}\bar{X}^{2lk}_{ab} = \left(\theta_{ab} -1 \right) \left( \mathcal{F}_{ab0\parallel}^{lk} \bar{u}_{{\parallel} a} +\mathcal{F}_{ab1\parallel}^{lk}\bar{u}_{{\perp} a} + \mathcal{F}_{abE}^{lk} \bar{T}_a\right) , \end{gather}

(3.87c)\begin{gather}\bar{X}^{3lk}_{ab} ={-} \frac{8 \nu_{ab} (\theta_{ab}-1)^2 \chi}{3 \sqrt{\rm \pi}(1 + \chi^2)^{1/2}}\left[ \mathcal{U}_{{\parallel}}^{lk} \bar{u}_{{\parallel} a} + \mathcal{U}_{{\perp}}^{lk} \bar{u}_{{\perp} a} + \frac{1}{1 + \chi^2} \mathcal{T}_{{\perp} }^{lk} \bar{T}_{a}\right] , \end{gather}

where we introduce

(3.88a)\begin{gather} \mathcal{R}_{ab}^{{\parallel} 1} ={-} \sum_{p=0}^\infty \sum_{j=0}^\infty \sum_{n=0}^{\infty} \sum_{s=0}^{n+j} \sum_{h=0}^{s + \lfloor p/2 \rfloor} \sum_{h_1=0}^h L_{hh_1}^{1} d_{njs}^0 \frac{\left( T^{{-}1} \right)^{1h}_{ps}}{\sqrt{2^p p!}}\frac{2 (1 + \chi^2) }{3 }\bar{\nu }_{ab}^{{\parallel} 3 + h_1} \mathcal{K}_{n}(b_a)n_a^{pj}, \end{gather}

(3.88b)\begin{gather}\mathcal{R}_{ab}^{{\parallel} 2} ={-} \sum_{p=0}^\infty \sum_{j=0}^\infty \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 d_{njs}^1 \frac{\left( T^{{-}1} \right)^{0h}_{ps}}{\sqrt{2^p p!}} \frac{ (1 + \chi^2)}{ (n+1)} \bar{\nu }_{ab}^{{\parallel} 2 + h_1} b_a \mathcal{K}_{n}(b_a) n_a^{pj},\end{gather}

(3.88c)\begin{gather}\mathcal{R}_{ab}^{{\parallel} 3} ={-} \frac{16 \nu_{ab}}{3 \chi^2 \sqrt{\rm \pi}} \sum_{p=0}^\infty \sum_{j=0}^\infty \sum_{n=0}^\infty \sum_{s=0}^{n+j} \sum_{h=0}^{s + \lfloor p/2 \rfloor} \sum_{h_1 =0}^h L_{hh_1}^0 d_{njs}^0 \frac{\left( T^{{-}1} \right)^{0h}_{ps}}{\sqrt{2^p p!}} \mathcal{K}_{n}(b_a) \nonumber\\ \times \left[ E_{ab}^{h_1} - \chi ( 1 + \chi^2) e_{ab}^{h_1+1}\right] n_a^{pj} \end{gather}

and

(3.89)\begin{gather} \mathcal{U}_{{\parallel}}^{lk} = \mathcal{K}_{k}(b_a) \frac{\delta_{l}^1 }{\sqrt{2}}, \end{gather}

(3.90)\begin{gather}\mathcal{U}_{{\perp}}^{lk} = \frac{b_a}{2} \left( \mathcal{K}_{k}(b_a) - \mathcal{K}_{k-1}(b_a)\right)\delta_l^0, \end{gather}

(3.91)\begin{gather}\mathcal{T}_{{\perp} }^{lk} =\left( \frac{\delta_l^2}{\sqrt{2}} - \delta_l^0 \right) \mathcal{K}_{k}+ \delta_{l}^0 \left[ (2k +1)\mathcal{K}_{k} - k \mathcal{K}_{k-1} - (k+1)\mathcal{K}_{k+1}\right]. \end{gather}

In (3.87b), we define

(3.92a)\begin{gather} \mathcal{F}_{ab 0 \parallel}^{lk} ={-} \frac{2 \left( 1 + \chi^2\right) }{3} \sum_{n=0}^\infty \sum_{s =0}^{n+k} \sum_{r =0}^{l + 2s} \sum_{t=0}^{s + \lfloor l/2\rfloor} \sum_{t_1=0}^t d_{nks}^0 L^1_{tt_1} \frac{\left( T^{{-}1}\right)^{1t}_{ls}}{\sqrt{2^l l!}} \delta_{r}^1 \mathcal{K}_{n}(b_a)\bar{ \nu}_{ab}^{{\parallel} 3 + t_1}, \end{gather}

(3.92b)\begin{gather}\mathcal{F}_{ab 1 \parallel}^{lk} ={-} 2 \left( 1 + \chi^2\right) \sum_{n=0}^\infty \sum_{s =0}^{n+k+1} \sum_{t=0}^{s + \lfloor l/2\rfloor} \sum_{t_1=0}^t d_{nks}^1 L_{tt_1}^{0} \frac{\left( T^{{-}1}\right)^{0t}_{ls}}{\sqrt{2^l l!}} \frac{b_a \mathcal{K}_{n}(b_a)}{(n+1)} \bar{\nu}_{ab}^{{\parallel} 2 + t_1}, \end{gather}

(3.92c)\begin{gather}\mathcal{F}_{abE}^{lk} ={-}\frac{8 \nu_{ab}}{\chi^2\sqrt{\rm \pi}} \sum_{n=0}^\infty \sum_{s =0}^{n+k} \sum_{t=0}^{s + \lfloor l/2\rfloor} \sum_{t_1=0}^t d_{nks}^0 \frac{\left( T^{{-}1}\right)^{0t}_{ls}}{\sqrt{2^l l!}} \mathcal{K}_{n}(b_a) L_{tt_1}^0 \left[ E_{ab}^{t_1} - \chi ( 1 + \chi^2) e_{ab}^{t_1+1}\right], \end{gather}

with

(3.93a)\begin{gather} \bar{u}_{{\parallel} a} = \sum_{j =0}^\infty \mathcal{K}_{j}(b_a) \frac{n_a^{1j}}{\sqrt{2}}, \end{gather}

(3.93b)\begin{gather}\bar{u}_{{\perp} a} = \frac{b_a}{2} \sum_{j=0}^\infty \mathcal{K}_{j}(b_a) \left( n_a^{0j} - n_a^{0j+1}\right), \end{gather}

(3.93c)\begin{gather}\bar{T}_{a} = \frac{2}{3} \sum_{j=0}^\infty \mathcal{K}_{j}(b_a) \left[ \frac{n_a^{2n}}{\sqrt{2}} + 2 jn_a^{0j} - j n_a^{0j-1} - (\,j+1)n_a^{0j+1} \right]. \end{gather}

In deriving (3.88) and (3.92), the closed expressions of $C_{ab}^{S}(f_{Ma} \boldsymbol {v} )$ and $C_{ab}^{S}( s_a^2 f_{Ma})$, given in (3.62), are used.

We now consider the Hermite–Laguerre expansion of the field component of the GK Sugama operator given in (3.78). Similarly to the test part of the collision operator, we multiply (3.78) by the Hermite–Laguerre basis and integrate the result over the velocity space in pitch-angle coordinates using the basis transformation from Hermite–Laguerre polynomials to Legendre and associated Laguerre polynomials, expressed in (3.32a) and (3.32b). This yields

(3.94)\begin{equation} \mathcal{C}_{ab}^{Flk} ={-} \bar{\mathcal{V}}_{{\parallel} ab}^{lk} \bar{V}_{{\parallel} ab} - \bar{\mathcal{V}}_{{\perp} ab}^{lk}\bar{V}_{{\perp} ab} - \bar{W}_{ab} \bar{\mathcal{W}}_{ab}^{lk}, \end{equation}

where we introduce

(3.95a)\begin{gather} \bar{\mathcal{V}}_{{\parallel} ab}^{lk} ={-} \frac{4 \theta_{ab}(1 + \chi^2)}{3 \sqrt{\rm \pi}} \sum_{n=0}^{\infty} \sum_{s =0}^{n+k}\frac{d_{nks}^0}{\sqrt{2^l l!}} \mathcal{K}_{n}(b_a) \mathcal{V}_{ab}^{ls}, \end{gather}

(3.95b)\begin{gather}\bar{\mathcal{V}}_{{\perp} ab}^{lk} = \sum_{n=0}^{\infty}\sum_{s=0}^{n + k +1}\frac{ d_{nks}^1}{\sqrt{2^l l!}} \frac{ b_a \mathcal{K}_{n}(b_a)}{ 2 (n+1)} \mathcal{I}_{ab}^{ls}, \end{gather}

(3.95c)\begin{gather}\bar{\mathcal{W}}_{ ab}^{lk} = \sum_{n=0}^\infty \sum_{s=0}^{n+k} \frac{d_{nks}^0}{\sqrt{2^l l!}} \mathcal{K}_{n}(b_a) \mathcal{W}_{ab}^{ls}, \end{gather}

with

(3.96)$$\begin{align} \mathcal{V}_{ab}^{lk} &= \frac{4 \nu_{ab}}{3 \sqrt{\rm \pi}} \sum_{g=0}^{l + 2k}\sum_{h =0}^{k + \lfloor l/2 \rfloor} \sum_{h_1 =0}^h L_{hh_1}^1 \left( T^{{-}1}\right)_{lk}^{1h} \left[ \frac{1}{ \chi^2 }\left(E_{ab}^{h_1}- \chi e_{ab}^{h_1 +1} \right) + \frac{2 \chi (\theta_{ab}-1) (3/2 + h_1)!}{3 \sqrt{\rm \pi}(1 + \chi^2)^{3/2}}\right] \nonumber\\ &\quad \times \left[ l > 0 \cup k > 0 \right],\end{align}$$

(3.97)$$\begin{align} \mathcal{I}_{ab}^{ls} & ={-} \frac{ 4 \nu_{ab}\theta_{ab} (1 + \chi^2)}{\chi^2\sqrt{\rm \pi}} \sum_{h =0}^{s +\lfloor l/2 \rfloor} \left( T^{{-}1}\right)^{0h}_{ls} \sum_{h_1 =0}^h L_{hh_1}^0 \left( E_{ab}^{h_1 -1} - \chi e_{ab}^{h_1}\right) \nonumber\\ &\quad - \frac{ \nu_{ab} \chi \theta_{ab}}{ 3 \sqrt{\rm \pi}} \frac{(\theta_{ab}-1)}{(1 + \chi^2)^{1/2}} \delta_{l}^0 \delta_s^0 , \end{align}$$

(3.98)$$\begin{align} \mathcal{W}_{ab}^{lk} & ={-} \frac{8\nu_{ab}\theta_{ab}}{\chi^2 \sqrt{\rm \pi}} \sum_{h=0}^{ k + \lfloor l/2\rfloor} \sum_{h_1=0}^h L_{hh_1}^0 \left( T^{{-}1}\right)_{lk}^{0h} \left(E_{ab}^{h_1} - \chi (1 + \chi^2) e_{ab}^{h_1 +1} \right) \nonumber\\ &\quad - \frac{8 \nu_{ab}\theta_{ab}}{3 \sqrt{\rm \pi}} \frac{ \chi (\theta_{ab}-1)}{(1 + \chi^2)^{3/2}} \left( 2 \delta_{l}^2 \delta_k^0 - \delta_l^0 \delta_k^1 \right). \end{align}$$

Finally, performing the velocity integrals in (3.79) using the gyromoment expansion for $h_b$ and the expressions of $C_{ba}^{T}( F_{Mb} v_\parallel )$ and $C_{ba}^{T}(F_{Mb} s_b^2)$, given in (3.64), yields

(3.99)\begin{gather} \bar{V}_{{\parallel} ab} ={-}\sigma \tau \theta_{ba}\frac{2 m_b}{ \gamma_{ab} } \left( 1 + 1/\chi\right) \sum_{p=0}^\infty \sum_{j=0}^\infty \sum_{n=0}^{\infty} \sum_{s =0}^{n+j} \frac{d_{njs}^0 }{\sqrt{2^p p!}} \mathcal{K}_{n}(b_b) n_b^{pj} \mathcal{V}_{ba}^{ps}, \end{gather}

(3.100)\begin{gather}\bar{V}_{{\perp} ab} = \frac{2m_b}{\chi\gamma_{ab}} \sum_{p=0}^\infty \sum_{j=0}^\infty \sum_{n=0}^{\infty} \sum_{s =0}^{n + j + 1}\frac{d^1_{njs}}{\sqrt{2^p p!}} \frac{b_b}{2(n+1)} \mathcal{K}_{n} (b_b)n_b^{pj}\mathcal{I}_{ba}^{ps}, \end{gather}

(3.101)\begin{gather}\bar{W}_{ab} = \frac{T_b}{\eta_{ab}} \sum_{p=0}^\infty \sum_{j=0}^\infty\sum_{n=0}^\infty \sum_{s =0}^{j + n} \frac{d_{njs}^0 }{\sqrt{2^p p!}} \mathcal{K}_{n}(b_b) n_b^{pj} \mathcal{W}_{ba}^{ps}. \end{gather}

The Hermite–Laguerre expansion of the test and the field components $\mathcal {C}_{ab}^{Tlk}$ and $\mathcal {C}_{ab}^{Flk}$, given in (3.80) and (3.94), respectively, can also be expressed in terms of $N_a^{pj}$ thanks to (2.17). Similarly to the linearized GK Coulomb operator, collisional FLR effects enter through the terms proportional to $\mathcal {K}_{n}(b_a)$.

3.2.2. Drift-kinetic limit of the Sugama collision operator

We now consider the DK limit of the GK Sugama collision operator. This limit can be deduced by neglecting FLR effects in (3.80) and (3.94), which appear through the kernel functions $\mathcal {K}_{n}(b_a)$. In addition, $f_a = {g}_a$ follows from (3.27) in the DK limit, such that $n_a^{pj} = N_a^{pj}$. Hence, we obtain the DK limit of the Hermite–Laguerre expansion of the test component,

(3.102)\begin{equation} \mathcal{C}_{ab}^{Tlk} = \mathcal{C}_{ab}^{Slk} + \sum_{n=1}^{3} X_{ab}^{nlk}, \end{equation}

where $\mathcal {C}_{ab}^{Slk}$ is given by (3.81) with $n_a^{pj} = N_a^{pj}$, and

(3.103a)\begin{gather} X_{ab}^{1lk} = (\theta_{ab}-1)\left[-(1 + \chi^2)\sqrt{2} \delta_l^1 \delta_k^0 \sum_{p,j} \frac{\textrm{I}_{ab1}^{ {\parallel} pj}}{\sqrt{2^p p!}}N_a^{pj}+ \frac{2}{3} \left( \frac{\delta_l^2 \delta_k^0}{\sqrt{2}}- \delta_{k}^1 \delta_l^0 \right) \sum_{p,j} \frac{\textrm{I}_{ab2}^{pj}}{\sqrt{2^p p!}} N_a^{pj} \right], \end{gather}

(3.103b)\begin{gather}X_{ab}^{2lk} = ( \theta_{ab}-1) \left[- \frac{\left( 1 + \chi^2\right)}{\sqrt{2^l l!}} \frac{ 2 u_{{\parallel} a}}{v_{Ta}} \textrm{I}_{ab1}^{ {\parallel} lk } + \frac{\delta T_a}{T_a} \frac{\textrm{I}_{ab2}^{lk}}{\sqrt{2^l l!}} \right], \end{gather}

(3.103c)\begin{gather}X_{ab}^{3lk} ={-} ( \theta_{ab}-1)^2 \frac{4 \nu_{ab}}{3 \sqrt{\rm \pi}} \left[ \frac{\sqrt{2}\chi}{(1 + \chi^2)^{1/2}} \frac{u_{{\parallel} a}}{v_{Ta}} \delta_l^1 \delta_k^0 + \frac{\delta T}{T_a} \frac{2 \chi}{ (1 + \chi^2)^{3/2}}\left( \frac{\delta_l^2 \delta_k^0}{\sqrt{2}} - \delta_l^0 \delta_k^1\right) \right], \end{gather}

where $u_{\parallel a} / v_{Ta}=N_a^{10}/ \sqrt {2}$ and $\delta T_a / T_a =(\sqrt {2} N_a^{20} - 2N_a^{01})/3$. In (3.103), we define

(3.104)\begin{gather} \textrm{I}_{ab1}^{ {\parallel} lk } = \frac{1}{3} \sum_{h=0}^{ k + \lfloor l/2 \rfloor } \sum_{h_1=0}^h L_{hh_1}^1 \left( T^{{-}1}\right)_{lk}^{1h} \bar{\nu}_{ab}^{{\parallel} h_1 + 3} \left[ l>0 \cup k>0\right], \end{gather}

(3.105)\begin{gather}\textrm{I}_{ab2}^{lk} ={-} \frac{8 \nu_{ab}}{\chi^2 \sqrt{\rm \pi} } \sum_{h=0}^{k + \lfloor l/2\rfloor}\sum_{h_1=0}^h L^{0}_{hh_1} \left( T^{{-}1}\right)_{lk}^{0h} \left[E_{ab}^{h_1}- \chi \left( 1 + \chi^2 \right)e_{ab}^{h_1+1} \right]. \end{gather}

Similar considerations lead to the DK limit of the Hermite–Laguerre gyromoment expansion of the field component, given in (3.94), that is

(3.106)\begin{gather}\mathcal{C}_{ab}^{Flk} ={-} \frac{2\theta_{ab} \theta_{ba}}{\chi \sigma} \frac{m_a}{ \gamma_{ab}}\frac{(1 + \chi^2)}{\sqrt{2^l l!}} \sum_{p=0}^\infty \sum_{j=0}^\infty\frac{(1 + 1/\chi^{2})}{\sqrt{2^p p!}} \mathcal{V}_{ab}^{lk} \mathcal{V}_{ba}^{pj} N_b^{pj} \nonumber\\ - \frac{1}{\sqrt{2^l l!}} \frac{1}{\tau} \frac{T_a}{\eta_{ab}} \sum_{p=0}^\infty \sum_{j=0}^\infty \frac{N_b^{pj}}{\sqrt{2^p p!}} \mathcal{W}_{ab}^{lk} \mathcal{W}_{ba}^{pj}. \end{gather}

4. Numerical implementations

Within the adopted Hermite–Laguerre gyromoment method, the numerical evaluation of the collision operators leads to the calculation of large sums of algebraic coefficients such as the ones appearing, for example, in (3.47). In the numerical implementation of the GK pseudospectral method proposed by (Mandell et al. Reference Mandell, Dorland and Landreman2018) where a GK Dougherty collision operator is considered, a closed analytical formula is used for the lowest-order coefficients from the products between Laguerre polynomials (such as $d_{hs_1}^{0}$ and $\bar {d}_{0kf}^{0}$ defined in (3.37) and (3.44), respectively). On the other hand, a Fourier dealiased scheme is employed to numerically compute the integral representation of the higher-order coefficients. This method can be generalized to obtain $\bar {d}_{nkf}^{m}$ and $d_{n'hs_1}^{m}$, but it is unpractical for the basis transformation coefficients $T^{lk}_{pjm}$ that appear in the operators considered here. This is because of the rapid increase in magnitude of these coefficients and the enhancement of the oscillatory nature of the basis functions (such as the Hermite polynomials) with the order of the polynomial basis. In fact, due to the presence of positive and negative numerical coefficients of different magnitude that eventually cancel each other, the Hermite–Laguerre formulation of the present collision models is highly sensitive to round-off errors. This becomes particularly pronounced as the order of the polynomials basis increases. As an example, figure 1 displays the magnitude (with sign) of a few of the basis transformation coefficients, in particular $T_{pqq}^{qq}$ defined in (3.32a). An increase of $40$ orders is observed by going from $p=1$ to polynomials of order $p =50$. We note that similar issues are observed in the numerical implementation of spectral methods based on generalized Laguerre polynomials in the energy variable $v^2$ (Belli & Candy Reference Belli and Candy2011). In addition, it should be mentioned that convergence tends to be relatively slow when using Laguerre polynomials in $v^2$ because they have very similar shapes for the bulk of the distribution (Landreman & Ernst Reference Landreman and Ernst2013). As a matter of fact, we find that the use of double and/or quadruple precision is not sufficient to carry out the evaluation of the coefficients required for the Coulomb and Sugama operators.

Figure 1.

Magnitude (with sign) of the coefficients of the basis transformation $T_{pqq}^{qq}$ in (3.32a) as a function of $p$ when $q = 1$, $2$ and $4$.

To avoid spurious unphysical effects associated with precision loss in the linearized collision operators, we evaluate the closed analytical expressions of the numerical coefficients and perform the sums numerically using a multiple precision arithmetic package. To avoid redundancy, the species-independent numerical coefficients that depend on the mass and temperature ratios of the colliding species ($\bar {\nu }_{*ab}^{Tpjt}$, $\bar {\nu }_{*ab}^{Fpjt}$ and $\bar {\nu }_{*ab}^{Spjd}$ defined in Appendices A and B) are computed once and stored. Then, given the normalized perpendicular wavenumber $b_a$, the mass and temperature ratios, the GK linearized collision operators can be obtained by performing linear combinations of the previously computed numerical coefficients.

Because of FLR effects, infinite sums arise in the GK Coulomb and Sugama operators from the expansion of the Bessel functions, $\textrm {J}_m$, in terms of associated Laguerre (see (2.12)). These sums appear in the Hermite–Laguerre projections, such as the ones in (3.47) and (3.88). For their numerical implementations, we truncate these sums at a finite $N > 0$. While a large number of terms in the FLR terms might a priori be needed for a good convergence, the nature of the kernel function, $\mathcal {K}_{n}$ defined in (2.13), ensures that only a small number of terms contribute in the sums at a given perpendicular wavenumber. This is because the kernel functions $\mathcal {K}_{n}$ decreases rapidly with $n$, as $\mathcal {K}_{n} \sim 1/ n!$ for large $n$ for all $b_a$, and $\mathcal {K}_{n} \sim b_a^{2n}/ n!$ when $b_a \ll 1$ (Frei et al. Reference Frei, Jorge and Ricci2020). Additionally, since the maximum of $\mathcal {K}_{n}$ occurs at $b_a = 2 \sqrt {n}$, as a rule of thumb, one can choose $N$ such that $N > b_a^2/4$. Figure 2 displays the relative error of the diagonals elements of the GK Coulomb operator corresponding to the gyromoments (blue lines) $(p,j) =(0,0)$ and (orange lines) $(p,j) = (1,1)$, for typical values of the normalized perpendicular wavenumber in the ion gyroscale when (solid lines) $b_a = 0.5$ and (dotted lines) $b_a = 1$, as a function of $N$. The relative error is computed with respect to the $N =4$ case. As observed, the error decreases linearly with $N$, and it is of the order of $\lesssim 0.1\,\%$ or smaller at $N =3$, when $b_a = 1$. A similar behaviour is found with the GK Sugama operator.

Figure 2.

Relative error of the diagonal elements, (blue lines) $(p,j) =(0,0)$ and (orange lines) $(p,j) = (1,1)$, of the GK Coulomb, when (solid) $b_a = 0.5$ and (dotted) $b_a = 1$, as a function of $N$. The relative error is measured with respect to the $N=4$.

Finally, we discuss the coupling between gyromoments introduced by the GK Coulomb collision operator. By mapping the gyromoments into a one-dimensional array, the GK Coulomb operator can be represented in matrix form, relating the projection of the collision operator with the gyromoments. Given $(P_a,J_a)< \infty$, we can construct the collisional matrix by defining the one-dimensional gyromoment array as $[N_{a}^{00}, N_{a}^{01},N_{a}^{02}, \ldots , N_{a}^{0J} , N_a^{10}, N_a^{11}, N_a^{12}, \dots , N_a^{P_a0}, \dots , N_a^{P_aJ_a}]^T$ and the one-dimensional index $\bar {l}(p,j) = (J_a +1) p + 1 + j$. Figure 3 displays the block matrix obtained from the GK Coulomb collision operator for like-species when $b_a = 0.5$ for $(p,j) \le (P_a,J_a)$. It can be seen that the $p$th Hermite blocks is coupled with the $l = p + 2n$ block (with $n$ a positive integer). We remark that the block structure of the matrix representing the Coulomb collision operator results from the vanishing basis transformation coefficients $T^{lk}_{pjm}$. Notice also the negative elements in the diagonal with a magnitude that increases with $\bar {l}$ at fixed $p$, yielding damping of higher-order gyromoments. A similar block matrix is obtained for the GK Sugama.

Figure 3.

Block matrix obtained from the GK Coulomb collision operator for like species when $b_a = 0.5$. We show, on the vertical axis, the one-dimensional index $\bar {l}(p,j) = (J_a +1) p + 1 + j$ and, on the horizontal axis, the one-dimensional index $\bar {k}(l,k) = (J_a +1) l + 1 + k$.

5. Numerical tests

This section presents the numerical tests and benchmarks of the implementation of the GK Coulomb and Sugama collision operators. First, we numerically investigate the linear growth rate of the ITG instability in § 5.1. For this test, we consider both the linearized GK Coulomb and Sugama collision operators and various perpendicular wavenumbers and levels of collisionality. As a second test, in § 5.2, we address the collisional damping of ZF. We benchmark our results against the DK linearized Coulomb and GK Sugama collision operators (Crandall et al. Reference Crandall, Jarema, Doerk, Pan, Merlo, Görler, Bañón Navarro, Told, Maurer and Jenko2020) currently available in the public release of the GK continuum code GENE (Jenko et al. Reference Jenko, Dorland, Kotschenreuther and Rogers2000), finding good agreement.

5.1. Ion-temperature gradient instability

Belonging to the class of ion-gyroscale gradient-driven drift-wave instability, the ITG mode is widely recognised as the main candidate to explain the level of anomalous ion heat transport observed experimentally in the tokamak core (Garbet et al. Reference Garbet, Mantica, Angioni, Asp, Baranov, Bourdelle, Budny, Crisanti, Cordey and Garzotti2004). As a matter of fact, the collisionless theory of the ITG mode in core conditions is well established and documented (see, e.g. Romanelli Reference Romanelli1989; Hahm & Tang Reference Hahm and Tang1989; Romanelli & Briguglio Reference Romanelli and Briguglio1990; Chen, Briguglio & Romanelli Reference Chen, Briguglio and Romanelli1991). On the other hand, the role of the ITG mode is less clear in the edge and in the SOL region due to the limitation of Braginskii-like fluid models that are most often applied in these regions (Hallatschek & Zeiler Reference Hallatschek and Zeiler2000; Mosetto et al. Reference Mosetto, Halpern, Jolliet, Loizu and Ricci2015). While a complete collisional theory of ITG modes using the gyromoment approach will be subject of a future publication, we consider here the slab growth rate of this mode, as a function of the collisionality and the perpendicular wavenumber, when the GK Coulomb and Sugama collision operators are used. We also perform a comparison with GENE.

We focus on a shearless slab configuration and measure the relative strength of the ion temperature and density gradients by the dimensionless parameter $\eta = L_N / L_{T_i}$, where $L_N = | \boldsymbol {\nabla }_\perp \ln N |^{-1}$ and $L_{T_i} = | \boldsymbol {\nabla }_\perp \ln T_i |^{-1}$. We assume that the electrons are adiabatic, and consider only like-species collisions, $C_i = C_{ii}$, neglecting the ion–electron collisions since they occur on a time scale larger by a factor $O(\sqrt {m_i / m_e})$ than the ion–ion collisions. We neglect magnetic gradient drifts and particle trapping, considering therefore the slab ITG branch. We define the normalized ion–ion collisionality by the parameter $\nu = L_N \nu _{ii} / c_s$, with $c_s^2 = T_e/ m_i$ the ion sound speed, and we impose $T_e = T_i$. We normalize the perpendicular wavenumber, $k_\perp$, to $\rho _s = c_s /\varOmega _i$ and the linear growth rate $\gamma$ to $c_s / L_N$. The parallel wavenumber is fixed at $k_\parallel = 0.1 / L_N$. We solve the gyromoment hierarchy, in (2.14), with $(P, J ) = (20,8)$ for the ion species by applying a closure by truncation such that $N_{i}^{pj} =0$ for $(p,j) > (P,J)$. We have carefully checked that the considered number of ion gyromoments provide a sufficient resolution of the velocity space for the numerical tests carried out in this section. We note the slab ITG mode relies mainly on parallel dynamics and thus emphasizes the faster convergence of the Hermite expansion in $v_\parallel$ relative to the Laguerre expansion.

The numerical results are shown in figure 4. We consider the full linearized GK Coulomb collision operator given in (3.46), and the GK Sugama collision operator, with test and field components defined in (3.80) and (3.94), respectively. The temperature gradient strength is fixed at $\eta =3$, and a scan in the perpendicular wavenumber $k_\perp$ is performed at different collisionalities. The GK Coulomb is represented by the solid blue lines, and the GK Sugama by the red solid lines, while their DK limits are shown for comparison by the dashed lines. Markers represent GENE results with the DK Coulomb and GK Sugama. First of all, we note that a good agreement is found between our results and GENE. We have also successfully compared the DK Sugama operator model with GENE (comparison not shown). This shows the ability of the gyromoment approach to accurately describe FLR collisional damping in advanced GK collision operators. From a more physical perspective, we observe that the GK Sugama yields an underestimate of the ITG growth rate (of the order of $15\ \%$ or less) compared with GK Coulomb for wavelengths smaller than the one corresponding to the peak growth rate. We also note, by comparing the DK operators, that FLR collisional effects provide a strong damping and that the DK Sugama follows closely the DK Coulomb, similar to results obtained for TEMs (Pan et al. Reference Pan, Ernst and Crandall2020, Reference Pan, Ernst and Hatch2021). In principle, comparison with the GK Coulomb operator in GENE (Pan et al. Reference Pan, Ernst and Crandall2020) is also possible when it becomes publicly available.

Figure 4.

Comparisons between the GK Coulomb (solid blue lines) and GK Sugama (solid red lines) operator models, given in (3.54), (3.80) and (3.94), and GENE (markers), as a function of the perpendicular wavenumber, $k_\perp$ for different collisionalities. The collisionless solution (dotted lines) is shown for comparison. A good agreement is found with the GENE operators. Here, $\eta = 3$.

To conclude our analysis of the ITG mode, we illustrate the effect of the collisional polarization terms in the GK Coulomb collision operator, expressed by $\sum _{i=1}^{2} (\mathcal {C}^{Tlk}_{ab\phi i} + \mathcal {C}^{Flk}_{ab\phi i} )$ appearing in (3.49). We remark that $\mathcal {C}^{Tlk}_{aa\phi 1} +\mathcal {C}^{Flk}_{aa\phi 1} =0$ for like species, such that only the $\mathcal {C}^{Tlk}_{aa\phi 2}$ and $\mathcal {C}^{Flk}_{aa\phi 2}$ terms in the polarization contribute in the GK linearized Coulomb operator. We perform a scan of the growth rate at a fixed collisionality ($\nu = 0.1$), as a function of the perpendicular wavenumber $k_\perp$ for different temperature gradient strengths $\eta$. The results are shown in figure 5 where the growth rates with the GK Coulomb that include and exclude the polarization terms are shown by the solid black and red lines, respectively. It is observed that the collisional polarization terms have a negligible effect on the ITG peak growth rate. On the other hand, the difference increases up to $15\,\%$ for short wavelengths.

Figure 5.

Effects of polarization given in (3.49) in the GK Coulomb collision operator on the ITG growth rate for different temperature gradients as the function of the perpendicular wavenumber $k_\perp$ when $\nu = 0.1$.

5.2. Collisional damping of ZF

Zonal flows can reduce the particle and heat exhaust from the core in magnetised plasma confinement devices (Hasegawa, Maclennan & Kodama Reference Hasegawa, Maclennan and Kodama1979; Hammett et al. Reference Hammett, Beer, Dorland, Cowley and Smith1993). These self-generated, primarily poloidal and axisymmetric, flows help to shear apart structures that propagate radially, caused by the growth of an instability driven by, e.g. a steep pressure gradient. In fact, ZF regulate and, ultimately, can suppress the turbulent transports driven by, for example, TEM and ITG modes (Candy & Waltz Reference Candy and Waltz2003; Ernst et al. Reference Ernst, Bonoli, Catto, Dorland, Fiore, Granetz, Greenwald, Hubbard, Porkolab and Redi2004).

The collisionless damping of self-generated ZF driven by ITG turbulence to a non-vanishing residual in the large radial wavelength limit has been addressed and calculated by Rosenbluth & Hinton (Reference Rosenbluth and Hinton1998), in a seminal work that considers a large aspect ratio and circular flux surface geometry. Later, including the effects of ion collisions modelled by a pitch-angle scattering operator, Hinton & Rosenbluth (Reference Hinton and Rosenbluth1999) showed that the ZF component of the electrostatic potential, $\phi _z(t)$, decays to a smaller value than the one predicted by the collisionless theory on long time scales. Hinton–Rosenbluth's analytical prediction,

(5.1)\begin{equation} \frac{\phi_z(t)}{\phi_z(0)} = \varpi \exp\left( \varpi^2\alpha^2 t\right) \left[ 1 - \text{erf} \left( \varpi \alpha t^{1/2}\right)\right], \end{equation}

has been the subject of numerous tests and benchmarks of GK codes (see, e.g. Idomura et al. Reference Idomura, Ida, Kano, Aiba and Tokuda2008; Merlo et al. Reference Merlo, Sauter, Brunner, Burckel, Camenen, Casson, Dorland, Fable, Görler and Jenko2016). In (5.1), we introduce $\varpi = 1/ (1 + 1.61 q^2 / \sqrt {\epsilon })$ and $\alpha = 3 {\rm \pi}q^2 \bar {\nu }^{1/2}/ (\epsilon \varLambda ^{3/2})$, where $\varLambda \simeq \ln [16 (\epsilon / \bar {\nu } t)^{1/2}]$ and $\bar {\nu } = 0.61 \nu _{ii}$. In addition, $\epsilon$ is the local inverse aspect ratio and $q$ is the safety factor. As an extension of the work by Hinton & Rosenbluth (Reference Hinton and Rosenbluth1999) to arbitrary collisionality, Xiao et al. (Reference Xiao, Catto and Molvig2007) retains the velocity dependence of the deflection frequency, $\nu _{ii}^D(v)$, appearing in the pitch-angle scattering operator, as well as a momentum restoring field component, in the limit of ZF with perpendicular wavelength larger than the ion poloidal Larmor radius. In fact, Xiao et al. (Reference Xiao, Catto and Molvig2007) shows that the time dependence of the collisional decay of the ZF residual is predicted by

(5.2)\begin{equation} \frac{\phi_z(t)}{\phi_z(0)} = \frac{\beta }{1 +\beta} \left[ 1 + \frac{1 - \varTheta}{\varTheta + \beta} \exp({- (1 + \beta) \varGamma \nu_{ii} t/ [\varTheta + \beta]})\right], \end{equation}

with $\beta = \epsilon ^2 /q^2$, $\varTheta = 1.635 \epsilon ^{3/2} + \epsilon ^2 /2 + 0.360 \epsilon ^{5/2}$ and $\varGamma = 0.4(1.46 \sqrt {\epsilon } + 1.32 \epsilon )$. Since the pitch-angle scattering is the dominant collisional damping process in a large aspect ratio tokamak (Hinton & Rosenbluth Reference Hinton and Rosenbluth1999), energy diffusion – neglected in (5.2) – is expected to play a subdominant role. Equation (5.2) provides an accurate analytical estimate to benchmark the ZF collisional damping using the DK Coulomb and DK Sugama collision operators.

As a test of the implementation of the collisional gyromoment hierarchy, given in (2.14), we consider a shearless, toroidal circular concentric flux surface geometry, with an initial zonal density perturbation radially varying with a wavevector $\boldsymbol {k} = k_r \boldsymbol {\nabla } r$. Density and temperature equilibrium gradients are neglected. A finite number of ion gyromoments, $(P,J)= (30,15)$, are evolved and electrons are assumed adiabatic. We consider a level of ion–ion neoclassical collisionality that is in the Pfirsch–Schlüter regime, i.e. $\nu _{i*} \simeq 3.1$. In this parameter regime, collisions smear out fine structures in $v_\parallel$, which develop at low collisionality because of passing ions and finite orbit width effects, and might require a larger number of gyromoments (Idomura et al. Reference Idomura, Ida, Kano, Aiba and Tokuda2008).

Figure 6 shows the ZF collisional damping predicted by the gyromoment approach (blue solid line) using the DK Coulomb collision operator and its comparison with GENE results (red line with cross markers). We also plot Hinton and Rosenbluth's, (5.1), and Xiao's, (5.2), predictions. An excellent match is observed between the gyromoment approach and GENE simulation. In addition, the collisional ZF damping agrees with the analytical prediction, given in (5.2). Ultimately, the Hermite–Laguerre spectrum at $t = 10 / \nu _{ii}$, in figure 7, reveals that the velocity-space dependence of the ion distribution function is well resolved. We remark that the decay of the spectrum is slower along the Hermite direction compared with the Laguerre direction. This is primarily because of the ballistic mode response of the passing ion population (Idomura et al. Reference Idomura, Ida, Kano, Aiba and Tokuda2008).

Figure 6.

Collisional ZF damping using (blue solid line) the DK Coulomb collision operator given in (3.54) and (red solid with markers) the GENE results. The Hinton–Rosenbluth collisional damping, (5.1), is shown by the black line, and the collisional damping predicted by Xiao, Catto & Molvig (Reference Xiao, Catto and Molvig2007), given in (5.2), is plotted by the solid green. A good agreement is found between GENE results and the gyromoment approach. Here, the parameters are $k_r = 0.05$, $\epsilon = 0.1$, $q = 1.4$ and $\nu _{i*} = 3.1$.

Figure 7.

Modulus of the normalized Hermite–Laguerre spectrum of the ZF residual at $\nu t =10$ taken from figure 6 on a logarithmic scale.

We now focus on the collisional ZF damping when the linearized GK Coulomb and GK Sugama collision operators are considered. The results are shown in figure 8, for (a) $k_r = 0.1$ and (b) $k_r =0.5$. The ZF damping by the GK Sugama operator is compared with GENE. Other parameters are the same as in figure 6. The GK Coulomb collision operator yields a smaller collisional damping of the ZF residual than the one predicted by the GK Sugama. The same qualitative result was obtained in Pan et al. (Reference Pan, Ernst and Crandall2020, Reference Pan, Ernst and Hatch2021), but in the low collisionality (banana) regime. We remark that the GK Sugama agrees well with GENE simulations. The observed difference between GK Coulomb and Sugama operators in the ZF residual collisional damping suggests that higher ion particle and heat fluxes can be expected in nonlinear simulations that use the GK Sugama operator rather than the GK Coulomb, particularly in the $\nu _{i*} \gtrsim 1$ Pfirsch–Schlüter regime. However, the opposite trend in TEM turbulence was observed in Pan et al. (Reference Pan, Ernst and Hatch2021).

Figure 8.

Collisional ZF damping using (blue solid line) the GK Coulomb collision operator, given in (3.46), (red solid) the GK Sugama, given in (3.80), and (red solid with markers) GENE when (a) $k_r = 0.1$ and (b) $k_r = 0.5$. Others parameters are the same in as figure 6.

Furthermore, we note that larger oscillations in the gyromoment approach are observed compared with GENE when $k_r$ increases. This is related to the large number of gyromoments required to accurately resolve the resonant passing particle dynamics and finite orbit width effects during the geodesic mode oscillation. The low resolution does not influence the collisional damping of the ZF residual.

6. Convergence test

The numerical tests in § 5 illustrate the ability of the gyromoment approach to model collisional effects using advanced linearized GK collision operators. Additionally, the direct comparisons with the GK continuum GENE code show that the Hermite–Laguerre approach is able to resolve the velocity dependence of the distribution functions in the presence of collisions, as shown in figure 7, with a finite number of gyromoments even when a simple closure by truncation is adopted.

As we now show, the number of gyromoments necessary to describe the distribution function decreases with the collisionality, as collisions provide a natural way to smear out any sharp velocity-space structures and damp the higher-order gyromoments in the Hermite–Laguerre spectrum. Thus, in the high collisional regime, only the lowest-order gyromoments are necessary, while, as the collisionality decreases, a larger number of gyromoments is necessary to resolve the essential kinetic effects in order to retrieve the collisionless solution (see, e.g. dotted black line in figure 4). We remark that the presence of kinetic effects related to, for example, finite magnetic gradient drifts, passing ions, particle trapping and finite orbit width, can highly affect the velocity-space structures, and ultimately deteriorates the convergence of the Hermite–Laguerre spectrum, when collisions are subdominant. These fine scale structures require a fine velocity-space resolution.

Contrary to velocity-space grid methods, the Hermite–Laguerre gyromoment approach features a spectral convergence rate. This convergence rate is faster than any algebraic schemes that characterize the grid methods (Boyd Reference Boyd2001), though some basis functions converge faster than others, depending on the velocity-space structure present in the distribution function (Landreman & Ernst Reference Landreman and Ernst2013). To illustrate the convergence properties of the gyromoment method, we compute the absolute error on the ITG growth rate as a function of the ion Hermite–Laguerre resolution $(P,J)$. Similarly, by using GENE, we consider convergence study with the number of velocity-space grid points, defined by $N_{v_\parallel }$ and $N_{\mu }$ (with $N_{v_\parallel }$ the number of grid points in the parallel velocity-space direction $v_\parallel$, and $N_\mu$ the number of grid points in the perpendicular velocity-space direction $\mu$) for different collisionality using the DK Coulomb collision operator. We first note that the difference between the most resolved growth rate evaluated by GENE and by the gyromoment method is less than $3\,\%$ in all cases. The results are shown in figure 9 ($k_\perp = 0.5$, $k_\parallel = 0.1$ and $\eta =3$). The convergence with respect to the parallel and perpendicular velocity direction refinement are plotted on the left-hand and on the right-hand panel, respectively.

Figure 9.

Absolute error of the ITG linear growth rate as a function of the number of gyromoments, $(P,J)$ (red), and the number of velocity-space grid points, $(N_{v_\parallel }, N_{\mu })$ used in GENE (black), for increasing collisionality regimes. We display, on the same horizontal axis, the gyromoment and GENE errors as a function of the discretization along the parallel direction given by $P$ for the gyromoment method and $N_{v_\parallel }$ for GENE (a), and as a function of the discretization along the perpendicular direction given by $J$ for the gyromoment method and $N_{\mu }$ for GENE (b), respectively. Collisions are modelled by the DK Coulomb operator. Here, $k_\perp = 0.5$, $k_\parallel = 0.1$ and $\eta = 3$.

As observed in figure 9, the gyromoment approach features a convergence rate in both velocity-space directions that is sensitive to the collisionality $\nu$ (as discussed above), while the velocity-space grid method, used in GENE, shows a convergence rate that is much less sensitive to $\nu$. In fact, while in both methods, fewer number of gyromoments and grid points in velocity-space are required in the high collisional regime as the perturbed distribution function does not feature sharp velocity-space structures, the effect of collision is much stronger on the number of gyromoments required than on grid points. While the increase of the convergence rate of the gyromoment method with $\nu$ is clearly visible in figure 9, we notice that the GENE error also decreases with the collisionality $\nu$, but is much less visible. We remark that, in the case of slab ITG mode, the perturbed distribution function does not feature sharp velocity-space structures which appear, e.g. in the ballistic mode response from passing particles and in the trapping of particles in velocity-space (Idomura et al. Reference Idomura, Ida, Kano, Aiba and Tokuda2008). However, preliminary studies show that the same conclusion on the convergence rate of the gyromoment approach can be made even if stronger kinetic drives are present. Thus, the convergence study we carried out suggests that the gyromoment approach offers a particularly efficient numerical scheme to simulate the plasma boundary, where the collisionality is larger (with respect to the strength of velocity-space structure drive) than in the core, compared with a finite difference. A careful analysis of the description of fine velocity-space structure using the gyromoment approach will be subject of a future publication.