The numerical implementation of the linearized gyrokinetic and drift-kinetic improved Sugama (IS) collision operators, recently introduced by Sugama *et al.* [Phys. Plasmas **26**, 102108 (2019)], is reported. The IS collision operator extends the validity of the widely used original Sugama (OS) operator [Sugama *et al.*, Phys. Plasmas **16**, 112503 (2009)] to the Pfirsch–Schlüter collisionality regime. Using a Hermite–Laguerre velocity–space decomposition of the perturbed gyrocenter distribution function that we refer to as the gyro-moment approach, the IS collision operator is written in a form of algebraic coefficients that depend on the mass and temperature ratios of the colliding species and perpendicular wavenumber. A comparison between the IS, OS, and Coulomb collision operators is performed, showing that the IS collision operator is able to approximate the Coulomb collision operator in the case of trapped electron mode in H-mode pedestal conditions better than the OS operator. In addition, the IS operator leads to a level of zonal flow residual which has an intermediate value between the Coulomb and the OS collision operators. The IS operator is also shown to predict a parallel electrical conductivity that approaches the one of the Coulomb operator within less than 1%, while the OS operator can underestimate the parallel electron current by at least 10%. Finally, closed analytical formulas of the lowest order gyro-moments of the IS, OS, and Coulomb operators are given, which are ready to use to describe the collisional effects in reduced gyro-moment fluid models.

## I. INTRODUCTION

While the core of fusion devices, such as tokamaks and stellarators, is sufficiently hot that plasma collisional processes are less important, the lower temperature in the boundary region enhances the role of collisions, calling for their accurate description. In fact, while the turbulent particle and heat transport in fusion devices is primarily anomalous, which is driven by small scale electromagnetic instabilities, collisions between particles can still play a crucial role since they may significantly affect the linear properties of these instabilities,^{1–5} their saturation mechanisms via zonal flow (ZF) generation,^{4–6} and the velocity–space structures of the particle distribution functions. In addition, collisions drive the neoclassical processes when turbulent transport is suppressed.^{7} An accurate description of collisions might also be important in the transport of heavy impurities, even in the core region, because of their large atomic number.

The relatively large value of the Coulomb logarithm ( $ ln \u2009 \Lambda \u2273 10$) in fusion devices, where small-angle deflections dominate, allows for the use of the Fokker–Planck collision operator^{8} to describe binary Coulomb collisions between particles. In this work, we refer to the Fokker–Planck collision operator as the Coulomb collision operator. The integrodifferential nature of the Coulomb collision operator makes its analytical and numerical treatment challenging. Hence, for practical applications, it is usually assumed that the distribution function is close to a Maxwellian allowing for the linearization of the Coulomb collision operator up to first order in the perturbed quantities. However, the numerical implementation and analytical treatment of the linearized Coulomb collision operator remains challenging. For this reason, approximated linearized collision operator models have been proposed in the literature.^{9–12}

One of the most widely used collision operator models, which we refer to as the linearized original Sugama (OS) collision operator in this work, is originally derived in Ref. 13. The OS collision operator is constructed from the linearized Coulomb collision operator to fulfill the conservation laws (particle, momentum, and energy), the entropy production criterion and the self-adjoint relations for arbitrary mass and temperatures ratios of the colliding species. This operator extends previous models^{11} to the case of colliding species with different temperatures. While the OS operator is implemented in numerous gyrokinetic (GK) codes and tested in neoclassical and turbulent studies,^{14,15} its deviation with respect to the Coulomb operator is expected to be enhanced when applied in the boundary plasma conditions. For instance, the OS produces a stronger collisional zonal flow damping^{5,6} and can therefore yield turbulent transport levels that can significantly differ from the ones obtained by the Coulomb collision operator.^{5} A recent study of the collisional effects of ion-temperature gradient (ITG) reports that the OS Sugama predicts a smaller ITG growth rate at perpendicular wavelength of the order of the ion gyroscale than the Coulomb collision operator.^{16} Additionally, previous collisional simulations of ZF demonstrate that the OS operator yields a stronger ZF damping than the Coulomb operator in both banana and Pfirsch–Schlüter regime.^{5,6}

To simulate highly collisional plasmas while still avoiding the use the Coulomb collision operator, Ref. 18 recently reported on the development of an operator that we refer to as the improved Sugama (IS) collision operator. The IS operator is designed to reproduce the same friction-flow relations (that we define below) of the linearized Coulomb collision operator by adding a correction term to the OS.^{18,19} While the IS operator has been successfully tested and implemented recently in neoclassical simulations using the GT5D code^{20} where like-species collisions are considered, no direct comparison between the IS, OS, and Coulomb collision operators has been reported yet on, e.g., microinstabilities and collisional zonal flow damping.

In this paper, we take advantage of recent analytical and numerical progress made in the development of GK collision operators based on a Hermite–Laguerre expansion,^{6} which we refer to as the gyro-moment approach, of the perturbed distribution function,^{21} and present the derivation and the expansion of the GK and drift-kinetic (DK) IS collision operators on the same basis. In particular, we leverage the gyro-moment expansion of the OS operator reported in Ref. 6, which is benchmarked with GENE.^{22} Despite that the gyro-moment expansion of the GK Coulomb is available in Ref. 6, simpler collision operator models are still important to develop accurate collisional descriptions of the plasma dynamics in the boundary region, yet simpler than the Coulomb operator. Using the Hermite–Laguerre approach, the integrodifferential nature of collision operator model is reduced to the evaluation of closed analytical expressions involving numerical coefficients that depend on the mass and temperature ratios of the colliding species and, when finite Larmor radius (FLR) terms are included, on the perpendicular wavenumber. The numerical implementation of the IS collision operator allows us to perform its comparison with the OS as well as the Coulomb collision operators on the study of instabilities and ZF damping for the first time. In particular, the linear properties of the trapped electron modes (TEMs) at steep pressure gradients, similar to H-mode conditions, are investigated and reveal that, indeed, the IS can approach better the Coulomb operator in the Pfirsch–Schlüter regime than the OS operator. Also, all operators yield results that agree within 10% at least for the parameters explored in this work. Nevertheless, larger deviations are expected in the case of heavy impurities^{23} with temperatures of the colliding species that can be significantly different. Additionally, we show that the IS yields a ZF damping intermediate between the OS and Coulomb collision operators in the Pfirsch–Schlüter regime. Finally, we evaluate the electrical Spitzer conductivity using the IS operator showing a good agreement with the Coulomb collision operator, while it is found that the OS operator underestimates the parallel electric current compared to the Coulomb operator by at least 10%. Taking advantage of the gyro-moment expansion, we explicitly evaluate the lowest order gyro-moments of the IS, OS, and Coulomb collision operators that can be used to model collisional effects in reduced gyro-moment models valid under the high-collisionality assumption.

The remainder of the present paper is organized as follows. In Sec. II, the IS collision operator is introduced. Then, in Sec. III, we derive the spherical harmonic expansion of the IS collision operator that allows us to evaluate the GK and DK limits of the same operator. In Sec. IV, we derive closed analytical expressions of the Braginksii matrices of the Coulomb and OS collision operators necessary for the evaluation of the correction terms added to the OS operator. Then, the gyro-moment method is detailed in Sec. V, where the Hermite–Laguerre expansion of the GK and DK IS collision operators are obtained analytically. In Sec. VI, numerical tests and comparisons are performed focusing on the TEM at steep pressure gradients, on the study of the collisional ZF damping in the Pfirsch–Schlüter regime, and on the Spitzer electrical conductivity. We conclude by discussing the results and future applications in Sec. VII. Appendixes A and B details the Coulomb and OS collision operators, and Appendix C reports on the analytical expressions of the lowest order gyro-moments of the Coulomb, OS, and IS collision operators, which are useful to derive reduced gyro-moment models for high collisional plasmas.

## II. IMPROVED SUGAMA COLLISION OPERATOR

We start by presenting the linearized IS collision operator following the notation and definitions of Ref. 18. We assume that the particle distribution function is perturbed with respect to a Maxwellian distribution of the particle of species *a*, $ f M a = f M a ( r , v ) = n a ( r ) / \pi 3 / 2 v T a ( r ) 3 / 2 e \u2212 v 2 / v T a ( r ) 2$, with $ n a ( r )$ the particle density, $ v T a ( r ) 2 = 2 T a ( r ) / m a$ the thermal particle velocity and $ z = ( r , v )$ the particle phase-space coordinates, being ** r** the particle position and $v$ the particle velocity such that $ v 2 = v \xb7 v$. The small-amplitude perturbation, $ f a = f a ( r , v )$, i.e., $ f a / f M a \u226a 1$, allows us to describe the collisions between species

*a*and

*b*by linearizing a nonlinear collision operator model $ C a b N L ( f a , f b )$. The linearized operator is denoted by $ C a b = C a b ( r , v ) = C a b ( f a , f b )$ and can be written as

where $ C a b T ( f a )$ and $ C a b F ( f b )$ are the test and field components of $ C a b ( f a , f b )$, and are defined from the nonlinear collision operator model, as $ C a b T ( f a ) = C a b N L ( f a , f M b )$ and $ C a b F ( f b ) = C a b N L ( f M a , f b )$.

Because the test and field components of the Coulomb collision operator, denoted by $ C a b L$ and defined in Appendix A, involve complex velocity–space derivatives of *f _{a}* [e.g., in $ C a b L T$ given in Eq. (A1)] and integrals

*f*[e.g., in $ C a b L F$ given in Eq. (A2)], approximated linearized collision operators have been proposed for implementation in numerical codes and analytical purposes in the past years.

_{b}^{9,11,13}Among these simplified models, the OS collision operator model,

^{13}that we denote as $ C a b S$, is widely used in present GK codes. The definitions of the test and field components of the OS collision operator, $ C a b S$, are introduced in Appendix A. The original Sugama operator, $ C a b S$, is derived from the linearized Coulomb collision operator, $ C a b L$, to conserve the three lowest order velocity moments of

*C*, i.e.,

_{ab}and satisfy the H-theorem and the adjointess relations even in the case of collisions between particles with temperatures $ T a \u2260 T b$, given by

and

*f*and

_{a}*g*are two arbitrary phase-space $ ( r , v )$ functions. The complete analytical derivation of $ C a b S$ is given in Ref. 13.

_{a}The difference between the OS and the Coulomb collision operators is expected to have a larger impact as the collisionality increases, in particular, in the case of particles with different mass and temperature.^{1} Hence, Ref. 17 proposes to improve the OS collision operator by adding a correction term to $ C a b S$, thus defining a new improved operator, which we refer to as the IS collision operator. The IS operator is designed such that it produces to same friction-flow relations,^{18} given by

than the Coulomb collision operator. In Eq. (5), $ s a 2 = v 2 / v T a 2$ is the energy coordinate and $ L i 3 / 2 ( x )$ (with $ i = 0 , 1 , 2 , \u2026$) is the associated Laguerre polynomial defined by

