2.1 Single particle motion

To derive a convenient equation of motion in the presence of a strong magnetic field $\boldsymbol{B}$ , we start with the Hamiltonian of a charged particle of species $a$ (Jackson Reference Jackson1998),

and its associated Lagrangian,

where $\boldsymbol{p}=m_{a}\boldsymbol{v}+q_{a}\,\boldsymbol{A}$ , $\boldsymbol{v}$ is the particle velocity, $\boldsymbol{A}$ is the magnetic vector potential, $\unicode[STIX]{x1D719}$ is the electrostatic potential, $m_{a}$ is the mass of the particle and $q_{a}$ its charge.

We now perform a coordinate transformation from the phase-space coordinates $\boldsymbol{z}=(\boldsymbol{x},\boldsymbol{v})$ to the guiding-centre coordinates $\boldsymbol{Z}=(\boldsymbol{R},v_{\Vert },\unicode[STIX]{x1D707},\unicode[STIX]{x1D703})$ by writing the particle velocity as (see, e.g. Littlejohn (Reference Littlejohn1983))

with $\boldsymbol{U}=\boldsymbol{v}_{E}(\boldsymbol{R})+v_{\Vert }\boldsymbol{b}(\boldsymbol{R})$ , $v_{\Vert }=\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{b}$ , $\boldsymbol{b}=\boldsymbol{B}/B$ and $\boldsymbol{v}_{E}=\boldsymbol{E}\times \boldsymbol{B}/B^{2}$ . The gyroangle $\unicode[STIX]{x1D703}=\tan ^{-1}[(\boldsymbol{v}-\boldsymbol{U})\boldsymbol{\cdot }\boldsymbol{e}_{2}/(\boldsymbol{v}-\boldsymbol{U})\boldsymbol{\cdot }\boldsymbol{e}_{1}]$ is introduced by defining the right-handed coordinate set $(\boldsymbol{e}_{1},\boldsymbol{e}_{2},\boldsymbol{b})$ , such that $\boldsymbol{c}=-\boldsymbol{a}\times \boldsymbol{b}=\boldsymbol{a}^{\prime }(\unicode[STIX]{x1D703})$ , with $\boldsymbol{a}=\cos \unicode[STIX]{x1D703}\boldsymbol{e}_{1}+\sin \unicode[STIX]{x1D703}\boldsymbol{e}_{2}$ . The decomposition in (2.3) allows us to isolate the high-frequency gyromotion, contained in the $v_{\bot }^{\prime }\boldsymbol{c}$ term, from the dominant guiding-centre velocity $\boldsymbol{U}$ . The adiabatic moment $\unicode[STIX]{x1D707}$ is defined as

whereas the guiding-centre position is

with $\unicode[STIX]{x1D70C}_{a}=\sqrt{2m_{a}\unicode[STIX]{x1D707}/(q_{a}^{2}B)}$ the Larmor radius. Incidentally, for the case of weakly varying magnetic fields, equation (2.5) describes the circular motion of a particle around its guiding centre $\boldsymbol{R}$ with radius $\unicode[STIX]{x1D70C}_{a}$ , i.e. $(\boldsymbol{x}-\boldsymbol{R})^{2}=\unicode[STIX]{x1D70C}_{a}^{2}$ .

As our goal is to develop a model that describes turbulent fluctuations occurring on a spatial scale longer than the sound Larmor radius $\unicode[STIX]{x1D70C}_{s}$ , and a time scale larger than the gyromotion one, we keep terms in the Lagrangian up to $O(\unicode[STIX]{x1D716})$ and order $T_{i}\lesssim T_{e}$ , which implies

We therefore expand the electromagnetic fields around $\boldsymbol{R}$ , to first order in $\unicode[STIX]{x1D716}$ , i.e.

and similarly for $\boldsymbol{A}$ . In the following, if not specified, the electromagnetic fields and potentials are evaluated at the guiding-centre position $\boldsymbol{R}$ , and we denote $\unicode[STIX]{x1D735}=\unicode[STIX]{x1D735}_{\boldsymbol{R}}$ . In addition, to take advantage of the difference between the turbulent and gyromotion time scales, we use the gyroaveraged Lagrangian $\langle L_{a}\rangle$ to evaluate the plasma particle motion, where the gyroaveraging operator $\langle \unicode[STIX]{x1D712}\rangle$ acting on a quantity $\unicode[STIX]{x1D712}(\unicode[STIX]{x1D703})$ is defined as

which is performed at fixed position $\boldsymbol{x}$ , as opposed to the gyrokinetic equation that can be obtained by gyroaveraging with $\boldsymbol{R}$ fixed (Hazeltine & Meiss Reference Hazeltine and Meiss2003).

To evaluate $\langle L_{a}\rangle$ we note that, with the expansion for $\unicode[STIX]{x1D719}$ and $\boldsymbol{A}$ , the Lagrangian in (2.2) can be expressed as $L_{a}=L_{0a}+L_{1a}+\tilde{L}_{a}$ where $L_{0a}$ is gyroangle independent,

$L_{1a}$ is proportional to $\unicode[STIX]{x1D70C}_{a}^{2}$ (and hence to $\unicode[STIX]{x1D707}$ ),

and the $\tilde{L}_{a}$ contribution contains the terms linearly proportional to $\cos \unicode[STIX]{x1D703}$ or $\sin \unicode[STIX]{x1D703}$ (Cary & Brizard Reference Cary and Brizard2009) which are not present in $\langle L_{a}\rangle$ as $\langle \tilde{L}_{a}\rangle =0$ .

