Numerical implementation of the improved Sugama collision operator using a moment approach

Frei, B. J.; Ernst, S.; Ricci, P.

doi:10.1063/5.0091244

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.⁷ 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 Λ ≳ 10$ ⁠) in fusion devices, where small-angle deflections dominate, allows for the use of the Fokker–Planck collision operator⁸ 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¹¹ 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.⁵ 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.¹⁶ 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²⁰ 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,⁶ which we refer to as the gyro-moment approach, of the perturbed distribution function,²¹ 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.²² 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²³ 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) / π^{3 / 2} v_{T a} {(r)}^{3 / 2} e^{- 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 \cdot v$ ⁠. The small-amplitude perturbation, $f_{a} = f_{a} (r, v)$ ⁠, i.e., $f_{a} / f_{M a} ≪ 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

C_{a b} (f_{a}, f_{b}) = C_{a b}^{T} (f_{a}) + C_{a b}^{F} (f_{b}),

(1)

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_b [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.^9,11,13 Among these simplified models, the OS collision operator model,¹³ 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_ab, i.e.,

\int d v C_{a b} (f_{a}, f_{b}) = 0,

(2a)

m_{a} \int d v v C_{a b}^{T} (f_{a}) = - m_{b} \int d v v C_{b a}^{F} (f_{a}),

(2b)

m_{a} \int d v v^{2} C_{a b}^{T} (f_{a}) = - m_{b} \int d v v^{2} C_{b a}^{F} (f_{a})

(2c)

and satisfy the H-theorem and the adjointess relations even in the case of collisions between particles with temperatures $T_{a} \neq T_{b}$ ⁠, given by

T_{a} \int d v \frac{f_{a}}{f_{M a}} C_{a b} (f_{a}, f_{b}) + T_{b} \int d v \frac{f_{b}}{f_{M b}} C_{b a} (f_{b}, f_{a}) \leq 0

(3)

and

(4)

\int d v \frac{f_{a}}{f_{M a}} C_{a b}^{T} (g_{a}) = \int d v \frac{g_{a}}{f_{M a}} C_{a b}^{T} (f_{a}),

(4a)

T_{a} \int d v \frac{f_{a}}{f_{M a}} C_{a b}^{F} (g_{b}) = T_{b} \int d v \frac{g_{b}}{f_{M b}} C_{b a}^{F} (f_{a}),

(4b)

respectively, where f_a and g_a are two arbitrary phase-space

(r, v)

functions. The complete analytical derivation of

C_{a b}^{S}

is given in Ref. 13.

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.¹ 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,¹⁸ given by

F_{a i} = {(- 1)}^{i - 1} \int d v m_{a} v L_{i - 1}^{3 / 2} (s_{a}^{2}) \sum_{b} C_{a b} (f_{a}, f_{b}),

(5)

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, \dots$ ⁠) is the associated Laguerre polynomial defined by

L_{i}^{α} (x) = \sum_{l = 0}^{i} L_{i l}^{α - 1 / 2} x^{l},

(6)

where $L_{i l}^{α - 1 / 2} = {(- 1)}^{l} (α + i)! / [(i - l)! (l + α)! l!]$ (with $α > - 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 $Δ C_{a b} (f_{a}, f_{b})$ to $C_{a b}^{S}$ (see Appendix A), i.e.,

C_{a b}^{I S} (f_{a}, f_{b}) = C_{a b}^{S} (f_{a}, f_{b}) + Δ C_{a b} (f_{a}, f_{b}),

(7)

where $Δ C_{a b}$ is defined as

Δ C_{a b} = Δ C_{a b}^{T} + Δ C_{a b}^{F},

(8)

with the test and field components of the correction term being

Δ C_{a b}^{T} = \sum_{ℓ = 0}^{L} \sum_{k = 0}^{K} \frac{m_{a}}{T_{a}} \frac{f_{M a}}{{\bar{τ}}_{a b}} c_{ℓ} Δ M_{a b}^{ℓ k} L_{ℓ}^{3 / 2} (s_{a}^{2}) v \cdot u_{a k} (f_{a})

(9)

and

Δ C_{a b}^{F} = \sum_{ℓ = 0}^{L} \sum_{k = 0}^{K} \frac{m_{a}}{T_{a}} \frac{f_{M a}}{{\bar{τ}}_{a b}} c_{ℓ} Δ N_{a b}^{ℓ k} L_{ℓ}^{3 / 2} (s_{a}^{2}) v \cdot u_{b k} (f_{b}),

(10)

respectively. Here, ${\bar{τ}}_{a b} = 3 \sqrt{π} / (4 ν_{a b})$ is the collisional time between the colliding species a and b, with ν_ab the associated collision frequency $ν_{a b} = 4 \sqrt{π} n_{b} q_{a}^{2} q_{b}^{2} ln Λ / [3 m_{a}^{1 / 2} T_{a}^{3 / 2}]$ (see Ref. 6). While the IS is derived in the limit $L \to \infty$ and $K \to \infty$ ⁠, 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 ≳ 2$ (see Ref. 19), we consider the cases of $L = K = 2$ ⁠, 5 and 10. Despite that neoclassical studies have revealed that the energy (⁠ $\sim v^{2}$ ⁠) expansion in Laguerre polynomials has slow convergence in the banana regime,²⁴ 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¹⁹

u_{a k} (f_{a}) = \frac{c_{k}}{n_{a}} \int d v f_{a} L_{k}^{3 / 2} (s_{a}^{2}) v,

(11)

with $c_{k} = 3 \cdot 2^{k} k! / (2 k + 3)!!$ ⁠. Finally, we introduce the correction Braginskii matrix elements, $Δ M_{a b}^{ℓ k}$ and $Δ N_{a b}^{ℓ k}$ ⁠, defined by

(12)

Δ M_{a b}^{ℓ k} = M_{a b}^{L ℓ k} - M_{a b}^{S ℓ k},

(12a)

Δ N_{a b}^{ℓ k} = N_{a b}^{L ℓ k} - N_{a b}^{S ℓ k},

(12b)

where the Braginskii matrices,

M_{a b}^{A ℓ k}

and

N_{a b}^{A ℓ k}

(being 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

(13)

\frac{n_{a}}{{\bar{τ}}_{a b}} M_{a b}^{A ℓ k} = \int d v v_{∥} L_{ℓ}^{3 / 2} (s_{a}^{2}) C_{a b}^{A T} (f_{M a} \frac{m_{a} v_{∥}}{T_{a}} L_{k}^{3 / 2} (s_{a}^{2})),

(13a)

\frac{n_{a}}{{\bar{τ}}_{a b}} N_{a b}^{A ℓ k} = \int d v v_{∥} L_{ℓ}^{3 / 2} (s_{a}^{2}) C_{a b}^{A F} (f_{M b} \frac{m_{b} v_{∥}}{T_{b}} L_{k}^{3 / 2} (s_{b}^{2})),

(13b)

with

v_{∥} = b \cdot v

the parallel component of the velocity along the magnetic field and

b = B / B

⁠.

The Braginskii matrices $M_{a b}^{A ℓ k}$ and $N_{a b}^{A ℓ 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

M_{a b}^{A 0 k} + \frac{T_{a} v_{T a}}{T_{b} v_{T b}} N_{b a}^{A 0 k} = 0 (k = 0, 1, 2, \dots) .

(14)

In addition, in the case of collisions between particle species with the same temperatures (T_a = T_b), 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

M_{a b}^{A ℓ k} = M_{a b}^{A k ℓ},

(15)

\frac{N_{a b}^{A ℓ k}}{T_{a} v_{T a}} = \frac{N_{b a}^{A k ℓ}}{T_{b} v_{T b}} (ℓ, k = 0, 1, 2, \dots),

(16)

which implies from Eq. (12) that

Δ M_{a b}^{ℓ k} = 0,

(17)

Δ N_{a b}^{ℓ k} = \frac{N_{a b}^{A 00} N_{a b}^{A ℓ k} - N_{a b}^{A ℓ 0} N_{a b}^{A 0 k}}{N_{a b}^{A 00}} (ℓ, k = 0, 1, 2, \dots),

(18)

Δ N_{a b}^{00} = Δ N_{a b}^{ℓ 0} = Δ N_{a b}^{0 k} = 0

(19)

for $ℓ, k = 1, 2, \dots$ in the case T_a = T_b. While the analytical expressions of the Braginskii matrices $M_{a}^{A ℓ k}$ and $N_{a}^{A ℓ k}$ up to $(ℓ, k) \leq 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 $(ℓ, k)$ in Sec. IV. This allows us to evaluate gyro-moment expansion of the IS for any (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

f_{a} (r, v) = f_{M a} \sum_{p, j} \frac{1}{σ_{j}^{p}} M_{a}^{p j} (r) \cdot Y^{p j} (s_{a})

(20)

with $s_{a} = v / v_{T a}$ and $σ_{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} \cdot e^{p m'} = {(- 1)}^{m} δ_{- m}^{m'}$ ⁠,²⁷ such that

Y^{p} (s_{a}) = s_{a}^{p} \sqrt{\frac{2 π^{3 / 2} p!}{2^{p} (p + 1 / 2)!}} \sum_{m = - p}^{p} Y_{p}^{m} (ξ, θ) e^{p m},

(21)

where $Y_{p}^{m} (ξ, θ)$ 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 - 1} = (e_{x} - i e_{y}) / \sqrt{2}, e^{10} = b$ and $e^{11} = - (e_{x} + i e_{y}) / \sqrt{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²⁷

e^{p m} = N_{p m} \sum_{n = 0}^{[(p - m) / 2]} a_{n}^{p m} {{(e^{11})}^{m + n} {(e^{1 - 1})}^{n} {(e^{10})}^{p - m - 2 n}}_{S},

(22)

where $N_{p m} = \sqrt{(p + m)! (p - m)! 2^{p - m} / (2 p)!}, a_{n}^{p m} = p! / [2^{n} n! (m + n)! (p - m - 2 n)!]$ ⁠, and ${\dots}_{S}$ denotes the symmetric part. A complete introduction to the spherical harmonic tensors can be found in Ref. 27.

The spherical harmonic basis in Eq. (20) satisfies the orthogonality relation²⁷

\begin{matrix} \frac{1}{π^{3 / 2} σ_{j}^{p}} \int d v e^{- v^{2}} L_{j}^{p + 1 / 2} (v^{2}) Y^{p'} (v) L_{j'}^{p' + 1 / 2} (v^{2}) Y^{p} (v) \cdot T^{p} = δ_{p p'} δ_{j j'} T^{p} \end{matrix}

(23)

with $T^{p}$ an arbitrary pth 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:

M_{a}^{p j} (r) = \frac{1}{n_{a}} \int d v f_{a} (r, v) Y^{p j} (s_{a}) .

(24)

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

\begin{matrix} u_{a k} (f_{a}) = \frac{c_{k} v_{T a}}{n_{a}} \int d v f_{a} L_{k}^{3 / 2} (s_{a}^{2}) Y^{1} (s_{a}) = c_{k} v_{T a} M_{a}^{1 k} (r) . \end{matrix}

(25)

Inserting Eq. (25) into Eqs. (9) and (10) yields the spherical harmonic expansion of the IS collision operator, which is useful to evaluate its GK limit.

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, μ, v_{∥}, θ, t)$ where R is the gyrocenter position, $μ = m v_{⊥}^{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.,

C_{a b}^{I S} = {⟨ C_{a b}^{I S} ⟩}_{R} = \int_{0}^{2 π} |_{R} \frac{d θ}{2 π} C_{a b}^{I S} (z (Z)),

(26)

where the integral over the gyroangle appearing in Eq. (26) is performed holding R constant, while collisions occur at the particle position r. 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, Z and z, respectively, can be written as $Z = z + δ z$ ⁠, where $δ z$ are functions of phase-space coordinates and perturbed fields and contain terms at all orders in the GK expansion parameter $ε \sim e ϕ / T_{e}$ (⁠ $ϕ$ being the small amplitude and small-scale electrostatic fluctuating potential^6,28) At the lowest order in the GK expansion, the coordinate transformation reduces to $δ v_{∥} = δ μ = δ θ = 0$ and $δ r ≃ - ρ_{a}$ ⁠. Hence, the IS collision operator can be gyro-averaged holding R constant in gyrocenter coordinates, such that $r = R (r, v) + ρ_{a} (r, v)$ ⁠.²¹

Focusing first on the test component of $Δ C_{a b}^{T}$ ⁠, we perform the gyro-average in Eq. (26) in gyrocenter coordinates using $r = R (r, v) + ρ_{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) = \int d k M_{a}^{1 k} (k) e^{i k \cdot R} e^{i k \cdot ρ_{a}}$ ⁠. For a single Fourier component, we derive

\begin{matrix} Δ C_{a b}^{T} = \sum_{ℓ = 0}^{L} \sum_{k = 0}^{K} \frac{m_{a}}{T_{a}} \frac{f_{M a}}{{\bar{τ}}_{a b}} c_{ℓ} c_{k} v_{T a} Δ M_{a b}^{ℓ k} L_{ℓ}^{3 / 2} (s_{a}^{2}) e^{i k \cdot R} \\ \times {⟨ e^{i k \cdot ρ_{a}} v ⟩}_{R} \cdot M_{a}^{1 k} (k) \\ = \sum_{ℓ = 0}^{L} \sum_{k = 0}^{K} \frac{m_{a}}{T_{a}} \frac{f_{M a}}{{\bar{τ}}_{a b}} c_{ℓ} c_{k} v_{T a} Δ M_{a b}^{ℓ k} L_{ℓ}^{3 / 2} (s_{a}^{2}) e^{i k \cdot R} \\ \times [v_{∥} J_{0 a} b \cdot M_{a}^{1 k} (k) + i v_{⊥} J_{1 a} e_{2} \cdot M_{a}^{1 k} (k)], \end{matrix}

(27)

where we introduce the basis vectors $(e_{1}, e_{2})$ ⁠, such that $v_{⊥} = v_{⊥} (e_{2} cos θ - e_{1} sin θ)$ ⁠. Additionally, $J_{0 a} = J_{0} (b_{a} \sqrt{x_{a}})$ and $J_{1 a} = J_{1} (b_{a} \sqrt{x_{a}})$ are the zeroth and first order Bessel functions (with $b_{a} = k_{⊥} ρ_{a}$ the normalized perpendicular wavenumber and $x_{a} = v_{⊥}^{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_a, and ${\bar{τ}}_{a b}$ ⁠) are assumed to vary on spatial scale lengths much larger that $ρ_{a}$ ⁠. Therefore, they are evaluated at R and are not affected by the gyro-average operator, in contrast to $M_{a}^{1 k} (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_a. The two gyrocenter distribution functions, h_a and g_a, are related by⁶

h_{a} (R, μ, v_{∥}) = g_{a} (R, μ, v_{∥}) + \frac{q_{a}}{T_{a}} F_{M a} (R, μ, v_{∥}) {⟨ ϕ ⟩}_{R}

(28)

in the electrostatic limit. Here, $F_{M a} = N_{a} (R) / π^{3 / 2} v_{T a}^{3} (R) e^{- v_{∥}^{2} / v_{T a}^{2} (R) - μ 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_a by the scalar invariance of the full particle and gyrocenter distribution functions, i.e.,

f_{a} (r, v) = g_{a} (R, μ, v_{∥}) + F_{M a} (R, μ, v_{∥}) - f_{M a} (r (Z), v (Z)) .

(29)

Using the pull-back operator $T$ ⁠, such that the functional forms of f_Ma and F_Ma are related by $f_{M a} = T F_{M a}$ ⁠,²¹ we derive that

f_{a} (r, v) = g_{a} (R, μ, v_{∥}) + (F_{M a} - T F_{M a}) (Z) = g_{a}^{g c} (Z) .

(30)

We remark that, while both g_a and h_a are gyrophase independent functions, $g_{a}^{g c}$ is gyrophase dependent via its arguments Z. An expression of the pull-back transformation $T$ can be obtained at the leading order in the GK expansion parameter $ε \sim e ϕ / T_{e}$ ⁠, yielding⁶

\begin{matrix} f_{a} (r, v) = g_{a}^{g c} (Z (z)) \\ = g_{a} (R (r, v), μ, v_{∥}) - \frac{q_{a}}{T_{a}} F_{M a} (ϕ (r) - {⟨ ϕ ⟩}_{R}) + O (ε^{2}) \\ = h_{a} (R (r, v), μ, v_{∥}) - \frac{q_{a}}{T_{a}} ϕ (r) F_{M a} + O (ε^{2}) \end{matrix}

(31)

being $r = R + ρ_{a} (μ, θ)$ ⁠. Using Eq. (31) allows us to finally express $M_{a}^{1 k}$ in terms of h_a,

\begin{matrix} v_{T a} M_{a}^{1 k} (r) = \frac{v_{T a}}{n_{a}} \int d v f_{a} (r, v) L_{k}^{3 / 2} (s_{a}^{2}) Y^{1} (s_{a}) \\ = \frac{1}{n_{a}} \int d v \int d r' δ (r' - r) f_{a} (r', v) L_{k}^{3 / 2} (s_{a}^{2}) v \\ = \frac{1}{N_{a}} \int d R d v_{∥} d μ d θ \frac{B}{m_{a}} δ (R + ρ_{a} - r) \\ \times h_{a} (R (r, v), μ, v_{∥}) L_{k}^{3 / 2} (s_{a}^{2}) v \\ = \frac{e^{i k \cdot r}}{N_{a}} \int d v_{∥} d μ d θ \frac{B}{m_{a}} L_{k}^{3 / 2} (s_{a}^{2}) \\ \times h_{a} (k, μ, v_{∥}) (b v_{∥} J_{0 a} - i v_{⊥} e_{2} J_{1 a}), \end{matrix}

(32)

where $N_{a} = \int d v_{∥} d μ d θ 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_{∥}^{*} / m_{a} = B (1 + v_{∥} b \cdot \nabla \times b / Ω_{a}) / m_{a} ≃ B / m_{a}$ ⁠, since the last term in $B_{∥}^{*}$ is of the order $v_{∥} b \cdot \nabla \times b / Ω_{a} \sim ρ_{a} / L_{B} ≪ 1$ (with L_B the typical equilibrium scale length of the magnetic field B). We remark that the contribution from the terms proportional to $ϕ$ ⁠, appearing in Eq. (31), is neglected in Eq. (32). In fact, while these terms are of the same order as h_a (i.e., they are order ε), they yield a small contribution in the collision operator.¹³ 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, $Δ C_{a b}^{T} = {⟨ Δ C_{a b}^{T} ⟩}_{R}$ ⁠, can be obtained in terms of the gyrocenter distribution function h_a. Focusing on the field component of $Δ C_{a b}$ ⁠, i.e., on $Δ C_{a b}^{F}$ ⁠, we remark that a similar derivation of its expression can be carried out as for $Δ 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 with b.

Finally, the GK IS collision operator can be expressed as follows:¹⁷

C_{a b}^{I S} = C_{a b}^{S} + Δ C_{a b}^{T} + Δ C_{a b}^{F},

(33)

where $C_{a b}^{S}$ is the OS GK collision operator, given in Ref. 13, and $Δ C_{a b}^{T}$ and $Δ C_{a b}^{F}$ are given by

(34)

\begin{matrix} Δ C_{a b}^{T} = \sum_{ℓ = 0}^{L} \sum_{k = 0}^{K} \frac{c_{ℓ}}{{\bar{τ}}_{a b}} \frac{m_{a} F_{M a}}{T_{a}} L_{ℓ}^{3 / 2} (s_{a}^{2}) Δ M_{a b}^{ℓ k} \\ \times ({\bar{u}}_{∥ a}^{k} (h_{a}) J_{0 a} v_{∥} + {\bar{u}}_{⊥ a}^{k} (h_{a}) J_{1 a} v_{⊥}), \end{matrix}

(34a)

\begin{matrix} Δ C_{a b}^{F} = \sum_{ℓ = 0}^{L} \sum_{k = 0}^{K} \frac{c_{ℓ}}{{\bar{τ}}_{a b}} \frac{m_{a} F_{M a}}{T_{a}} L_{ℓ}^{3 / 2} (s_{a}^{2}) Δ N_{a b}^{ℓ k} \\ \times ({\bar{u}}_{∥ b}^{k} (h_{b}) J_{0 a} v_{∥} + {\bar{u}}_{⊥ b}^{k} (h_{b}) J_{1 a} v_{⊥}), \end{matrix}

(34b)

where we introduce the quantities

(35)

{\bar{u}}_{∥ s}^{k} (h_{s}) = \frac{c_{k}}{n_{s}} \int d v_{∥} d μ d θ \frac{B}{m_{s}} L_{k}^{3 / 2} (s_{s}^{2}) h_{s} J_{0 s} v_{∥},

(35a)

{\bar{u}}_{⊥ s}^{k} (h_{s}) = \frac{c_{k}}{n_{s}} \int d v_{∥} d μ d θ \frac{B}{m_{s}} L_{k}^{3 / 2} (s_{s}^{2}) h_{s} J_{1 s} v_{⊥} .

(35b)

The GK IS collision operator can be obtained by adding to the GK OS collision operator,

C_{a b}^{S}

⁠, the terms in Eqs. and . In Appendix B, we discuss the GK formulation of the Coulomb and OS operators that we use to perform the numerical tests in Sec. VI. In particular, we note that the GK operators considered in this work are derived from the linearized collision operators applied to the perturbed particle distribution function f_a, following the derivation of the full-F nonlinear GK Coulomb collision operator in Ref. 26 (see Appendix B).

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 ≃ R$ ⁠. Hence, the DK IS collision operator is derived in the zero gyroradius limit of the GK IS operator approximating $J_{0 s} ≃ 1, J_{1 s} = 0$ and $f_{s} ≃ 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 $Δ C_{a b}$ are

(36)

Δ C_{a b}^{T} = \sum_{ℓ = 0}^{L} \sum_{k = 0}^{K} \frac{m_{a}}{T_{a}} \frac{F_{M a}}{{\bar{τ}}_{a b}} c_{ℓ} c_{k} v_{T a} L_{ℓ}^{3 / 2} (s_{a}^{2}) Δ M_{a b}^{ℓ k} v_{∥} b \cdot M_{a}^{1 k},

(36a)

Δ C_{a b}^{F} = \sum_{ℓ = 0}^{L} \sum_{k = 0}^{K} \frac{m_{a}}{T_{a}} \frac{F_{M a}}{{\bar{τ}}_{a b}} c_{ℓ} c_{k} v_{T b} L_{ℓ}^{3 / 2} (s_{a}^{2}) Δ N_{a b}^{ℓ k} v_{∥} b \cdot M_{b}^{1 k},

(36b)

respectively. In Eq. , the spherical particle moments,

M_{s}^{1 k}

(⁠

s = a, b

⁠), are expressed in terms of g_s, such that

b \cdot M_{s}^{1 k} = \frac{1}{N_{s}} \int d μ d v_{∥} d θ \frac{B}{m_{s}} g_{s} L_{k}^{3 / 2} (s_{s}^{2}) \frac{v_{∥}}{v_{T s}} .

(37)

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 $Δ M_{a b}^{ℓ k}$ and $Δ N_{a b}^{ℓ k}$ ⁠, in Sec. IV.

IV. BRAGINKSII MATRICES

In order to evaluate the correction term $Δ C_{a b}$ given in Eq. (8), analytical expressions for the Braginskii matrices, $M_{a b}^{A ℓ k}$ and $N_{a b}^{A ℓ k}$ ⁠, associated with the Coulomb and OS collision operators are derived. This extends the evaluation of the $(ℓ, k) \leq 2$ Braginksii matrices reported in Ref. 17 to arbitrary $(ℓ, 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 ℓ k}$ and $N_{a b}^{L ℓ 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

(38)

C_{a b}^{L T} (f_{a}) = C_{a b}^{L T} (r, v) = \sum_{p = 0}^{\infty} \sum_{j = 0}^{\infty} M_{a}^{p j} (r) \cdot C_{a b}^{L T} (\frac{f_{M a}}{σ_{j}^{p}} Y^{p j} (s_{a})),

(38a)

C_{a b}^{L F} (f_{b}) = C_{a b}^{L F} (r, v) = \sum_{p = 0}^{\infty} \sum_{j = 0}^{\infty} M_{b}^{p j} (r) \cdot C_{a b}^{L F} (\frac{f_{M b}}{σ_{j}^{p}} Y^{p j} (s_{b}))

(38b)

with

C_{a b}^{L T} (\frac{f_{M a}}{σ_{j}^{p}} Y^{p j} (s_{a})) = \frac{f_{M a}}{σ_{j}^{p}} Y^{p} (\hat{v}) ν_{a b}^{Tpj} (v),

(39a)

C_{a b}^{L F} (\frac{f_{M b}}{σ_{j}^{p}} Y^{p j} (s_{b})) = \frac{f_{M a}}{σ_{j}^{p}} Y^{p} (\hat{v}) ν_{a b}^{Fpj} (v),

(39b)

where the expressions of the test and field speed functions, $ν_{a b}^{Tpj} (v)$ and $ν_{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 ℓ 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_{∥} L_{k}^{3 / 2} (s_{a}^{2}) / T_{a}$ ⁠, i.e.,

C_{a b}^{L T} (f_{M a} \frac{m_{a} v_{∥}}{T_{a}} L_{k}^{3 / 2} (s_{a}^{2})) = \frac{2 f_{M a}}{v_{T a}} \frac{v_{∥}}{v} ν_{a b}^{T 1 k} (v) .

(40)

Then, the matrix element $M_{a b}^{L ℓ 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 $ν_{a b}^{T 1 k} (v)$ ⁠. It yields

M_{a b}^{L ℓ k} = \sum_{l = 0}^{ℓ} \frac{2}{3} \frac{{\bar{τ}}_{a b}}{n_{a}} L_{ℓ l}^{1} {\bar{ν}}_{* a b}^{T 1 k l} .

(41)

A similar derivation can be carried out to evaluate $N_{a b}^{L ℓ k}$ [see Eq. (13b)] by using the expression in Eq. (39b), i.e.,

N_{a b}^{L ℓ k} = \sum_{l = 0}^{ℓ} \frac{2}{3} \frac{{\bar{τ}}_{a b}}{n_{a}} χ_{a b} L_{ℓ l}^{1} {\bar{ν}}_{* a b}^{F 1 k l}

(42)

with $χ_{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 ${\bar{ν}}_{* a b}^{T 1 k l}$ and ${\bar{ν}}_{* 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, $Δ M_{a b}^{ℓ k}$ and $Δ N_{a b}^{ℓ 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 ℓ k}$ and $N_{a b}^{S ℓ k}$ ⁠, appearing in Eq. (12). For $f_{a} = f_{M a} m_{a} v_{∥} 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

C_{a b}^{S T} (f_{M a} L_{k}^{3 / 2} (s_{a}^{2}) \frac{m_{a} v_{∥}}{T_{a}}) = \frac{2}{v_{T a}} f_{M a} \frac{v_{∥}}{v} ν_{a b}^{S 1 k} (v) + \sum_{i = 1}^{3} X_{a b}^{i},

(43)

where the quantities $X_{a b}^{i}$ are defined by

X_{a b}^{1} = - \frac{16}{3 \sqrt{π}} (1 + χ_{a b}^{2}) (θ_{a b} - 1) f_{M a} \frac{m_{a} v_{∥}}{T_{a}} \sum_{l = 0}^{k} L_{k l}^{1} {\bar{ν}}_{a b}^{∥ l + 3},

(44a)

X_{a b}^{2} = - 2 (1 + χ_{a b}^{2}) (θ_{a b} - 1) \frac{m_{a}}{T_{a}} u_{∥ a}^{k} v_{∥} f_{M a} ν_{a b}^{∥} (v) s_{a}^{2},

(44b)

X_{a b}^{3} = - \frac{2}{{\bar{τ}}_{a b}} f_{M a} \frac{χ_{a b} {(θ_{a b} - 1)}^{2}}{\sqrt{1 + χ_{a b}^{2}}} \frac{m_{a}}{T_{a}} v_{∥} u_{∥ a}^{k},

(44c)

with $ν_{a b}^{S 1 k} (v)$ the velocity-dependent speed function (whose expression is given in Ref. 6), being ${\bar{ν}}_{a b}^{∥ k} = \int_{0}^{\infty} d s_{a} s_{a}^{2 k} ν_{a b}^{∥} (v) e^{- s_{a}^{2}}$ (with $ν_{a b}^{∥} (v) = 2 ν_{a b} [erf (s_{b}) - s_{b} erf' (s_{b})] / (2 s_{b}^{2} s_{a}^{3})$ the velocity dependent energy diffusion frequency, and $u_{∥ a}^{k} = 4 / (3 \sqrt{π}) (k + 3 / 2)! / k! δ_{k}^{0}$ ⁠. In deriving Eq. (43), we remark that the terms proportional to $T_{a} \int d v f_{a} (2 s_{a}^{2} / 3 - 1) / n_{a}$ vanish exactly when applied to $f_{a} = f_{M a} m_{a} v_{∥} L_{k}^{3 / 2} (s_{a}^{2}) / T_{a}$ because of the velocity integration over the pitch-angle variable $v_{∥} / v$ ⁠.

Similarly, the field component of the OS collision for $f_{b} = f_{M b} m_{b} v_{∥} L_{k}^{3 / 2} (s_{b}^{2}) / T_{b}$ yields

\begin{matrix} C_{a b}^{S F} (f_{M b} L_{k}^{3 / 2} (s_{b}^{2}) \frac{m_{b} v_{∥}}{T_{a}}) = \frac{2 θ_{a b}}{{\bar{τ}}_{a b}} (1 + χ_{a b}^{2}) f_{M a} \frac{m_{a} v_{∥}}{T_{a}} [\frac{3 \sqrt{π}}{2} \frac{Φ (s_{b})}{s_{a}} + \frac{χ_{a b} (θ_{a b} - 1)}{{(1 + χ_{a b}^{2})}^{3 / 2}}] V_{a b} (f_{M b} L_{k}^{3 / 2} (s_{b}^{2}) \frac{m_{b} v_{∥}}{T_{b}}), \end{matrix}

(45)

where

\begin{matrix} V_{a b} (f_{M b} L_{k}^{3 / 2} (s_{b}^{2}) \frac{m_{b} v_{∥}}{T_{b}}) = - \frac{θ_{b a}}{{\bar{τ}}_{b a}} (1 + χ_{b a}^{2}) \frac{m_{b}}{γ_{a b}} \int d v L_{k}^{3 / 2} (s_{b}^{2}) \frac{m_{b} v_{∥}^{2}}{T_{b}} f_{M b} [\frac{3 \sqrt{π}}{2} \frac{Φ (s_{a})}{s_{b}} + \frac{χ_{b a} (θ_{b a} - 1)}{{(1 + χ_{b a}^{2})}^{3 / 2}}], \end{matrix}

(46)

with $Φ (x) = [erf (x) - x erf' (x)] / (2 x^{2})$ ⁠. The velocity integral in Eq. (46) can be performed analytically using Eq. (6) and leads to

V_{a b} (f_{M b} L_{k}^{3 / 2} (s_{b}^{2}) \frac{m_{b} v_{∥}}{T_{b}}) = \frac{m_{b}}{γ_{a b}} V_{a b}^{l},

(47)

with the following definition:

\begin{matrix} V_{a b}^{l} = - 2 θ_{a b} (1 + χ_{a b}^{2}) \frac{n_{b}}{{\bar{τ}}_{a b}} [\frac{1}{2} \frac{χ_{a b} (θ_{a b} - 1)}{{(1 + χ_{a b}^{2})}^{3 / 2}} u_{∥ a}^{l} \\ + \sum_{m = 0}^{l} L_{l m}^{1} (\frac{1}{χ_{a b}^{2}} E_{a b}^{m} - \frac{1}{χ_{a b}} e_{a b}^{m + 1})], \end{matrix}

(48)

where we introduce $e_{b a}^{k} = \int_{0}^{\infty} d s_{b} s_{b}^{2 k} erf' (s_{a}) e^{- s_{b}^{2}}$ and $E_{b a}^{k} = \int_{0}^{\infty} d s_{b} s_{b}^{2 k + 1} erf (s_{a}) e^{- 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 ℓ k}$ and $N_{a b}^{S ℓ 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

M_{a b}^{S ℓ k} = \sum_{n = 1}^{3} M_{abn}^{S ℓ k},

(49)

where we introduce the quantities

M_{a b 1}^{S ℓ k} = \frac{{\bar{τ}}_{a b}}{n_{a}} \sum_{l = 0}^{ℓ} \frac{2}{3} L_{ℓ l}^{1} {\bar{ν}}_{* a b}^{S 1 k l},

(50a)

\begin{matrix} M_{a b 2}^{S ℓ k} = - \frac{16}{3 \sqrt{π}} {\bar{τ}}_{a b} (θ_{a b} - 1) (1 + χ_{a b}^{2}) \\ \times [u_{∥ a}^{k} \sum_{l = 0}^{ℓ} L_{ℓ l}^{1} {\bar{ν}}_{a b}^{∥ l + 3} + u_{∥ a}^{ℓ} \sum_{l = 0}^{k} L_{k l}^{1} {\bar{ν}}_{a b}^{∥ l + 3}], \end{matrix}

(50b)

M_{a b 3}^{S ℓ k} = - \frac{2 χ_{a b} {(θ_{a b} - 1)}^{2}}{\sqrt{1 + χ_{a b}^{2}}} u_{∥ a}^{k} u_{∥ a}^{ℓ} .

(50c)

In Eq. (50a), the analytical expression of the speed integrated function, ${\bar{ν}}_{* a b}^{S 1 k l}$ ⁠, is reported in Ref. 6. Similarly for $N_{a b}^{S ℓ j}$ ⁠, using Eq. (45) and employing the expansion of the associated Laguerre polynomials in Eq. (6) yields

N_{a b}^{S ℓ k} = - \frac{2 {\bar{τ}}_{a b}}{n_{a}} \frac{m_{b}}{γ_{a b}} V_{b a}^{k} V_{a b}^{ℓ} .

(51)

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 $Δ M_{a b}^{ℓ k}$ and $Δ N_{a b}^{ℓ k}$ for arbitrary $(ℓ, 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.²¹ 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.¹⁶

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 $Δ C_{a b}^{T}$ and $Δ 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_a [see Eq. (28)]. Therefore, we expand the collision operator in terms of gyro-moments of h_a. More precisely, h_a is written as

h_{a} = \sum_{p = 0}^{\infty} \sum_{j = 0}^{\infty} n_{a}^{p j} \frac{H_{p} (s_{∥ a}) L_{j} (x_{a})}{\sqrt{2^{p} p!}} F_{M a}

(52)

with $s_{∥ a} = v_{∥} / v_{T a}$ ⁠. In Eq. (52), the Hermite and Laguerre polynomials, H_p and L_j, are defined via their Rodrigues's formulas $H_{p} (x) = {(- 1)}^{p} e^{x^{2}} d^{p} (e^{- x^{2}}) / d x^{p}$ and $L_{j} (x) = e^{x} / j! d^{j} (e^{- x} x^{j}) / d x^{j}$ ⁠, and are orthogonal over the intervals, $[- \infty, \infty]$ weighted by $e^{- x^{2}}$ ⁠, and $[0, + \infty]$ weighted by $e^{- x}$ ⁠, respectively, such that

\int_{- \infty}^{\infty} d x H_{p} (x) H_{p'} (x) e^{- x^{2}} = 2^{p} p! \sqrt{π} δ_{p}^{p'},

(53)

\int_{0}^{\infty} d x L_{j} (x) L_{j'} (x) e^{- x} = δ_{j}^{j'} .

(54)

Because of the orthogonality relations in Eq. (53), the non-adiabatic gyro-moments of h_a, $n_{a}^{p j}$ ⁠, are defined by

n_{a}^{p j} = \frac{1}{N_{a}} \int d μ d v_{∥} d θ \frac{B}{m_{a}} h_{a} \frac{H_{p} (s_{∥ a}) L_{j} (x_{a})}{\sqrt{2^{p} p!}},

(55)

where N_a is the gyrocenter density and F_Ma 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_{∥}^{*} / m_{a}$ ⁠,^21,29 such that

n_{a}^{* p j} = n_{a}^{p j} + \frac{v_{T a}}{\sqrt{2}} \frac{b \cdot \nabla \times b}{Ω_{a}} (\sqrt{p + 1} n_{a}^{p + 1 j} + \sqrt{p} n_{a}^{p - 1 j}) .

(56)

We detail the projection of the GK IS collision operator in Sec. V A, and obtain the DK limit of the same operator in Sec. V B in terms of $n_{a}^{p j}$ ⁠.

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

C_{a b}^{ISpj} = C_{a b}^{Spj} + Δ C_{a b}^{Tpj} + Δ C_{a b}^{Fpj},

(57)

where the Hermite–Laguerre projection of $C_{a b}^{I S}$ is defined by

C_{a b}^{ISpj} = \frac{1}{N_{a}} \int d μ d v_{∥} d θ \frac{B}{m_{a}} C_{a b}^{I S} \frac{H_{p} (s_{∥ a}) L_{j} (x_{a})}{\sqrt{2^{p} p!}},

(58)

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

Δ C_{a b}^{Tpj} = \frac{1}{N_{a}} \int d μ d v_{∥} d θ \frac{B}{m_{a}} Δ C_{a b}^{T} \frac{H_{p} (s_{∥ a}) L_{j} (x_{a})}{\sqrt{2^{p} p!}}

(59)

and

Δ C_{a b}^{Fpj} = \frac{1}{N_{a}} \int d μ d v_{∥} d θ \frac{B}{m_{a}} Δ C_{a b}^{F} \frac{H_{p} (s_{∥ a}) L_{j} (x_{a})}{\sqrt{2^{p} p!}} .

(60)

Since the expression of $C_{a b}^{Spj}$ is reported in Ref. 6, we focus here on the gyro-moment expansion of $Δ C_{a b}^{T}$ and $Δ 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 $Δ C_{a b}^{Tpj}$ ⁠. As an initial step, we express ${\bar{u}}_{∥ s}^{k}$ and ${\bar{u}}_{⊥ 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 ${\bar{u}}_{∥ s}^{k}$ and ${\bar{u}}_{⊥ s}^{k}$ yields

(61)

{\bar{u}}_{∥ s}^{k} = v_{T s} c_{k} \sum_{p = 0}^{\infty} \sum_{j = 0}^{\infty} n_{a}^{p j} I_{∥ s}^{pjk},

(61a)

{\bar{u}}_{⊥ s}^{k} = v_{T s} c_{k} \sum_{p = 0}^{\infty} \sum_{j = 0}^{\infty} n_{a}^{p j} I_{⊥ s}^{pjk},

(61b)

where we introduce the velocity integrals

I_{∥ s}^{pjk} = \frac{1}{N_{a}} \int d μ d v_{∥} d θ \frac{B}{m_{a}} \frac{H_{p} (s_{∥ s}) L_{j} (x_{s})}{\sqrt{2^{p} p!}} F_{M s} L_{k}^{3 / 2} (s_{s}^{2}) J_{0 s} s_{∥ s}

(62)

and

I_{⊥ s}^{pjk} = \frac{1}{N_{a}} \int d μ d v_{∥} d θ \frac{B}{m_{a}} \frac{H_{p} (s_{∥ s}) L_{j} (x_{s})}{\sqrt{2^{p} p!}} F_{M s} L_{k}^{3 / 2} (s_{s}^{2}) J_{1 s} \sqrt{x_{s}} .

(63)

To analytically evaluate the velocity integrals in $I_{∥ s}^{pjk}$ and $I_{⊥ 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:³⁰

J_{m} (b_{s} \sqrt{x_{s}}) = {(\frac{b_{s} \sqrt{x_{s}}}{2})}^{m} \sum_{n = 0}^{\infty} \frac{n! K_{n} (b_{s})}{(n + m)!} L_{n}^{m} (x_{s})

(64)

with the nth-order kernel function $K_{n} (b_{s}) = {(b_{s} / 2)}^{2 n} e^{- b_{s}^{2} / 4} / n!$ describing FLR effects. Then, using Eq. (64) with

L_{n}^{m} (x) L_{j} (x) x^{m} = \sum_{f = 0}^{n + m + j} d_{njf}^{m} L_{f} (x),

(65)

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_{∥ s}^{pjk}$ can be computed,

\begin{matrix} I_{∥ s}^{pjk} = \frac{2}{3 \sqrt{π}} \sum_{n = 0}^{\infty} \sum_{f = 0}^{n + j} {(T^{- 1})}_{p f}^{1 k 0} \frac{K_{n} (b_{s})}{\sqrt{2^{p} p!}} d_{njf}^{0} \frac{(k + 3 / 2)!}{k!} \\ \times [(p \geq 1 \cup f \geq 1) \cap (f + ⌊ p / 2 ⌋ \geq k)], \end{matrix}

(66)

with $[\cdot]$ the Iverson bracket (⁠ $[A] = 1$ if A is true, and 0 otherwise). Finally, the velocity integral contained in $I_{⊥ s}^{pjk}$ can be evaluated similarly to the one in Eq. (66) and yields

\begin{matrix} I_{⊥ s}^{pjk} = \sum_{n = 0}^{\infty} \sum_{f = 0}^{n + j + 1} \sum_{r = 0}^{f + ⌊ p / 2 ⌋} \sum_{q = 0}^{k} \sum_{r_{1} = 0}^{r} \frac{{(T^{- 1})}_{p f}^{0 r 0}}{\sqrt{π}} \frac{b_{s} K_{n} (b_{s})}{(n + 1) \sqrt{2^{p} p!}} \\ \times d_{njf}^{1} L_{k q}^{1} L_{r r_{1}}^{0} (1 / 2 + r_{1} + q)! . \end{matrix}

(67)

We remark that the expression of the numerical coefficients ${(T^{- 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 $Δ C_{a b}^{Tpj}$ obtained by projecting Eq. (9) onto the Hermite–Laguerre basis. In Eq. (59), one recognizes the quantities $I_{∥ s}^{pjk}$ and $I_{⊥ s}^{pjk}$ defined in Eqs. (66) and (67). Thus, using their definitions, the gyro-moment expansion of the test component of the correction terms, $Δ C_{a b}^{Tpj}$ ⁠, is deduced,

Δ C_{a b}^{Tpj} = \sum_{ℓ = 0}^{L} \sum_{k = 0}^{K} \frac{2 c_{ℓ}}{{\bar{τ}}_{a b}} Δ M_{a b}^{ℓ k} (\frac{{\bar{u}}_{∥ a}^{k}}{v_{T a}} I_{∥ a}^{p j ℓ} + \frac{{\bar{u}}_{⊥ a}^{k}}{v_{T a}} I_{⊥ a}^{p j ℓ}) .

(68)

We remark that $Δ 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_{∥ a}^{k}$ and $u_{⊥ 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, $Δ C_{a b}^{Fpj}$ ⁠, defined in Eq. (60), and inverting the species role between a and b in $u_{∥ s}^{k}$ and $u_{⊥ s}^{k}$ yields

Δ C_{a b}^{Fpj} = \sum_{ℓ = 0}^{L} \sum_{k = 0}^{K} \frac{2 c_{ℓ}}{{\bar{τ}}_{a b}} Δ N_{a b}^{ℓ k} (\frac{{\bar{u}}_{∥ b}^{k}}{v_{T a}} I_{∥ a}^{p j ℓ} + \frac{{\bar{u}}_{⊥ b}^{k}}{v_{T a}} I_{⊥ a}^{p j ℓ}),

(69)

where ${\bar{u}}_{∥ b}^{k}$ and ${\bar{u}}_{⊥ 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

Δ C_{a b}^{T} = \sum_{p = 0}^{\infty} \sum_{j = 0}^{\infty} Δ C_{a b}^{Tpj} \frac{H_{p} (s_{∥ a}) L_{j} (x_{a})}{\sqrt{2^{p} p!}} .

(70)

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:

\begin{matrix} Δ C_{a b}^{Tpj} = \sum_{ℓ = 0}^{min (L, j + ⌊ p / 2 ⌋)} \sum_{k = 0}^{K} \frac{{(T^{- 1})}_{p j}^{1 ℓ}}{\sqrt{2^{p} p!}} \frac{(ℓ + 3 / 2)!}{ℓ!} \frac{4 c_{ℓ} c_{k}}{3 \sqrt{π} {\bar{τ}}_{a b}} Δ M_{a b}^{ℓ k} \\ \times b \cdot M_{a}^{1 k} [p \geq 1 \cup j \geq 1] \end{matrix}

(71)

and

\begin{matrix} Δ C_{a b}^{Fpj} = \sum_{ℓ = 0}^{min (L, j + ⌊ p / 2 ⌋)} \sum_{k = 0}^{K} \frac{{(T^{- 1})}_{p j}^{1 ℓ}}{\sqrt{2^{p} p!}} \frac{(ℓ + 3 / 2)!}{ℓ!} \frac{4 c_{ℓ} c_{k}}{3 \sqrt{π} {\bar{τ}}_{a b}} \frac{1}{χ} Δ N_{a b}^{ℓ k} \\ b \cdot M_{b}^{1 k} [p \geq 1 \cup j \geq 1] . \end{matrix}

(72)

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} = \int d θdμd v_{∥} B g_{a} H_{p} L_{j} / m_{a} \sqrt{2^{p} p!}$ [see Eq. (55)], by

b \cdot M_{s}^{1 k} = \sum_{g = 0}^{1 + 2 k} \sum_{h = 0}^{k} T_{1 k}^{g h} \sqrt{2^{g} g!} N_{s}^{g h},

(73)

where we use the fact that $f_{a} ≃ 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 ≃ 3$ terms in $Δ 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 ℓ k}$ and $N_{a b}^{A ℓ 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^{- 50}$ ⁠, illustrating the robustness of the numerical framework used in this work. Additionally, the Braginskii matrices, $M_{a a}^{A ℓ k}$ and $N_{a a}^{A ℓ 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_b case), the Braginskii matrices $M_{a a}^{L ℓ k}$ and $M_{a a}^{S ℓ k}$ are equal as shown in the left panel of Fig. 2, implying $Δ M_{a b}^{ℓ k} = 0$ ⁠. On the other hand, the difference in the Braginskii matrices associated with the field components, i.e., $N_{a a}^{A ℓ 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.

FIG. 1.

View large Download slide

Verification of the momentum conservation, Eq. (14), obtained by using the closed expressions of the Braginksii matrices (numerically computed using 50 significant digits) as a function of k for the Coulomb (A = L, solid lines) and OS (A = S, dotted lines) collision operators. Collisions between electrons and ions with T_e = T_i and $m_{e} / m_{i} = 0.0027$ are considered.

FIG. 2.

View large Download slide

Braginskii matrices of (top) Coulomb and (bottom) OS collision operators for like-species (a = b) associated with the test (left) and the field (right) components, respectively, defined by $M_{a a}^{A ℓ k}$ and $N_{a a}^{A ℓ k}$ ⁠, as a function of $(ℓ, k)$ ⁠.

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_a in the DK limit) at $(p, j) = (P, J)$ ⁠, therefore, assuming that higher-order gyro-moments, p > P and j > J, vanish. Given (P, J), the gyro-moment expansion of the collision operator A can be written as

C_{a b}^{Apj} = \sum_{p' = 0}^{P} \sum_{j' = 0}^{J} C_{abp' j'}^{ATpj} n_{a}^{p' j'} + \sum_{p' = 0}^{P} \sum_{j' = 0}^{J} C_{abp' j'}^{AFpj} n_{b}^{p' j'} .

(74)

The coefficients $C_{abp' j'}^{ATpj}$ and $C_{abp' j'}^{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

\begin{matrix} C_{a b}^{ISpj} = \sum_{p' = 0}^{P} \sum_{j' = 0}^{J} (C_{abp' j'}^{STpj} + Δ C_{abp' j'}^{Tpj}) n_{a}^{p' j'} \\ + \sum_{p' = 0}^{P} \sum_{j' = 0}^{J} (C_{abp' j'}^{SFpj} + Δ C_{abp' j'}^{Fpj}) n_{b}^{p' j'} . \end{matrix}

(75)

In the case of the GK IS, the explicit expressions of the coefficients, $Δ C_{abp' j'}^{Tpj}$ and $Δ C_{abp' j'}^{Fpj}$ ⁠, can be obtained by inserting the analytical expressions of ${\bar{u}}_{∥ s}^{k}$ and ${\bar{u}}_{⊥ s}^{k}$ ⁠, given in Eq. (61), into Eqs. (68) and (69), respectively. In the case of the DK IS, the coefficients, $Δ C_{abp' j'}^{Tpj}$ and $Δ C_{abp' j'}^{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' j'}^{ATpj}$ and $C_{abp' j'}^{AFpj}$ ⁠, associated with the DK Coulomb and OS collision operators as well as $Δ C_{abp' j'}^{Tpj}$ and $Δ C_{abp' j'}^{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 $Δ C_{abp' j'}^{Tpj}$ and $Δ C_{abp' j'}^{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 $\bar{l} (p, j) = (J + 1) p + j + 1$ (where p and j run from 0 to P and J, respectively) and, similarly, the column index $\bar{l} (p', j')$ ⁠. 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,⁶ a consequence of vanishing polynomial basis coefficients ${(T^{- 1})}_{p j}^{lkm}$ ⁠.

FIG. 3.

View large Download slide

Matrix representation of the coupling between gyro-moments associated with Eqs. (74) and (75). The IS (left), OS (center), and their relative difference associated with the gyro-moment expansion of the correction term $Δ C_{a b}$ in Eq. (7) (right) are represented using the row index $\bar{l} (p, j) = (J + 1) p + j + 1$ and column index $\bar{l} (p', j') = (J + 1) p' + j' + 1$ ⁠. The matrix elements are obtained by evaluating numerically the analytical expressions derived in Sec. V. Both the GK (top) and DK (bottom) limits of the IS and OS operator are shown. The colorbars associated with the differences of the OS and IS operator are adjusted for better visualization. Here, we consider $(P, J) = (6, 3)$ gyro-moments and $L = K = 5$ and $k_{⊥} = 0.5$ ⁠.

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),³⁶ 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,

\begin{matrix} \frac{\partial}{\partial t} g_{a} + i ω_{B a} h_{a} + v_{| |} b \cdot \nabla h_{a} - \frac{μ}{m_{a}} (b \cdot \nabla B) \frac{\partial}{\partial v_{| |}} h_{a} \\ - i ω_{T a}^{*} {⟨ χ_{a} (r) ⟩}_{R} F_{M a} = \sum_{b} C_{a b}^{A}, \end{matrix}

(76)

where $χ_{a} (r) = ϕ (r) - v_{∥} ψ (r)$ ⁠, with $ϕ$ 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,

\sum_{a} \frac{q_{a}^{2}}{N_{a} T_{a}} (1 - Γ_{0} (a_{a})) ϕ (r) = \sum_{a} q_{a} \int d μ d v_{∥} d θ \frac{B}{m_{a}} J_{0} (b_{a} \sqrt{x_{a}}) g_{a}

(77)

and

(\frac{k_{⊥}^{2}}{4 π} + \sum_{a} \frac{q_{a}^{2} N_{a}}{m_{a}} Γ_{0} (a_{a})) ψ (r) = \sum_{a} q_{a} \int d μ d v_{∥} d θ \frac{B}{m_{a}} J_{0} (b_{a} \sqrt{x_{a}}) v_{∥} g_{a},

(78)

respectively. In Eqs. (76)–(78), we introduce the magnetic drift frequency $ω_{B a} = v_{T a}^{2} (x_{a} + 2 s_{∥ a}^{2}) ω_{B} / (2 Ω_{a})$ with $ω_{B} = (b \times \nabla ln B) \cdot k$ ⁠, the diamagnetic frequency $ω_{T a}^{*} = [ω_{N} + ω_{T_{a}} (x_{a} + s_{∥ a}^{2} - 3 / 2)]$ with $ω_{N} = b \times \nabla ln N \cdot k / B$ and $ω_{T_{a}} = b \times \nabla ln T_{a} \cdot k / B, a_{a} = b_{a}^{2} / 2$ ⁠, and $Γ_{0} (x) = I_{0} (x) e^{- x}$ (with I₀ 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 - α$ model (with α = 0).³⁷ In the local flux-tube approach, we assume constant radial density and temperature gradients, $L_{N}^{- 1}$ and $L_{T a}^{- 1}$ ⁠, with values $R_{N} = R_{0} / L_{N} = R_{0} / L_{T i} = R_{0} / L_{T e} = 20$ ⁠, where R₀ is the tokamak major radius. Electromagnetic effects are introduced with $β_{e} = 8 π 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 $ε = 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 ≳ 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_x spectrum around k_x = 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.

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 ω_r are plotted as a function of k_y 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-k_y and high-k_y modes. We also notice that the low-k_y and the high-k_y modes change continuously from electron to the ion diamagnetic directions as k_y increases. While the mode peaking at $k_{y} = 0.25$ (k_y 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 (⁠ $ω_{r} > 0$ ⁠), we identity it here as a TEM with a conventional mode structure. Indeed, despite $ω_{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} ≳ 0.5$ with an unconventional ballooning mode structures with secondary peaks located near $χ = \pm π / 2$ (χ is the ballooning angle) away from the outboard midplane (where the mode at $k_{y} \sim 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_y TEM and to the TEM peaking near $k_{y} = 0.6$ ⁠, with an unconventional mode structure, as the high-k_y 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.

FIG. 4.

View large Download slide

Collisionless growth rate γ (red lines) and frequency ω_r (blue lines) as a function of the binormal wavenumber k_y plotted on the same axis. The colored lines are the results obtained by the gyro-moment approach using $(P, J) = (24, 10)$ ⁠, and the black cross markers are the collisionless results using the GENE eigensolver. A positive mode frequency (⁠ $ω_{r} > 0$ ⁠) corresponds to the mode propagating in the ion diamagnetic direction and a negative mode frequency (⁠ $ω_{r} < 0$ ⁠) to the electron diamagnetic direction. Here, $T_{i} / T_{e} = 1$ ⁠.

FIG. 5.

View large Download slide

Modulus of the electrostatic ballooning eigenmode function, $ϕ (χ)$ [normalized to $ϕ (χ = 0)$ ], as a function of the ballooning angle χ corresponding to the case of the low-k_y TEM at $k_{y} = 0.25$ (blue solid line) and to the high-k_y TEM developing at $k_{y} = 0.6$ (red solid line).

We now investigate the collisionality dependence of both modes. While the electron collisionality is expected to be mainly in the banana regime with $ν_{e}^{*} = \sqrt{2} q ν / ε^{3 / 2} ≲ 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, $ν_{e}^{*} \sim 1 / ε^{3 / 2} ≫ 1$ ⁠, at the bottom of the pedestal.³⁶ 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 $ν_{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_y 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-k_y TEM is weakly affected by collisions when $ν_{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 $ν_{e}^{*}$ ⁠, while the DK IS is able to correct the DK OS and to approach the DK Coulomb as 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.

FIG. 6.

View large Download slide

The low-k_y TEM growth rate γ (left) and mode frequency ω_r (right) as a function of $ν_{e}^{*}$ obtained using the DK Coulomb (blue markers), the DK OS (red markers), the DK IS with K = 2 (solid magenta line), the DK IS with K = 5 (solid green line), and the DK IS with K = 10 collision operators (solid light blue line). We note that the DK IS with K = 5 (light green line without marker) overlaps with the K = 10 case (light blue line without marker). The results with the DK Coulomb and DK OS operators using $(P, J) = (16, 8)$ gyro-moments and 6 gyro-moments (6GM) with the analytical expressions reported in Appendix C are shown by the circle and cross markers, respectively. The insets focus on the Pfirsch–Schlüter regime, $20 < ν_{e}^{*} < 80$ ⁠. Here, the parameters are the same as in Fig. 4 at $k_{y} = 0.25$ with $T_{i} / T_{e} = 1$ ⁠.

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_N (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 γ, $σ (γ) = (γ - γ_{C}) / γ_{C}$ ⁠, where γ is obtained using the OS or IS collision operator models and γ_C is the one obtained using the DK Coulomb operator (see Fig. 7). The same definition $σ (ω_{r}) = (ω_{r} - ω_{r C}) / ω_{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 γ and ω_r (compared to the DK Coulomb) at low collisionality (with a peak near $ν_{e}^{*} \sim 1$ ⁠), and that difference changes sign in the Pfirsch–Schlüter regime, when $ν_{e}^{*} ≳ 10$ ⁠, where the DK OS operator overestimates γ and ω_r. 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 $ν_{e}^{*}$ (see the red areas in the left panels in Fig. 8), the correction terms in the IS operator reduce $σ (γ)$ [and $σ (ω_{r})$ ] below $2 %$ for all density gradients at high collisionalities. More precisely, the agreement between the DK IS and DK Coulomb improves with 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 $σ (γ)$ in the case K = 2 are reduced in the case of K = 5 (and K = 10) when $ν_{e}^{*} ≳ 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 ≳ 5$ ⁠. We notice that, in general, $σ (γ) ≲ σ (ω_{r})$ for all operators. Finally, we remark that, in both the DK OS and DK IS operators, the deviations in γ and ω_r peak near $ν_{e}^{*} \sim 1$ and increase at lower gradients with $σ (γ) \sim 5 %$ and $σ (ω_{r}) \sim 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_N, is mainly attributed to the decrease in γ and ω_r, as they are of the order of the diamagnetic frequency.

FIG. 7.

View large Download slide

Low-k_y TEM growth rate γ (left), and real mode frequency ω_r (right) as a function of the normalized density gradient, R_N, and collisionality, $ν_{e}^{*}$ ⁠, using the DK Coulomb collision operator.

FIG. 8.

View large Download slide

Signed relative difference of the growth rate, $σ (γ)$ (top), and frequency, $σ (ω_{r})$ (bottom), as a function of the density gradient R_N and electron collisionality $ν_{e}^{*}$ for the OS and IS (K = 2, K = 5 and K = 10) operators (from left to right, respectively). The results using the DK Coulomb, used as reference, are shown in Fig. 7. The other parameters are the same as in Fig. 4, with $k_{y} = 0.25$ and $T_{i} / T_{e} = 1$ ⁠.

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 $σ (γ)$ plotted as a function of the temperature ratio $T_{i} / T_{e}$ and shown in Fig. 9 in the Pfirsch–Schlüter regime when $ν_{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 $σ (γ) ≲ 1 %$ ⁠. The same observations can be made for the real mode frequency, ω_r.

FIG. 9.

View large Download slide

Signed relative difference of the growth rate γ, $σ (γ)$ ⁠, predicted by the OS (red cross) and IS with K = 2 (magenta), K = 5 (green), and K = 10 (light blue) operators. The dashed back line represents perfect agreement with the DK Coulomb collision operator. The same parameters as in Fig. 6 are used, except for $ν_{e}^{*} = 50$ ⁠.

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.¹⁶ We evaluate numerically these sums in $I_{∥ s}^{pjk}$ and $I_{⊥ 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_y TEM growth rate, γ, and mode frequency, ω_r, 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 $ν_{e}^{*}$ using the GK Coulomb, GK OS and GK IS with different values of 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_y 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^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 ≳ 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, ω_r, predicted by the GK Coulomb is well retrieved by the GK IS operators independently for 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_e, the test components of the GK OS and GK Coulomb are equivalent [see Eqs. (A1) and (A3)] since all the terms $Δ M_{a b}^{ℓ k}$ in $Δ 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 $Δ 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-k_y 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-k_y 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 %$ ⁠.

FIG. 10.

View large Download slide

High-k_y TEM growth rate γ (left) and the real mode frequency ω_r (right) as a function of $ν_{e}^{*}$ obtained using the GK Coulomb (blue circle markers), the GK OS (red circle markers), the GK IS with K = 2 (solid magenta line), the GK IS with K = 5 (solid green line), and the GK IS with K = 10 (solid light blue line) collision operators. For comparison, the predictions of the DK Coulomb operator (blue cross markers) are also shown. The parameters are the same as in Fig. 4 at $k_{y} = 0.6$ with $T_{i} / T_{e} = 1$ ⁠.

FIG. 11.

View large Download slide

Same as Fig. 10, but using the field component in the DK limit (DK-F) in all operators.

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.⁴¹ 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,⁴² the later study was extended to the assess the collisional ZF damping in the banana regime when $ν_{i}^{*} ≲ 1$ (with $ν_{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.⁴³ A refinement of the exponential decay in Ref. 43 of the ZF residual prediction $R_{z} (\infty) = ϕ_{z} (\infty) / ϕ_{z} (0)$ (with $ϕ_{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,

R_{z} (\infty) \to ς = \frac{ε^{2}}{q^{2}} \frac{1}{(1 + ε^{2} / q^{2})} .

(79)

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 $ν_{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) = ϕ_{z} (t) / ϕ_{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 ≃ 3$ is necessary for the DK IS to converge also in this case.

FIG. 12.

View large Download slide

Collisional ZF damping obtained using the DK Coulomb (blue line), DK OS (red line) and DK IS with $L = K = 2, 5, 10$ (magenta, green and light blue lines, respectively) collision operators in the Pfirsch–Shlüter regime $ν_{i}^{*} = 3.13$ ⁠. The collisionless (see Ref. 43) and collisional [see Eq. (79)] long time predictions, $ϖ = 1 / (1 + 1.6 q / \sqrt{ε})$ and ς, are plotted by the dashed black and blue lines, respectively. Here, $k_{x} = 0.05, q = 1.4$ ⁠, and $ε = 0.1$ ⁠.

FIG. 13.

View large Download slide

Collisional ZF damping obtained using the GK Coulomb (blue lines), GK OS (red lines) and GK IS collision operators with K = 2 (magenta lines) and K = 5 (green lines). The damping of an initial density perturbation with a radial wavenumber $k_{x} = 0.1$ (left) and $k_{x} = 0.2$ (right) are shown. The analytical collisionless and collisional predictions, ϖ and ς, are plotted for comparison. It is observed that the GK IS collision operator provides better approximation to the GK Coulomb operator. The parameters and the number of gyro-moments are the same as in Fig. 12.

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¹³ 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 ν ≲ 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 ν = 5$ after the damping of the GAM oscillations, for the case $k_{x} = 0.2$ ⁠. We plot $| g_{i} |$ as a function of $s_{∥_{i}}$ (at x_i = 0) and x_i (at $s_{∥ 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_{∥} | ≲ v_{T i}$ more strongly than the GK Coulomb operator. A similar observation can be made when $v_{⊥} ≲ v_{T i}$ ⁠, as shown in the right panel of Fig. 14.

FIG. 14.

View large Download slide

Modulus of the perturbed ion gyrocenter distribution function $| g_{i} |$ (normalized to its maximum value), at time $t ν = 5$ after the damping of the GAM oscillations and at the outboard midplane in Fig. 13, plotted as a function of $s_{∥ i}$ at x_i = 0 (left) and as a function of x_i at $s_{∥ i} = 0$ (right). The case of a Maxwellian distribution function $e^{- s_{∥ i}^{2} - x_{i}}$ is shown for comparison (dashed black lines).

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

\partial_{t} f_{e} - \frac{e}{m_{e}} E \partial_{v_{∥}} f_{e} = C_{e e}^{A T} (f_{e}) + C_{e e}^{A F} (f_{e}) + C_{e i}^{A T} (f_{e})

(80)

with $v_{∥} = v \cdot 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.,

\partial_{t} N_{e}^{p j} + \frac{e}{v_{T e} m_{e}} \sqrt{2 p} E N_{e}^{p - 1 j} = C_{e e}^{ATpj} + C_{e e}^{AFpj} + C_{e i}^{ATpj} .

(81)

We remark that the term associated with electric force, proportional to $\sqrt{p} E N_{e}^{p - 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_{∥ e} = - e N_{e} u_{e}$ (with $u_{e} = \int d v v_{∥} f_{e} / N_{e} = N_{e}^{10} v_{T e} / \sqrt{2}$ the parallel electron velocity), is established resulting from an applied electric field of normalized amplitude $e E / [\sqrt{m_{e} T_{e}} ν_{e e}] = 10^{- 3}$ (see Fig. 15). From the saturated current, the electrical conductivity, $σ_{∥ e} = j_{∥ e} / E$ ⁠, can be computed.

FIG. 15.

View large Download slide

Normalized saturated electron parallel velocity, $u_{e} / v_{T e}$ ⁠, as a function of $t ν_{e e}$ obtained using the Coulomb (blue lines), OS (red line), and IS with K = 2 (magenta line) and K = 5 (green dashed line) operators. Here, we use $(P, J) = (20, 5)$ gyro-moments.

In Fig. 16, we compare our numerical estimates of $σ_{∥ e}$ with the analytical prediction of the Spitzer conductivity,^25,45 given by $σ_{∥ e} = 16 T_{e}^{3 / 2} / [{(2 π)}^{3 / 2} \sqrt{m_{e}} Z^{2} e^{2} ln Λ])$ obtained in the large ion ionization Z limit. We note that in the $Z \to \infty$ limit, the electron and ion collisions can be modeled by the pitch-angle scattering collision operator, $C_{e i} (f_{e}) ≃ - ν_{e i}^{D} (v) L^{2} f_{e}$ [where $ν_{e i}^{D} (v)$ is the velocity-dependent deflection collision frequency, defined below Eq. (A4)], and the collisions between electrons can be neglected since $ν_{e e} / ν_{e i} \sim 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 \to \infty$ 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,⁴⁶ 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.

FIG. 16.

View large Download slide

Normalized electrical conductivity, $σ_{∥ e}$ ⁠, as a function of the ion ionization degree Z obtained using the same operators [and number of $(P, J) = (20, 5)$ gyro-moments] as in Fig. 15. The numerical results are compared with the Spitzer analytical conductivity⁴⁵ as shown by the black markers.

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.¹⁷ 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 ℓ k}$ and $N_{a}^{A ℓ k}$ ⁠, associated with the Coulomb and OS collision operators are obtained for arbitrary $(ℓ, k)$ using the spherical harmonic expansion.⁶ 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 $(ℓ, 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, ξ, θ)$ (with v the modulus of v and $ξ = v_{∥} / v$ the pitch angle) by⁸

\begin{matrix} C_{a b}^{L T} (f_{a}) = \frac{m_{a} ν_{a b} v_{T a}}{n_{b}} [(\frac{2 \partial_{v} G (f_{M b})}{v^{2}} + (1 - \frac{m_{a}}{m_{b}}) \partial_{v} H (f_{M b})) \partial_{v} f_{a} \\ - \frac{1}{v^{3}} \partial_{v} G (f_{M b}) L^{2} f_{a} + \partial_{v}^{2} G (f_{M b}) \partial_{v}^{2} f_{a} + \frac{m_{a}}{m_{b}} 8 π f_{M b} f_{a}] \end{matrix}

(A1)

and

\begin{matrix} C_{a b}^{L F} (f_{b}) = \frac{2 ν_{a b} v_{T a} f_{M a}}{n_{b}} [2 s_{a}^{2} \partial_{v}^{2} G (f_{b}) - H (f_{b}) \\ - (1 - \frac{m_{a}}{m_{b}}) v \partial_{v} H (f_{b}) + \frac{m_{a}}{m_{b}} 4 π v_{T a}^{2} f_{b}], \end{matrix}

(A2)

respectively. In Eq. (A1), we introduce the Rosenbluth potentials $H (f) = 2 \int d v' f (v') / u$ and $G (f) = \int d v' f (v') u$ (with $u = | v - v' |$ ⁠) and the operator $L^{2} f = \partial_{v} (v^{2} \partial_{v} f) - v^{2} \nabla_{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' (s_{b})] / 2$ and $H (f_{M b}) = 2 n_{b} erf (s_{b}) / v$ [with the error function, $erf (x) = 2 \int_{0}^{x} d s e^{- s^{2}} / \sqrt{π}$ ⁠, and its derivative, $erf' (x) = 2 e^{- x^{2}} / \sqrt{π}$ ]. 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 r 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})$ ].

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¹³

C_{a b}^{S T} (f_{a}) = C_{a b}^{0} (f_{a}) + \sum_{n = 1}^{3} X_{a b}^{n},

(A3)

where

C_{a b}^{0} (f_{a}) = - ν_{a b}^{D} (v) L^{2} f_{a} + \frac{1}{v^{2}} \partial_{v} [ν_{a b}^{∥} (v) v^{4} F_{a M} \partial_{v} (\frac{f_{a}}{f_{a M}})]

(A4)

and

\begin{matrix} X_{a b}^{1} = 2 (θ_{a b} - 1) f_{M a} [\frac{m_{a}}{T_{a}} \frac{1}{n_{a}} \int d v' v' C_{a b}^{0} (f_{a}) \cdot v \\ + (s_{a}^{2} - \frac{3}{2}) \frac{1}{n_{a}} \int d v \frac{2}{3} s_{a}^{2} C_{a b}^{0} (f_{a})], \end{matrix}

(A5)

X_{a b}^{2} = 2 (θ_{a b} - 1) [\frac{m_{a}}{T_{a}} u_{a} (f_{a}) \cdot C_{a b}^{0} (f_{M a} v) + \frac{δ T_{a} (f_{a})}{T_{a}} C_{a b}^{0} (f_{M a} s_{a}^{2})],

(A6)

\begin{matrix} X_{a b}^{3} = - \frac{8 ν_{a b}}{3 \sqrt{π}} f_{M a} \frac{χ {(θ_{a b} - 1)}^{2}}{{(1 + χ_{a b}^{2})}^{1 / 2}} [\frac{m_{a}}{T_{a}} v \cdot u_{a} (f_{a}) \\ + \frac{δ T_{a} (f_{a})}{T_{a}} (s_{a}^{2} - \frac{3}{2}) \frac{2}{(1 + χ_{a b}^{2})}] \end{matrix}

(A7)

and by

C_{a b}^{S F} (f_{b}) = - \frac{m_{a} V_{a b} (f_{b})}{T_{a}} \cdot C_{a b}^{S T} (f_{M a} v) - W_{a b} (f_{b}) C_{a b}^{S T} (f_{M a} s_{a}^{2}),

(A8)

respectively. In Eq. (A3), we introduce the velocity-dependent collision frequencies, $ν_{a b}^{D} (v) = ν_{a b} [erf (s_{b}) - Φ (s_{b})] / s_{a}^{3}$ and $ν_{a b}^{∥} (v) = 2 ν_{a b} Φ (s_{b}) / s_{a}^{3}$ [with $Φ (x) = [erf (x) - x erf' (x)] / (2 x^{2})$ being the Chandrasekhar function], the fluid quantities, $u_{a} (f_{a}) = \int d v v f_{a} / n_{a}, δ T_{a} (f_{a}) = T_{a} \int d v f_{a} (2 s_{a}^{2} / 3 - 1) / n_{a}$ ⁠, the coefficients $θ_{a b} = \sqrt{(T_{a} / T_{b} + χ_{a b}^{2}) / (1 + χ_{a b}^{2})}, τ = T_{a} / T_{b}, χ_{a b} = v_{T a} / v_{T b} = \sqrt{τ / σ}$ ⁠, and $σ = m_{a} / m_{b}$ ⁠. Finally, in Eq. (A8), we have

(A9)

V_{a b} (f_{b}) = \frac{m_{b}}{γ_{a b}} \int d v \frac{f_{b}}{f_{M b}} C_{b a}^{S T} (f_{M b} v),

(A9a)

W_{a b} (f_{b}) = \frac{T_{b}}{η_{a b}} \int d v \frac{f_{b}}{f_{M b}} C_{b a}^{S T} (f_{M b} s_{b}^{2})

(A9b)

with

γ_{a b} = - (n_{a} m_{a} χ_{a b}) (T_{a} / T_{b} + χ_{a b}^{2}) / [{\bar{τ}}_{a b} {(1 + χ_{a b}^{2})}^{3 / 2}]

and

η_{a b} = - 3 χ_{a b} T_{a} (T_{a} / T_{b} + χ_{a b}^{2}) / [{\bar{τ}}_{a b} {(1 + χ_{a b}^{2})}^{5 / 2}]

⁠. The expressions of

V_{a b} (f_{b})

and

W_{a b} (f_{b})

given in Eqs. are chosen such that particle, momentum, and energy conservations are satisfied [see Eq. ]. The closed analytical expressions of

C_{a b}^{0} (f_{M a} v), C_{a b}^{0} (f_{M a} s_{a}^{2})

⁠, and

C_{a b}^{S T} (f_{M a} v)

with

C_{a b}^{S T} (f_{M a} s_{a}^{2})

can be found in Ref. 13. We remark that

θ_{a b} = 1

when T_a = T_b 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. . Hence, the OS operator differs from the Coulomb collision operator only by its field component, given in Eq. , when T_a = T_b. 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. and , respectively, act on the perturbed particle distribution function 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.,

C_{a b}^{T} = {⟨ (T^{- 1} C_{a b}^{T}) (Z) ⟩}_{R} = \int_{0}^{2 π} |_{R} \frac{d θ}{2 π} C_{a b}^{T} (z (Z))

(B1)

and

C_{a b}^{F} = {⟨ (T^{- 1} C_{a b}^{F}) (Z) ⟩}_{R} = \int_{0}^{2 π} |_{R} \frac{d θ}{2 π} C_{a b}^{F} (z (Z)),

(B2)

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^{- 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 Z, i.e., $f (z (Z)) = (T^{- 1} f) (Z) = f (T^{- 1} Z)$ ⁠, being $z (Z) = T^{- 1} 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

\begin{matrix} M_{a}^{p j} (r) = \frac{1}{N_{a}} \int d R d v_{∥} d μ d θ \frac{B}{m_{a}} δ (R + ρ_{a} - r) \\ \times h_{a} (R (r, v), μ, v_{∥}) Y^{p} (s) L_{j}^{p + 1 / 2} (s_{a}^{2}) \end{matrix},

(B3)

with $R (r, v) ≃ r - ρ_{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^{- 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

\begin{matrix} C_{a b}^{T} (z (Z)) = (T^{- 1} C_{a b}^{T}) (Z) = C_{a b}^{T} (T^{- 1} Z) ≃ C_{a b}^{T} (R + ρ_{a} (μ, θ), μ, v_{∥}, θ), \end{matrix}

(B4)

where we use the lowest order gyrocenter coordinate transformations, in particular $r (Z) ≃ R + ρ_{a} (μ, θ)$ ⁠. 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 r appearing in $M_{a}^{p j}$ is transformed to $R + ρ_{a} (μ, θ)$ according to Eq. (B4), yields

\begin{matrix} C_{a b}^{L T} = \sum_{p = 0}^{\infty} \sum_{j = 0}^{\infty} \sum_{m = - p}^{p} \sum_{n = - \infty}^{\infty} \frac{f_{M a}}{σ_{j}^{p}} ν_{a b}^{Tpj} (v) e^{i k \cdot R} M_{a}^{p j} (k) \cdot e^{p m} \\ \times i^{n} J_{n} (b_{a} \sqrt{x_{a}}) \sqrt{\frac{2 π^{3 / 2} p!}{2^{p} (p + 1 / 2)!}} {⟨ e^{i n θ} Y_{p}^{m} (ξ, θ) ⟩}_{R}, \end{matrix}

(B5)

where Eq. (21) and the Jacobi–Anger identity, $e^{i k \cdot ρ_{a}} = \sum_{n} i^{n} J_{n} (b_{a} \sqrt{x_{a}}) e^{i n θ}$ ⁠, 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 $ρ_{a} (μ, θ)$ and the θ-dependence of $Y_{p}^{m} (ξ, θ)$ ⁠, 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} \sqrt{x_{a}}) = \sum_{j = 0}^{\infty} c_{j} {(b_{a} \sqrt{x_{a}} / 2)}^{2 j + n}$ ⁠, with $c_{j} = {(- 1)}^{j} / [j! (n + j)!]$ ⁠. These FLR terms ultimately damp small scale fluctuations,¹⁶ 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⁶) 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, μ, v_{∥}, θ)$ ⁠, as follows:²⁶

\frac{\partial}{\partial v} |_{r} f_{a} (r, v) = \sum_{p = 0}^{\infty} \sum_{j = 0}^{\infty} \sum_{ℓ = 0}^{j} \frac{L_{j ℓ}^{p}}{σ_{j}^{p}} \cdot M_{a}^{p j} (r (Z)) \frac{\partial}{\partial v} |_{r} (s_{a}^{2 ℓ} Y^{p} (s_{a}) f_{M a}),

(B6)

with $r (Z) ≃ R + ρ_{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, μ, v_{∥}, θ)$ ⁠, and in terms of h_a by using the chain rule while holding r constant, i.e.,

\begin{matrix} \frac{\partial}{\partial v} |_{r} f_{a} (r, v) = \frac{\partial Z}{\partial v} |_{r} \frac{\partial}{\partial Z} f_{a} (r (Z), v (Z)) ≃ \frac{\partial}{\partial v} |_{R} h_{a} - \frac{\partial ρ_{a}}{\partial v} |_{r} \cdot \nabla h_{a} \end{matrix}

(B7)

with $f_{a} ≃ h_{a}$ [see Eq. (31)], and by approximating

\begin{matrix} G (f_{b} (r, v)) ≃ G (h_{b} (R (r, v), μ, v_{∥})) \\ = \int d v' | v' - v | h_{b} (r - ρ_{b} (θ', μ'), μ', v_{∥}^{'}) \end{matrix}

(B8)

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 $\sim ν_{a b} k_{⊥}^{2} ρ_{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_a 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

M_{a}^{p j} (r (Z)) ≃ M_{a}^{p j} (R) + ρ_{a} \cdot \nabla M_{a}^{p j} (R) .

(B9)

Using Eq. (B9) into Eq. (B6), the fact that $ρ_{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

\begin{matrix} \frac{\partial}{\partial v} |_{r} f_{a} (r, v) ≃ \frac{\partial}{\partial v} |_{r} f_{a} (R + ρ_{a}, v) - \frac{\partial ρ_{a}}{\partial v} |_{r} \cdot \nabla f_{a} (R, v) \\ ≃ \frac{\partial}{\partial v} |_{R} h_{a} - \frac{\partial ρ_{a}}{\partial v} |_{r} \cdot \nabla h_{a} . \end{matrix}

(B10)

In Eq. (B10), we also use the scalar invariance of the velocity–space derivative, i.e., $\partial_{v} |_{r} f_{a} (R + ρ_{a}, v) ≃ \partial_{v} |_{R} h_{a} (R, μ, v_{∥})$ ⁠, and $f_{a} ≃ 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_{⊥} ρ_{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 $\sim ν_{a b} k_{⊥}^{2} ρ_{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 $ρ_{a}^{2} k_{⊥}^{2}$ ⁠, depends on the perpendicular energy coordinate $v_{⊥}^{2}$ ⁠. Hence, at large $k_{⊥}^{2} v_{⊥}^{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_{⊥}^{2}$ and $k_{⊥}^{2}$ ⁠, respectively. Numerical tests reveal that truncating these sums at $j \sim 12$ and $n \sim 8$ is sufficient in the cases discussed herein where $k_{⊥} ρ_{i} ≲ 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., $λ_{mfp} / L_{∥} ≪ 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 $λ_{mfp} / L_{∥}$ ⁠. On the other hand, the gyro-moments associated with the non-Maxwellian component of g_a, i.e., the parallel and perpendicular heat fluxes (3, 0) and (1, 1), respectively, are first order in $λ_{mfp} / L_{∥}$ ⁠.

By projecting the GK Boltzmann equation on the Hermite–Laguerre basis²¹ [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_{⊥}$ ⁠). Using the closed analytical formulas of the DK IS, given in Eqs. (71) and (72), the non-vanishing terms $Δ C_{abp' j'}^{Tpj}$ in the gyro-moment expansion of the test component $Δ C_{a b}^{Tpj}$ ⁠, are given by

Δ C_{a b 30}^{T 10} = \frac{4 τ^{3 / 2} ν_{a b}}{5 {(σ + τ)}^{5 / 2}} \sqrt{\frac{6}{π}} [τ (- \sqrt{\frac{(σ + 1) τ}{σ + τ}} + σ + 1) - σ \sqrt{\frac{(σ + 1) τ}{σ + τ}}],

(C1a)

Δ C_{a b 11}^{T 10} = \frac{8 τ^{3 / 2} ν_{a b}}{5 \sqrt{π} {(σ + τ)}^{5 / 2}} [σ \sqrt{\frac{(σ + 1) τ}{σ + τ}} + τ (\sqrt{\frac{(σ + 1) τ}{σ + τ}} - σ - 1)],

(C1b)

Δ C_{a b 10}^{T 30} = - \frac{4 \sqrt{τ} ν_{a b}}{5 {(σ + τ)}^{5 / 2}} \sqrt{\frac{2}{3 π}} [10 σ^{2} (τ - 1) + 3 \sqrt{\frac{(σ + 1) τ^{5}}{σ + τ}} + σ τ (3 \sqrt{\frac{(σ + 1) τ}{σ + τ}} + τ - 4) - 3 τ^{2}],

(C1c)

Δ C_{a b 30}^{T 30} = - \frac{12 ν_{a b} σ (τ - 1) \sqrt{τ} (10 σ^{2} - 2 σ τ + 3 τ^{2})}{25 \sqrt{π} {(σ + τ)}^{7 / 2}},

(C1d)

Δ C_{a b 11}^{T 30} = \frac{4 ν_{a b} \sqrt{\frac{6}{π}} σ (τ - 1) \sqrt{τ} (10 σ^{2} - 2 σ τ + 3 τ^{2})}{25 {(σ + τ)}^{7 / 2}},

(C1e)

Δ C_{a b 10}^{T 11} = \frac{8 ν_{a b} \sqrt{τ}}{15 \sqrt{π} {(σ + τ)}^{5 / 2}} [10 σ^{2} (τ - 1) + 3 \sqrt{\frac{(σ + 1) τ^{5}}{σ + τ}} + σ τ (3 \sqrt{\frac{(σ + 1) τ}{σ + τ}} + τ - 4) - 3 τ^{2}],

(C1f)

Δ C_{a b 30}^{T 11} = Δ C_{a b 11}^{T 30},

(C1g)

Δ C_{a b 11}^{T 11} = - \frac{8 ν_{a b} σ (τ - 1) \sqrt{τ} (10 σ^{2} - 2 σ τ + 3 τ^{2})}{25 \sqrt{π} {(σ + τ)}^{7 / 2}},

(C1h)

where $τ = T_{a} / T_{b}, σ = m_{a} / m_{b}$ are the temperature and mass ratios of the colliding species, respectively. Similarly, the non-vanishing terms $Δ C_{a b, p' j'}^{Fpj}$ in the gyro-moment expansion of the field component $Δ C_{a b}^{Fpj}$ [see Eq. (74)] are

Δ C_{a b 30}^{F 10} = - \frac{4 \sqrt{\frac{6}{π}} σ^{3 / 2} τ (- \sqrt{(σ + 1) (σ + τ)} + σ + 1)}{5 {(σ + τ)}^{5 / 2}},

(C2a)

Δ C_{a b 11}^{F 10} = \frac{8 σ \sqrt{σ τ} (- \sqrt{(σ + 1) τ^{3}} + σ (\sqrt{τ (σ + τ)} - \sqrt{(σ + 1) τ}) + \sqrt{τ (σ + τ)})}{5 \sqrt{π} {(σ + τ)}^{3}},

(C2b)

Δ C_{a b 10}^{F 30} = \frac{4 \sqrt{\frac{6}{π}} \sqrt{σ} τ (\sqrt{(σ + 1) τ^{3}} + τ (- \sqrt{σ + τ}) + σ (- 3 τ \sqrt{σ + τ} + \sqrt{(σ + 1) τ} + 2 \sqrt{σ + τ}))}{5 {(σ + τ)}^{3}},

(C2c)

Δ C_{a b 30}^{F 30} = \frac{36 σ^{3 / 2} τ (σ (5 τ - \sqrt{τ} - 2) - (\sqrt{τ} - 3) τ)}{25 \sqrt{π} {(σ + τ)}^{7 / 2}},

(C2d)

Δ C_{a b 11}^{F 30} = \frac{12 \sqrt{\frac{6}{π}} σ^{3 / 2} τ (σ (- 5 τ + \sqrt{τ} + 2) + (\sqrt{τ} - 3) τ)}{25 {(σ + τ)}^{7 / 2}},

(C2e)

Δ C_{a b 10}^{F 11} = - \frac{8 \sqrt{σ τ} (σ (- 3 \sqrt{τ^{3} (σ + τ)} + \sqrt{σ + 1} τ + 2 \sqrt{τ (σ + τ)}) - \sqrt{τ^{3} (σ + τ)} + \sqrt{σ + 1} τ^{2})}{5 \sqrt{π} {(σ + τ)}^{3}},

(C2f)

Δ C_{a b 30}^{F 11} = Δ C_{a b 11}^{F 30},

(C2g)

Δ C_{a b 11}^{F 11} = \frac{24 σ^{3 / 2} τ (σ (5 τ - \sqrt{τ} - 2) - (\sqrt{τ} - 3) τ)}{25 \sqrt{π} {(σ + τ)}^{7 / 2}} .

(C2h)

We notice that, from the conservation of particle in Eq. (2a), it follows $Δ C_{a b, p' j'}^{T 00} = Δ C_{a b, p' j'}^{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, $Δ C_{aap' j'}^{Tpj} = 0$ ⁠, while the non-vanishing terms of the field component $Δ C_{aap' j'}^{Fpj}$ are given by

Δ C_{a a 30}^{F 30} = \frac{9}{25} \sqrt{\frac{2}{π}},

(C3a)

Δ C_{a a 11}^{F 30} = - \frac{6}{25} \sqrt{\frac{3}{π}},

(C3b)

Δ C_{a a 11}^{F 11} = \frac{6}{25} \sqrt{\frac{2}{π}} .

(C3c)

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

C_{a b 10}^{S T 10} = - \frac{8 ν_{a b} (σ + 1)}{3 \sqrt{π}} {(\frac{τ}{σ + τ})}^{3 / 2},

(C4a)

C_{a b 30}^{S T 10} = \frac{4 ν_{a b}}{5} \sqrt{\frac{6}{π}} \frac{\sqrt{σ + 1} τ^{2}}{{(σ + τ)}^{2}},

(C4b)

C_{a b 11}^{S T 10} = - \frac{8 ν_{a b} \sqrt{σ + 1} τ^{2}}{5 \sqrt{π} {(σ + τ)}^{2}},

(C4c)

C_{a b 20}^{S T 20} = - \frac{16 ν_{a b} \sqrt{τ} (5 σ^{2} (τ + 2) + 21 σ τ + 6 τ^{2})}{45 \sqrt{π} {(σ + τ)}^{5 / 2}},

(C4d)

C_{a b 01}^{S T 20} = - \frac{16 ν_{a b} \sqrt{\frac{2}{π}} \sqrt{τ} (- 5 σ^{2} (τ - 1) + 3 σ τ + 3 τ^{2})}{45 {(σ + τ)}^{5 / 2}},

(C4e)

C_{a b 01}^{S T 01} = - \frac{16 ν_{a b} \sqrt{τ} (5 σ^{2} + 2 (5 σ + 9) σ τ + 3 τ^{2})}{45 \sqrt{π} {(σ + τ)}^{5 / 2}},

(C4f)

C_{a b 30}^{S T 30} = - \frac{4 ν_{a b} \sqrt{τ} (70 σ^{2} + 56 σ τ + 31 τ^{2})}{35 \sqrt{π} {(σ + τ)}^{5 / 2}},

(C4g)

C_{a b 11}^{S T 30} = - \frac{4 ν_{a b} \sqrt{\frac{2}{3 π}} τ^{3 / 2} (28 σ + τ)}{35 {(σ + τ)}^{5 / 2}},

(C4h)

C_{a b 11}^{S T 11} = - \frac{8 ν_{a b} \sqrt{τ} (105 σ^{2} + 98 σ τ + 47 τ^{2})}{105 \sqrt{π} {(σ + τ)}^{5 / 2}}

(C4i)

with $C_{ablk}^{STpj} = C_{abpj}^{STlk}$ ⁠, and

C_{a b 10}^{S F 10} = \frac{8 ν_{a b} σ (σ + 1) τ}{3 \sqrt{π} \sqrt{σ {(σ + τ)}^{3}}},

(C5a)

C_{a b 30}^{S F 10} = - \frac{4 ν_{a b} \sqrt{\frac{6}{π}} \sqrt{σ^{3} (σ + 1)} τ}{5 {(σ + τ)}^{2}},

(C5b)

C_{a b 11}^{S F 10} = \frac{8 ν_{a b} \sqrt{σ^{3} (σ + 1)} τ}{5 \sqrt{π} {(σ + τ)}^{2}},

(C5c)

C_{a b 20}^{S F 20} = \frac{16 ν_{a b} σ (σ + 1) τ^{3 / 2}}{9 \sqrt{π} {(σ + τ)}^{5 / 2}},

(C5d)

C_{a b 01}^{S F 20} = - \frac{16 ν_{a b} \sqrt{\frac{2}{π}} σ (σ + 1) τ^{3 / 2}}{9 {(σ + τ)}^{5 / 2}},

(C5e)

C_{a b 20}^{S F 01} = - \frac{16 ν_{a b} \sqrt{\frac{2}{π}} σ (σ + 1) τ^{3 / 2}}{9 {(σ + τ)}^{5 / 2}},

(C5f)

C_{a b 01}^{S F 01} = \frac{32 ν_{a b} σ (σ + 1) τ^{3 / 2}}{9 \sqrt{π} {(σ + τ)}^{5 / 2}},

(C5g)

C_{a b 10}^{S F 30} = - \frac{4 ν_{a b} \sqrt{\frac{6}{π}} \sqrt{σ (σ + 1) τ^{3}}}{5 {(σ + τ)}^{2}},

(C5h)

C_{a b 30}^{S F 30} = \frac{36 ν_{a b} {(σ τ)}^{3 / 2}}{25 \sqrt{π} {(σ + τ)}^{5 / 2}},

(C5i)

C_{a b 11}^{S F 30} = - \frac{12 ν_{a b} \sqrt{\frac{6}{π}} {(σ τ)}^{3 / 2}}{25 {(σ + τ)}^{5 / 2}},

(C5j)

C_{a b 10}^{S F 11} = \frac{8 ν_{a b} \sqrt{σ (σ + 1) τ^{3}}}{5 \sqrt{π} {(σ + τ)}^{2}},

(C5k)

C_{a b 30}^{S F 11} = - \frac{12 ν_{a b} \sqrt{\frac{6}{π}} {(σ τ)}^{3 / 2}}{25 {(σ + τ)}^{5 / 2}},

(C5l)

C_{a b 11}^{S F 11} = \frac{24 ν_{a b} {(σ τ)}^{3 / 2}}{25 \sqrt{π} {(σ + τ)}^{5 / 2}}

(C5m)

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

C_{a a 20}^{S 20} = - \frac{64}{45} \sqrt{\frac{2}{π}},

(C6a)

C_{a a 01}^{S 20} = - \frac{64}{45} \sqrt{\frac{1}{π}},

(C6b)

C_{a a 01}^{S 01} = - \frac{32}{45} \sqrt{\frac{2}{π}},

(C6c)

C_{a a 30}^{S 30} = - \frac{361}{175} \sqrt{\frac{2}{π}},

(C6d)

C_{a a 11}^{S 30} = - \frac{208}{175} \sqrt{\frac{1}{3 π}},

(C6e)

C_{a a 11}^{S 11} = - \frac{1187}{525} \sqrt{\frac{2}{π}} .

(C6f)

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

C_{a b 10}^{L T 10} = - \frac{8 ν_{a b} (σ + 1)}{3 \sqrt{π}} {(\frac{τ}{σ + τ})}^{3 / 2},

(C7a)

C_{a b 30}^{L T 10} = \frac{4}{5} \sqrt{\frac{6}{π}} (σ + 1) {(\frac{τ}{σ + τ})}^{5 / 2},

(C7b)

C_{a b 11}^{L T 10} = - \frac{8 ν_{a b} (σ + 1)}{5 \sqrt{π} {(\frac{τ}{σ + τ})}^{5 / 2}},

(C7c)

C_{a b 00}^{L T 20} = - \frac{8 ν_{a b} \sqrt{\frac{2}{π}} (τ - 1) \sqrt{\frac{τ}{σ}}}{3 {(\frac{σ + τ}{σ})}^{3 / 2}},

(C7d)

C_{a b 20}^{L T 20} = - \frac{ν_{a b} 8 {(\frac{τ}{σ})}^{3 / 2} (σ (10 σ + τ + 13) + 4 τ)}{15 \sqrt{π} σ {(\frac{σ + τ}{σ})}^{5 / 2}},

(C7e)

C_{a b 01}^{L T 20} = \frac{8 ν_{a b} \sqrt{\frac{2}{π}} {(\frac{τ}{σ})}^{3 / 2} (- 3 σ τ + σ - 2 τ)}{15 σ {(\frac{σ + τ}{σ})}^{5 / 2}},

(C7f)

C_{a b 00}^{L T 01} = \frac{16 ν_{a b} (τ - 1) \sqrt{\frac{τ}{σ}}}{3 \sqrt{π} {(\frac{σ + τ}{σ})}^{3 / 2}},

(C7g)

C_{a b 20}^{L T 01} = \frac{8 \sqrt{\frac{2}{π}} {(\frac{τ}{σ})}^{3 / 2} (- 3 σ τ + σ - 2 τ)}{15 σ {(\frac{σ + τ}{σ})}^{5 / 2}},

(C7h)

C_{a b 01}^{L T 01} = - \frac{16 {(\frac{τ}{σ})}^{3 / 2} (σ (5 σ - τ + 7) + τ)}{15 \sqrt{π} σ {(\frac{σ + τ}{σ})}^{5 / 2}},

(C7i)

C_{a b 10}^{L T 30} = - \frac{4 ν_{a b}}{5 σ^{2}} \sqrt{\frac{2}{3 π}} \sqrt{\frac{τ}{σ}} {(\frac{σ}{σ + τ})}^{5 / 2} (10 σ^{2} (τ - 1) + σ (τ - 4) τ - 3 τ^{2}),

(C7j)

C_{a b 30}^{L T 30} = - \frac{4 ν_{a b} {(\frac{τ}{σ + τ})}^{3 / 2} (14 σ^{2} (τ + 8) + 70 σ^{3} + σ τ (19 τ + 68) + 31 τ^{2})}{35 \sqrt{π} {(σ + τ)}^{2}},

(C7k)

C_{a b 11}^{L T 30} = \frac{4 ν_{a b} \sqrt{\frac{2}{3 π}} {(\frac{τ}{σ + τ})}^{3 / 2} (σ^{2} (14 - 42 τ) + σ τ (3 τ - 32) - τ^{2})}{35 {(σ + τ)}^{2}},

(C7l)

C_{a b 10}^{L T 11} = \frac{8 ν_{a b} \sqrt{\frac{τ}{σ}} (10 σ^{2} (τ - 1) + σ (τ - 4) τ - 3 τ^{2})}{15 \sqrt{π} σ^{2} {(\frac{σ + τ}{σ})}^{5 / 2}},

(C7m)

C_{a b 30}^{L T 11} = \frac{4 ν_{a b} \sqrt{\frac{2}{3 π}} {(\frac{τ}{σ + τ})}^{3 / 2} (σ^{2} (14 - 42 τ) + σ τ (3 τ - 32) - τ^{2})}{35 {(σ + τ)}^{2}},

(C7n)

C_{a b 11}^{L T 11} = - \frac{8 ν_{a b} {(\frac{τ}{σ + τ})}^{3 / 2} (7 (15 σ + 23) σ^{2} + (27 σ + 47) τ^{2} + 2 (21 σ + 59) σ τ)}{105 \sqrt{π} {(σ + τ)}^{2}},

(C7o)

and

C_{a b 10}^{L F 10} = \frac{8 ν_{a b} (σ + 1) τ}{3 \sqrt{π} σ} {(\frac{σ}{σ + τ})}^{3 / 2},

(C8a)

C_{a b 30}^{L F 10} = - \frac{4 ν_{a b} \sqrt{\frac{6}{π}} (σ + 1) τ}{5 σ} {(\frac{σ}{σ + τ})}^{5 / 2},

(C8b)

C_{a b 11}^{L F 10} = \frac{8 ν_{a b} (σ + 1) τ}{5 \sqrt{π} σ} {(\frac{σ}{σ + τ})}^{5 / 2},

(C8c)

C_{a b 00}^{L F 20} = - \frac{8 ν_{a b} \sqrt{\frac{2}{π}} (τ - 1) \sqrt{\frac{τ}{σ}}}{3 {(\frac{σ + τ}{σ})}^{3 / 2}},

(C8d)

C_{a b 20}^{L F 20} = \frac{8 ν_{a b} \sqrt{\frac{τ}{σ}} (σ (3 τ - 1) + 2 τ)}{5 \sqrt{π} σ {(\frac{σ + τ}{σ})}^{5 / 2}},

(C8e)

C_{a b 01}^{L F 20} = \frac{8 ν_{a b} \sqrt{\frac{2}{π}} \sqrt{\frac{τ}{σ}} (- 3 σ τ + σ - 2 τ)}{15 σ {(\frac{σ + τ}{σ})}^{5 / 2}},

(C8f)

C_{a b 00}^{L F 01} = \frac{16 (τ - 1) \sqrt{\frac{τ}{σ}}}{3 \sqrt{π} {(\frac{σ + τ}{σ})}^{3 / 2}},

(C8g)

C_{a b 20}^{L F 01} = \frac{8 \sqrt{\frac{2}{π}} \sqrt{\frac{τ}{σ}} (- 3 σ τ + σ - 2 τ)}{15 σ {(\frac{σ + τ}{σ})}^{5 / 2}},

(C8h)

C_{a b 01}^{L F 01} = \frac{32 \sqrt{\frac{τ}{σ}} (σ (3 τ - 1) + 2 τ)}{15 \sqrt{π} σ {(\frac{σ + τ}{σ})}^{5 / 2}},

(C8i)

C_{a b 10}^{L F 30} = - \frac{4 ν_{a b} \sqrt{\frac{6}{π}} τ (σ (3 τ - 2) + τ)}{5 σ^{2} {(\frac{σ + τ}{σ})}^{5 / 2}},

(C8j)

C_{a b 30}^{L F 30} = \frac{12 ν_{a b} τ (σ (5 τ - 2) + 3 τ)}{7 \sqrt{π} {(σ + τ)}^{2} {(\frac{σ + τ}{σ})}^{3 / 2}},

(C8k)

C_{a b 11}^{L F 30} = - \frac{12 ν_{a b} \sqrt{\frac{6}{π}} τ (σ (5 τ - 2) + 3 τ)}{35 {(σ + τ)}^{2} {(\frac{σ + τ}{σ})}^{3 / 2}},

(C8l)

C_{a b 10}^{L F 11} = \frac{8 ν_{a b} τ (σ (3 τ - 2) + τ)}{5 \sqrt{π} σ^{2} {(\frac{σ + τ}{σ})}^{5 / 2}},

(C8m)

C_{a b 30}^{L F 11} = - \frac{12 ν_{a b} \sqrt{\frac{6}{π}} τ (σ (5 τ - 2) + 3 τ)}{35 {(σ + τ)}^{2} {(\frac{σ + τ}{σ})}^{3 / 2}},

(C8n)

C_{a b 11}^{L F 11} = \frac{48 ν_{a b} τ (σ (5 τ - 2) + 3 τ)}{35 \sqrt{π} {(σ + τ)}^{2} {(\frac{σ + τ}{σ})}^{3 / 2}},

(C8o)

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

C_{a b 20}^{L 20} = - \frac{16}{15} \sqrt{\frac{2}{π}},

(C9a)

C_{a b 01}^{L 20} = - \frac{16}{15} \sqrt{\frac{1}{π}},

(C9b)

C_{a b 01}^{L 01} = - \frac{8}{15} \sqrt{\frac{2}{π}},

(C9c)

C_{a b 30}^{L 30} = - \frac{8}{5} \sqrt{\frac{2}{π}},

(C9d)

C_{a b 11}^{L 30} = - \frac{8}{5} \sqrt{\frac{1}{3 π}},

(C9e)

C_{a b 11}^{L 11} = - \frac{28}{15} \sqrt{\frac{2}{π}} .

(C9f)

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.

References

1.

E. A.

Belli

and

J.

Candy

, “

Implications of advanced collision operators for gyrokinetic simulation

,”

Plasma Phys. Controlled Fusion

59

,

045005

(

2017

).

https://doi.org/10.1088/1361-6587/aa5c94

Google Scholar

Crossref

2.

M.

Barnes

,

I. G.

Abel

,

W.

Dorland

,

D. R.

Ernst

,

G. W.

Hammett

,

P.

Ricci

,

B. N.

Rogers

,

A. A.

Schekochihin

, and

T.

Tatsuno

, “

Linearized model Fokker–Planck collision operators for gyrokinetic simulations. II. Numerical implementation and tests

,”

Phys. Plasmas

16

,

072107

(

2009

).

https://doi.org/10.1063/1.3155085

Google Scholar

Crossref

3.

P.

Manas

,

Y.

Camenen

,

S.

Benkadda

,

W. A.

Hornsby

, and

A. G.

Peeters

, “

Enhanced stabilisation of trapped electron modes by collisional energy scattering in tokamaks

,”

Phys. Plasmas

22

,

062302

(

2015

).

https://doi.org/10.1063/1.4922754

Google Scholar

Crossref

4.

Q.

Pan

,

D. R.

Ernst

, and

P.

Crandall

, “

First implementation of gyrokinetic exact linearized Landau collision operator and comparison with models

,”

Phys. Plasmas

27

,

042307

(

2020

).

https://doi.org/10.1063/1.5143374

Google Scholar

Crossref

5.

Q.

Pan

,

D. R.

Ernst

, and

D. R.

Hatch

, “

Importance of gyrokinetic exact Fokker-Planck collisions in fusion plasma turbulence

,”

Phys. Rev. E

103

,

L051202

(

2021

).

https://doi.org/10.1103/PhysRevE.103.L051202

Google Scholar

Crossref

PubMed

6.

B. J.

Frei

,

J.

Ball

,

A. C. D.

Hoffmann

,

R.

Jorge

,

P.

Ricci

, and

L.

Stenger

, “

Development of advanced linearized gyrokinetic collision operators using a moment approach

,”

J. Plasma Phys.

87

,

905870501

(

2021

).

https://doi.org/10.1017/S0022377821000830

Google Scholar

Crossref

7.

E. A.

Belli

and

J.

Candy

, “

Full linearized Fokker–Planck collisions in neoclassical transport simulations

,”

Plasma Phys. Controlled Fusion

54

,

015015

(

2011

).

https://doi.org/10.1088/0741-3335/54/1/015015

Google Scholar

Crossref

8.

M. N.

Rosenbluth

,

W. M.

MacDonald

, and

D. L.

Judd

, “

Fokker-Planck equation for an inverse-square force

,”

Phys. Rev.

107

,

1

(

1957

).

Google Scholar

Crossref

9.

J. P.

Dougherty

, “

Model Fokker-Planck equation for a plasma and its solution

,”

Phys. Fluids

7

,

1788

(

1964

).

https://doi.org/10.1063/1.2746779

Google Scholar

Crossref

10.

S. P.

Hirshman

and

D. J.

Sigmar

, “

Approximate Fokker-Planck collision operator for transport theory applications

,”

Phys. Fluids

19

,

1532

(

1976

).

https://doi.org/10.1063/1.861356

Google Scholar

Crossref

11.

I. G.

Abel

,

M.

Barnes

,

S. C.

Cowley

,

W.

Dorland

, and

A. A.

Schekochihin

, “

Linearized model Fokker-Planck collision operators for gyrokinetic simulations. I. Theory

,”

Phys. Plasmas

15

,

122509

(

2008

).

https://doi.org/10.1063/1.3046067

Google Scholar

Crossref

12.

M.

Francisquez

,

J.

Juno

,

A.

Hakim

,

G. W.

Hammett

, and

D. R.

Ernst

, “

Improved multispecies Dougherty collisions

,” arXiv:2109.10381 (

2021

).

Google Scholar

13.

H.

Sugama

,

T.-H.

Watanabe

, and

M.

Nunami

, “

Linearized model collision operators for multiple ion species plasmas and gyrokinetic entropy balance equations

,”

Phys. Plasmas

16

,

112503

(

2009

).

https://doi.org/10.1063/1.3257907

Google Scholar

Crossref

14.

M.

Nunami

,

M.

Nakata

,

T.-H.

Watanabe

, and

H.

Sugama

, “

Development of linearized collision operator for multiple ion species in gyrokinetic flux-tube simulations

,”

Plasma Fusion Res.

10

,

1403058

(

2015

).

https://doi.org/10.1585/pfr.10.1403058

Google Scholar

Crossref

15.

M.

Nakata

,

M.

Nunami

,

T.-H.

Watanabe

, and

H.

Sugama

, “

Improved collision operator for plasma kinetic simulations with multi-species ions and electrons

,”

Comput. Phys. Commun.

197

,

61

(

2015

).

https://doi.org/10.1016/j.cpc.2015.08.007

Google Scholar

Crossref

16.

B. J.

Frei

,

A. C. D.

Hoffmann

, and

P.

Ricci

, “

Local gyrokinetic collisional theory of the ion-temperature gradient mode

,”

J. Plasma Phys.

88

,

905880304

(

2022

).

https://doi.org/10.1017/S0022377822000344

Google Scholar

Crossref

17.

H.

Sugama

,

S.

Matsuoka

,

S.

Satake

,

M.

Nunami

, and

T.-H.

Watanabe

, “

Improved linearized model collision operator for the highly collisional regime

,”

Phys. Plasmas

26

,

102108

(

2019

).

https://doi.org/10.1063/1.5115440

Google Scholar

Crossref

18.

S. P.

Hirshman

and

D. J.

Sigmar

, “

Neoclassical transport of impurities in tokamak plasmas

,”

Nucl. Fusion

21

,

1079

(

1981

).

https://doi.org/10.1088/0029-5515/21/9/003

Google Scholar

Crossref

19.

M.

Honda

, “

Impact of higher-order flows in the moment equations on Pfirsch-Schlüter friction coefficients

,”

Phys. Plasmas

21

,

092508

(

2014

).

https://doi.org/10.1063/1.4895795

Google Scholar

Crossref

20.

S.

Matsuoka

,

H.

Sugama

, and

Y.

Idomura

, “

Neoclassical transport simulations with an improved model collision operator

,”

Phys. Plasmas

28

,

064501

(

2021

).

https://doi.org/10.1063/5.0047204

Google Scholar

Crossref

21.

B. J.

Frei

,

R.

Jorge

, and

P.

Ricci

, “

A gyrokinetic model for the plasma periphery of tokamak devices

,”

J. Plasma Phys.

86

,

905860205

(

2020

).

https://doi.org/10.1017/S0022377820000100

Google Scholar

Crossref

22.

F.

Jenko

,

W.

Dorland

,

M.

Kotschenreuther

, and

B. N.

Rogers

, “

Electron temperature gradient driven turbulence

,”

Phys. Plasmas

7

,

1904

(

2000

).

https://doi.org/10.1063/1.874014

Google Scholar

Crossref

23.

F. J.

Casson

,

C.

Angioni

,

E. A.

Belli

,

R.

Bilato

,

P.

Mantica

,

T.

Odstrcil

,

T.

Pütterich

,

M.

Valisa

,

L.

Garzotti

,

C.

Giroud

et al, “

Theoretical description of heavy impurity transport and its application to the modelling of tungsten in jet and ASDEX upgrade

,”

Plasma Phys. Controlled Fusion

57

,

014031

(

2014

).

https://doi.org/10.1088/0741-3335/57/1/014031

Google Scholar

Crossref

24.

M.

Landreman

and

D. R.

Ernst

, “

New velocity-space discretization for continuum kinetic calculations and Fokker–Planck collisions

,”

J. Comput. Phys.

243

,

130

(

2013

).

https://doi.org/10.1016/j.jcp.2013.02.041

Google Scholar

Crossref

25.

P.

Helander

and

D. J.

Sigmar

,

Collisional Transport in Magnetized Plasmas

(

Cambridge University Press

,

2005

), Vol.

4

.

Google Scholar

26.

R.

Jorge

,

B. J.

Frei

, and

P.

Ricci

, “

Nonlinear gyrokinetic Coulomb collision operator

,”

J. Plasma Phys.

85

,

905850604

(

2019

).

https://doi.org/10.1017/S0022377819000734

Google Scholar

Crossref

27.

R. F.

Snider

,

Irreducible Cartesian Tensors

, De Gruyter Studies in Mathematical Physics (

De Gruyter

,

2017

), No. 1.

Google Scholar

Crossref

28.

A. J.

Brizard

and

T.

Hahm

, “

Foundations of nonlinear gyrokinetic theory

,”

Rev. Mod. Phys.

79

,

421

(

2007

).

https://doi.org/10.1103/RevModPhys.79.421

Google Scholar

Crossref

29.

R.

Jorge

,

P.

Ricci

, and

N. F.

Loureiro

, “

A drift-kinetic analytical model for scrape-off layer plasma dynamics at arbitrary collisionality

,”

J. Plasma Phys.

83

,

905830606

(

2017

).

https://doi.org/10.1017/S002237781700085X

Google Scholar

Crossref

30.

I. S.

Gradshteyn

and

I. M.

Ryzhik

,

Table of Integrals, Series, and Products

(

Academic Press

,

2014

).

Google Scholar

31.

D. P.

Fulton

,

Z.

Lin

,

I.

Holod

, and

Y.

Xiao

, “

Microturbulence in DIII-D tokamak pedestal. I. Electrostatic instabilities

,”

Phys. Plasmas

21

,

042110

(

2014

).

https://doi.org/10.1063/1.4871387

Google Scholar

Crossref

32.

M.

Kotschenreuther

,

D. R.

Hatch

,

S.

Mahajan

,

P.

Valanju

,

L.

Zheng

, and

X.

Liu

, “

Pedestal transport in H-mode plasmas for fusion gain

,”

Nucl. Fusion

57

,

064001

(

2017

).

https://doi.org/10.1088/1741-4326/aa6416

Google Scholar

Crossref

33.

M. J.

Pueschel

,

D. R.

Hatch

,

D. R.

Ernst

,

W.

Guttenfelder

,

P. W.

Terry

,

J.

Citrin

, and

J. W.

Connor

, “

On microinstabilities and turbulence in steep-gradient regions of fusion devices

,”

Plasma Phys. Controlled Fusion

61

,

034002

(

2019

).

https://doi.org/10.1088/1361-6587/aaf8c1

Google Scholar

Crossref

34.

H.-S.

Xie

and

B.

Li

, “

Global theory to understand toroidal drift waves in steep gradient

,”

Phys. Plasmas

23

,

082513

(

2016

).

https://doi.org/10.1063/1.4960964

Google Scholar

Crossref

35.

M.

Han

,

Z.-X.

Wang

,

J.

Dong

, and

H.

Du

, “

Multiple ion temperature gradient driven modes in transport barriers

,”

Nucl. Fusion

57

,

046019

(

2017

).

https://doi.org/10.1088/1741-4326/aa5d02

Google Scholar

Crossref

36.

D. M.

Thomas

,

A. W.

Leonard

,

T. H.

Osborne

,

R. J.

Groebner

,

W. P.

West

, and

K. H.

Burrell

, “

The effect of plasma collisionality on pedestal current density formation in DIII-D

,”

Plasma Phys. Controlled Fusion

48

,

A183

(

2006

).

https://doi.org/10.1088/0741-3335/48/5A/S17

Google Scholar

Crossref

37.

X.

Lapillonne

,

S.

Brunner

,

T.

Dannert

,

S.

Jolliet

,

A.

Marinoni

,

L.

Villard

,

T.

Görler

,

F.

Jenko

, and

F.

Merz

, “

Clarifications to the limitations of the s-α equilibrium model for gyrokinetic computations of turbulence

,”

Phys. Plasmas

16

,

032308

(

2009

).

https://doi.org/10.1063/1.3096710

Google Scholar

Crossref

38.

B.

Coppi

,

S.

Migliuolo

, and

Y.-K.

Pu

, “

Candidate mode for electron thermal energy transport in multi-keV plasmas

,”

Phys. Fluids B

2

,

2322

(

1990

).

https://doi.org/10.1063/1.859497

Google Scholar

Crossref

39.

D. R.

Ernst

,

K.

Zeller

,

N.

Basse

,

L.

Lin

,

M.

Porkolab

,

W.

Dorland

, and

A.

Long

, “

New developments in trapped electron mode turbulence

,”

APS Div. Plasma Phys. Meet. Abstr.

50

,

235

(

2005

).

Google Scholar

40.

E.

Wang

,

X.

Xu

,

J.

Candy

,

R. J.

Groebner

,

P. B.

Snyder

,

Y.

Chen

,

S. E.

Parker

,

W.

Wan

,

G.

Lu

, and

J. Q.

Dong

, “

Linear gyrokinetic analysis of a DIII-D H-mode pedestal near the ideal ballooning threshold

,”

Nucl. Fusion

52

,

103015

(

2012

).

https://doi.org/10.1088/0029-5515/52/10/103015

Google Scholar

Crossref

41.

P. H.

Diamond

,

S. I.

Itoh

,

K.

Itoh

, and

T. S.

Hahm

, “

Zonal flows in plasma—A review

,”

Plasma Phys. Controlled Fusion

47

,

R35

(

2005

).

https://doi.org/10.1088/0741-3335/47/5/R01

Google Scholar

Crossref

42.

M. N.

Rosenbluth

and

F. L.

Hinton

, “

Poloidal flow driven by ion-temperature-gradient turbulence in tokamaks

,”

Phys. Rev. Lett.

80

,

724

(

1998

).

https://doi.org/10.1103/PhysRevLett.80.724

Google Scholar

Crossref

43.

F. L.

Hinton

and

M. N.

Rosenbluth

, “

Dynamics of axisymmetric (E × B) and poloidal flows in tokamaks

,”

Plasma Phys. Controlled Fusion

41

,

A653

(

1999

).

https://doi.org/10.1088/0741-3335/41/3A/059

Google Scholar

Crossref

44.

Y.

Xiao

,

P. J.

Catto

, and

K.

Molvig

, “

Collisional damping for ion temperature gradient mode driven zonal flow

,”

Phys. Plasmas

14

,

032302

(

2007

).

https://doi.org/10.1063/1.2536297

Google Scholar

Crossref

45.

L.

Spitzer

and

R.

Härm

, “

Transport phenomena in a completely ionized gas

,”

Phys. Rev.

89

,

977

(

1953

).

https://doi.org/10.1103/PhysRev.89.977

Google Scholar

Crossref

46.

M.

Giacomin

,

A.

Pau

,

P.

Ricci

,

O.

Sauter

,

T.

Eich

,

J.

Contributors

et al, “

First-principles density limit scaling in tokamaks based on edge turbulent transport and implications for ITER

,”

Phys. Rev. Lett.

128

,

185003

(

2022

).

https://doi.org/10.1103/PhysRevLett.128.185003

Google Scholar

Crossref

PubMed

47.

B.

Li

and

D. R.

Ernst

, “

Gyrokinetic Fokker-Planck collision operator

,”

Phys. Rev. Lett.

106

,

195002

(

2011

).

https://doi.org/10.1103/PhysRevLett.106.195002

Google Scholar

Crossref

PubMed

48.

P.

Crandall

,

D.

Jarema

,

H.

Doerk

,

Q.

Pan

,

G.

Merlo

,

T.

Görler

,

A. B.

Navarro

,

D.

Told

,

M.

Maurer

, and

F.

Jenko

, “

Multi-species collisions for delta-f gyrokinetic simulations: Implementation and verification with GENE

,”

Comput. Phys. Commun.

255

,

107360

(

2020

).

https://doi.org/10.1016/j.cpc.2020.107360

Google Scholar

Crossref

2022

Author(s)

Numerical implementation of the improved Sugama collision operator using a moment approach

I. INTRODUCTION

II. IMPROVED SUGAMA COLLISION OPERATOR

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

A. Spherical harmonic expansion of the improved Sugama collision operator

B. Gyrokinetic improved Sugama collision operator

C. Drift-kinetic improved Sugama collision operator

IV. BRAGINKSII MATRICES

A. Braginskii matrix of the Coulomb collision operator

B. Braginskii matrix of the original Sugama collision operator

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

A. Expansion of the GK IS collision operator

B. Expansion of the DK IS collision operator

VI. NUMERICAL TESTS AND COMPARISON BETWEEN COLLISION OPERATORS

A. Numerical implementation

B. Trapped electron mode in steep pressure gradient conditions

C. Collisional zonal flow damping

D. Spitzer electrical conductivity

VII. CONCLUSION

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: LINEARIZED COULOMB AND ORIGINAL SUGAMA COLLISION OPERATORS

APPENDIX B: GK FORMULATION OF COLLISION OPERATORS

APPENDIX C: LOWEST ORDER GYRO-MOMENT ANALYTICAL EXPRESSIONS

References

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Numerical implementation of the improved Sugama collision operator using a moment approach

I. INTRODUCTION

II. IMPROVED SUGAMA COLLISION OPERATOR

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

A. Spherical harmonic expansion of the improved Sugama collision operator

B. Gyrokinetic improved Sugama collision operator

C. Drift-kinetic improved Sugama collision operator

IV. BRAGINKSII MATRICES

A. Braginskii matrix of the Coulomb collision operator

B. Braginskii matrix of the original Sugama collision operator

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

A. Expansion of the GK IS collision operator

B. Expansion of the DK IS collision operator

VI. NUMERICAL TESTS AND COMPARISON BETWEEN COLLISION OPERATORS

A. Numerical implementation

B. Trapped electron mode in steep pressure gradient conditions

C. Collisional zonal flow damping

D. Spitzer electrical conductivity

VII. CONCLUSION

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: LINEARIZED COULOMB AND ORIGINAL SUGAMA COLLISION OPERATORS

APPENDIX B: GK FORMULATION OF COLLISION OPERATORS

APPENDIX C: LOWEST ORDER GYRO-MOMENT ANALYTICAL EXPRESSIONS

References

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only