where $ L i l \alpha \u2212 1 / 2 = ( \u2212 1 ) l ( \alpha + i ) ! / [ ( i \u2212 l ) ! ( l + \alpha ) ! l ! ]$ (with $ \alpha > \u2212 1$). The derivation of the correction term added to $ C a b S$ can be found in Ref. 17 and we report the results here.

The linearized IS collision operator, denoted by $ C a b I S$, is obtained by adding the correction term $ \Delta C a b ( f a , f b )$ to $ C a b S$ (see Appendix A), i.e.,

where $ \Delta C a b$ is defined as

with the test and field components of the correction term being

and

respectively. Here, $ \tau \xaf a b = 3 \pi / ( 4 \nu a b )$ is the collisional time between the colliding species *a* and *b*, with *ν _{ab}* the associated collision frequency $ \nu a b = 4 \pi n b q a 2 q b 2 \u2009 ln \u2009 \Lambda / [ 3 m a 1 / 2 T a 3 / 2 ]$ (see Ref. 6). While the IS is derived in the limit $ L \u2192 \u221e$ and $ K \u2192 \u221e$, here we consider (

*L*,

*K*) to be two positive integers that we choose equal, i.e.,

*L*=

*K*. Since previous neoclassical transport studies suggest that accurate friction coefficients in the Pfirsch–Schlüter regime require that $ L = K \u2273 2$ (see Ref. 19), we consider the cases of $ L = K = 2$, 5 and 10. Despite that neoclassical studies have revealed that the energy ( $ \u223c v 2$) expansion in Laguerre polynomials has slow convergence in the banana regime,

^{24}we show here that small

*L*and

*K*are required at high collisionality in the test cases considered in Sec. VI. We note that the first neoclassical studies reported in Ref. 20 are performed with the like-species IS using only $ L = K = 1$, 2.

The quantities $ u a k ( f a )$ appearing in Eqs. (9) and (10) are defined as the flow vectors and are expressed by^{19}

with $ c k = 3 \xb7 2 k k ! / ( 2 k + 3 ) ! !$. Finally, we introduce the correction Braginskii matrix elements, $ \Delta M a b \u2113 k$ and $ \Delta N a b \u2113 k$, defined by

*A*=

*L*for the Coulomb and

*A*=

*S*for the OS operators) are obtained from the test and field components of the operator collision model $ C a b A T$ and $ C a b A F$, respectively, and are defined by

^{19,25}

The Braginskii matrices $ M a b A \u2113 k$ and $ N a b A \u2113 k$ satisfy a set of relations stemming from the conservation laws and symmetries of the collision operator. In particular, from the momentum conservation law in Eq. (2b), one obtains that

In addition, in the case of collisions between particle species with the same temperatures (*T _{a}* =

*T*), the Braginskii matrices admit symmetry properties because of the self-adjoint relations of $ C a b A T$ and $ C a b A F$ given in Eq. (4). From Eqs. (13) and (4), one obtains that

_{b}which implies from Eq. (12) that

for $ \u2113 , k = 1 , 2 , \u2026$ in the case *T _{a}* =

*T*. While the analytical expressions of the Braginskii matrices $ M a A \u2113 k$ and $ N a A \u2113 k$ up to $ ( \u2113 , k ) \u2264 2$ can be found in Ref. 17, we provide a new derivation of these coefficients here that allows us to extend them to arbitrary order $ ( \u2113 , k )$ in Sec. IV. This allows us to evaluate gyro-moment expansion of the IS for any (

_{b}*L*,

*K*). Using these analytical expressions, we demonstrate numerically in Sec. VI that the relations and symmetry properties of the Braginskii matrices are satisfied.

## III. SPHERICAL HARMONIC EXPANSION AND GYRO-AVERAGE OF THE IMPROVED SUGAMA OPERATOR

We expand the IS collision operator in terms of spherical harmonic particle moments in Sec. III A.^{6,26} This allow us to evaluate its gyro-average and to derive the GK and DK limits of the IS collision operators in Secs. III B and III C, respectively. We notice that, while the GK form of the IS collision operator is presented in Ref. 17, here we follow a different methodology, which is based in the spherical harmonic technique used to obtain the GK and DK Coulomb collision operators in Ref. 6. We note that the GK formulations of the Coulomb and of the OS collision operators used in this work are reported in Appendix A. Finally, we note that the spherical harmonic expansion is particularly useful in deriving the expressions of the Braginksii matrix elements in Sec. IV.

### A. Spherical harmonic expansion of the improved Sugama collision operator

The perturbed particle distribution function, $ f a ( r , v )$, is expanded in the spherical harmonic basis according to^{6,26}

with $ s a = v / v T a$ and $ \sigma j p = p ! ( p + j + 1 / 2 ) ! / [ 2 p ( p + 1 / 2 ) ! j ! ]$. In Eq. (20), the spherical harmonic basis is defined by $ Y p j ( s a ) = Y p ( s a ) L j p + 1 / 2 ( s a 2 )$, where we introduce the spherical harmonic tensors of order *p*, $ Y p$, and the associated Laguerre polynomials, $ L j p + 1 / 2 ( s a 2 )$, defined in Eq. (6). The tensor $ Y p$ can be explicitly defined by introducing the spherical harmonic basis, $ e p m$, which satisfies the orthogonality relation $ e p m \xb7 e p m \u2032 = ( \u2212 1 ) m \delta \u2212 m m \u2032$,^{27} such that

where $ Y p m ( \xi , \theta )$ are the scalar harmonic functions. The spherical harmonic basis, $ e 1 m$, can be related to the orthogonal Cartesian velocity–space basis $ ( e x , e y , e z )$ with $ e z = b$, by $ e 1 \u2212 1 = ( e x \u2212 i e y ) / 2 , \u2009 e 10 = b$ and $ e 11 = \u2212 ( e x + i e y ) / 2$ for *p* = 1, while $ e p m$ with *p* > 1 can be expressed in terms of $ e 1 m$ by using the closed formula^{27}

where $ N p m = ( p + m ) ! ( p \u2212 m ) ! 2 p \u2212 m / ( 2 p ) ! , \u2009 a n p m = p ! / [ 2 n n ! ( m + n ) ! ( p \u2212 m \u2212 2 n ) ! ]$, and $ { \u2026} S$ denotes the symmetric part. A complete introduction to the spherical harmonic tensors can be found in Ref. 27.

with $ T p$ an arbitrary *p*th order tensor. From the orthogonality relation in Eq. (23), it follows that the spherical harmonic particle moments $ M a p j ( r )$ are defined as follows:

Here, we emphasize that the spherical harmonic particle moments depend on the particle position ** r** only.

We now use the definition of the spherical harmonic particle moments, $ M a p j$ given Eq. (24), and relate them to the flow vectors expressed in Eq. (11), to obtain the spherical harmonic expansion of the IS collision operator. Therefore, from Eq. (11) and noticing that $ s a = Y 1 ( s a )$, we derive

### B. Gyrokinetic improved Sugama collision operator

We now consider the GK limit of the IS collision operator where the fast particle gyro-motion is analytically averaged out. Contrary to the IS collision operator, defined on the particle phase-space $ z = ( r , v )$ [see Eq. (7)], the GK IS collision operator, which we denote by $ C a b I S$, is defined on the gyrocenter phase-space coordinates $ Z = ( R , \mu , v \u2225 , \theta , t )$ where ** R** is the gyrocenter position, $ \mu = m v \u22a5 2 / [ 2 B ]$ is the magnetic moment and

*θ*is the gyroangle. More precisely, $ C a b I S$, is obtained by performing the gyro-average of the IS collision operator $ C a b I S$, i.e.,

where the integral over the gyroangle appearing in Eq. (26) is performed holding ** R** constant, while collisions occur at the particle position

**. We detail the transformation given in Eq. (26) in Appendix B that we use to obtain, for instance, the GK Coulomb operator. In general, the coordinate transformation that relates the gyrocenter and particles coordinates,**

*r***and**

*Z***, respectively, can be written as $ Z = z + \delta z$, where $ \delta z$ are functions of phase-space coordinates and perturbed fields and contain terms at all orders in the GK expansion parameter $ \epsilon \u223c e \varphi / T e$ ( $\varphi $ being the small amplitude and small-scale electrostatic fluctuating potential**

*z*^{6,28}) At the lowest order in the GK expansion, the coordinate transformation reduces to $ \delta v \u2225 = \delta \mu = \delta \theta = 0$ and $ \delta r \u2243 \u2212 \rho a$. Hence, the IS collision operator can be gyro-averaged holding

**constant in gyrocenter coordinates, such that $ r = R ( r , v ) + \rho a ( r , v )$.**

*R*^{21}

Focusing first on the test component of $ \Delta C a b T$, we perform the gyro-average in Eq. (26) in gyrocenter coordinates using $ r = R ( r , v ) + \rho a ( r , v )$ in the spatial argument of $ M a 1 k ( r )$, such that in Fourier space it yields $ M a 1 k ( r ) = \u222b d k M a 1 k ( k ) e i k \xb7 R e i k \xb7 \rho a$. For a single Fourier component, we derive

where we introduce the basis vectors $ ( e 1 , e 2 )$, such that $ v \u22a5 = v \u22a5 ( e 2 \u2009 cos \u2009 \theta \u2212 e 1 \u2009 sin \u2009 \theta )$. Additionally, $ J 0 a = J 0 ( b a x a )$ and $ J 1 a = J 1 ( b a x a )$ are the zeroth and first order Bessel functions (with $ b a = k \u22a5 \rho a$ the normalized perpendicular wavenumber and $ x a = v \u22a5 2 / v T a 2$), resulting from the presence of FLR effects in the IS collision operator. We remark that the equilibrium quantities in Eq. (27) (e.g., *f _{Ma}*,

*T*and $ \tau \xaf a b$) are assumed to vary on spatial scale lengths much larger that $ \rho a$. Therefore, they are evaluated at

_{a},**and are not affected by the gyro-average operator, in contrast to $ M a 1 k ( r )$.**

*R*We now relate the spherical harmonic moments $ M a 1 k$ of the perturbed particle to the gyrocenter perturbed distribution functions. More precisely, we express the $ M a 1 k$ in terms of the non-adiabatic part, *h _{a}*, of the perturbed gyrocenter distribution function

*g*. The two gyrocenter distribution functions,

_{a}*h*and

_{a}*g*, are related by

_{a}^{6}

in the electrostatic limit. Here, $ F M a = N a ( R ) / \pi 3 / 2 v T a 3 ( R ) e \u2212 v \u2225 2 / v T a 2 ( R ) \u2212 \mu B ( R ) / T a ( R )$ is the gyrocenter Maxwellian distribution function. The perturbed particle distribution function *f _{a}* is related to the perturbed gyrocenter distribution function

*g*by the scalar invariance of the full particle and gyrocenter distribution functions, i.e.,

_{a}Using the pull-back operator $T$, such that the functional forms of *f _{Ma}* and

*F*are related by $ f M a = T F M a$,

_{Ma}^{21}we derive that