We note that $\langle L_{1a}\rangle$ can be simplified since $\langle (\boldsymbol{a}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{A}\boldsymbol{\cdot }\boldsymbol{c}\rangle =-\boldsymbol{b}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\times \boldsymbol{A})/2$ , and $\langle (\boldsymbol{a}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{A}\boldsymbol{\cdot }\boldsymbol{a}\rangle =\unicode[STIX]{x1D735}_{\bot }\boldsymbol{\cdot }\boldsymbol{A}/2$ . Subtracting the total derivative $-q_{a}\,\text{d}/\text{d}t(\unicode[STIX]{x1D70C}_{a}^{2}\unicode[STIX]{x1D735}_{\bot }\boldsymbol{A})/4$ from $\langle L_{a}\rangle$ , which does not alter the resulting equations of motion, we redefine the gyroaveraged Lagrangian as

We now order the terms appearing in $\langle L_{a}\rangle$ . As imposed by the Bohm sheath conditions (Stangeby Reference Stangeby2000), both electrons and ions stream along the field lines with parallel velocities comparable to the sound speed $c_{s}=\sqrt{T_{e}/m_{i}}$ in the SOL. The Bohm boundary conditions at the sheath also set the electrostatic potential $e\unicode[STIX]{x1D719}\sim \unicode[STIX]{x1D6EC}T_{e}$ across the SOL, where $\unicode[STIX]{x1D6EC}=\log \sqrt{m_{i}/(m_{e}2\unicode[STIX]{x03C0})}\simeq 3$ . Therefore, we keep the $m_{a}v_{E}^{2}/2$ term in the Lagrangian in (2.11), as to take into account the presence of the numerically large factor $\unicode[STIX]{x1D6EC}^{2}$ in $v_{E}^{2}\sim \unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D6EC}^{2}c_{s}^{2}$ .

By neglecting the higher-order terms in (2.11), i.e. $-(\unicode[STIX]{x1D70C}_{a}^{2}/4)\text{d}[\unicode[STIX]{x1D735}_{\bot }\boldsymbol{\cdot }(q_{a}\,\boldsymbol{A})]/\text{d}t$ , the expression for the gyroaveraged Lagrangian describing SOL single particle dynamics, up to $O(\unicode[STIX]{x1D716})$ , can be written as

where $q_{a}\unicode[STIX]{x1D719}^{\ast }=q_{a}\unicode[STIX]{x1D719}+m_{a}v_{E}^{2}/2+\unicode[STIX]{x1D707}B$ , and $q_{a}\boldsymbol{A}^{\ast }=q_{a}\boldsymbol{A}+m_{a}v_{\Vert }\boldsymbol{b}+m_{a}\boldsymbol{v}_{E}$ . The Euler–Lagrange equations applied to the Lagrangian in (2.12) for the coordinates $\unicode[STIX]{x1D703}$ , $v_{\Vert }$ , and $\unicode[STIX]{x1D707}$ , yield, respectively, $\dot{\unicode[STIX]{x1D707}}=0$ , $v_{\Vert }=\boldsymbol{b}\boldsymbol{\cdot }\dot{\boldsymbol{R}}$ , and $\dot{\unicode[STIX]{x1D703}}=\unicode[STIX]{x1D6FA}_{a}$ . For the $\boldsymbol{R}$ coordinate, we obtain

where the relation $(\unicode[STIX]{x1D735}\boldsymbol{A}-(\unicode[STIX]{x1D735}\boldsymbol{A})^{\text{T}})\boldsymbol{\cdot }\dot{\boldsymbol{R}}=\dot{\boldsymbol{R}}\times (\unicode[STIX]{x1D735}\times \boldsymbol{A})$ has been used, and we defined $\boldsymbol{E}^{\ast }=-\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}^{\ast }-\unicode[STIX]{x2202}_{t}\boldsymbol{A}^{\ast }$ , and $\boldsymbol{B}^{\ast }=\unicode[STIX]{x1D735}\times \boldsymbol{A}^{\ast }$ , with the parallel component of $\boldsymbol{B}^{\ast }$ given by

By projecting (2.13) along $\boldsymbol{B}^{\ast }$ , we derive $m\dot{v}_{\Vert }B_{\Vert }^{\ast }=e\boldsymbol{E}^{\ast }\boldsymbol{\cdot }\boldsymbol{B}^{\ast }$ , while crossing with $\boldsymbol{b}$ yields the guiding-centre velocity $\dot{\boldsymbol{R}}B_{\Vert }^{\ast }=v_{\Vert }\boldsymbol{B}^{\ast }+\boldsymbol{E}^{\ast }\times \boldsymbol{B}/B$ . Using the expressions for the fields $\boldsymbol{E}^{\ast }$ and $\boldsymbol{B}^{\ast }$ , we obtain

and

In (2.15), (2.16), in addition to the time derivatives of the phase-space coordinates $\dot{\boldsymbol{R}},\dot{v}_{\Vert }$ , that only have an explicit time dependence, we define the total derivative $\text{d}/\text{d}t$ of a field $\unicode[STIX]{x1D719}(\boldsymbol{R},t)$ that has an explicit time and $\boldsymbol{R}$ dependence as $\text{d}\unicode[STIX]{x1D719}/\text{d}t\equiv \unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D719}+\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}$ . The ${\mathcal{A}}$ term represents the higher-order nonlinear terms in $\dot{v}_{\Vert }$ that ensure phase-space conservation properties (Cary & Brizard Reference Cary and Brizard2009), and it is given by

The guiding-centre equations of motion (2.15) and (2.16) satisfy the energy, $E_{gc}=q_{a}\unicode[STIX]{x1D719}^{\ast }+m_{a}v_{\Vert }^{2}/2$ (Cary & Brizard Reference Cary and Brizard2009), and momentum, $\boldsymbol{P}_{gc}=e\boldsymbol{A}^{\ast }$ (Cary & Brizard Reference Cary and Brizard2009), conservation laws, i.e.

and