We remark that, while both *g _{a}* and

*h*are gyrophase independent functions, $ g a g c$ is gyrophase dependent via its arguments

_{a}**. An expression of the pull-back transformation $T$ can be obtained at the leading order in the GK expansion parameter $ \epsilon \u223c e \varphi / T e$, yielding**

*Z*^{6}

being $ r = R + \rho a ( \mu , \theta )$. Using Eq. (31) allows us to finally express $ M a 1 k$ in terms of *h _{a}*,

where $ N a = \u222b d v \u2225 d \mu d \theta F M a$ is the equilibrium gyrocenter density. In the third line of Eq. (32), the lowest order guiding-center contribution to the gyro-center phase-space volume element is neglected such that $ B \u2225 * / m a = B ( 1 + v \u2225 b \xb7 \u2207 \xd7 b / \Omega a ) / m a \u2243 B / m a$, since the last term in $ B \u2225 *$ is of the order $ v \u2225 b \xb7 \u2207 \xd7 b / \Omega a \u223c \rho a / L B \u226a 1$ (with *L _{B}* the typical equilibrium scale length of the magnetic field

*B*). We remark that the contribution from the terms proportional to $\varphi $, appearing in Eq. (31), is neglected in Eq. (32). In fact, while these terms are of the same order as

*h*(i.e., they are order

_{a}*ε*), they yield a small contribution in the collision operator.

^{13}We neglect them here, but notice that their contributions to the Coulomb collision operator are included in Ref. 6 and have little effects at the gyroradius scale. With Eqs. (32) and (27), the GK test component of the correction term, $ \Delta C a b T = \u27e8 \Delta C a b T \u27e9 R$, can be obtained in terms of the gyrocenter distribution function

*h*. Focusing on the field component of $ \Delta C a b$, i.e., on $ \Delta C a b F$, we remark that a similar derivation of its expression can be carried out as for $ \Delta C a b T$. In particular, the expression of the spherical harmonic moment $ M b 1 k$, appearing in Eq. (10), is obtained with Eq. (32) having replaced

_{a}*a*with

*b*.

Finally, the GK IS collision operator can be expressed as follows:^{17}

where $ C a b S$ is the OS GK collision operator, given in Ref. 13, and $ \Delta C a b T$ and $ \Delta C a b F$ are given by

*f*, following the derivation of the full-F nonlinear GK Coulomb collision operator in Ref. 26 (see Appendix B).

_{a}### C. Drift-kinetic improved Sugama collision operator

We now derive the DK IS operator from the GK IS collision operator, given in Eq. (33), by neglecting the difference between the particle and gyrocenter position, such that $ r \u2243 R$. Hence, the DK IS collision operator is derived in the zero gyroradius limit of the GK IS operator approximating $ J 0 s \u2243 1 , \u2009 J 1 s = 0$ and $ f s \u2243 g s$ [see Eq. (28)]. The DK OS is given in Ref. 6, and the DK limits of the GK test and field components of $ \Delta C a b$ are

*g*, such that

_{s}The GK and DK expressions of the IS collision operator, given in Eqs. (34) and (36), can be implemented in continuum GK codes using a discretization scheme in velocity–space (e.g., based on a finite volume approach or a finite difference scheme) to evaluate numerically the velocity-integrals. In this work, we use a gyro-moment approach to carry out these velocity integrals and implement these operators numerically. For this purpose, we derive closed analytical expressions of the Braginksii matrices, required to evaluate the quantities $ \Delta M a b \u2113 k$ and $ \Delta N a b \u2113 k$, in Sec. IV.

## IV. BRAGINKSII MATRICES

In order to evaluate the correction term $ \Delta C a b$ given in Eq. (8), analytical expressions for the Braginskii matrices, $ M a b A \u2113 k$ and $ N a b A \u2113 k$, associated with the Coulomb and OS collision operators are derived. This extends the evaluation of the $ ( \u2113 , k ) \u2264 2$ Braginksii matrices reported in Ref. 17 to arbitrary $ ( \u2113 , k )$. For this calculation, we leverage the spherical harmonic expansions of the Coulomb and OS collision operators presented in Ref. 6. More precisely, we use the spherical harmonic expansion of *f _{a}* (see Sec. III A) to obtain the Braginskii matrices of the Coulomb and OS collision operators in Secs. IV A and IV B, respectively.

### A. Braginskii matrix of the Coulomb collision operator

We first derive the Braginskii matrix associated with the Coulomb collision operator, namely, $ M a b L \u2113 k$ and $ N a b L \u2113 k$, appearing in Eq. (12). For this purpose, we use the expansion of the perturbed particle distribution function, given in Eq. (20), to obtain the spherical harmonic expansion of the test and field components of the Coulomb collision operator, $ C a b L T$ and $ C a b L F$, derived in Ref. 6. These expressions, reported here, are given by

where the expressions of the test and field speed functions, $ \nu a b Tpj ( v )$ and $ \nu a b Fpj ( v )$, can be found in Appendix A of Ref. 6. We first use the expression in Eq. (39a) to evaluate $ M a b L \u2113 k$ [see Eq. (13a)]. We remark that the GK Coulomb collision operator, which we compare with the GK IS operator in Sec. VI, is derived in Ref. 6 by gyro-averaging Eq. (38) according to Eq. (26).

By noticing that $ Y 1 ( s a ) = v / v T a$ and by considering the component parallel to ** b** of Eq. (39a), we obtain the test part of the Coulomb collision operator evaluated with $ f M a m a v \u2225 L k 3 / 2 ( s a 2 ) / T a$, i.e.,

Then, the matrix element $ M a b L \u2113 k$, defined in Eq. (13a), can be computed by expanding the associated Laguerre polynomial using Eq. (6) with *p* = 1 and by performing the velocity integrals over the speed function $ \nu a b T 1 k ( v )$. It yields

A similar derivation can be carried out to evaluate $ N a b L \u2113 k$ [see Eq. (13b)] by using the expression in Eq. (39b), i.e.,

with $ \chi a b = v T a / v T b$ the ratio between the thermal velocities. The closed analytical expressions of test and field speed integrated functions $ \nu \xaf * a b T 1 k l$ and $ \nu \xaf * a b F 1 k l$ appearing in Eqs. (41) and (42), respectively, are given in Ref. 6. Equations (41) and (42) allows for the evaluation of the terms associated with the Coulomb collision operator appearing in the correction matrix elements, $ \Delta M a b \u2113 k$ and $ \Delta N a b \u2113 k$ defined in Eq. (12). They are evaluated in terms of mass and temperature ratios of the colliding species.

### B. Braginskii matrix of the original Sugama collision operator

We now evaluate the Braginskii matrix associated with the OS collision operator, namely, $ M a b S \u2113 k$ and $ N a b S \u2113 k$, appearing in Eq. (12). For $ f a = f M a m a v \u2225 L k 3 / 2 ( s a 2 ) / T a$, using the spherical harmonic expansion of the OS operator in Ref. 6, the test component of the OS collision operator yields

where the quantities $ X a b i$ are defined by

with $ \nu a b S 1 k ( v )$ the velocity-dependent speed function (whose expression is given in Ref. 6), being $ \nu \xaf a b \u2225 k = \u222b 0 \u221e d s a s a 2 k \nu a b \u2225 ( v ) e \u2212 s a 2$ (with $ \nu a b \u2225 ( v ) = 2 \nu a b [ erf ( s b ) \u2212 s b erf \u2032 ( s b ) ] / ( 2 s b 2 s a 3 )$ the velocity dependent energy diffusion frequency, and $ u \u2225 a k = 4 / ( 3 \pi ) ( k + 3 / 2 ) ! / k ! \delta k 0$. In deriving Eq. (43), we remark that the terms proportional to $ T a \u222b d v f a ( 2 s a 2 / 3 \u2212 1 ) / n a$ vanish exactly when applied to $ f a = f M a m a v \u2225 L k 3 / 2 ( s a 2 ) / T a$ because of the velocity integration over the pitch-angle variable $ v \u2225 / v$.

Similarly, the field component of the OS collision for $ f b = f M b m b v \u2225 L k 3 / 2 ( s b 2 ) / T b$ yields

where

with $ \Phi ( x ) = [ erf ( x ) \u2212 x erf \u2032 ( x ) ] / ( 2 x 2 )$. The velocity integral in Eq. (46) can be performed analytically using Eq. (6) and leads to

with the following definition:

where we introduce $ e b a k = \u222b 0 \u221e d s b s b 2 k erf \u2032 ( s a ) e \u2212 s b 2$ and $ E b a k = \u222b 0 \u221e d s b s b 2 k + 1 erf ( s a ) e \u2212 s b 2$. The test and field components of the OS collision operator, given in Eqs. (43) and (45), are now in a suitable form to evaluate the analytical expressions of $ M a b S \u2113 k$ and $ N a b S \u2113 k$ defined by Eq. (13).

Starting with the Braginksii matrix element associated with $ C a b T S$, the velocity integral in $ M a b Sjk$, given in Eq. (13a), is evaluated using the series expansion of the associated Laguerre polynomials, Eq. (6). Thus, we derive

where we introduce the quantities

In Eq. (50a), the analytical expression of the speed integrated function, $ \nu \xaf * a b S 1 k l$, is reported in Ref. 6. Similarly for $ N a b S \u2113 j$, using Eq. (45) and employing the expansion of the associated Laguerre polynomials in Eq. (6) yields

The Braginskii matrix elements associated with the Coulomb collision operator, given in Eqs. (41) and (42), and the ones associated with the OS operator, given in Eqs. (49) and (51), allow us to obtain the correction Braginskii matrix elements $ \Delta M a b \u2113 k$ and $ \Delta N a b \u2113 k$ for arbitrary $ ( \u2113 , k )$.

## V. GYRO-MOMENT EXPANSION OF THE IMPROVED SUGAMA COLLISION OPERATOR

We now project the GK and DK IS collision operators onto a Hermite–Laguerre polynomial basis, a technique that we refer to as the gyro-moment approach. Previous works^{6,16,29} demonstrate the advantage of the gyro-moment approach in modeling the plasma dynamics in the boundary region, where the time evolution of the gyro-moments is obtained by projecting the GK Boltzmann equation onto the Hermite–Laguerre basis yielding an infinite set of fluid-like equations.^{21} At high collisionality, high-order gyro-moments are damped such that only the lowest order ones are sufficient to evolve the dynamics. As a consequence, in these conditions, the gyro-moment hierarchy can be reduced to a fluid model where collisional effects are obtained using the Hermite–Laguerre expansion of advanced collision operators at the lowest order in the ratio between the particle mean-free-path to the parallel scale length. In Appendix C, we evaluate the lowest order gyro-moments of the DK IS, OS and Coulomb collision operators that enter in the evolution equations of the lowest order gyro-moments associated with fluid quantities. Ultimately, this allows us to compare analytically the fundamental differences between the IS, OS and the Coulomb collision operators. In addition, the closed analytical expressions reported in Appendix C can be used to derive high-collisional closures of the gyro-moment hierarchy.^{16}

Since the gyro-moment expansion of the OS GK and DK Sugama collision operator is obtained in Ref. 6, we focus here on the projections of the corrections $ \Delta C a b T$ and $ \Delta C a b F$, given in Eq. (34). The GK IS collision operator $ C a b I S$ is formulated in terms of moments of *h _{a}*, the non-adiabatic part of the perturbed gyrocenter distribution function

*g*[see Eq. (28)]. Therefore, we expand the collision operator in terms of gyro-moments of

_{a}*h*. More precisely,

_{a}*h*is written as

_{a}with $ s \u2225 a = v \u2225 / v T a$. In Eq. (52), the Hermite and Laguerre polynomials, *H _{p}* and

*L*, are defined via their Rodrigues's formulas $ H p ( x ) = ( \u2212 1 ) p e x 2 d p ( e \u2212 x 2 ) / d x p$ and $ L j ( x ) = e x / j ! d j ( e \u2212 x x j ) / d x j$, and are orthogonal over the intervals, $ [ \u2212 \u221e , \u221e ]$ weighted by $ e \u2212 x 2$, and $ [ 0 , + \u221e ]$ weighted by $ e \u2212 x$, respectively, such that

_{j}Because of the orthogonality relations in Eq. (53), the non-adiabatic gyro-moments of *h _{a}*, $ n a p j$, are defined by

where *N _{a}* is the gyrocenter density and

*F*the gyrocenter Maxwellian distribution function. We remark that the velocity-dependence contained in the gyrocenter phase-space Jacobian [neglected in Eq. (32)] can be retained in the gyro-moment expansion, carried out below, by introducing the modified gyro-moments $ n a * p j$, defined by Eq. (55) having replaced $ B / m a$ by $ B \u2225 * / m a$,

_{Ma}^{21,29}such that

### A. Expansion of the GK IS collision operator

We first derive the gyro-moment expansion of the GK IS collision operator in Eq. (33), which is

where the Hermite–Laguerre projection of $ C a b I S$ is defined by

and similar definitions are used for the remaining terms in Eq. (57), in particular,

and

Since the expression of $ C a b Spj$ is reported in Ref. 6, we focus here on the gyro-moment expansion of $ \Delta C a b T$ and $ \Delta C a b F$, defined in Eqs. (59) and (60), respectively. First, we derive the Hermite–Laguerre projection of the test component of the correction term $ \Delta C a b Tpj$. As an initial step, we express $ u \xaf \u2225 s k$ and $ u \xaf \u22a5 s k$, defined in Eq. (35) and appearing in Eq. (34), in terms of the non-adiabatic gyro-moments $ n a p j$. Injecting the expansion of *h _{a}* into $ u \xaf \u2225 s k$ and $ u \xaf \u22a5 s k$ yields

and

To analytically evaluate the velocity integrals in $ I \u2225 s pjk$ and $ I \u22a5 s pjk$, we expand the Bessel functions $ J 0 s$ and $ J 1 s$ in terms of associated Laguerre polynomials, $ L n m ( x s )$ [see Eq. (6)], as follows:^{30}

with the *n*th-order kernel function $ K n ( b s ) = ( b s / 2 ) 2 n e \u2212 b s 2 / 4 / n !$ describing FLR effects. Then, using Eq. (64) with

where *L _{f}* is the Laguerre polynomial defined in Eq. (6) with

*α*= 0 and $ d njf m$ are numerical coefficients whose closed analytical expressions are given in Ref. 6, the velocity integral in $ I \u2225 s pjk$ can be computed,

with $ [ \xb7 ]$ the Iverson bracket ( $ [ A ] = 1$ if *A* is true, and 0 otherwise). Finally, the velocity integral contained in $ I \u22a5 s pjk$ can be evaluated similarly to the one in Eq. (66) and yields

We remark that the expression of the numerical coefficients $ ( T \u2212 1 ) p j lkm$ in Eqs. (66) and (67) can be found in Ref. 26 and arise from the basis transformation from Hermite–Laguerre to associated Legendre–Laguerre polynomials.

We now have all elements necessary to focus on the evaluation of $ \Delta C a b Tpj$ obtained by projecting Eq. (9) onto the Hermite–Laguerre basis. In Eq. (59), one recognizes the quantities $ I \u2225 s pjk$ and $ I \u22a5 s pjk$ defined in Eqs. (66) and (67). Thus, using their definitions, the gyro-moment expansion of the test component of the correction terms, $ \Delta C a b Tpj$, is deduced,

We remark that $ \Delta C a b Tpj$ can then be expressed explicitly in terms of the non-adiabatic gyro-moments $ n a l k$ appearing in the definitions of $ u \u2225 a k$ and $ u \u22a5 a k$ given in Eq. (61).

Carrying out the same derivation yielding Eq. (68) for the gyro-moment expansion of the field component of the correction term, $ \Delta C a b Fpj$, defined in Eq. (60), and inverting the species role between *a* and *b* in $ u \u2225 s k$ and $ u \u22a5 s k$ yields

where $ u \xaf \u2225 b k$ and $ u \xaf \u22a5 b k$ can be expressed in terms of the non-adiabatic gyro-moments of *h _{b}* by using Eq. (61). With the gyro-moment expansion of the test and field components of the correction term, given in Eqs. (68) and (69), the GK IS collision operator, Eq. (57), is expressed in terms of the non-adiabatic gyro-moments $ n a p j$. We remark that, given the gyro-moment expansion of the GK IS collision operator in Eqs. (68) and (69), the test and field components of this operator, Eq. (34), can be recovered by applying the inverse transformation of Eqs. (59) and (60), respectively, e.g., for the test component

### B. Expansion of the DK IS collision operator

We now consider the gyro-moment expansion of the GK IS collision operator, Eq. (57), in the DK limit. The gyro-moment expansion of the DK IS collision operator can be obtained by projecting Eq. (36) or by taking explicitly the zeroth-limit of the gyro-moment expansion of the GK IS collision operator Eq. (57). In both cases, it yields the following Hermite–Laguerre projections of the test and field components of the correction term:

and

In Eqs. (71) and (72), the spherical harmonic moments $ M a 1 k$, defined by Eq. (20) and related to the flow vectors in Eq. (25), are expressed as a function of the gyro-moments of the perturbed distribution function *g _{a}*, namely, $ N a p j = \u222b d \theta d\mu d v \u2225 B g a H p L j / m a 2 p p !$ [see Eq. (55)], by

where we use the fact that $ f a \u2243 g a$ in the DK limit. The basis transformation coefficients, $ T 1 k g h$, are the DK basis transformation coefficients defined and derived in Ref. 29. With the gyro-moment expansion of the DK OS collision operator derived in Ref. 6, the DK IS collision operator follows by adding Eqs. (71) and (72) to the former. Finally, we remark that the DK IS operator can be recovered by applying the same inverse transformation as the one given in Eq. (70) by using Eqs. (71) and (72) instead.

## VI. NUMERICAL TESTS AND COMPARISON BETWEEN COLLISION OPERATORS

Using the gyro-moment expansion of the IS collision operator, we perform the first numerical tests and comparisons between the IS, OS, and the Coulomb collision operators. The discussion of the numerical results is organized as follows. First, we discuss the numerical implementation of the closed analytical formulas appearing in the IS operator, in particular of the Braginksii matrices in Sec. VI A. We show that the numerical implementation satisfies the conservation laws and associated symmetry properties of Braginskii matrices to machine precision, regardless of the values of (*L*, *J*). Then, as a first application of the IS collision operator using the gyro-moment approach, we investigate the collisionality dependence of TEM that develops at steep pressure gradients, such as those in H-mode pedestal, and compare the IS with the OS and Coulomb collision operators in Sec. VI B. The Coulomb and OS collision operators and their GK limits used in the present comparisons are detailed in Appendix B. In Sec. VI C, we perform tests to study the collisional ZF damping and compare the numerical results with analytical predictions. Finally, in Sec. VI D, we compare the parallel electrical Spitzer conductivities predicted by the OS, IS, and Coulomb operators. The numerical tests show that the IS collision operator approaches better the Coulomb operator than the OS collision operator in the Pfirsch–Schlüter regime. In all cases investigated, $ L = K \u2243 3$ terms in $ \Delta C a b$ are required for convergence.

### A. Numerical implementation

In the present work, the need of the numerical integration for the evaluation of the velocity integrals appearing in the definitions of $ M a b A \u2113 k$ and $ N a b A \u2113 k$, given in Eq. (13), and of the velocity integrals involved in the gyro-moment expansion of the IS operator, Eqs. (62) and (63), is removed by using the closed analytical expressions given, for instance, in Eqs. (66) and (67) to evaluate the velocity integrals contained in Eqs. (62) and (63). We note that the numerical error stemming from the numerical integration of the velocity integrals become arbitrary large as the order of polynomials increases, and thus alter the accuracy of the Hermite–Laguerre projection for high-order gyro-moments. However, the evaluation of a large number of numerical coefficients is still required in the gyro-moment approach, which is affected by cancelation and round-off errors. Therefore, the analytical expressions in the present work are evaluated using an arbitrary-precision arithmetic software.

In Fig. 1, we verify the momentum conservation law, given in Eq. (14), using the analytical expressions derived in Sec. IV using the arbitrary-precision arithmetic software with 50 significant digits (numerical experiments show that fewer digits could lead to round off numerical errors in the Braginskii matrices and, more generally, in the gyro-moment expansions) as a function of *k* for both the Coulomb and OS collision operators. The momentum conservation law is shown to be satisfied within an error of less than $ 10 \u2212 50$, illustrating the robustness of the numerical framework used in this work. Additionally, the Braginskii matrices, $ M a a A \u2113 k$ and $ N a a A \u2113 k$, of the Coulomb and OS collision operators are shown in Fig. 2 for the case of like-species collisions. Since the test components of both operators [see Eqs. (A1) and (A3)] are equal in the case of like-species collisions (and, more generally, in the *T _{a}* =

*T*case), the Braginskii matrices $ M a a L \u2113 k$ and $ M a a S \u2113 k$ are equal as shown in the left panel of Fig. 2, implying $ \Delta M a b \u2113 k = 0$. On the other hand, the difference in the Braginskii matrices associated with the field components, i.e., $ N a a A \u2113 k$ (right panel of Fig. 2), arises due to the difference in the field component of the OS with respect to the Coulomb collision operator.

_{b}The gyro-moments approach allows us to investigate the coupling between gyro-moments induced by the IS and OS collision operators. This can be done by truncating the Hermite–Laguerre expansion of *h _{a}*, see Eq. (52) (or

*g*in the DK limit) at $ ( p , j ) = ( P , J )$, therefore, assuming that higher-order gyro-moments,

_{a}*p*>

*P*and

*j*>

*J*, vanish. Given (

*P*,

*J*), the gyro-moment expansion of the collision operator

*A*can be written as

The coefficients $ C abp \u2032 j \u2032 ATpj$ and $ C abp \u2032 j \u2032 AFpj$ associated with the Coulomb (*A* = *L*) and OS (*A* = *S*) collision operators can be obtained from Ref. 6. In the case of the IS collision operator (*A* = *IS*), Eq. (74) becomes

In the case of the GK IS, the explicit expressions of the coefficients, $ \Delta C abp \u2032 j \u2032 Tpj$ and $ \Delta C abp \u2032 j \u2032 Fpj$, can be obtained by inserting the analytical expressions of $ u \xaf \u2225 s k$ and $ u \xaf \u22a5 s k$, given in Eq. (61), into Eqs. (68) and (69), respectively. In the case of the DK IS, the coefficients, $ \Delta C abp \u2032 j \u2032 Tpj$ and $ \Delta C abp \u2032 j \u2032 Fpj$, are obtained by using Eq. (73) (given *L* = *K*) into Eqs. (71) and (72), respectively. We note that Appendix C reports the closed analytical expressions of the lowest order coefficients, $ C abp \u2032 j \u2032 ATpj$ and $ C abp \u2032 j \u2032 AFpj$, associated with the DK Coulomb and OS collision operators as well as $ \Delta C abp \u2032 j \u2032 Tpj$ and $ \Delta C abp \u2032 j \u2032 Fpj$ in the high-collisional case where the gyro-moments is truncated by neglecting the moments with $ p + 2 j > 3$. In this case, the analytical expressions of $ \Delta C abp \u2032 j \u2032 Tpj$ and $ \Delta C abp \u2032 j \u2032 Fpj$ become independent of *K* (and *L*) when $ L = K > 1$ [see Eq. (73)].

For a numerical implementation, Eq. (74) can be recast in a matrix form by introducing the one-dimensional row index $ l \xaf ( p , j ) = ( J + 1 ) p + j + 1$ (where *p* and *j* run from 0 to *P* and *J*, respectively) and, similarly, the column index $ l \xaf ( p \u2032 , j \u2032 )$. This formulation allows us also to illustrate the coupling between the gyro-moments associated with the GK IS and GK OS operators, as well as their difference for the case of like-species (see Fig. 3). The same analysis is carried out also for the DK operators in Fig. 3. We observe that, because of the self-adjoint relations of like-species collisions [see Eq. (4)], the gyro-moment matrices are symmetric. We also observe the block structure of the GK IS and OS,^{6} a consequence of vanishing polynomial basis coefficients $ ( T \u2212 1 ) p j lkm$.

### B. Trapped electron mode in steep pressure gradient conditions

The study of microinstabilities appearing in steep pressure gradient conditions have gained large interest in the past years because of their role in determining the turbulent transport in H-mode pedestals.^{31–33} The linear properties of microinstabilities at steep pressure gradients can significantly differ from the one at weaker gradients, typically found in the core. For example, unconventional ballooning mode structures can be encountered if the pressure gradients are above a certain linear threshold with the location of the largest mode amplitude being shifted from the outboard midplane position, in contrast to the conventional mode structure found at lower gradients.^{34,35} Because the pedestal is characterized by a wide range of collisionalities ranging from the low-collisionality banana (at the top of the pedestal) to the high-collisionality Pfisch–Schlüter regime (at the bottom of the pedestal and in the scrape-off layer),^{36} an accurate collision operator is necessary for the proper description and interaction of these modes. Thus, we compare the properties of steep pressure gradient TEM, when the IS, the OS, and the Coulomb collision operators are used.

To carry out this numerical investigation, a linear flux-tube code using the gyro-moment approach has been implemented to solve the linearized electromagnetic GK Boltzmann equation, which we obtain from the full-F GK equation in Ref. 21,

where $ \chi a ( r ) = \varphi ( r ) \u2212 v \u2225 \psi ( r )$, with $\varphi $ the perturbed electrostatic potential and *ψ* the perturbed magnetic vector potential, being evaluated at ** r**, and given by the self-consistent GK quasineutrality and GK Ampère's law,

and

respectively. In Eqs. (76)–(78), we introduce the magnetic drift frequency $ \omega B a = v T a 2 ( x a + 2 s \u2225 a 2 ) \omega B / ( 2 \Omega a )$ with $ \omega B = ( b \xd7 \u2207 \u2009 ln \u2009 B ) \xb7 k$, the diamagnetic frequency $ \omega T a * = [ \omega N + \omega T a ( x a + s \u2225 a 2 \u2212 3 / 2 ) ]$ with $ \omega N = b \xd7 \u2207 \u2009 ln \u2009 N \xb7 k / B$ and $ \omega T a = b \xd7 \u2207 \u2009 ln \u2009 T a \xb7 k / B , \u2009 a a = b a 2 / 2$, and $ \Gamma 0 ( x ) = I 0 ( x ) e \u2212 x$ (with *I*_{0} the modified Bessel function). On the right-hand side of Eq. (76) is the linearized GK collision operator composed of the sum of the GK test ( $ C a b A T$) and field ( $ C a b A F$) components, such that $ C a b A = C a b A T + C a b A F$, which are defined in Appendix B.

While a detailed description of the resolution of Eq. (76) using the gyro-moment approach will be the subject to a future publication. Here, we mention that we assume concentric, circular and closed magnetic flux surfaces using the $ s \u2212 \alpha $ model (with *α* = 0).^{37} In the local flux-tube approach, we assume constant radial density and temperature gradients, $ L N \u2212 1$ and $ L T a \u2212 1$, with values $ R N = R 0 / L N = R 0 / L T i = R 0 / L T e = 20$, where *R*_{0} is the tokamak major radius. Electromagnetic effects are introduced with $ \beta e = 8 \pi P e 2 / B 0 2 = 0.01 %$. For numerical reasons, we use a larger electron to ion mass ratio $ m e / m i = 0.0027$. The local safety factor *q* and the inverse aspect ratio *ε* are fixed at *q* = 2.7 and $ \epsilon = 0.18$, while a magnetic shear *s* = 0.5 is used. We remark that the value of magnetic shear, smaller than typical edge values ( $ s \u2273 3$), is chosen to limit the additional computational cost related to the rapid increase in the radial wavenumber, *k _{x}*, with

*s*along the parallel direction. Additionally, we center the

*k*spectrum around

_{x}*k*= 0. Collisional effects are introduced by using the gyro-moment expansion of the IS operator derived in this work, and the gyro-moment expansions of the Coulomb and OS Sugama operators are reported in Ref. 6.

_{x}Given these parameters, we identify two branches of unstable modes developing at different values of the binormal wavenumber *k _{y}* in the collisionless limit. This is shown in Fig. 4 where the growth rate

*γ*and the real mode frequency

*ω*are plotted as a function of

_{r}*k*using $ ( P , J ) = ( 20 , 10 )$ gyro-moments. We remark that the collisionless GENE simulations are in excellent agreement with the gyro-moment approach, verifying its validity in the collisionless regime. A discontinuous frequency jump from positive to negative values indicates a mode transition, separating a branch of low-

_{y}*k*and high-

_{y}*k*modes. We also notice that the low-

_{y}*k*and the high-

_{y}*k*modes change continuously from electron to the ion diamagnetic directions as

_{y}*k*increases. While the mode peaking at $ k y = 0.25$ (

_{y}*k*is normalized to the ion sound Larmor radius) is normally identified as the ITG mode at lower gradients because of its propagation along the ion diamagnetic direction ( $ \omega r > 0$), we identity it here as a TEM with a conventional mode structure. Indeed, despite $ \omega r > 0$, the mode persists if the ion and electron temperature gradients are removed from the system, while it is stabilized if the electrons are assumed adiabatic. The discontinuity observed in Fig. 4 is due to a transition to a TEM developing at $ k y \u2273 0.5$ with an unconventional ballooning mode structures with secondary peaks located near $ \chi = \xb1 \pi / 2$ (

_{y}*χ*is the ballooning angle) away from the outboard midplane (where the mode at $ k y \u223c 0.2$ peaks), as shown in Fig. 5. Hence, we refer to the TEM peaking at $ k y = 0.25$, with a conventional mode structure, as the low-

*k*TEM and to the TEM peaking near $ k y = 0.6$, with an unconventional mode structure, as the high-

_{y}*k*TEM. We remark that the TEM modes identified in Fig. 4 have very similar features to the ones found in, e.g., Refs. 38–40.

_{y}We now investigate the collisionality dependence of both modes. While the electron collisionality is expected to be mainly in the banana regime with $ \nu e * = 2 q \nu / \epsilon 3 / 2 \u2272 1$ in the middle of the pedestal of present and future tokamak devices (with *ν* the ion–ion collision frequency normalized to $ c s / R 0$, see Ref. 6), the temperature drop yields collisionalities than can be in the Pfirsch–Schlüter regime, $ \nu e * \u223c 1 / \epsilon 3 / 2 \u226b 1$, at the bottom of the pedestal.^{36} Focusing first on the low-*k _{y}* TEM, we notice that since it develops at low binormal wavenumber ( $ k y = 0.25$), the DK limits of the collision operators are considered (numerical tests show that FLR effects at $ k y = 0.25$ change the value of the growth rate

*γ*by less than $ 1 %$ at the level of collisionalities explored). The growth rate and the real mode frequency obtained by using the DK Coulomb, DK OS and DK IS using

*K*= 2, 5, and 10 terms in Eq. (36) are shown as a function of $ \nu e *$ in Fig. 6 in the case of $ T i / T e = 1$. Here, we consider the case $ ( P , J ) = ( 16 , 8 )$ and the case where only six gyro-moments (6GM) are retained. For the 6GM model, the closed analytical expressions of the collision operators reported in Appendix C are used to evaluate the collisional terms. We first observe that the low-

*k*TEM is destabilized by collisions in the Pfirsh–Schlüter regime where both the growth rate and frequency increase with collisionality. We remark that the low-

_{y}*k*TEM is weakly affected by collisions when $ \nu e * < 1$, with collisions that have a stabilizing effects on the mode when $ R N < 20$. Second, it is remarkable that the deviation (of the order or smaller than $ 10 %$) between the DK OS and DK Coulomb increases with $ \nu e *$, while the DK IS is able to correct the DK OS and to approach the DK Coulomb as

_{y}*L*and

*K*increase. This is better shown in the insets of Fig. 6 where a good agreement between the DK IS and Coulomb operator is observed. Third, the 6GM model (cross markers) approach the results of the full calculations, i.e., with $ ( P , J ) = ( 16 , 8 )$, as the collisionality increases, while they deviate from each other at low collisionality because of kinetic effects that are not resolved in the 6GM model. In addition, we note the good agreement between the DK IS and Coulomb in the 6GM model in the Pfirsch–Schlüter regime, demonstrating the robustness of the present numerical results.

To analyze the relative difference of the IS and OS operators with respect to the DK Coulomb, we scan the low-*k _{y}* TEM mode over the density gradient

*R*(with the temperature gradients, $ R 0 / L T i = R 0 / L T e = R N$) for different collisionalities and compute the signed relative difference of the growth rate

_{N}*γ*, $ \sigma ( \gamma ) = ( \gamma \u2212 \gamma C ) / \gamma C$, where

*γ*is obtained using the OS or IS collision operator models and

*γ*is the one obtained using the DK Coulomb operator (see Fig. 7). The same definition $ \sigma ( \omega r ) = ( \omega r \u2212 \omega r C ) / \omega r C$ is used for the real mode frequency. We plot the results in Fig. 8. First, we observe that the DK OS operator underestimates both

_{C}*γ*and

*ω*(compared to the DK Coulomb) at low collisionality (with a peak near $ \nu e * \u223c 1$), and that difference changes sign in the Pfirsch–Schlüter regime, when $ \nu e * \u2273 10$, where the DK OS operator overestimates

_{r}*γ*and

*ω*. A difference of the order of $ 5 %$ is found in growth rate and of the order of $ 10 %$ in the frequency. While these deviations increase with $ \nu e *$ (see the red areas in the left panels in Fig. 8), the correction terms in the IS operator reduce $ \sigma ( \gamma )$ [and $ \sigma ( \omega r )$] below $ 2 %$ for all density gradients at high collisionalities. More precisely, the agreement between the DK IS and DK Coulomb improves with

_{r}*K*in the Pfirsch–Schlüter regime. In fact, the deviations from the DK Coulomb operator observed by the presence of the red area in $ \sigma ( \gamma )$ in the case

*K*= 2 are reduced in the case of

*K*= 5 (and

*K*= 10) when $ \nu e * \u2273 10$. Additionally, the small differences observed between the DK IS with

*K*= 10 and

*K*= 5 show that the results of the IS operator are converged when $ K \u2273 5$. We notice that, in general, $ \sigma ( \gamma ) \u2272 \sigma ( \omega r )$ for all operators. Finally, we remark that, in both the DK OS and DK IS operators, the deviations in

*γ*and

*ω*peak near $ \nu e * \u223c 1$ and increase at lower gradients with $ \sigma ( \gamma ) \u223c 5 %$ and $ \sigma ( \omega r ) \u223c 10 %$. The larger deviations are explained by the fact that the OS and, hence, the IS operators are based on a truncated moment approximation of the Coulomb collision operator (see, e.g., the field component) and that the effects of high-order moments in the collision operator models can no longer be ignored at this intermediate level of collisionality. On the other hand, in the Pfirsch–Schlüter regime, the contribution of high-order moments becomes small as being damped by collisions. The increase in their relative differences with decreasing density gradient,

_{r}*R*, is mainly attributed to the decrease in

_{N}*γ*and

*ω*, as they are of the order of the diamagnetic frequency.

_{r}The deviation between the DK OS and DK Coulomb operators depends on the temperature ratio of the colliding species (as well as on the mass ratio). In order to study the impact of the ion to electron temperature ratio, we consider $ \sigma ( \gamma )$ plotted as a function of the temperature ratio $ T i / T e$ and shown in Fig. 9 in the Pfirsch–Schlüter regime when $ \nu e * = 50$. It is confirmed that the correction terms [see Eq. (7)] enable the IS operator to approximate the Coulomb collision operator better than the OS operator as $ T i / T e$ increases, with $ \sigma ( \gamma ) \u2272 1 %$. The same observations can be made for the real mode frequency, *ω _{r}*.

We now turn to the collisionality dependence of the high-*k _{y}* TEM mode developing near $ k y = 0.6$ (see Fig. 4) using the GK collision operators. In particular, we use the GK Coulomb, GK OS, and GK IS operators, using the spherical harmonic expansion detailed in Appendix B. As the perpendicular wavenumber in the argument of the Bessel functions increases, an increasingly large number of terms in the infinite sums arising from the expansion of the Bessel functions, Eq. (64), is required for convergence.

^{16}We evaluate numerically these sums in $ I \u2225 s pjk$ and $ I \u22a5 s pjk$, given in Eqs. (66) and (67), by truncating them at

*n*= 6 (we have verified the convergence of our results). The collisionality dependence of the high-

*k*TEM growth rate,

_{y}*γ*, and mode frequency,

*ω*, obtained by using the same parameters in Fig. 4 (for $ k y = 0.6$) are shown in Fig. 10, as a function of the electron collisionality $ \nu e *$ using the GK Coulomb, GK OS and GK IS with different values of

_{r}*L*=

*K*. We also show the results of the DK Coulomb to illustrate the effects of FLR terms in the collision operators. First, we observe that the high-

*k*TEM is stabilized by the GK operators compared to the DK Coulomb (and other DK operators) because of the presence of the FLR terms in the former. Second, it is remarkable that, while the GK OS operator is able to capture the trend of the growth rate and of the mode frequency observed with the GK Coulomb, it yields systematically a smaller growth rate at all collisionalities. We remark that, while a direct comparison between the GK operators implemented in the GENE code is outside of the scope of the present work, similar observations are made in TEM simulations at weaker gradients

_{y}^{4,5}based on the GENE code. Third, and finally, we observe that the GK IS yields a growth rate similar (yet smaller) to the GK OS collision operator. In particular, the GK IS with

*K*= 2 reproduces the same growth rate as the GK OS, while larger values of

*K*, i.e.,

*K*= 5 and 10 (convergence is achieved with $ L = K \u2273 3$), produce a smaller growth rate than the GK OS within $ 5 %$ (we have carefully verified that our numerical results are converged by increasing the number of gyro-moments and the number of points in the parallel direction in the simulations). On the other hand, the real mode frequency,

*ω*, predicted by the GK Coulomb is well retrieved by the GK IS operators independently for

_{r}*K*, while it is underestimated by the GK OS (see Fig. 10). We remark that, in the DK limit, the IS operator yields very similar growth rates and real mode frequencies than the DK Coulomb at high collisionalities. To understand these observations, we first note that, in the case considered here where

*T*=

_{i}*T*, the test components of the GK OS and GK Coulomb are equivalent [see Eqs. (A1) and (A3)] since all the terms $ \Delta M a b \u2113 k$ in $ \Delta C a b T$ [see Eq. (12)] vanish exactly. This implies that the GK OS and GK IS differ only by the GK correction terms in their field components, i.e., by $ \Delta C a b F$ given in Eq. (34b). To illustrate the contribution from the FLR terms, we repeat the simulations in Fig. 10 considering the DK limits of the field components in all operators, but retain the GK test components being equivalent to the one of the GK Coulomb operator. The results are displayed in Fig. 11 and show that the IS operator yields a growth rate larger and closer to the GK Coulomb compared than the GK OS in the absence of FLR terms in the field components, while (not shown) the real mode frequency agrees between the GK IS and GK Coulomb operators. We remark that the high-

_{e}*k*TEM is strongly damped at high collisionality if the FLR terms in the field component are neglected. Overall, the GK IS and GK OS yield a similar collisionality dependence of the high-

_{y}*k*TEM with a good agreement in the mode frequency between the GK OS and GK Coulomb, despite that the growth rate predicted by the GK IS and GK OS differ from the GK Coulomb within $ 10 %$.

_{y}### C. Collisional zonal flow damping

Axisymmetric, poloidal zonal flows (ZFs) are believed to be among the key physical mechanisms at play in the L–H mode transition by ultimately, regulating the level of turbulent transport.^{41} It is, therefore, of primary importance to test and compare the effect of collisions, modeled by using the IS, OS and Coulomb collision operators, on their dynamics.

While the damping to a residual level of the ZFs has been originally studied in the collisionless regime,^{42} the later study was extended to the assess the collisional ZF damping in the banana regime when $ \nu i * \u2272 1$ (with $ \nu i *$ the ion–ion collisionality) for radial wavelengths much longer than the ion poloidal gyroradius, assuming adiabatic electrons, and using a pitch-angle scattering operator mimicking the collisional drag of energetic ions.^{43} A refinement of the exponential decay in Ref. 43 of the ZF residual prediction $ R z ( \u221e ) = \varphi z ( \u221e ) / \varphi z ( 0 )$ (with $ \varphi z ( t )$ the flux-surface averaged potential) with a momentum restoring pitch-angle scattering operator was later derived in Ref. 44 for long-wavelength modes,

Even if Eq. (79) does not include energy diffusion as well as the effects of GK terms in the collision operator, it still provides a good estimate asymptotic ZF residual predicted by the DK IS, as shown below.

For our tests, we consider only ion–ion collisions and, by including the IS operator, we extend the study of Ref. 6, where the differences between the GK OS and GK Coulomb operators in the collisional ZF damping are illustrated, finding a stronger damping by the former. We focus on the Pfirsch–Schlüter regime with $ \nu i * = 3.13$, with convergence being achieved with $ ( P , J ) = ( 24 , 10 )$ gyro-moments. The collisional time traces of the ZF residual, $ R z ( t ) = \varphi z ( t ) / \varphi z ( 0 )$, obtained for the DK IS (with $ K = L = 2 , 5 , 10$), DK OS, and DK Coulomb collision operators are shown in Fig. 12 for $ k x = 0.05$ and in Fig. 13 for the $ k x = 0.1$ and $ k x = 0.2$ using the GK operators. Starting from the case of small radial wavenumber $ k x = 0.05$, we first observe that all DK operators agree with the analytical long time prediction given in Eq. (79), despite the absence of energy diffusion in the latter. Consistently with Ref. 6, the DK Sugama yields a stronger damping of the ZF than the DK Coulomb. The addition of the correction terms to the OS operator allows the DK IS to better approximate the DK Coulomb operator, yielding a weaker damping of the ZF. We remark that only a small difference between the $ L = K = 2 , 5$ cases is noticeable showing that $ L = K \u2243 3$ is necessary for the DK IS to converge also in this case.

We now consider the collisional ZF damping at larger *k _{x}* values using the GK IS, GK OS, and GK Coulomb collision operators, $ k x = 0.1$ and $ k x = 0.2$, and plot the results in Fig. 13. The same collisionality of the $ k x = 0.05$ case is used. Only the

*K*= 2, 5 cases are considered for the GK IS for simplicity, since convergence is achieved with these parameters. First, consistently with Ref. 6, it is observed that the GK OS produces a stronger ZF damping with respect to the GK Coulomb. We note that a similar observation based on the results of the GENE code is reported in Refs. 4 and 5 (the collisional damping predicted by the gyro-moment method and the GENE code using the GK OS

^{13}is successfully benchmarked in Ref. 6). Second, the GK IS collision operator provides a better approximation to the GK Coulomb collision operator than the GK OS operator. This is particularly true in the early phase of the damping, i.e., $ t \nu \u2272 10$. Although the GK IS still yields a better approximation than the GK OS, larger GK IS departures from the GK Coulomb are observed at later times when the GAM oscillations are completely suppressed. For larger collisionalities and smaller scale lengths, the same observations about the rapid ZF decay made at lower collisionalities hold.

Finally, we investigate the effects of the GK collision operators on the ion velocity–space distribution function. We consider the modulus of the perturbed ion distribution function, $ | g i |$, obtained by using Eq. (52), at time $ t \nu = 5$ after the damping of the GAM oscillations, for the case $ k x = 0.2$. We plot $ | g i |$ as a function of $ s \u2225 i$ (at *x _{i}* = 0) and

*x*(at $ s \u2225 i = 0$) in Fig. 14, respectively, for the different GK collision operators. It is observed that the GK IS yields similar velocity–space structures along both the parallel and perpendicular directions than the GK OS operator with the distribution function being depleted in the region of the velocity–space $ | v \u2225 | \u2272 v T i$ more strongly than the GK Coulomb operator. A similar observation can be made when $ v \u22a5 \u2272 v T i$, as shown in the right panel of Fig. 14.

_{i}### D. Spitzer electrical conductivity

As a final numerical test, we consider the evaluation of the Spitzer electrical conductivity. This is obtained from the stationary current resulting from the balance between the collisional drag and a constant electric force. We consider an unmagnetized fully ionized plasma with a fixed background of ions with electrical charge $ q i Z$ (with *Z* the ion ionization degree), subject to an electric field, $ E = E e z$, with *E* a constant amplitude. Hence, the kinetic equation describing the evolution of the electron distribution function, *f _{e}*, is given by

with $ v \u2225 = v \xb7 e z$. Collisional effects between electrons ( $ C e e A T$ and $ C e e A F$) and the stationary background of ions ( $ C e i A T$) are modeled using the Coulomb, OS, and IS collision operators. Projecting Eq. (81) onto the Hermite–Laguerre basis yields the evolution equation of the gyro-moments of *f _{e}* denoted by $ N e p j$, i.e.,

We remark that the term associated with electric force, proportional to $ p E N e p \u2212 1 j$, vanishes when *p* = 0. We evolve the gyro-moment hierarchy, given in Eq. (81), with $ ( P , J ) = ( 20 , 5 )$ gyro-moments in time until a stationary electron current, $ j \u2225 e = \u2212 e N e u e$ (with $ u e = \u222b d v v \u2225 f e / N e = N e 10 v T e / 2$ the parallel electron velocity), is established resulting from an applied electric field of normalized amplitude $ e E / [ m e T e \nu e e ] = 10 \u2212 3$ (see Fig. 15). From the saturated current, the electrical conductivity, $ \sigma \u2225 e = j \u2225 e / E$, can be computed.

In Fig. 16, we compare our numerical estimates of $ \sigma \u2225 e$ with the analytical prediction of the Spitzer conductivity,^{25,45} given by $ \sigma \u2225 e = 16 T e 3 / 2 / [ ( 2 \pi ) 3 / 2 m e Z 2 e 2 \u2009 ln \u2009 \Lambda ] )$ obtained in the large ion ionization *Z* limit. We note that in the $ Z \u2192 \u221e$ limit, the electron and ion collisions can be modeled by the pitch-angle scattering collision operator, $ C e i ( f e ) \u2243 \u2212 \nu e i D ( v ) L 2 f e$ [where $ \nu e i D ( v )$ is the velocity-dependent deflection collision frequency, defined below Eq. (A4)], and the collisions between electrons can be neglected since $ \nu e e / \nu e i \u223c 1 / Z 2$. First, we observe that the OS operator produces an electrical conductivity that is approximately $ 10 %$ smaller than the one obtained by the Coulomb operator at low *Z*, and that the difference between the predictions of the IS and Coulomb operators is less than $ 1 %$ for all values of *Z*. While the deviations in the electrical conductivity decrease with *Z* because the contribution of the pitch-angle scattering that dominates with *Z* in all operators, the deviation at low *Z* observed between OS and Coulomb is associated with the different field components of these operators, i.e., $ C e e L F$ and $ C e e S F$, given in Eqs. (A2) and (A8), respectively, since the test components of these operators are equal in the case of like-species collisions. Second, we observe that all operators converge to the analytical $ Z \u2192 \u221e$ Spitzer conductivity as *Z* becomes large, within less than $ 8 %$ for *Z* = 10, providing a verification of the numerical implementation. Given the importance of the plasma resistivity in setting the level of turbulent transport in the scrape-off-layer,^{46} the present test shows that the OS operator underestimates the parallel current (see Fig. 15), which might have a significant effect on boundary turbulent simulations.