In addition, we note that using (2.15), (2.16) and Maxwell’s equations, a conservation equation for $B_{\Vert }^{\ast }$ can be derived:

2.2 The guiding-centre Boltzmann equation

The Boltzmann equation for the evolution of the distribution function $f_{a}(\boldsymbol{x},\boldsymbol{v})$ of the particles in $(\boldsymbol{x},\boldsymbol{v})$ coordinates is

where $C(f_{a})=\sum _{b}C(f_{a},f_{b})=\sum _{b}C_{ab}$ is the collision operator. Because $f_{a}$ can significantly deviate from a Maxwellian distribution function in the SOL (Battaglia et al. Reference Battaglia, Burrell, Chang, Ku, Degrassie and Grierson2014), we consider the bilinear Coulomb operator $C_{ab}$ (Balescu Reference Balescu1988), to model collisions between particles of species $a$ and $b$

with

and

the Rosenbluth potentials satisfying $\unicode[STIX]{x1D6FB}_{v}^{2}G_{b}=H_{b}$ . In (2.22) we introduced $L_{ab}=q_{a}^{2}q_{b}^{2}\unicode[STIX]{x1D706}/(4\unicode[STIX]{x03C0}\unicode[STIX]{x1D716}_{0}^{2}m_{a}^{2})=\unicode[STIX]{x1D708}_{ab}v_{tha}^{3}/n_{b}$ , where $\unicode[STIX]{x1D706}$ is the Coulomb logarithm, $\unicode[STIX]{x1D708}_{ab}$ the collision frequency between species $a$ and $b$ and $v_{tha}^{2}=2T_{a}/m_{a}$ .

Taking advantage of the small electron to ion mass ratio, the collision operator between unlike species can be simplified (see, e.g. Balescu (Reference Balescu1988), Helander & Sigmar (Reference Helander and Sigmar2005)). The electron–ion collisions are modelled by using the operator $C_{ei}(f_{e})=C_{ei}^{0}+C_{ei}^{1}$ , where $C_{ei}^{0}$ is the Lorentz pitch-angle scattering operator

and $C_{ei}^{1}$ the momentum-conserving term

with $\boldsymbol{c}_{e}=(\boldsymbol{v}-\boldsymbol{u}_{e})/v_{the}$ the normalized difference between the electron particle $\boldsymbol{v}$ and fluid $\boldsymbol{u}_{e}$ velocities, and $\boldsymbol{u}_{i}$ the ion fluid velocity.

Ion–electron collisions are modelled with the operator

where $\boldsymbol{R}_{ei}=\int m_{e}\boldsymbol{v}C_{ei}\,\text{d}\boldsymbol{v}$ is the electron–ion friction force.

We take advantage of (1.2) to order the electron collision frequency $\unicode[STIX]{x1D708}_{e}$ and the ion collision frequency $\unicode[STIX]{x1D708}_{i}$ as

where we used the relation $\unicode[STIX]{x1D708}_{i}\sim \sqrt{m_{e}/m_{i}}(T_{e}/T_{i})^{3/2}\unicode[STIX]{x1D708}_{e}$ . The orderings in (2.7) and (2.28) yield the lower bound in (1.3) for the ion to electron temperature ratio. We now express the particle distribution function $f_{a}$ in terms of the guiding-centre coordinates by defining $F_{a}$ , a function of guiding-centre coordinates, as

Using the chain rule to rewrite (2.21) in guiding-centre coordinates, we obtain

where $\dot{\boldsymbol{R}}$ and $\dot{v}_{\Vert }$ are given by (2.15) and (2.16) respectively, $\dot{\unicode[STIX]{x1D703}}=\unicode[STIX]{x1D6FA}_{a}$ , and $\dot{\unicode[STIX]{x1D707}}=0$ . Equation (2.30) can be simplified by applying the gyroaveraging operator. This results in the drift-kinetic equation

We now write (2.31) in a form useful to take gyrofluid moments of the form $\int \langle F_{a}\rangle B\,\text{d}v_{\Vert }\,\text{d}\unicode[STIX]{x1D707}\,\text{d}\unicode[STIX]{x1D703}$ (see § 4). Using the conservation law in (2.20) for $B_{\Vert }^{\ast }$ , we can write the guiding-centre Boltzmann equation in conservative form as

Moreover, in order to relate the gyrofluid moments $\int \langle F_{a}\rangle B\,\text{d}v_{\Vert }\,\text{d}\unicode[STIX]{x1D707}\,\text{d}\unicode[STIX]{x1D703}$ with the usual fluid moments $\int f_{a}\,\text{d}^{3}v$ , we estimate the order of magnitude of the gyrophase dependent part of the distribution function $\tilde{F}_{a}=F_{a}-\langle F_{a}\rangle$ where $\langle F_{a}\rangle$ obeys (2.31). The equation for the evolution of $\tilde{F}_{a}$ is obtained by subtracting (2.31) from the Boltzmann equation, (2.30), that is

Using the orderings in (1.2) and (2.28), as well as $\unicode[STIX]{x2202}_{t}\sim \dot{\boldsymbol{R}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\sim \dot{v}_{\Vert }\unicode[STIX]{x2202}_{v_{\Vert }}\sim \unicode[STIX]{x1D716}\unicode[STIX]{x1D6FA}_{i}$ and $\unicode[STIX]{x1D6FA}_{a}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}\sim \unicode[STIX]{x1D6FA}_{a}$ , the comparison of the leading-order term on the left-hand side of (2.33) with the right-hand side of the same equation imply the following ordering for $\tilde{F}_{e}$

and $\tilde{F}_{i}$

To evaluate the leading-order term of $\tilde{F}_{a}$ , we expand the collision operator $C(F_{a})=C_{0}(\langle F_{a}\rangle )+\unicode[STIX]{x1D716}C_{1}(F_{a})+\cdots \,$ , such that