## VII. CONCLUSION

In this work, the gyro-moment method has been applied to implement the recently developed improved Sugama (IS) collision operator in both the gyrokinetic (GK) and drift-kinetic (DK) regimes.^{17} Designed to extend the validity of the original Sugama (OS) collision operator to the Pfirsch–Schlüter regime, the Hermite–Laguerre expansion of the perturbed distribution function allows expressing the IS collision operator as a linear combination of gyro-moments with coefficients that are analytical functions of the mass and temperature ratios of the colliding species and, in the GK limit, of the perpendicular wavenumber. Analytical expressions of the Braginskii matrices, $ M a b A \u2113 k$ and $ N a A \u2113 k$, associated with the Coulomb and OS collision operators are obtained for arbitrary $ ( \u2113 , k )$ using the spherical harmonic expansion.^{6} This allows the evaluations of the correction terms that are added to the OS operator yielding the IS collision operators for arbitrary gyro-moments.

We describe the numerical implementation of the IS collision operator. This is based on an arbitrary precision arithmetic library to avoid the numerical loss of precision and round-off errors when evaluating the Braginskii matrices. We demonstrated that the conservation laws (particle, momentum, and energy) are satisfied at arbitrary order in $ ( \u2113 , k )$. The IS collision operator is tested and compared with the OS and Coulomb collision operator in edge conditions at steep pressure gradients and high collisionality. Three test cases are considered, which involve the evaluation of the TEM growth rate and frequency at steep gradients, the collisional damping of ZF and, finally, the evaluation of the Spitzer electrical conductivity. The analysis of the linear properties of a conventional TEM, developing at long wavelengths, reveals that the IS is able to approximate the TEM growth rate and frequency predicted by the DK Coulomb operator better than the OS, particularly in the Pfirsch–Schlüter regime. At small perpendicular wavelengths and high collisionality, the study an unconventional TEM shows that the GK IS and GK OS essentially yields the same FLR damping, within a $ 10 %$ of difference. The collisional ZF damping is also explored and compared with analytical results. The analysis shows that the IS operator yields an intermediate value of the ZF residual, between the weaker value produced by the OS operator and the larger ZF residual predicted by the Coulomb collision operator. Finally, the electrical conductivity predicted by the IS operator is found to be in good agreement with the one produced by the Coulomb operator (within less than $ 1 %$), while deviations larger than $ 10 %$ are obtained when using the OS operator. In conclusion, the IS and OS produce linear results that differ by around $ 10 %$ in all cases explored in the present work. Nonlinear turbulent simulations are required to investigate the impact of these differences on the saturated turbulent state, in particular, near the nonlinear stability threshold.

We demonstrate that the computational cost of the gyro-moment approach decreases with collisionality, as illustrated in the case of TEM at steep gradients (see Fig. 6). In a future publication, we show that the gyro-moment approach is particularly efficient in steep pressure gradient conditions, where the number of gyro-moments can be significantly reduced compared to weaker gradients (typically found in the core) even at low collisionality. This is particular relevant to H-mode pedestal applications and offers an ideal framework to construct reduced fluid-like models to explore turbulent transport in the boundary of fusion devices. As an example, the analytical expression derived in this work allow us to evaluate the lowest order gyro-moment terms of the IS, OS and Coulomb collision operators in Appendix C. Finally, we remark that the expressions presented in this work can be easily extended to the study of multicomponent plasmas with species of different mass and temperatures.

## ACKNOWLEDGMENTS

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

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Baptiste Jimmy Frei:** Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). **Samuel Ernst:** Formal analysis (equal); Software (equal). **Paolo Ricci:** Supervision (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: LINEARIZED COULOMB AND ORIGINAL SUGAMA COLLISION OPERATORS

The linearized particle Coulomb collision operator, denoted by $ C a b L$, is defined as the sum of its test and field components, $ C a b L T$ and $ C a b L F$, respectively, given in pitch-angle coordinates $ ( r , v , \xi , \theta )$ (with *v* the modulus of ** v** and $ \xi = v \u2225 / v$ the pitch angle) by

^{8}

and

respectively. In Eq. (A1), we introduce the Rosenbluth potentials $ H ( f ) = 2 \u222b d v \u2032 f ( v \u2032 ) / u$ and $ G ( f ) = \u222b d v \u2032 f ( v \u2032 ) u$ (with $ u = | v \u2212 v \u2032 |$) and the operator $ L 2 f = \u2202 v ( v 2 \u2202 v f ) \u2212 v 2 \u2207 v 2 f$. When evaluated with a Maxwellian distribution function, analytical expressions of the Rosenbluth potentials can be obtained and are given by $ G ( f M b ) = n b v T b [ ( 1 + 2 s b 2 ) erf ( s b ) / s b + erf \u2032 ( s b ) ] / 2$ and $ H ( f M b ) = 2 n b erf ( s b ) / v$ [with the error function, $ erf ( x ) = 2 \u222b 0 x d s e \u2212 s 2 / \pi $, and its derivative, $ erf \u2032 ( x ) = 2 e \u2212 x 2 / \pi $]. We remark that the velocity–space derivatives and integrals, appearing in Eqs. (A1) and (A2), act on the perturbed particle distribution function *f _{a}* and are evaluated holding

**constant, such that $ C a b L T ( f a ) = C a b L T ( f a ( r , v ) ) = C a b L T ( r , v )$ [and the same for $ C a b L F ( f b )$].**

*r*The original particle Sugama (OS) collision operator, denoted by $ C a b S$, is expressed by the sum of its test and field components, $ C a b S T$ and $ C a b S F$, given by^{13}

where

and

and by

respectively. In Eq. (A3), we introduce the velocity-dependent collision frequencies, $ \nu a b D ( v ) = \nu a b [ erf ( s b ) \u2212 \Phi ( s b ) ] / s a 3$ and $ \nu a b \u2225 ( v ) = 2 \nu a b \Phi ( s b ) / s a 3$ [with $ \Phi ( x ) = [ erf ( x ) \u2212 x erf \u2032 ( x ) ] / ( 2 x 2 )$ being the Chandrasekhar function], the fluid quantities, $ u a ( f a ) = \u222b d v v f a / n a , \u2009 \delta T a ( f a ) = T a \u222b d v f a ( 2 s a 2 / 3 \u2212 1 ) / n a$, the coefficients $ \theta a b = ( T a / T b + \chi a b 2 ) / ( 1 + \chi a b 2 ) , \u2009 \tau = T a / T b , \u2009 \chi a b = v T a / v T b = \tau / \sigma $, and $ \sigma = m a / m b$. Finally, in Eq. (A8), we have

*T*=

_{a}*T*such that the test component of the OS operator reduces to $ C a b S T ( f a ) = C a b 0$, with $ C a b 0$ being the test component with the linearized Coulomb collision operator given in Eq. (A1). Hence, the OS operator differs from the Coulomb collision operator only by its field component, given in Eq. (A8), when

_{b}*T*=

_{a}*T*. Similarly to the linearized Coulomb collision operator, we remark that the velocity–space derivatives and integrals contained in $ C a b 0$ and $ C a b S F$ given in Eqs. (A4) and (A8), respectively, act on the perturbed particle distribution function

_{b}*f*.

_{a}### APPENDIX B: GK FORMULATION OF COLLISION OPERATORS

We now discuss the GK formulations of the GK Coulomb and GK OS operators based on the spherical harmonic technique used in this work. The GK linearized collision operator associated with the operator *C _{ab}* [see Eq. (1)], $ C a b$, is given by the sum of the GK test and field components, i.e.,

and

where we note that $ C a b T = C a b T ( z ) = C a b T ( f a ( z ) , f M b )$ and, similarly, $ C a b F = C a b F ( z ) = C a b F ( f M a , f b ( z ) )$. In Eqs. (B1) and (B2), we define the push-forward operator, $ T \u2212 1$ (inverse of $T$) that allows us to express a scalar function, such as $ f ( z )$, defined on the particle phase-space ** z**, in terms of the gyrocenter phase-space coordinates

**, i.e., $ f ( z ( Z ) ) = ( T \u2212 1 f ) ( Z ) = f ( T \u2212 1 Z )$, being $ z ( Z ) = T \u2212 1 Z$.**

*Z*In the spherical harmonic expansion, $ C a b T ( r , v )$ and $ C a b F ( r , v )$ are expressed as linear combinations of particle fluid moments, $ M a p j ( r )$, that are written in terms of *h _{a}*, i.e.,

^{6,26,29}

with $ R ( r , v ) \u2243 r \u2212 \rho a$. We notice that Eq. (B3) reduces to Eq. (32) with *p* = 1. The analytical expression of the $ M a p j$ moments in terms of gyro-moments, $ n a p j$, is detailed in Ref. 6. Using Eq. (B3) into the spherical harmonic expansion of the collision operator, Eq. (38), allows us to obtain the linearized collision operator (Coulomb or OS operators) in the particle phase-space $ z = ( r , v )$ in terms of the gyrocenter distribution function *h _{a}*.

In order to evaluate the gyro-average in Eqs. (B1) and (B2), the collision operator is first transformed to gyrocenter coordinates using the push-forward operator $ T \u2212 1$. Because the collision operator is a scalar phase-space function defined on the particle phase-space, i.e., $ C a b T = C a b T ( z )$, it transforms as

where we use the lowest order gyrocenter coordinate transformations, in particular $ r ( Z ) \u2243 R + \rho a ( \mu , \theta )$. Equation (B4) allows the expression of the linearized collision operator in gyrocenter phase-space coordinates ** Z**, as needed to perform analytically the gyro-average. For instance, evaluating the gyro-average of the test component of the Coulomb operator for a single Fourier component using Eq. (38), where the particle position

**appearing in $ M a p j$ is transformed to $ R + \rho a ( \mu , \theta )$ according to Eq. (B4), yields**

*r*where Eq. (21) and the Jacobi–Anger identity, $ e i k \xb7 \rho a = \u2211 n i n J n ( b a x a ) e i n \theta $, are used. A similar derivation for the GK field component can be performed. Finally, the gyro-moment expansion of the GK Coulomb collision operator that we use in this work is derived by projecting Eq. (B5) onto the Hermite–Laguerre basis, as detailed in Ref. 6. The gyro-averaging produces FLR terms, associated with the spatial shift $ \rho a ( \mu , \theta )$ and the *θ*-dependence of $ Y p m ( \xi , \theta )$, in both the test and field components of the linearized GK collision operators, that appear through Bessel functions of all orders, i.e., $ J n ( b a x a ) = \u2211 j = 0 \u221e c j ( b a x a / 2 ) 2 j + n$, with $ c j = ( \u2212 1 ) j / [ j ! ( n + j ) ! ]$. These FLR terms ultimately damp small scale fluctuations,^{16} and vanish when the DK limit is considered.

We remark that the spherical harmonic expansion allows us to evaluate the velocity derivatives (and integrals contained in the Rosenbluth potentials^{6}) exactly and express them as linear combination of $ M a p j ( r )$ given in Eq. (B3). For instance, the velocity derivative contained in the test component is computed using the spherical harmonic expansion, Eq. (20), and evaluated in gyrocenter coordinates, $ Z = ( R , \mu , v \u2225 , \theta )$, as follows:^{26}

with $ r ( Z ) \u2243 R + \rho a$. A similar expression as Eq. (B6) can be derived for the second order velocity derivatives contained in Eq. (A1).

On the other hand, in previous formulations^{13,47} and numerical implementations^{4,48} of GK collision operators, the velocity derivatives and integrals contained in the particle collision operators are evaluated in terms of the lowest order gyrocenter coordinates, $ Z = ( R , \mu , v \u2225 , \theta )$, and in terms of *h _{a}* by using the chain rule while holding

**constant, i.e.,**

*r*with $ f a \u2243 h a$ [see Eq. (31)], and by approximating

in the Rosenbluth potential *G* (and similarly in $ H ( f b )$). Using the transformation in Eq. (B7) to express the second-order velocity derivatives in Eq. (A1) yields FLR terms of the order of $ \u223c \nu a b k \u22a5 2 \rho a 2 h a$ in the test component when gyro-averaged, while the transformation in Eq. (B8) produces FLR terms proportional to Bessel functions *J _{n}* in the field components.

^{13,47}However, despite the difference in the analytical treatment of the velocity derivatives using the spherical harmonic expansion [see Eq. (B6)] and the lowest order gyrocenter coordinates [see Eq. (B7)], both approaches contain the same lowest order FLR terms, proportional to the spatial gradients of

*h*if proper approximations are applied to the spherical harmonic approach. In fact, Eq. (B7) can be recovered from Eq. (B6) by Taylor expanding the spatial dependence of particle spherical moments, such that

_{a}Using Eq. (B9) into Eq. (B6), the fact that $ \rho a$ is a velocity-dependent function and the spherical harmonic expansion of $ f a ( R , v )$ [i.e., Eq. (20) with $r$ replaced by ** R**], we derive that

In Eq. (B10), we also use the scalar invariance of the velocity–space derivative, i.e., $ \u2202 v | r f a ( R + \rho a , v ) \u2243 \u2202 v | R h a ( R , \mu , v \u2225 )$, and $ f a \u2243 h a$ in the second term. Equation (B10) shows that the velocity–space derivative evaluated using the spherical moments in Eq. (B6) agrees with the transformation in Eq. (B7) at the lowest order in $ k \u22a5 \rho a$. Using the first-order derivative, given in Eq. (B10), to express the second-order derivatives in the GK test component yields a FLR term of the order of $ \u223c \nu a b k \u22a5 2 \rho a 2 n a p j$ when gyro-averaged, similar to the transformation in Eq. (B7).

We finally remark that the last term in Eq. (B10), proportional to $ \rho a 2 k \u22a5 2$, depends on the perpendicular energy coordinate $ v \u22a5 2$. Hence, at large $ k \u22a5 2 v \u22a5 2$, a large number of Laguerre gyro-moments is required to capture accurately these FLR effects. In these cases, a careful truncation of the sum over *j* in Eq. (38) and in the Bessel function expansions [e.g., the sum over *n* in Eq. (64)] contained in the GK collision operators must be performed to obtain a well-behaved numerical results at larger $ v \u22a5 2$ and $ k \u22a5 2$, respectively. Numerical tests reveal that truncating these sums at $ j \u223c 12$ and $ n \u223c 8$ is sufficient in the cases discussed herein where $ k \u22a5 \rho i \u2272 1$.

### APPENDIX C: LOWEST ORDER GYRO-MOMENT ANALYTICAL EXPRESSIONS

At high collisionality, the number of gyro-moments necessary to describe the perturbed distribution function *g _{a}* is reduced since higher-order gyro-moments are strongly damped by collisions. In this regime, the perturbed distribution function is well approximated by a perturbed Maxwellian, with its perturbation that has a relative amplitude of the order of the ratio of the particle mean-free path to the typical parallel scale length, i.e., $ \lambda mfp / L \u2225 \u226a 1$. The Maxwellian gyro-moments (

*p*,

*j*), i.e., the gyrocenter density (0, 0), the parallel gyrocenter velocity (1, 0), the parallel and perpendicular temperatures, (2, 0) and (0, 1), respectively, are leading order in $ \lambda mfp / L \u2225$. On the other hand, the gyro-moments associated with the non-Maxwellian component of

*g*, i.e., the parallel and perpendicular heat fluxes (3, 0) and (1, 1), respectively, are first order in $ \lambda mfp / L \u2225$.

_{a}By projecting the GK Boltzmann equation on the Hermite–Laguerre basis^{21} [or the linearized GK Boltzmann equation given Eq. (76)], the gyro-moment expansion of the IS collision operator, presented in this work, allows us to evaluate explicitly the collisional terms that enter in the evolution equations of the lowest order gyro-moment enumerated above. We consider the DK IS collision operator (FLR effects yield complicated coefficients that rely on sums with the number of significant terms that depends on $ k \u22a5$). Using the closed analytical formulas of the DK IS, given in Eqs. (71) and (72), the non-vanishing terms $ \Delta C abp \u2032 j \u2032 Tpj$ in the gyro-moment expansion of the test component $ \Delta C a b Tpj$, are given by

where $ \tau = T a / T b , \u2009 \sigma = m a / m b$ are the temperature and mass ratios of the colliding species, respectively. Similarly, the non-vanishing terms $ \Delta C a b , p \u2032 j \u2032 Fpj$ in the gyro-moment expansion of the field component $ \Delta C a b Fpj$ [see Eq. (74)] are

We notice that, from the conservation of particle in Eq. (2a), it follows $ \Delta C a b , p \u2032 j \u2032 T 00 = \Delta C a b , p \u2032 j \u2032 F 00$. Finally, we give the non-vanishing terms in the case of like-species collision. Because the test component of the OS operator reduces to the Coulomb collision operator when *a* = *b*, $ \Delta C aap \u2032 j \u2032 Tpj = 0$, while the non-vanishing terms of the field component $ \Delta C aap \u2032 j \u2032 Fpj$ are given by

We also provide the non-vanishing lowest order terms in the gyro-moment expansion of the OS collision operator, which are given by

with $ C ablk STpj = C abpj STlk$, and

for the test and field components, respectively. In the case of like-species collision, the non-vanishing terms are

Finally and for completeness, we report the non-vanishing coefficients of the DK Coulomb collision operator of the test $ C a b L T$ and the field $ C a b L F$ components, respectively, i.e.,

and

respectively. Finally, the non-vanishing terms for like-species collisions are

The lowest order gyro-moments of the IS, OS, and Coulomb reported above can be used to obtain reduced-fluid models to study the plasma dynamics in the Pfirsch–Schlüter regime.