The relation in (2.36) can be further simplified by expanding the $\unicode[STIX]{x1D703}$ dependence of $F_{a}$ in Fourier harmonics,

so that for $m=0$ we have $\langle F_{a}\rangle =F_{0a}$ , and similarly for $C_{0}(\langle F_{a}\rangle )$

We can then write (2.36) as

for $m\not =0$ .

3.1 Guiding-centre moment expansion of $\langle F_{a}\rangle$

To take advantage of the anisotropy introduced by a strong magnetic field, and efficiently treat the left-hand side of (2.32) where the parallel and perpendicular directions appear decoupled, we express $\langle F_{a}\rangle$ by using a Hermite polynomial basis expansion for the parallel velocity coordinate (Grad Reference Grad1949; Armstrong Reference Armstrong1967; Grant & Feix Reference Grant and Feix1967; Ng, Bhattacharjee & Skiff Reference Ng, Bhattacharjee and Skiff1999; Zocco & Schekochihin Reference Zocco and Schekochihin2011; Loureiro, Schekochihin & Zocco Reference Loureiro, Schekochihin and Zocco2013; Parker & Dellar Reference Parker and Dellar2015; Schekochihin et al. Reference Schekochihin, Parker, Highcock, Dellar, Dorland and Hammett2016; Tassi Reference Tassi2016) and a Laguerre polynomial basis for the perpendicular velocity coordinate (Omotani et al. Reference Omotani, Dudson, Havlíčková and Umansky2015; Zocco et al. Reference Zocco, Loureiro, Dickinson, Numata and Roach2015; Mandell et al. Reference Mandell, Dorland and Landreman2017). More precisely, we use the following expansion

where the physicists’ Hermite polynomials $H_{p}$ of order $p$ are defined by Abramowicz et al. (Reference Abramowicz, Czerny, Lasota and Szuszkiewicz1988)

and normalized via

and the Laguerre polynomials $L_{j}$ of order $j$ are defined by Abramowicz et al. (Reference Abramowicz, Czerny, Lasota and Szuszkiewicz1988)

which are orthonormal with respect to the weight $\text{e}^{-x}$

Because of the orthogonality of the Hermite–Laguerre basis, the coefficients $N_{a}^{pj}$ of the expansion in (3.1) are

and correspond to the guiding-centre moments of $\langle F_{a}\rangle$ .

In (3.1), the shifted bi-Maxwellian is introduced

where $s_{\Vert a}$ and $s_{\bot a}$ are the normalized parallel and perpendicular shifted velocities respectively, defined by

and

which provide an efficient representation of the distribution function in both the weak $(u_{\Vert a}\ll v_{tha})$ and strong flow ( $u_{\Vert a}\sim v_{tha}$ ) regimes.

The guiding-centre density $N_{a}$ , appearing in (3.7), the guiding-centre fluid velocity $u_{\Vert a}$ , in (3.8) and the guiding-centre parallel $T_{\Vert a}=P_{\Vert a}/N_{a}$ and perpendicular $T_{\bot a}=P_{\bot a}/N_{a}$ temperatures in (3.8) , (3.9) are defined as $N_{a}=\Vert 1\Vert _{a}$ , $N_{a}u_{\Vert a}=\Vert v_{\Vert }\Vert _{a}$ , $P_{\Vert a}=m_{a}\Vert (v_{\Vert }-u_{\Vert a})^{2}\Vert _{a}$ and $P_{\bot a}=\Vert \unicode[STIX]{x1D707}B\Vert _{a}$ , where

The definitions of $N_{a}$ , $u_{\Vert a}$ , $P_{\Vert a}$ and $P_{\bot a}$ imply that $N_{a}^{00}=1,~N_{a}^{10}=0,~N_{a}^{20}=0,~N_{a}^{01}=0$ respectively. Later, we will consider the parallel and perpendicular heat flux, defined as

which are related to the coefficients $N_{a}^{30},N_{a}^{11}$ by

3.2 Fluid moment expansion of the collision operator

A polynomial expansion of the collision operators in (2.22) was carried out in Ji & Held (Reference Ji and Held2006), and later extended to effectively take into account finite fluid velocity and unlike-species collisions in Ji & Held (Reference Ji and Held2008). This allowed expressing $C_{ab}$ as products of fluid moments of $f_{a}$ and $f_{b}$ . We summarize here the main steps of Ji & Held (Reference Ji and Held2006, Reference Ji and Held2008).

Similarly to (3.1), the particle distribution function $f_{a}$ is expanded as

where $f_{aM}=n_{a}\exp (-c_{a}^{2})/(\unicode[STIX]{x03C0}^{3/2}v_{tha}^{3})$ is a shifted Maxwell–Boltzmann distribution function, and $\boldsymbol{c}_{a}$ the shifted velocity defined as $\boldsymbol{c}_{a}=(\boldsymbol{v}-\boldsymbol{U}_{a})/v_{tha}$ , with $\boldsymbol{U}_{a}$ the fluid velocity. The fluid variables $n_{a},\boldsymbol{U}_{a},$ and $T_{a}$ are defined as the usual moments of the particle distribution function $f_{a}$ , i.e. $n_{a}=\int f_{a}\,\text{d}^{3}\boldsymbol{v},n_{a}\boldsymbol{U}_{a}=\int f_{a}\boldsymbol{v}\,\text{d}^{3}v,n_{a}T_{a}=\int mf_{a}(\boldsymbol{v}-\boldsymbol{u}_{a})^{2}\,\text{d}^{3}v/3$ .

The tensors $\unicode[STIX]{x1D64B}_{a}^{lk}(\boldsymbol{c}_{a})=\unicode[STIX]{x1D64B}^{l}(\boldsymbol{c}_{a})L_{k}^{l+1/2}(c_{a}^{2})$ constitute an orthogonal basis, where $\unicode[STIX]{x1D64B}^{l}(\boldsymbol{c}_{a})$ is the symmetric and traceless tensor

with $\unicode[STIX]{x1D644}$ denoting the identity matrix, $\{\unicode[STIX]{x1D63C}^{i}\}$ denoting the symmetrization of the tensor $\unicode[STIX]{x1D63C}^{i}$ , $\lfloor l/2\rfloor$ denoting the largest integer less than or equal to $l/2$ , and the coefficients $d_{i}^{l}$ and $S_{i}^{l}$ defined by

and

The tensor $\unicode[STIX]{x1D64B}^{l}(\boldsymbol{c}_{a})$ is normalized via

with $\unicode[STIX]{x1D70E}_{l}=l!/[2^{l}(l+1/2)!]$ . We note that the tensor $\unicode[STIX]{x1D63C}^{i}$ is formed by $i$ multiplications of the $\unicode[STIX]{x1D63C}$ elements (e.g. if $\unicode[STIX]{x1D63C}$ is a rank-2 tensor, $\unicode[STIX]{x1D63C}^{3}\equiv \unicode[STIX]{x1D63C}\unicode[STIX]{x1D63C}\unicode[STIX]{x1D63C}$ , which in index notation can be written as $(\unicode[STIX]{x1D608}^{3})_{ijlkmn}=\unicode[STIX]{x1D608}_{ij}\unicode[STIX]{x1D608}_{lk}\unicode[STIX]{x1D608}_{mn}$ ).

In the expansion in (3.13), $L_{k}^{l+1/2}(x)$ are the associated Laguerre polynomials

normalized via

with $\unicode[STIX]{x1D706}_{k}^{l}=(l+k+1/2)!/k!$ and $L_{km}^{l}=[(-1)^{m}(l+k+1/2)!]/[(k-m)!(l+m+1/2)!m!]$ . The $\unicode[STIX]{x1D70E}_{k}^{l}=\unicode[STIX]{x1D70E}_{l}\unicode[STIX]{x1D706}_{k}^{l}$ term is a normalization factor from the orthogonality relations in (3.17), (3.19).

Finally, the coefficients of the expansion in (3.13) $\boldsymbol{M}_{a}^{lk}$ are

which correspond to the moments of $f_{a}$ due to the orthogonality relations in (3.17), (3.19).

By using the expansion in (3.13) in the collision operator in (2.22), a closed form for $C_{ab}$ in terms of products of $\boldsymbol{M}_{a}^{lk}$ can be obtained. For like-species collisions it reads

with

where $\hat{\boldsymbol{c}_{a}}=\boldsymbol{c}_{a}/c_{a}$ , $\boldsymbol{\cdot }^{n}$ is the $n$ -fold inner product (e.g. for the matrix $\unicode[STIX]{x1D63C}=\unicode[STIX]{x1D608}_{ij}$ , $(\unicode[STIX]{x1D63C}\boldsymbol{\cdot }^{1}\unicode[STIX]{x1D63C})_{ij}=\sum _{k}\unicode[STIX]{x1D608}_{ki}\unicode[STIX]{x1D608}_{kj}$ ) and $\overline{\unicode[STIX]{x1D63C}}$ the traceless symmetrization of $\unicode[STIX]{x1D63C}$ (e.g. $\overline{\unicode[STIX]{x1D63C}}=(\unicode[STIX]{x1D608}_{ij}+\unicode[STIX]{x1D608}_{ji})/2-\unicode[STIX]{x1D6FF}_{ij}\sum _{k}\unicode[STIX]{x1D608}_{kk}/3$ ). We refer the reader to Ji & Held (Reference Ji and Held2009) for the explicit form of the $\unicode[STIX]{x1D708}_{\ast abu}^{lm,nr}$ coefficients. The expansion of the unlike-species collisions is reported in Ji & Held (Reference Ji and Held2008).

3.3 Guiding-centre moment expansion of the collision operators

In order to apply the gyroaveraging operator to the like-species collision operator $C_{aa}$ in (3.21), we expand the fluid moments as $\boldsymbol{M}_{a}^{lk}=\boldsymbol{M}_{a0}^{lk}+\unicode[STIX]{x1D716}\boldsymbol{M}_{a1}^{lk}+\cdots \,,$ aiming at representing the collision operator up to $O(\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D708}}\unicode[STIX]{x1D716})$ . An analytical expression for the leading order $\boldsymbol{M}_{a0}^{lk}$ in terms of guiding-centre moments $N_{a}^{pj}$ can be obtained as follows. By splitting $f_{a}=\langle f_{a}\rangle +\tilde{f}_{a}$ when evaluating the fluid moments $\boldsymbol{M}_{a}^{lk}$ according to (3.20), we obtain

where the Dirac delta function was introduced to convert the velocity integral into an $(\boldsymbol{x},\boldsymbol{v})$ integral that encompasses the full phase space. Since the volume element in phase space can be written as $\text{d}^{3}\boldsymbol{x}\,\text{d}^{3}\,\boldsymbol{v}=(B_{\Vert }^{\ast }/m)\,\text{d}\boldsymbol{R}\,\text{d}v_{\Vert }\,\text{d}\unicode[STIX]{x1D707}\,\text{d}\unicode[STIX]{x1D703}$ (Cary & Brizard Reference Cary and Brizard2009), and defining $\boldsymbol{x}^{\prime }=\boldsymbol{R}+\unicode[STIX]{x1D70C}_{a}\boldsymbol{a}$ , we can write the fluid moments in (3.23) as

where $\langle f_{a}\rangle$ and $\tilde{f}_{a}$ in (3.23) are written in terms of guiding-centre coordinates using (2.29). Neglecting the higher order $\unicode[STIX]{x1D746}_{a}$ and $\tilde{F}_{a}$ terms, the leading-order fluid moments $\boldsymbol{M}_{a0}^{lk}$ are given by

The $\unicode[STIX]{x1D703}$ integration can be performed by making use of the gyroaveraging formula of the $\unicode[STIX]{x1D64B}^{l}$ tensor (Ji & Held Reference Ji and Held2009)

where $\unicode[STIX]{x1D709}_{a}=\boldsymbol{c}_{a}\boldsymbol{\cdot }\boldsymbol{b}/c_{a}$ is the pitch-angle velocity coordinate, and $P_{l}$ is a Legendre polynomial defined by

and normalized via

yielding

Finally, we use the basis transformation

with the inverse

to obtain an expression in terms of the Hermite–Laguerre basis. A numerical evaluation of $T_{alk}^{pj}$ and $(T_{a}^{-1})_{pj}^{lk}$ was carried out in Omotani et al. (Reference Omotani, Dudson, Havlíčková and Umansky2015). Instead, in appendix A, we derive the analytic expressions of both $T_{alk}^{pj}$ and $(T_{a}^{-1})_{pj}^{lk}$ .

Using the definition of guiding-centre moments $N_{a}^{pj}$ in (3.6), the leading-order fluid moment $\boldsymbol{M}_{a0}^{lk}$ is then given by

where we define

The leading-order part $C_{aa0}$ of the collision operator $C_{aa}$ can be calculated by approximating $M_{a}^{lk}$ appearing in (3.22) with $\boldsymbol{M}_{a0}^{lk}$ . For the ions, the largest contribution to $\boldsymbol{M}_{i}^{lk}-\boldsymbol{M}_{i0}^{lk}$ is of order $\unicode[STIX]{x1D716}$ and it is given by the $\unicode[STIX]{x1D70C}_{i}$ appearing in (3.24) (the $\tilde{F}_{i}$ correction is smaller since $\tilde{F}_{i}\lesssim \unicode[STIX]{x1D716}^{2}\langle F_{i}\rangle$ , see (2.35)). Therefore, by using the ordering in (2.28), the largest correction to $C_{ii0}$ is $O(\sqrt{m_{e}/m_{i}}\unicode[STIX]{x1D716}\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D708}})$ . The correction to $C_{ee0}$ is of the same order. It follows that we can approximate $C_{aa}$ appearing in (3.22) with $C_{aa0}$ to represent the collision operator up to $O(\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D708}}\unicode[STIX]{x1D716})$ .

As an aside, we note that the relationship between the guiding-centre and fluid moments in (3.32) provides, for the indices $(l,k)=(0,0)$ ,

while, for $(l,k)=(0,1)$ ,

Moreover, the $(l,k)=(2,0)$ moment provides a relationship useful to express the viscosity tensor $\unicode[STIX]{x1D72B}_{a}=\int (\boldsymbol{c}_{a}\boldsymbol{c}_{a}-c_{a}^{2}\unicode[STIX]{x1D644})f_{a}\,\text{d}\boldsymbol{v}$

and $(l,k)=(1,1)$ gives

with $\boldsymbol{q}_{a}$ the heat flux density $\boldsymbol{q}_{a}=m\int \boldsymbol{c}_{a}c_{a}^{2}f_{a}\,\text{d}\boldsymbol{v}/2$ .

In order to express the Boltzmann equation, (2.32), in terms of the guiding-centre moments $N_{a}^{pj}$ , we evaluate the guiding-centre moments of $\langle C_{aa}\rangle$ , namely

By using the gyroaveraging property (3.26) of $\unicode[STIX]{x1D64B}^{l}(\boldsymbol{c}_{a})$ in the like-species operator in (3.21), (3.22) (with $\boldsymbol{M}_{a}^{lk}=\boldsymbol{M}_{a0}^{lk}$ ), and the relation (3.32) between $\boldsymbol{M}_{a0}^{lk}$ and $N_{a}^{pj}$ , the gyroaveraged collision operator coefficients $\langle c(f_{a}^{lkm},f_{a}^{nqr})\rangle$ are given by

with ${\mathcal{P}}_{i+u}^{l,n}=\boldsymbol{P}^{l+n-2(i+u)}\boldsymbol{\cdot }\overline{\boldsymbol{P}^{l}\boldsymbol{\cdot }^{i+u}\unicode[STIX]{x1D64B}^{n}}$ .

Using the basis transformation of (3.31) to express $H_{p}(s_{\Vert a})L_{j}(s_{\bot a}^{2})$ in (3.38) in terms of $c_{a}^{l}P_{l}(\unicode[STIX]{x1D709}_{a})L_{k}^{l+1/2}(c_{a}^{2})$ , and performing the resulting integral, we obtain

with $C_{\ast aabu}^{jw,lm,nr}=\int \,\text{d}\boldsymbol{v}\,c_{a}^{2w+j}f_{Ma}\unicode[STIX]{x1D708}_{\ast aau}^{lm,nr}$ .

We now turn to the electron–ion collision operator, $C_{ei}=C_{ei}^{0}+C_{ei}^{1}$ , with $C_{ei}^{0}$ given by (2.25) and $C_{ei}^{1}$ given by (2.26). As the basis $L_{k}^{l+1/2}\unicode[STIX]{x1D64B}^{l}(\boldsymbol{c}_{a})$ is an eigenfunction of the Lorentz pitch-angle scattering operator $C_{ei}^{0}$ with eigenvalue $-l(l+1)$ (Ji & Held Reference Ji and Held2008), we write $C_{ei}^{0}$ as

Similarly to like-species collisions, we approximate $\boldsymbol{M}_{e}^{lk}\simeq \boldsymbol{M}_{e0}^{lk}$ in (3.41), representing $C_{ei}^{0}$ accurately up to $O(\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D708}}\unicode[STIX]{x1D716})$ . Using the basis transformation in (3.31) and the gyroaverage property of $\unicode[STIX]{x1D64B}^{l}(\boldsymbol{c}_{a})$ in (3.26), we take guiding-centre moments of $C_{ei}$ of the form (3.38), and obtain

where the $A_{ei}$ coefficients are given by

Finally, for the ion–electron collision operator, $C_{ie}$ , we neglect $O(\sqrt{m_{e}/m_{i}}\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D708}}\unicode[STIX]{x1D716})$ corrections by approximating $F_{i}\simeq \langle F_{i}\rangle$ , and use the transformation in (2.3) to convert the $C_{ie}$ operator in (2.27) to guiding-centre variables, yielding

By evaluating $\boldsymbol{R}_{ei}$ at the guiding-centre position $\boldsymbol{R}$ , we write $\boldsymbol{R}_{ei}\boldsymbol{\cdot }\boldsymbol{b}=N_{e}m_{e}v_{th\Vert e}C_{ei}^{10}/\sqrt{2}+O(\sqrt{m_{e}/m_{i}}\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D708}}\unicode[STIX]{x1D716})$ and gyroaverage (3.44), yielding

where we used $c_{\Vert i}^{2}=s_{\Vert i}^{2}T_{\Vert i}/T_{i}$ . Taking guiding-centre moments of the form (3.38) of $\langle C_{ie}\rangle$ in (3.45), we obtain

with

6.1 Fluid equations

We first look at the $(p,j)=(0,0)$ case of (4.9). Noting that $\overline{N}_{a}^{00}=0$ and $C_{ab}^{00}=0$ , we obtain

Evaluating $\Vert \dot{\boldsymbol{R}}\Vert _{a}^{\ast pj}$ in (4.11) and ${\mathcal{F}}_{a}^{pj}$ in (4.10), for $(p,j)=(0,0)$ , equation (6.1) yields the continuity equation

The upper-convective derivative $\text{d}_{a}^{0}/\text{d}t$ , defined by

is related to the guiding-centre fluid velocity $\boldsymbol{u}_{0a}$

and it differs from the lower-convective derivative $\text{d}_{0a}/\text{d}t$ in (4.5) by the addition of the last two terms in (6.4). The vorticity $\unicode[STIX]{x1D6FB}_{\bot }^{2}\unicode[STIX]{x1D719}$ is related to the $\boldsymbol{E}\times \boldsymbol{B}$ drift by

and it appears in (6.2) due to the difference between $\overline{N}_{a}^{00}$ and $N_{a}^{00}$ [see (4.7)]. To derive (6.2), we use the low- $\unicode[STIX]{x1D6FD}$ limit expression for $\boldsymbol{b}\times \boldsymbol{k}\simeq (\boldsymbol{b}\times \unicode[STIX]{x1D735}B)/B$ and neglect $u_{\Vert a}\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \boldsymbol{b}/\unicode[STIX]{x1D6FA}_{a}$ as

therefore keeping up to $O(\unicode[STIX]{x1D716}^{2})$ terms (namely the $\unicode[STIX]{x1D6FB}_{\bot }^{2}\unicode[STIX]{x1D719}$ term in (6.5)).

The parallel momentum equation is obtained by setting $(p,j)=(1,0)$ in (4.9), yielding

with

The expression for $C_{ab}^{10}$ is given in appendix D, as well as all the $C_{ab}^{pj}$ coefficients relevant for the present fluid model. The left-hand side of (6.7) describes the convection of $u_{\Vert a}$ , while the first term in the right-hand side is related to pressure and heat flux gradients, the second term to resistivity (collisional effects), the third term consists of high-order terms kept to ensure phase-space conservation properties, and the last term is the parallel fluid acceleration, namely due to parallel electric fields, mirror force and inertia.

The parallel and perpendicular temperature equations are obtained by setting $(p,j)=(2,0)$ and $(0,1)$ respectively in (4.9). This yields for the parallel temperature

where

and for the perpendicular temperature

with

The equations for the evolution of the parallel $Q_{\Vert a}$ and perpendicular $Q_{\bot a}$ heat fluxes are obtained by setting $(p,j)=(3,0)$ and $(1,1)$ respectively in (4.9), yielding

and

In (6.13) and (6.14) we neglected the higher-order moments with respect to $N^{30}$ and $N^{11}$ , an approximation that we will scrutinize in the next section. Equations (6.2)–(6.14) constitute a closed set of six coupled nonlinear partial differential equations for both the hydrodynamical variables $n_{a},u_{\Vert a},T_{\Vert a},T_{\bot a}$ and the kinetic variables $Q_{\Vert a}$ and $Q_{\bot a}$ .

With respect to previous delta-F (Brizard Reference Brizard1992; Dorland & Hammett Reference Dorland and Hammett1993) and full-F gyrofluid models (Madsen Reference Madsen2013), our fluid model, equations (6.2)–(6.14), while neglecting $k_{\bot }\unicode[STIX]{x1D70C}_{i}\sim 1$ effects, includes the velocity contributions from the $B_{\Vert }^{\ast }$ denominator in the equations of motion (2.15) and (2.16) and includes the effects of full Coulomb collisions up to order $\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D708}}\unicode[STIX]{x1D716}$ . Also, due to the choice of basis functions with shifted velocity arguments $H_{p}(s_{\Vert a})$ instead of $H_{p}(v_{\Vert }/v_{tha})$ , we obtain a set of equations that can efficiently describe both weak flow $(u_{\Vert a}\ll v_{tha})$ and strong flow $(u_{\Vert a}\sim v_{tha})$ regimes.

6.2 High-collisionality regime

We now consider the high-collisionality regime, where the characteristic fluctuation frequency of the hydrodynamical variables $\unicode[STIX]{x1D714}$

is much smaller than the collision frequency $\unicode[STIX]{x1D708}_{a}\simeq \unicode[STIX]{x1D708}_{aa}$ , that is

where the mean free path $\unicode[STIX]{x1D706}_{mfpa}$ in (6.16) is defined as

Equation (6.16) describes the so-called linear transport regime (Balescu Reference Balescu1988). In this case, the distribution function can be expanded around a Maxwell–Boltzmann equilibrium, according to the Chapman–Enskog asymptotic closure scheme (Chapman Reference Chapman1962) and, to first order in $\unicode[STIX]{x1D6FF}_{a}$ , we have

According to (6.18), all moments $N_{a}^{pj}$ in the Hermite–Laguerre expansion (3.1) with $(p,j)\not =(0,0)$ are order $\unicode[STIX]{x1D6FF}_{a}$ . Since $Q_{\Vert a}$ and $Q_{\bot a}$ are determined at first order in $\unicode[STIX]{x1D6FF}_{a}$ only by the moments $(p,j)=(0,0),(3,0),(1,1)$ , the truncation of Sec. (6.1), i.e. neglecting $(p,j)\not =(0,0),(3,0),(1,1)$ is justified. For a more detailed discussion on this topic see Balescu (Reference Balescu1988). Moreover, in the linear regime, a relationship between the hydrodynamical and kinetic variables can be obtained along the lines of the semi-collisional closure. This allows us to express $Q_{\Vert a}$ and $Q_{\bot a}$ as a function of $N_{a},u_{\Vert a},T_{\Vert a}$ and $T_{\bot a}$ , therefore reducing the number of equations. We now derive this functional relationship.

We consider (6.13)–(6.14) in the linear regime, and neglect the polarization terms that are proportional to $\unicode[STIX]{x1D6FB}_{\bot }^{2}\unicode[STIX]{x1D719}/(\unicode[STIX]{x1D6FA}_{a}B)$ . This yields $\sqrt{3/2}\sum _{b}C_{ab}^{30}/v_{th\Vert a}\simeq R_{\Vert a}$ and $\sum _{b}C_{ab}^{11}/(\sqrt{2}v_{th\Vert a})\simeq R_{\bot a}$ , with $R_{\Vert a}$ and $R_{\bot a}$ given by

since $\text{d}_{a}^{0}/\text{d}t\sim \text{d}_{0a}/\text{d}t\sim \unicode[STIX]{x1D714}$ and $(\text{d}^{0}Q_{\Vert ,\bot }/\text{d}t)/Q_{\Vert ,\bot a}\sim \unicode[STIX]{x1D6FF}_{a}^{2}\unicode[STIX]{x1D708}_{a}$ . We compute the guiding-centre moments of the collision operator $C_{ab}^{30}$ and $C_{ab}^{11}$ by truncating the series for the like-species collision operator in (3.40) at $(l,k,n,q)=(2,1,2,1)$ . The resulting $C_{ab}^{pj}$ coefficients are presented in appendix D.

With the expression of $C_{ab}^{30}$ and $C_{ab}^{11}$ , we can solve for $Q_{\Vert a}$ and $Q_{\bot a}$ . In the regime $(T_{\Vert a}-T_{\bot a})/T_{a}\sim \unicode[STIX]{x1D6FF}$ , at lowest order, we obtain for the electron species

and

Analogous expressions are obtained for the ion species.

Equations (6.2), (6.7), (6.9) and (6.11), with $Q_{\Vert a}$ and $Q_{\bot a}$ given by (6.21) and (6.22) are valid in the high-collisionality regime, and can be compared with the drift-reduced Braginskii equations in Zeiler et al. (Reference Zeiler, Drake and Rogers1997). We first rewrite the continuity equation, equation (6.2), in the form

where we expand the convective derivative $\text{d}^{0}a/\text{d}t$ using (6.3) and (6.4), and neglect polarization terms proportional to the electron mass $m_{e}$ . By noting that the diamagnetic drift $v_{de}$ can be written as

and by considering the isotropic regime $T_{\Vert e}\sim T_{\bot e}\sim T_{e}$ , we obtain

which corresponds to the continuity equation in the drift-reduced Braginskii model in Zeiler et al. (Reference Zeiler, Drake and Rogers1997). In that model, the polarization equation is obtained by subtracting both electron and ion continuity equations, using Poisson’s equation $n_{e}\simeq n_{i}$ with $n_{e}$ and $n_{i}$ the electron and ion particle densities respectively, and neglecting the electron to ion mass ratio. Applying the same procedure to the present fluid model, we obtain

In (6.26), the first three terms, which are not present in the drift-reduced Braginskii model, correspond to the difference between ion guiding-centre density $N_{i}$ and particle density $n_{i}$ , proportional to both $\unicode[STIX]{x1D6FB}_{\bot }^{2}\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D6FB}_{\bot }^{2}P_{i}$ . The parallel momentum and temperature equations, equation (6.7) and (6.9), with respect to (Zeiler et al. Reference Zeiler, Drake and Rogers1997), contain the higher-order term ${\mathcal{A}}$ that ensures phase-space conservation, mirror force terms proportional to $(\unicode[STIX]{x1D6FB}_{\Vert }B)/B$ and polarization terms proportional to $\unicode[STIX]{x1D6FB}_{\bot }^{2}\unicode[STIX]{x1D719}/(\unicode[STIX]{x1D6FA}_{a}B)$ due to the difference between guiding-centre and particle fluid quantities. This set of fluid equations constitute an improvement over the drift-reduced Braginskii model. With respect to the original Braginskii equations (Braginskii Reference Braginskii1965), they include the nonlinear terms that arise when retaining full Coulomb collisions, and the effect of ion–electron collisions.