1 Introduction
The plasma periphery, which encompasses the edge and the scrape-off layer regions, plays a central role in determining the overall performance of a fusion device, as it regulates the overall plasma confinement, it controls the plasma–wall interactions, it is responsible for power exhaust and it governs the plasma refuelling and the removal of fusion ashes (Ricci Reference Ricci2015). Understanding the plasma dynamics in the periphery is therefore crucial for the success of the whole fusion program (Connor et al. Reference Connor, Hastie, Wilson and Miller1998).
While the plasma dynamics in the scrape-off layer has been described mainly using drift-reduced fluid models valid at low frequencies compared to the ion cyclotron frequency,
$\unicode[STIX]{x1D714}\ll \unicode[STIX]{x1D6FA}_{i}$
, and in the limits
$k_{\Vert }\unicode[STIX]{x1D706}_{mfpe}\ll 1$
and
$k_{\bot }\unicode[STIX]{x1D70C}_{i}\ll 1$
, i.e. short electron mean free path in comparison to the parallel wavelength and long perpendicular wavelengths with respect to the Larmor radius (Dudson et al.
Reference Dudson, Umansky, Xu, Snyder and Wilson2009; Tamain et al.
Reference Tamain, Ghendrih, Tsitrone, Sarazin, Garbet, Grandgirard, Gunn, Serre, Ciraolo and Chiavassa2009; Ricci et al.
Reference Ricci, Halpern, Jolliet, Loizu, Mosetto, Fasoli, Furno and Theiler2012; Halpern et al.
Reference Halpern, Ricci, Jolliet, Loizu, Morales, Mosetto, Musil, Riva, Tran and Wersal2016; Stegmeir et al.
Reference Stegmeir, Coster, Maj, Hallatschek and Lackner2016; Paruta et al.
Reference Paruta, Ricci, Riva, Wersal, Beadle and Frei2018; Zhu, Francisquez & Rogers Reference Zhu, Francisquez and Rogers2018), these approximations are often marginal near the separatrix and inside it, i.e. in the edge region. In fact, even though turbulence is still dominated by low-frequency fluctuations, the plasma in the edge is hotter and less collisional than in the scrape-off layer and the use of a fluid model becomes questionable. Moreover, in the edge region, small scale
$k_{\bot }\unicode[STIX]{x1D70C}_{i}\sim 1$
fluctuations are important (Hahm, Wang & Madsen Reference Hahm, Wang and Madsen2009). This is especially relevant in the high-temperature tokamak H-mode regime (Zweben et al.
Reference Zweben, Boedo, Grulke, Hidalgo, LaBombard, Maqueda, Scarin and Terry2007), the regime of operation relevant for ITER and future devices. Despite recent progress (Chang et al.
Reference Chang, Ku, Tynan, Hager, Churchill, Cziegler, Greenwald, Hubbard and Hughes2017; Shi et al.
Reference Shi, Hammett, Stoltzfus-Dueck and Hakim2017; Pan et al.
Reference Pan, Told, Shi, Hammett and Jenko2018), overcoming the limitation of the drift-reduced fluid models in the description of the tokamak periphery region by using a gyrokinetic model valid at
$k_{\bot }\unicode[STIX]{x1D70C}_{i}\sim 1$
has proven to be exceptionally demanding. Among the numerous challenges, the effort is undermined by the lack of a proper collisional gyrokinetic model for the periphery.
In fact, with respect to the core, due to the lower temperature values and associated high collisionality, the use of a gyrokinetic model to simulate the plasma dynamics in the tokamak periphery requires an accurate collision operator. This is necessary as collisions set the level of neoclassical transport and strongly influence the turbulent dynamics by affecting the linear growth rate and nonlinear evolution of turbulent modes (Barnes et al. Reference Barnes, Abel, Dorland, Ernst, Hammett, Ricci, Rogers, Schekochihin and Tatsuno2009).
Since the first formulations of the gyrokinetic theory, there have been significant research efforts to take collisions into account (Catto & Tsang Reference Catto and Tsang1977; McCoy, Mirin & Killeen Reference McCoy, Mirin and Killeen1981; Start Reference Start2002; Brizard Reference Brizard2004; Abel et al. Reference Abel, Barnes, Cowley, Dorland and Schekochihin2008b ; Barnes et al. Reference Barnes, Abel, Dorland, Ernst, Hammett, Ricci, Rogers, Schekochihin and Tatsuno2009; Li & Ernst Reference Li and Ernst2011; Dorf et al. Reference Dorf, Cohen, Compton, Dorr, Rognlien, Angus, Krasheninnikov, Colella, Martin and Mccorquodale2012; Estève et al. Reference Estève, Garbet, Sarazin, Grandgirard, Cartier-Michaud, Dif-Pradalier, Ghendrih, Latu and Norscini2015; Hakim et al. Reference Hakim, Francisquez, Juno and Hammet2019; Pan & Ernst Reference Pan and Ernst2019). The first effort devoted to a gyrokinetic collision operator can be traced back to the work of Catto & Tsang (Reference Catto and Tsang1977), later improved by Abel et al. (Reference Abel, Barnes, Cowley, Dorland, Hammett, Schekochihin, Tatsuno, Sauter, Garbet and Sindoni2008a ) by adding the necessary terms needed to ensure non-negative entropy production. The result of this effort is a linearized gyrokinetic collision operator that contains pitch-angle scattering effects and retains important conservation properties. A linearized gyrokinetic Coulomb collision operator derived from first principles was then presented in Li & Ernst (Reference Li and Ernst2011) and Madsen (Reference Madsen2013). However, as turbulence in the tokamak periphery is essentially nonlinear, the relative level of fluctuations in this region being of order unity (Scott Reference Scott2002), and the level of collisions is not sufficient for a local thermalization, the distribution function may significantly deviate from a local Maxwellian distribution (Tskhakaya Reference Tskhakaya2012). Therefore, a nonlinear formulation of the gyrokinetic Coulomb collision operator is crucial to adequately describe the dynamics in the periphery. Only recently, several theoretical studies have emerged in order to derive non-linearized collisional gyrokinetic operators that keep conservation laws in their differential form. In particular, we mention the recent Poisson bracket formulations of the full nonlinear Coulomb collision operator (Brizard Reference Brizard2004; Burby, Brizard & Qin Reference Burby, Brizard and Qin2015; Sugama, Watanabe & Nunami Reference Sugama, Watanabe and Nunami2015). While the formulation of these operators represent significant progress, the presence of a six-dimensional phase-space integral in these expressions makes their numerical implementation still extremely difficult.
In this work, the Coulomb gyrokinetic collision operator is derived in a form that can be efficiently implemented in numerical simulation codes as it involves only integrals over the two gyrokinetic velocity coordinates. The derivation of the full Coulomb collision operator is based on a multipole expansion of the Rosenbluth potentials. This allows us to write the Coulomb collision operator in terms of moments of the distribution function and to apply the gyroaverage operator to the resulting expansion. The Coulomb collision operator is then expressed in terms of two-dimensional velocity integrals of the distribution function. We show that the gyroangle dependence of the expansion coefficients, given in terms of scalar spherical harmonics, allows for analytical gyroaveraging integrations at arbitrary values of the perpendicular wavevector. Furthermore, motivated by recent work based on a pseudo-spectral approach to the gyrokinetic equation (Mandell, Dorland & Landreman Reference Mandell, Dorland and Landreman2018; Frei, Jorge & Ricci Reference Frei, Jorge and Ricci2019), the collision operator is projected onto a Hermite–Laguerre polynomial basis, and is expressed in terms of moments of the distribution function on the same basis. The set of moment-hierarchy equations can then be rigorously closed by using systematic techniques (such as the semi-collisional closure (Zocco & Schekochihin Reference Zocco and Schekochihin2011; Jorge, Ricci & Loureiro Reference Jorge, Ricci and Loureiro2017)) without requiring ad hoc truncation of infinite series.
This paper is organized as follows. Section 2 derives the gyrokinetic equation and § 3 presents the multipole expansion of the Coulomb collision operator. In § 4, the Coulomb operator is ported to a gyrocentre coordinate system, while § 5 projects the collision operator onto a Hermite–Laguerre polynomial basis to obtain a closed-form expression of the Coulomb collision operator in terms of the moments of the distribution function on the same basis. The gyrokinetic collision operator for unlike species is presented in § 6 using an expansion based on the smallness of the electron-to-ion mass ratio. The conclusions follow.
2 Gyrokinetic model
The evolution of the distribution function
$f_{a}=f_{a}(\boldsymbol{x},\boldsymbol{v})$
is given by the Boltzmann equation
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}f_{a}}{\unicode[STIX]{x2202}t}+\dot{\boldsymbol{x}}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}f_{a}}{\unicode[STIX]{x2202}\boldsymbol{x}}+\dot{\boldsymbol{v}}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}f_{a}}{\unicode[STIX]{x2202}\boldsymbol{v}}=\mathop{\sum }_{b}C(f_{a},f_{b}),\end{eqnarray}$$
with
$C(f_{a},f_{b})$
the Coulomb (also known as Landau) collision operator (Landau Reference Landau1936; Rosenbluth, MacDonald & Judd Reference Rosenbluth, MacDonald and Judd1957). This is an operator of the Fokker–Planck type, derived from first principles, and valid in the common case where small-angle Coulomb collisions dominate. Its expression is given by
$$\begin{eqnarray}C(f_{a},f_{b})=L_{ab}\mathop{\sum }_{j,k}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}v_{k}}\left[\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}v_{j}}\left(\frac{\unicode[STIX]{x2202}^{2}G_{b}}{\unicode[STIX]{x2202}v_{k}\unicode[STIX]{x2202}v_{j}}f_{a}\right)-2\left(1+\frac{m_{a}}{m_{b}}\right)\frac{\unicode[STIX]{x2202}H_{b}}{\unicode[STIX]{x2202}v_{k}}f_{a}\right],\end{eqnarray}$$
with
$$\begin{eqnarray}H_{b}=\int \frac{f_{b}(\boldsymbol{v}^{\prime })}{|\boldsymbol{v}-\boldsymbol{v}^{\prime }|}\,\text{d}\boldsymbol{v}^{\prime },\end{eqnarray}$$
and
$$\begin{eqnarray}G_{b}=\int f_{b}(\boldsymbol{v}^{\prime })|\boldsymbol{v}-\boldsymbol{v}^{\prime }|\,\text{d}\boldsymbol{v}^{\prime },\end{eqnarray}$$
the Rosenbluth potentials satisfying
$\unicode[STIX]{x1D6FB}_{v}^{2}G_{b}=2H_{b}$
and
$\unicode[STIX]{x1D6FB}_{v}^{2}H_{b}=-4\unicode[STIX]{x03C0}f_{b}$
(Rosenbluth et al.
Reference Rosenbluth, MacDonald and Judd1957). In (2.2),
$L_{ab}=q_{a}^{2}q_{b}^{2}\unicode[STIX]{x1D706}_{ab}/(8\unicode[STIX]{x03C0}\unicode[STIX]{x1D716}_{0}^{2}m_{a}^{2})=\unicode[STIX]{x1D708}_{ab}v_{tha}^{3}/n_{b}$
is introduced, where
$\unicode[STIX]{x1D706}_{ab}$
and
$\unicode[STIX]{x1D708}_{ab}$
are the Coulomb logarithm and the collision frequency between species
$a$
and
$b$
respectively,
$v_{tha}=\sqrt{2T_{a}/m_{a}}$
is the thermal speed and
$q_{a}$
and
$m_{a}$
are the charge and the mass of particles of species
$a$
,
$a=e,i$
.
In the present paper we consider a plasma with properties that satisfy the gyrokinetic ordering (Brizard & Mishchenko Reference Brizard and Mishchenko2009; Frei et al.
Reference Frei, Jorge and Ricci2019). More precisely, denoting typical turbulent frequencies as
$\unicode[STIX]{x1D714}\sim |\unicode[STIX]{x2202}_{t}\log n|\sim |\unicode[STIX]{x2202}_{t}\log T_{e}|$
with
$n_{e}$
and
$T_{e}$
the electron density and temperature respectively, and typical wavenumbers
$\boldsymbol{k}=k_{\Vert }\boldsymbol{b}+\boldsymbol{k}_{\bot }$
, with
$\boldsymbol{k}\sim |\unicode[STIX]{x1D735}\log n_{e}|\sim |\unicode[STIX]{x1D735}\log T_{e}|$
and
$\boldsymbol{b}=\boldsymbol{B}/B$
the magnetic field unit vector, we assume
$$\begin{eqnarray}\unicode[STIX]{x1D716}\sim \frac{|\boldsymbol{v}_{\boldsymbol{E}}|}{c_{s}}\sim \frac{k_{\Vert }}{k_{\bot }}\ll 1,\end{eqnarray}$$
where
$c_{s}=\sqrt{T_{e}/m_{i}}$
is the sound speed,
$\unicode[STIX]{x1D70C}_{s}=c_{s}/\unicode[STIX]{x1D6FA}_{i}$
the sound Larmor radius,
$\unicode[STIX]{x1D6FA}_{i}=eB/m_{i}$
the ion gyrofrequency and
$\boldsymbol{v}_{\boldsymbol{E}}=\boldsymbol{E}\times \boldsymbol{B}/B^{2}$
the
$\boldsymbol{E}\times \boldsymbol{B}$
drift velocity with
$\boldsymbol{E}=-\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}-\unicode[STIX]{x2202}_{t}\boldsymbol{A}$
the electric field. As strong radial electric fields are known to play a role in the tokamak edge (particularly in the H-mode pedestal), and large-scale fluctuations are the ones at play in the tokamak scrape-off layer, we split the electrostatic potential as
$\unicode[STIX]{x1D719}=\unicode[STIX]{x1D719}_{0}+\unicode[STIX]{x1D719}_{1}$
(Dimits, LoDestro & Dubin Reference Dimits, LoDestro and Dubin1992; Qin et al.
Reference Qin, Cohen, Nevins and Xu2006; Frei et al.
Reference Frei, Jorge and Ricci2019), i.e. into a possibly large-scale drift-kinetic component,
$\unicode[STIX]{x1D719}_{0}$
, satisfying
$$\begin{eqnarray}\frac{e\unicode[STIX]{x1D719}_{0}}{T_{e}}\sim 1,\end{eqnarray}$$
and its small-amplitude gyrokinetic,
$\unicode[STIX]{x1D719}_{1}$
, component
$$\begin{eqnarray}\frac{\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x1D719}_{0}}\sim \unicode[STIX]{x1D716}_{\unicode[STIX]{x1D6FF}}\ll 1.\end{eqnarray}$$
A similar decomposition into large- and small-scale fluctuations is applied to the magnetic vector potential
$\boldsymbol{A}=\boldsymbol{A}_{0}+\boldsymbol{A}_{1}$
, with
$|\boldsymbol{A}_{1}|/|\boldsymbol{A}_{0}|\sim \unicode[STIX]{x1D716}_{\unicode[STIX]{x1D6FF}}$
. Both
$\unicode[STIX]{x1D719}_{0}$
and
$\unicode[STIX]{x1D719}_{1}$
are assumed to yield a similar contribution to the total electric field
$$\begin{eqnarray}\boldsymbol{E}\sim \unicode[STIX]{x1D735}_{\bot }\unicode[STIX]{x1D719}_{0}\sim \unicode[STIX]{x1D735}_{\bot }\unicode[STIX]{x1D719}_{1}.\end{eqnarray}$$
Therefore, by ordering typical gradient lengths of
$\unicode[STIX]{x1D719}_{1}$
to be comparable to
$\unicode[STIX]{x1D70C}_{s}$
, we set
$$\begin{eqnarray}\unicode[STIX]{x1D70C}_{s}\left|\frac{\unicode[STIX]{x1D735}_{\bot }\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x1D719}_{1}}\right|\sim 1,\end{eqnarray}$$
which, using (2.7) and (2.8), constrains typical gradient lengths of
$\unicode[STIX]{x1D719}_{0}$
to be much larger than
$\unicode[STIX]{x1D70C}_{s}$
, as
$$\begin{eqnarray}\unicode[STIX]{x1D70C}_{s}\left|\frac{\unicode[STIX]{x1D735}_{\bot }\unicode[STIX]{x1D719}_{0}}{\unicode[STIX]{x1D719}_{0}}\right|\sim \unicode[STIX]{x1D716}_{\unicode[STIX]{x1D6FF}}.\end{eqnarray}$$
In the following, we set
$\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D6FF}}\sim \unicode[STIX]{x1D716}$
, which, using (2.5), yields
$$\begin{eqnarray}k_{\bot }\unicode[STIX]{x1D70C}_{s}\frac{e\unicode[STIX]{x1D719}_{0}}{T_{e}}\sim \frac{e\unicode[STIX]{x1D719}_{1}}{T_{e}}\sim \unicode[STIX]{x1D716}.\end{eqnarray}$$
The scale length
$L_{B}\sim R_{0}$
of the equilibrium magnetic field (with
$R_{0}$
the major radius of the tokamak device), is ordered by the small parameter
$\unicode[STIX]{x1D716}_{B}\sim \unicode[STIX]{x1D70C}_{s}/L_{B}$
. We note that the collision operator developed here is valid for both
$\unicode[STIX]{x1D716}_{B}\sim \unicode[STIX]{x1D716}^{2}$
and
$\unicode[STIX]{x1D716}_{B}\sim \unicode[STIX]{x1D716}^{3}$
, with the second case being more of interest for the periphery since the plasma temperature is lower than the tokamak core (Hahm et al.
Reference Hahm, Wang and Madsen2009). Finally, the collision frequency is ordered as
$$\begin{eqnarray}\frac{\unicode[STIX]{x1D708}_{i}}{\unicode[STIX]{x1D6FA}_{i}}\sim \unicode[STIX]{x1D716}_{\unicode[STIX]{x1D708}}\sim \unicode[STIX]{x1D716}^{2},\end{eqnarray}$$
with
$\unicode[STIX]{x1D708}_{i}=\unicode[STIX]{x1D708}_{ii}$
the ion–ion collision frequency. For
$T_{e}\sim T_{i}$
, the ordering in (2.12) implies that
$k_{\Vert }\unicode[STIX]{x1D706}_{mfpe}\sim k_{\Vert }\unicode[STIX]{x1D706}_{mfpi}\sim k_{\bot }\unicode[STIX]{x1D70C}_{s}/\unicode[STIX]{x1D716}$
, with
$\unicode[STIX]{x1D706}_{mfpa}=v_{tha}/\unicode[STIX]{x1D708}_{a}$
.
By taking advantage of the ordering in (2.5)–(2.7) and (2.12), the gyrokinetic model effectively removes the fast time scale associated with the cyclotron motion and reduces the dimensionality of the kinetic equation from six phase-space variables, i.e. (
$\boldsymbol{x},\boldsymbol{v}$
), to five. While linear and nonlinear gyrokinetic equations of motion were originally derived using recursive techniques (Rutherford & Frieman Reference Rutherford and Frieman1968; Taylor & Hastie Reference Taylor and Hastie1968; Catto Reference Catto1978), more recent derivations of the gyrokinetic equation based on Hamiltonian Lie perturbation theory (Cary Reference Cary1981) ensure the existence of phase-space volume and magnetic moment conservation laws (Hahm Reference Hahm1988; Brizard & Hahm Reference Brizard and Hahm2007; Hahm et al.
Reference Hahm, Wang and Madsen2009; Frei et al.
Reference Frei, Jorge and Ricci2019). The Hamiltonian derivations are carried out in two steps. In the first step, small-scale electromagnetic fluctuations with perpendicular wavelengths comparable to the particle Larmor radius (
$\unicode[STIX]{x1D719}_{1}$
and
$A_{\Vert 1}$
) are neglected (Cary & Brizard Reference Cary and Brizard2009). Within this approximation, the coordinate transformation from particle phase-space coordinates (
$\boldsymbol{x},\boldsymbol{v}$
) to guiding-centre coordinates
$\boldsymbol{Z}=(\boldsymbol{R},v_{\Vert },\unicode[STIX]{x1D707},\unicode[STIX]{x1D703})$
is derived, where
$\boldsymbol{R}$
is the guiding centre,
$v_{\Vert }$
the parallel velocity,
$\unicode[STIX]{x1D707}$
the magnetic moment and
$\unicode[STIX]{x1D703}$
the gyroangle. The second step introduces small-scale and small-amplitude electromagnetic fluctuations,
$\unicode[STIX]{x1D719}_{1}$
and
$\boldsymbol{A}_{1}$
. For this purpose, a gyrocentre coordinate system
$\overline{\boldsymbol{Z}}=(\overline{\boldsymbol{R}},\overline{v}_{\Vert },\overline{\unicode[STIX]{x1D707}},\overline{\unicode[STIX]{x1D703}})$
is constructed perturbatively from the guiding-centre coordinates
$\boldsymbol{Z}$
via a transformation
$T$
of the form
$$\begin{eqnarray}\overline{\boldsymbol{Z}}=T\boldsymbol{Z}=\boldsymbol{Z}+\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D6FF}}\boldsymbol{Z}_{1}+\cdots \,,\end{eqnarray}$$
where
$\boldsymbol{Z}_{1}$
contains terms proportional to
$\unicode[STIX]{x1D719}_{1}$
and
$\boldsymbol{A}_{1}$
, such that
$\overline{\unicode[STIX]{x1D707}}=T\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}+\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D6FF}}\unicode[STIX]{x1D707}_{1}+\cdots \,$
, remains an adiabatic invariant [see, e.g., Brizard & Hahm (Reference Brizard and Hahm2007)]. This allows us to reduce the number of phase-space variables in the kinetic Boltzmann equation describing the evolution of the particle distribution function from six to five, simplifying the analytical and numerical treatment of magnetized plasma systems.
More precisely, in order to simplify and remove the gyroangle dependence of the Boltzmann equation, the distribution function
$f_{a}(\boldsymbol{x},\boldsymbol{v})$
is first expressed in terms of the guiding-centre coordinates
$\boldsymbol{Z}$
by defining the guiding-centre distribution function
$F_{a}(\boldsymbol{Z})$
as
$$\begin{eqnarray}F_{a}(\boldsymbol{Z})=f_{a}(\boldsymbol{x}(\boldsymbol{Z}),\boldsymbol{v}(\boldsymbol{Z})).\end{eqnarray}$$
The coordinate transformation
$\boldsymbol{v}(\boldsymbol{Z})$
is given by
$$\begin{eqnarray}\boldsymbol{v}=v_{\Vert }\boldsymbol{b}+\boldsymbol{v}_{E0}+v_{\bot }^{\prime }(\cos \unicode[STIX]{x1D703}\boldsymbol{e}_{1}+\sin \unicode[STIX]{x1D703}\boldsymbol{e}_{2})=v[\cos \unicode[STIX]{x1D711}\boldsymbol{b}+\sin \unicode[STIX]{x1D711}(\cos \unicode[STIX]{x1D703}\boldsymbol{e}_{1}+\sin \unicode[STIX]{x1D703}\boldsymbol{e}_{2})],\end{eqnarray}$$
with
$(\boldsymbol{b},\boldsymbol{e}_{1},\boldsymbol{e}_{2})$
a fixed right-handed coordinate set,
$v_{E0}=-\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{0}$
the drift-kinetic
$\boldsymbol{E}\times \boldsymbol{B}$
drift,
$\cos \unicode[STIX]{x1D711}=v_{\Vert }/v$
the cosine of the pitch angle and
$\unicode[STIX]{x1D703}$
the gyroangle. The magnetic moment
$\unicode[STIX]{x1D707}$
is defined as
$$\begin{eqnarray}\unicode[STIX]{x1D707}=\frac{m_{a}v_{\bot }^{\prime 2}}{2B},\end{eqnarray}$$
whereas the particle position
$\boldsymbol{x}(\boldsymbol{Z})$
is written as
$$\begin{eqnarray}\boldsymbol{x}=\boldsymbol{R}+\unicode[STIX]{x1D746}_{a},\end{eqnarray}$$
with
$$\begin{eqnarray}\unicode[STIX]{x1D746}_{a}=\unicode[STIX]{x1D746}_{a}(\boldsymbol{R},\unicode[STIX]{x1D707},\unicode[STIX]{x1D703})=\sqrt{\frac{2m_{a}\unicode[STIX]{x1D707}}{q_{a}^{2}B}}(-\sin \unicode[STIX]{x1D703}\boldsymbol{e}_{1}+\cos \unicode[STIX]{x1D703}\boldsymbol{e}_{2})\end{eqnarray}$$
the Larmor radius and
$\boldsymbol{R}$
the guiding centre of the particle. The Jacobian of the guiding-centre transformation of (2.15)–(2.17) is given by
$B_{\Vert }^{\ast }/m_{a}=(B/m_{a})(1+\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \boldsymbol{v}_{E}/\unicode[STIX]{x1D6FA}_{a}+v_{\Vert }\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \boldsymbol{b}/\unicode[STIX]{x1D6FA}_{a})$
(Cary & Brizard Reference Cary and Brizard2009). We note that, in a weak-flow regime where
$\boldsymbol{v}_{E}$
is absent from (2.15), the calculation that follows remains valid, except for the Jacobian of the guiding-centre transformation
$B_{\Vert }^{\ast }$
, which should be replaced by
$B_{\Vert }^{\ast }=(B/m_{a})(1+v_{\Vert }\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \boldsymbol{b}/\unicode[STIX]{x1D6FA}_{a})$
.
To account for the small-scale fluctuations and magnetic inhomogeneity, the gyrokinetic distribution function
$\overline{F}_{a}(\overline{\boldsymbol{Z}})$
is then defined as (Brizard & Hahm Reference Brizard and Hahm2007)
$$\begin{eqnarray}\overline{F_{a}}(\overline{\boldsymbol{Z}})=F_{a}(\boldsymbol{Z})=T\overline{F}_{a}(\boldsymbol{Z}),\end{eqnarray}$$
with the coordinate transformation between
$\overline{\boldsymbol{Z}}$
and
$\boldsymbol{Z}$
given perturbatively by (2.13). Indeed, using the chain rule to rewrite the Boltzmann equation, equation (2.1), in terms of gyrocentre
$\overline{\boldsymbol{Z}}$
coordinates, we obtain
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\overline{F_{a}}}{\unicode[STIX]{x2202}t}+\dot{\overline{\boldsymbol{Z}}}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}\overline{F_{a}}}{\unicode[STIX]{x2202}\overline{\boldsymbol{Z}}}=\mathop{\sum }_{b}C(\overline{F_{a}},\overline{F_{b}}).\end{eqnarray}$$
We now introduce the gyroaverage operator
$\langle \cdots \rangle _{\overline{\boldsymbol{R}}}$
defined by
$$\begin{eqnarray}\langle \unicode[STIX]{x1D712}\rangle _{\overline{\boldsymbol{R}}}=\frac{1}{2\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}\unicode[STIX]{x1D712}(\overline{\boldsymbol{Z}})\,\text{d}\overline{\unicode[STIX]{x1D703}},\end{eqnarray}$$
where all gyrocentre coordinates
$\overline{\boldsymbol{Z}}$
but
$\overline{\unicode[STIX]{x1D703}}$
are kept fixed during the integration. By applying the gyroaverage operator to (2.20), the gyrokinetic equation is obtained
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\langle \overline{F_{a}}\rangle _{\overline{\boldsymbol{R}}}+\left\langle \dot{\overline{\boldsymbol{Z}}}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}\overline{F_{a}}}{\unicode[STIX]{x2202}\overline{\boldsymbol{Z}}}\right\rangle _{\overline{\boldsymbol{R}}}=\mathop{\sum }_{b}\langle C(\overline{F_{a}},\overline{F_{b}})\rangle _{\overline{\boldsymbol{R}}}.\end{eqnarray}$$
Equation (2.22) can be further simplified by noting that, in the gyrokinetic framework, the transformation in (2.13) is constructed in such a way that the gyrocentre equations of motion, i.e. the equations that determine
$\dot{\overline{\boldsymbol{Z}}}=(\dot{\overline{\boldsymbol{R}}},\dot{\overline{v_{\Vert }}},\dot{\overline{\unicode[STIX]{x1D707}}},\dot{\overline{\unicode[STIX]{x1D703}}})$
, are gyroangle independent and that
$\overline{\unicode[STIX]{x1D707}}$
is an adiabatic invariant satisfying
$\dot{\overline{\unicode[STIX]{x1D707}}}=0$
(Brizard & Hahm Reference Brizard and Hahm2007). Therefore, the gyrokinetic equation in (2.22) can be written as
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\langle \overline{F_{a}}\rangle _{\overline{\boldsymbol{R}}}+\dot{\overline{\boldsymbol{R}}}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}\langle \overline{F_{a}}\rangle _{\overline{\boldsymbol{R}}}}{\unicode[STIX]{x2202}\overline{\boldsymbol{R}}}+\dot{\overline{v_{\Vert }}}\frac{\unicode[STIX]{x2202}\langle \overline{F_{a}}\rangle _{\overline{\boldsymbol{R}}}}{\unicode[STIX]{x2202}\overline{v_{\Vert }}}=\mathop{\sum }_{b}\langle C(\overline{F_{a}},\overline{F_{b}})\rangle _{\overline{\boldsymbol{R}}}.\end{eqnarray}$$
In order to further simplify (2.23), we estimate the order of magnitude of the gyrophase-dependent part of the distribution function
$\widetilde{\overline{F_{a}}}=\overline{F_{a}}-\langle \overline{F_{a}}\rangle _{\overline{\boldsymbol{R}}}$
, where
$\overline{F_{a}}$
obeys (2.20) and
$\langle \overline{F_{a}}\rangle _{\overline{\boldsymbol{R}}}$
obeys (2.23). For this purpose, we note that the equation for the evolution of
$\widetilde{\overline{F_{a}}}=\overline{F_{a}}-\langle \overline{F_{a}}\rangle _{\overline{\boldsymbol{R}}}$
can be obtained by subtracting (2.23) from (2.20), yielding
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\widetilde{\overline{F_{a}}}}{\unicode[STIX]{x2202}t}+\dot{\overline{\boldsymbol{R}}}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}\widetilde{\overline{F_{a}}}}{\unicode[STIX]{x2202}\overline{\boldsymbol{R}}}+\dot{\overline{v_{\Vert }}}\frac{\unicode[STIX]{x2202}\widetilde{\overline{F_{a}}}}{\unicode[STIX]{x2202}\overline{v_{\Vert }}}+\dot{\overline{\unicode[STIX]{x1D703}}}\frac{\unicode[STIX]{x2202}\widetilde{\overline{F_{a}}}}{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D703}}}=\mathop{\sum }_{b}C(\overline{F_{a}},\overline{F_{b}})-\langle C(\overline{F_{a}},\overline{F_{b}})\rangle _{\overline{\boldsymbol{R}}}.\end{eqnarray}$$
To lowest order,
$\dot{\overline{\unicode[STIX]{x1D703}}}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}\widetilde{\overline{F_{a}}}\sim \unicode[STIX]{x1D6FA}_{a}\widetilde{\overline{F_{a}}}$
and
$\unicode[STIX]{x2202}_{t}\sim \dot{\overline{\boldsymbol{R}}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\overline{\boldsymbol{R}}}\sim \dot{\overline{v_{\Vert }}}\unicode[STIX]{x2202}_{\overline{v_{\Vert }}}\sim \unicode[STIX]{x1D716}\unicode[STIX]{x1D6FA}_{i}$
. Therefore, the leading-order estimate of (2.24) gives
$$\begin{eqnarray}\widetilde{\overline{F_{a}}}\simeq \frac{1}{\unicode[STIX]{x1D6FA}_{a}}\mathop{\sum }_{b}\int _{0}^{\overline{\unicode[STIX]{x1D703}}}[C(\langle \overline{F_{a}}\rangle _{\overline{\boldsymbol{R}}},\langle \overline{F_{b}}\rangle _{\overline{\boldsymbol{R}}})-\langle C(\langle \overline{F_{a}}\rangle _{\overline{\boldsymbol{R}}},\langle \overline{F_{b}}\rangle _{\overline{\boldsymbol{R}}})\rangle _{\overline{\boldsymbol{R}}}]\,\text{d}\overline{\unicode[STIX]{x1D703}}^{\prime }.\end{eqnarray}$$
Using the fact that
$C(\overline{F_{a}},\overline{F_{b}})\sim \unicode[STIX]{x1D708}_{a}\overline{F}_{a}$
, together with (2.12), and expanding
$\overline{F_{a}}$
as
$\overline{F_{a}}=\langle \overline{F_{a}}\rangle _{\overline{\boldsymbol{R}}}+\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D708}}\overline{F_{a1}}+\cdots$
, we find that
$$\begin{eqnarray}\frac{\widetilde{\overline{F_{e}}}}{\langle \overline{F_{e}}\rangle _{\overline{\boldsymbol{R}}}}\sim \left(\frac{T_{i}}{T_{e}}\right)^{3/2}\sqrt{\frac{m_{e}}{m_{i}}}\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D708}}\sim \left(\frac{T_{i}}{T_{e}}\right)^{3/2}\sqrt{\frac{m_{e}}{m_{i}}}\unicode[STIX]{x1D716}^{2},\end{eqnarray}$$
and
$$\begin{eqnarray}\frac{\widetilde{\overline{F}_{i}}}{\langle \overline{F}_{i}\rangle _{\overline{\boldsymbol{R}}}}\sim \frac{\unicode[STIX]{x1D708}_{i}}{\unicode[STIX]{x1D6FA}_{i}}\sim \unicode[STIX]{x1D716}_{\unicode[STIX]{x1D708}}\sim \unicode[STIX]{x1D716}^{2},\end{eqnarray}$$
showing that, up to second order in
$\unicode[STIX]{x1D716}$
, the gyroangle dependence of the distribution function can be neglected in (2.23). We remark that a similar estimate for the gyrophase dependent part of the guiding-centre distribution function was found in Jorge et al. (Reference Jorge, Ricci and Loureiro2017).
We now evaluate the magnitude of the collisional term in (2.24). Using the expansion
$C(\overline{F_{a}},\overline{F_{b}})=C_{0}(\overline{F_{a}},\overline{F_{b}})+\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D6FF}}C_{1}(\overline{F_{a}},\overline{F_{b}})+\cdots \!$
with
$C_{0}(\overline{F_{a}},\overline{F_{b}})\sim \unicode[STIX]{x1D708}_{a}\overline{F_{a}}$
, and noting that the first-order gyrocentre transformation
$\boldsymbol{Z}_{1}=\boldsymbol{Z}-\overline{\boldsymbol{Z}}+O(\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D6FF}}^{2})$
in (2.13) is mass dependent, i.e.
$\overline{\boldsymbol{Z}_{1e}}\sim \sqrt{m_{e}/m_{i}}\overline{\boldsymbol{Z}_{1i}}$
(see, e.g., Brizard & Hahm (Reference Brizard and Hahm2007)), the magnitude of the Coulomb collision operator for electrons can be estimated as
$$\begin{eqnarray}C(\overline{F}_{e},\overline{F_{b}})\sim \unicode[STIX]{x1D708}_{e}\overline{F}_{e}\sim \sqrt{\frac{m_{i}}{m_{e}}}\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D708}}\unicode[STIX]{x1D6FA}_{i}\overline{F}_{e}+O(\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D708}}\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D6FF}}\unicode[STIX]{x1D6FA}_{i}\overline{F}_{e})\sim \sqrt{\frac{m_{i}}{m_{e}}}\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D6FA}_{i}\overline{F}_{e}+O(\unicode[STIX]{x1D716}^{3}\unicode[STIX]{x1D6FA}_{i}\overline{F}_{e}).\end{eqnarray}$$
Thus, the third-order term in the expansion in (2.28) does not contain a
$\sqrt{m_{i}/m_{e}}$
factor, in contrast to its lower-order counterpart. A similar argument holds for the ions, yielding
$$\begin{eqnarray}C(\overline{F}_{i},\overline{F_{b}})\sim \unicode[STIX]{x1D708}_{i}\overline{F}_{i}\sim \unicode[STIX]{x1D716}_{\unicode[STIX]{x1D708}}\unicode[STIX]{x1D6FA}_{i}\overline{F}_{i}+O(\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D708}}\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D6FF}}\unicode[STIX]{x1D6FA}_{i}\overline{F}_{i})\sim \unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D6FA}_{i}\overline{F}_{i}+O(\unicode[STIX]{x1D716}^{3}\unicode[STIX]{x1D6FA}_{i}\overline{F}_{i}).\end{eqnarray}$$
Equations (2.28) and (2.29) show that the lowest-order collision operator
$C_{0}(\overline{F_{a}},\overline{F_{b}})$
is, in fact,
$O(\unicode[STIX]{x1D716}^{2})$
. Therefore, the gyrokinetic equation valid up to second order in
$\unicode[STIX]{x1D716}$
, considered in a large number of edge gyrokinetic models (see, e.g., Qin et al. (Reference Qin, Cohen, Nevins and Xu2006, Reference Qin, Cohen, Nevins and Xu2007), Hahm et al. (Reference Hahm, Wang and Madsen2009) and Frei et al. (Reference Frei, Jorge and Ricci2019)), can thus be written as
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\langle \overline{F_{a}}\rangle _{\overline{\boldsymbol{R}}}+\dot{\overline{\boldsymbol{R}}}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\overline{\boldsymbol{R}}}\langle \overline{F_{a}}\rangle _{\overline{\boldsymbol{R}}}+\dot{\overline{v_{\Vert }}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\overline{v_{\Vert }}}\langle \overline{F_{a}}\rangle _{\overline{\boldsymbol{R}}}=\mathop{\sum }_{b}\langle C_{0}(\langle \overline{F_{a}}\rangle _{\overline{\boldsymbol{R}}},\langle \overline{F_{b}}\rangle _{\overline{\boldsymbol{R}}})\rangle _{\overline{\boldsymbol{R}}},\end{eqnarray}$$
assuming
$\dot{\overline{\boldsymbol{R}}}$
and
$\dot{\overline{v_{\Vert }}}$
to be at least
$O(\unicode[STIX]{x1D716}^{2})$
accurate. We note that, although only the lowest order in
$\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D6FF}}$
collision operator
$C_{0}(\langle \overline{F_{a}}\rangle _{\overline{\boldsymbol{R}}},\langle \overline{F_{b}}\rangle _{\overline{\boldsymbol{R}}})$
is used in (2.30), all orders in
$k_{\bot }\unicode[STIX]{x1D70C}_{s}$
are kept.
3 Multipole expansion of the Coulomb collision operator
The goal of this section is to find a suitable basis to expand
$f_{a}$
such that the Coulomb operator in (2.2) can be cast as a function of moments of
$f_{a}$
. This first step considerably simplifies the derivation of the gyrokinetic collision operator. We start by noting that the Rosenbluth potential
$H_{b}$
in (2.3) is analogous to the expression of the electrostatic potential due to a charge distribution, a similarity already noted by Rosenbluth et al. (Reference Rosenbluth, MacDonald and Judd1957). This fact allows us to make use of known electrostatic expansion techniques (Jackson Reference Jackson1998) to perform a multipole expansion of the Rosenbluth potentials. We first Taylor expand the factor
$1/|\boldsymbol{v}-\boldsymbol{v}^{\prime }|$
in (2.3) around
$\boldsymbol{v}=0$
if
$v\leqslant v^{\prime }$
or around
$\boldsymbol{v}^{\prime }=0$
if
$v>v^{\prime }$
, yielding
$$\begin{eqnarray}\frac{1}{|\boldsymbol{v}-\boldsymbol{v}^{\prime }|}=\left\{\begin{array}{@{}ll@{}}\displaystyle \mathop{\sum }_{l=0}^{\infty }{\displaystyle \frac{(-\boldsymbol{v}^{\prime })^{l}}{l!}}\boldsymbol{\cdot }{\displaystyle \frac{\unicode[STIX]{x2202}^{l}}{\unicode[STIX]{x2202}\boldsymbol{v}^{l}}}\left({\displaystyle \frac{1}{v}}\right),\quad & v^{\prime }\leqslant v,\\ \displaystyle \mathop{\sum }_{l=0}^{\infty }{\displaystyle \frac{(-\boldsymbol{v})^{l}}{l!}}\boldsymbol{\cdot }{\displaystyle \frac{\unicode[STIX]{x2202}^{l}}{\unicode[STIX]{x2202}(\boldsymbol{v}^{\prime })^{l}}}\left({\displaystyle \frac{1}{v^{\prime }}}\right),\quad & v<v^{\prime },\end{array}\right.\end{eqnarray}$$
where we used the identity
$\unicode[STIX]{x2202}_{\boldsymbol{v}}(1/|\boldsymbol{v}-\boldsymbol{v}^{\prime }|)_{v=0}=-\unicode[STIX]{x2202}_{\boldsymbol{v}^{\prime }}(1/v^{\prime })$
and where we denote the inner product between all the
$l$
indices of two
$l$
-rank tensors,
$\unicode[STIX]{x1D64F}_{1}^{l}$
and
$\unicode[STIX]{x1D64F}_{2}^{l}$
, as
$\unicode[STIX]{x1D64F}_{1}^{l}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}_{2}^{l}$
. Both
$v\leqslant v^{\prime }$
and
$v>v^{\prime }$
cases are included in order to take into account the fact that
$f_{b}(\boldsymbol{v}^{\prime })$
is, in general, finite over the entire velocity space
$\boldsymbol{v}^{\prime }$
. Denoting
$\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})$
the spherical harmonic tensor (Weinert Reference Weinert1980)
$$\begin{eqnarray}\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})=\frac{(-1)^{l}v^{2l+1}}{(2l-1)!!}\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{v}}\right)^{l}\frac{1}{v},\end{eqnarray}$$
we obtain the following expression for
$H_{b}$
$$\begin{eqnarray}H_{b}=2\mathop{\sum }_{l=0}^{\infty }\frac{(2l-1)!!}{l!}\left(\int _{v>v^{\prime }}f_{b}(\boldsymbol{v}^{\prime })\frac{(\boldsymbol{v}^{\prime })^{l}}{v^{2l+1}}\boldsymbol{\cdot }\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})\,\text{d}\boldsymbol{v}^{\prime }+\int _{v^{\prime }\geqslant v}f_{b}(\boldsymbol{v}^{\prime })\frac{(\boldsymbol{v})^{l}}{(v^{\prime })^{2l+1}}\boldsymbol{\cdot }\unicode[STIX]{x1D654}^{l}(\boldsymbol{v}^{\prime })\,\text{d}\boldsymbol{v}^{\prime }\right).\end{eqnarray}$$
In order to simplify (3.3), we note that the tensor
$\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})=Y_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\cdots \unicode[STIX]{x1D6FE}}^{l}(\boldsymbol{v})$
is symmetric and totally traceless, i.e. traceless between any combination of two of its indices. Symmetry arises from the fact that any couple of indices in
$Y_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\ldots \unicode[STIX]{x1D6FE}}^{l}(\boldsymbol{v})$
are interchangeable as the velocity derivatives commute for
$v\not =0$
. The totally traceless feature (i.e.
$\sum _{\unicode[STIX]{x1D6FC}}Y_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FC}\ldots \unicode[STIX]{x1D6FE}}^{l}(\boldsymbol{v})=0$
and the same for any other pair of indices) stems from the fact that the contraction between any two indices in
$Y_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\cdots \unicode[STIX]{x1D6FE}}^{l}(\boldsymbol{v})$
leads to the multiplicative factor
$\unicode[STIX]{x1D6FB}_{\boldsymbol{v}}^{2}=\unicode[STIX]{x2202}_{\boldsymbol{v}}\boldsymbol{\cdot }\unicode[STIX]{x2202}_{\boldsymbol{v}}(1/v)$
, which vanishes for
$v\not =0$
(we note that
$Y^{l}(\boldsymbol{v})$
vanishes in the limit
$\boldsymbol{v}=0$
). Furthermore, by defining the tensor
$(\unicode[STIX]{x1D66B})_{TS}^{l}$
as the traceless symmetric counterpart of
$(\boldsymbol{v})^{l}$
(e.g.
$(\unicode[STIX]{x1D66B})_{TS}^{2}=\boldsymbol{v}\boldsymbol{v}-\boldsymbol{I}v^{2}/3$
with
$\boldsymbol{I}$
the identity matrix), we replace the tensors
$(\boldsymbol{v}^{\prime })^{l}$
and
$(\boldsymbol{v})^{l}$
in (3.3) by their traceless symmetric counterparts
$(\boldsymbol{v}^{\prime })_{TS}^{l}$
and
$(\unicode[STIX]{x1D66B})_{TS}^{l}$
respectively
$$\begin{eqnarray}H_{b}=2\mathop{\sum }_{l=0}^{\infty }\frac{(2l-1)!!}{l!}\left(\int _{v>v^{\prime }}f_{b}(\boldsymbol{v}^{\prime })\frac{(\boldsymbol{v}^{\prime })_{TS}^{l}}{v^{2l+1}}\boldsymbol{\cdot }\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})\,\text{d}\boldsymbol{v}^{\prime }+\int _{v^{\prime }\geqslant v}f_{b}(\boldsymbol{v}^{\prime })\frac{(\unicode[STIX]{x1D66B})_{TS}^{l}}{(v^{\prime })^{2l+1}}\boldsymbol{\cdot }\unicode[STIX]{x1D654}^{l}(\boldsymbol{v}^{\prime })\,\text{d}\boldsymbol{v}^{\prime }\right),\end{eqnarray}$$
as they differ only by terms proportional to the identity matrix that vanish when summed with
$\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})$
and
$\unicode[STIX]{x1D654}^{l}(\boldsymbol{v}^{\prime })$
. In fact, for the
$l=2$
case we have
$(\boldsymbol{v}^{2}-(\unicode[STIX]{x1D66B})_{TS}^{2})\boldsymbol{\cdot }\unicode[STIX]{x1D654}^{2}(\boldsymbol{v})=(v^{2}/3)\boldsymbol{I}\boldsymbol{\cdot }\unicode[STIX]{x1D654}^{2}(\boldsymbol{v})=(v^{2}/3)\sum _{\unicode[STIX]{x1D6FC}}Y_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FC}}^{2}=0$
, and similarly for
$l>2$
. In addition, following the convention in Snider (Reference Snider2017), scalars (
$l=0$
) and vectors (
$l=1$
) are considered to be traceless symmetric quantities. Finally, we relate the tensors
$(\unicode[STIX]{x1D66B})_{TS}^{l}$
and
$\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})$
. For
$l=0$
and
$l=1$
, we have
$\unicode[STIX]{x1D654}^{0}(\boldsymbol{v}^{\prime })=1=(\boldsymbol{v}^{\prime })_{TS}^{0}$
and
$\unicode[STIX]{x1D654}^{1}(\boldsymbol{v}^{\prime })=\boldsymbol{v}^{\prime }=(\boldsymbol{v}^{\prime })_{TS}^{1}$
. For
$l=2$
, equation (3.2) gives
$$\begin{eqnarray}\unicode[STIX]{x1D654}^{2}(\boldsymbol{v}^{\prime })=\boldsymbol{v}^{\prime }\boldsymbol{v}^{\prime }-\frac{v^{\prime 2}}{3}\boldsymbol{I}=(\boldsymbol{v}^{\prime })_{TS}^{2}.\end{eqnarray}$$
The results obtained for
$l=0,1$
and 2 can be generalized, i.e.
$(\boldsymbol{v}^{\prime })_{TS}^{l}=\unicode[STIX]{x1D654}^{l}(\boldsymbol{v}^{\prime })$
as proved by induction (Weinert Reference Weinert1980). The Rosenbluth potential
$H_{b}$
can therefore be written as
$$\begin{eqnarray}H_{b}=2\mathop{\sum }_{l=0}^{\infty }\frac{(2l-1)!!}{l!}\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})\boldsymbol{\cdot }\left[\frac{1}{(v^{2})^{l+1/2}}\int _{v^{\prime }<v}f_{b}(\boldsymbol{v}^{\prime })\unicode[STIX]{x1D654}^{l}(\boldsymbol{v}^{\prime })d\boldsymbol{v}^{\prime }+\int _{v^{\prime }\geqslant v}f_{b}(\boldsymbol{v}^{\prime })\frac{\unicode[STIX]{x1D654}^{l}(\boldsymbol{v}^{\prime })}{[(v^{\prime })^{2}]^{l+1/2}}\,\text{d}\boldsymbol{v}^{\prime }\right].\end{eqnarray}$$
The first term in (3.6) can be regarded as the potential due to the charge distribution
$f_{b}(\boldsymbol{v}^{\prime })$
inside a sphere of radius
$v$
, while the second term is the potential due to a finite charge distribution
$f_{b}(\boldsymbol{v}^{\prime })$
at
$v^{\prime }\geqslant v$
.
We now look for an expansion of
$f_{b}$
that allows us to perform the integrals in (3.6) analytically by writing
$H_{b}$
as a sum of velocity moments of
$f_{b}$
. We consider the basis functions (Hirshman & Sigmar Reference Hirshman and Sigmar1976)
$$\begin{eqnarray}\unicode[STIX]{x1D654}^{lk}(\boldsymbol{v})=\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})L_{k}^{l+1/2}(v^{2}),\end{eqnarray}$$
with
$L_{k}^{l+1/2}(v)$
an associated Laguerre polynomial (Abramowitz, Stegun & Miller Reference Abramowitz, Stegun and Miller1965), i.e.
$$\begin{eqnarray}L_{k}^{l+1/2}(v)=\mathop{\sum }_{m=0}^{k}L_{km}^{l}v^{m},\end{eqnarray}$$
where
$$\begin{eqnarray}L_{km}^{l}=\frac{(-1)^{m}(l+k+1/2)!}{(k-m)!(l+m+1/2)!m!}.\end{eqnarray}$$
The basis
$\unicode[STIX]{x1D654}^{lk}(\boldsymbol{v})$
is orthogonal, with the orthogonality relation given by (Banach & Piekarski Reference Banach and Piekarski1989; Snider Reference Snider2017)
$$\begin{eqnarray}\int \text{e}^{-v^{2}}\unicode[STIX]{x1D654}^{l^{\prime }k^{\prime }}(\boldsymbol{v})\unicode[STIX]{x1D654}^{lk}(\boldsymbol{v})\,\text{d}\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}^{lk}=\unicode[STIX]{x1D6FF}_{ll^{\prime }}\unicode[STIX]{x1D6FF}_{kk^{\prime }}\unicode[STIX]{x03C0}^{3/2}\unicode[STIX]{x1D70E}_{k}^{l}\unicode[STIX]{x1D64F}^{lk},\end{eqnarray}$$
with
$\unicode[STIX]{x1D64F}^{lk}$
an arbitrary symmetric and traceless tensor, and
$\unicode[STIX]{x1D70E}_{k}^{l}$
the normalization constant
$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{k}^{l}=\frac{l!(l+k+1/2)!}{2^{l}(l+1/2)!k!}.\end{eqnarray}$$
A proof that
$\unicode[STIX]{x1D654}^{lk}(\boldsymbol{v})$
is a complete basis, i.e. that each
$l$
and
$k$
element of
$\unicode[STIX]{x1D654}^{lk}(\boldsymbol{v})$
is linearly independent and that a linear combination of its elements spans any smooth function
$f(\boldsymbol{v})$
, can be found in Banach & Piekarski (Reference Banach and Piekarski1989), where the equivalence between Grad’s moment expansion in tensorial Hermite polynomials (which forms a complete basis) and
$\unicode[STIX]{x1D654}^{lk}(\boldsymbol{v})$
is shown. We then write
$f_{b}$
as
$$\begin{eqnarray}f_{b}=f_{Mb}\mathop{\sum }_{l,k=0}^{\infty }\unicode[STIX]{x1D654}^{lk}\left(\frac{\boldsymbol{v}}{v_{thb}}\right)\boldsymbol{\cdot }\frac{\boldsymbol{M}_{b}^{lk}}{\unicode[STIX]{x1D70E}_{k}^{l}},\end{eqnarray}$$
with
$f_{Mb}$
a Maxwellian distribution function
$$\begin{eqnarray}f_{Mb}=\frac{n_{b}}{v_{thb}^{3}\unicode[STIX]{x03C0}^{3/2}}\text{e}^{-v^{2}/v_{thb}^{2}}.\end{eqnarray}$$
According to (3.10), with
$\unicode[STIX]{x1D64F}^{lk}=\boldsymbol{M}_{b}^{lk}$
, the coefficients
$\boldsymbol{M}_{b}^{lk}$
are obtained by taking velocity moments of
$f_{b}$
of the form
$$\begin{eqnarray}\boldsymbol{M}_{b}^{lk}=\frac{1}{n_{b}}\int f_{b}(\boldsymbol{v})\unicode[STIX]{x1D654}^{lk}\left(\frac{\boldsymbol{v}}{v_{thb}}\right)\text{d}\boldsymbol{v}.\end{eqnarray}$$
Finally, we note that (3.12) allows us to retain only the
$l=k=0$
moment when the plasma is in thermal equilibrium.
Plugging the expansion for
$f_{b}$
given by (3.12) into (3.6), we obtain the following expression
$$\begin{eqnarray}\displaystyle H_{b} & = & \displaystyle \frac{n_{b}}{v_{thb}\unicode[STIX]{x03C0}^{3/2}}\mathop{\sum }_{l,l^{\prime },k}\frac{(2l-1)!!}{l!\unicode[STIX]{x1D70E}_{k}^{l^{\prime }}}\nonumber\\ \displaystyle & & \displaystyle \times \left(\frac{\unicode[STIX]{x1D654}^{l}(\hat{v})}{x_{b}^{(l+1)/2}}\boldsymbol{\cdot }\int _{0}^{x_{b}}\text{e}^{-x}L_{k}^{l^{\prime }+1/2}(x)x^{(l+l^{\prime }+1)/2}\,\text{d}x\int \unicode[STIX]{x1D654}^{l}(\hat{v}^{\prime })\unicode[STIX]{x1D654}^{l^{\prime }}(\hat{v}^{\prime })\,\text{d}\unicode[STIX]{x1D70E}^{\prime }\boldsymbol{\cdot }\boldsymbol{M}_{b}^{l^{\prime }k}\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,x_{b}^{l/2}\unicode[STIX]{x1D654}^{l}(\hat{v})\boldsymbol{\cdot }\int _{x_{b}}^{\infty }\text{e}^{-x}L_{k}^{l+1/2}(x)\,\text{d}x\int \unicode[STIX]{x1D654}^{l}(\hat{v}^{\prime })\unicode[STIX]{x1D654}^{l^{\prime }}(\hat{v}^{\prime })\,\text{d}\unicode[STIX]{x1D70E}^{\prime }\boldsymbol{\cdot }\boldsymbol{M}_{b}^{l^{\prime }k}\right),\end{eqnarray}$$
where we define the normalized velocity
$x_{b}=v^{2}/v_{thb}^{2}$
, the solid angle
$\unicode[STIX]{x1D70E}$
such that
$\text{d}\boldsymbol{v}=v^{2}\,\text{d}v\,\text{d}\unicode[STIX]{x1D70E}$
and use the relation
$\unicode[STIX]{x1D654}^{l}(v)=v^{l}\unicode[STIX]{x1D654}^{l}(\hat{v})$
with
$\boldsymbol{v}=v\hat{v}$
(Weinert Reference Weinert1980). Applying the orthogonality relation of (3.10) for
$k=0$
, and expanding the associated Laguerre polynomials using (3.8), we write
$H_{b}$
as
$$\begin{eqnarray}\displaystyle H_{b} & = & \displaystyle \frac{2n_{b}}{v_{thb}}\mathop{\sum }_{l,k}\mathop{\sum }_{m=0}^{k}\frac{L_{km}^{l}}{\unicode[STIX]{x1D70E}_{k}^{l}}\frac{\unicode[STIX]{x1D654}^{l}(\hat{v})\boldsymbol{\cdot }\boldsymbol{M}_{b}^{lk}}{2l+1}\nonumber\\ \displaystyle & & \displaystyle \times \,\frac{1}{\sqrt{\unicode[STIX]{x03C0}}}\left(\frac{1}{x_{b}^{(l+1)/2}}\int _{0}^{x_{b}}\text{e}^{-x}x^{m+l+1/2}\,\text{d}x+x_{b}^{l/2}\int _{x_{b}}^{\infty }\text{e}^{-x}x^{m}\,\text{d}x\right),\end{eqnarray}$$
where the identity
$$\begin{eqnarray}\frac{(2l-1)!!}{2^{l}(l+1/2)!}=\frac{2}{\sqrt{\unicode[STIX]{x03C0}}}\frac{1}{2l+1},\end{eqnarray}$$
is used to simplify (3.16).
We note that the expression of
$H_{b}$
in (3.16) corresponds to the one in Ji & Held (Reference Ji and Held2006), having replaced the
$\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})$
tensors by the
$\unicode[STIX]{x1D64B}^{l}(\boldsymbol{v})$
tensors defined by the recursion relation (see equation (14) of Ji & Held (Reference Ji and Held2006))
$$\begin{eqnarray}\unicode[STIX]{x1D64B}^{l+1}(\boldsymbol{v})=\boldsymbol{v}\unicode[STIX]{x1D64B}^{l}(\boldsymbol{v})-\frac{v^{2}}{2l+1}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{v}}\unicode[STIX]{x1D64B}^{l}(\boldsymbol{v}),\end{eqnarray}$$
with
$\unicode[STIX]{x1D64B}^{0}(\boldsymbol{v})=1$
and
$\unicode[STIX]{x1D64B}^{1}(\boldsymbol{v})=\boldsymbol{v}$
. We can indeed prove that
$\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})=\unicode[STIX]{x1D64B}^{l}(\boldsymbol{v})$
by deriving the tensor
$\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})$
using (3.2). This yields
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{v}}\unicode[STIX]{x1D654}^{l}(\boldsymbol{v}) & = & \displaystyle \frac{(-1)^{l}}{(2l-1)!!}\left[(2l+1)v^{2l-1}\boldsymbol{v}\frac{\unicode[STIX]{x2202}^{l}}{\unicode[STIX]{x2202}\boldsymbol{v}^{l}}\frac{1}{v}+v^{2l+1}\frac{\unicode[STIX]{x2202}^{l+1}}{\unicode[STIX]{x2202}\boldsymbol{v}^{2l+1}}\frac{1}{v}\right]\nonumber\\ \displaystyle & = & \displaystyle \frac{2l+1}{v^{2}}\left[\boldsymbol{v}\frac{v^{2l+1}(-1)^{l}}{(2l-1)!!}\frac{\unicode[STIX]{x2202}^{l}}{\unicode[STIX]{x2202}\boldsymbol{v}^{l}}\frac{1}{v}-\frac{(-1)^{l+1}v^{2(l+1)+1}}{(2l+1)!!}\frac{\unicode[STIX]{x2202}^{l+1}}{\unicode[STIX]{x2202}\boldsymbol{v}^{l+1}}\frac{1}{v}\right]\nonumber\\ \displaystyle & = & \displaystyle \frac{2l+1}{v^{2}}[\boldsymbol{v}\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})-\unicode[STIX]{x1D654}^{l+1}(\boldsymbol{v})],\end{eqnarray}$$
which is the same recursion relation present in (3.18). Since
$\unicode[STIX]{x1D654}^{0}(\boldsymbol{v})=\unicode[STIX]{x1D64B}^{0}(\boldsymbol{v})$
and
$\unicode[STIX]{x1D654}^{1}(\boldsymbol{v})=\unicode[STIX]{x1D64B}^{1}(\boldsymbol{v})$
, the proof is complete.
The integrals in (3.16) can be put in terms of upper
$$\begin{eqnarray}I_{+}^{k}=\frac{1}{\sqrt{\unicode[STIX]{x03C0}}}\int _{0}^{x_{b}}\,\text{d}x\text{e}^{-x}x^{(k-1)/2},\end{eqnarray}$$
and lower
$$\begin{eqnarray}I_{-}^{k}=\frac{1}{\sqrt{\unicode[STIX]{x03C0}}}\int _{x_{b}}^{\infty }\,\text{d}x\text{e}^{-x}x^{(k-1)/2},\end{eqnarray}$$
incomplete gamma functions (Abramowitz et al. Reference Abramowitz, Stegun and Miller1965), yielding
$$\begin{eqnarray}H_{b}=\frac{2n_{b}}{v_{thb}}\mathop{\sum }_{l,k}\mathop{\sum }_{m=0}^{k}\frac{L_{km}^{l}}{\unicode[STIX]{x1D70E}_{k}^{l}}\frac{\unicode[STIX]{x1D654}^{l}(\hat{v})\boldsymbol{\cdot }\boldsymbol{M}_{b}^{lk}}{2l+1}\left(\frac{I_{+}^{2l+2m+2}}{x_{b}^{(l+1)/2}}+x_{b}^{l/2}I_{-}^{2m+1}\right).\end{eqnarray}$$
A procedure similar to the one used to obtain (3.22) can be followed for the second Rosenbluth potential
$G_{b}$
by expanding the distribution function
$f_{b}$
appearing in
$G_{b}$
according to (3.12), therefore obtaining
$$\begin{eqnarray}\displaystyle G_{b} & = & \displaystyle \frac{2n_{b}}{v_{thb}}\mathop{\sum }_{l,k}\mathop{\sum }_{m=0}^{k}\frac{L_{km}^{l}}{\unicode[STIX]{x1D70E}_{k}^{l}}\frac{\unicode[STIX]{x1D654}^{l}(\hat{v})\boldsymbol{\cdot }\boldsymbol{M}_{b}^{lk}}{2l+1}\left[\frac{1}{2l+3}\left(\frac{I_{+}^{2l+2m+4}}{x_{b}^{(l+1)/2}}+x_{b}^{l/2+1}I_{-}^{2m+1}\right)\right.\nonumber\\ \displaystyle & & \displaystyle \left.-\frac{1}{2l-1}\left(\frac{I_{+}^{2l+2m+2}}{x_{b}^{(l-1)/2}}+x_{b}^{l/2}I_{-}^{2m+3}\right)\right].\end{eqnarray}$$
Having derived a closed-form expression for the Rosenbluth potentials, we now turn to the full Coulomb collision operator. We first note that, although the Rosenbluth potentials
$H_{b}$
and
$G_{b}$
are linear functions of
$f_{b}$
, the Coulomb collision operator is, in fact, bilinear in
$f_{a}$
and
$f_{b}$
. In order to rewrite the Coulomb collision operator in (2.2) in terms of a single spherical harmonic tensor
$\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})$
, we make use of the following identity between symmetric traceless tensors (Ji & Held Reference Ji and Held2009)
$$\begin{eqnarray}[\unicode[STIX]{x1D654}^{l-u}(\hat{v})\boldsymbol{\cdot }\boldsymbol{M}_{a}^{lk}]\boldsymbol{\cdot }^{u}[\unicode[STIX]{x1D654}^{n-u}(\hat{v})\boldsymbol{\cdot }\boldsymbol{M}_{a}^{nk}]=\mathop{\sum }_{j=0}^{\min (l,n)-u}d_{j}^{l-u,n-u}\unicode[STIX]{x1D654}^{l+n-2(j+u)}(\hat{v})\boldsymbol{\cdot }(\boldsymbol{M}_{a}^{lk}\boldsymbol{\cdot }^{j+u}\boldsymbol{M}_{b}^{nq})_{TS},\end{eqnarray}$$
where
$\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}$
). The
$d_{j}^{l,n}$
coefficient can be written in terms of the
$t_{j}^{l,n}$
coefficient
$$\begin{eqnarray}t_{j}^{l,n}=\frac{l!n!(-2)^{j}(2l+2n-2j)!(l+n)!}{(2l+2n)!j!(l-j)!(n-j)!(l+n-j)!},\end{eqnarray}$$
as
$$\begin{eqnarray}d_{j}^{l,n}=\mathop{\sum }_{j_{k}|\mathop{\sum }_{k=1}^{h}j_{k}=j}(-1)^{h}\mathop{\prod }_{k=1}^{h}t_{j_{k}}^{l-\mathop{\sum }_{g=1}^{k-1}j_{g},n-\mathop{\sum }_{g=1}^{k-1}j_{g}}.\end{eqnarray}$$
In (3.26), the summation range involves the integer partitions of
$j$
, i.e. the decomposition of
$j$
into different sums of
$h$
positive integers, here labelled as
$j_{k}$
, with
$k$
ranging from
$1$
to
$h$
(e.g. for
$j=3$
, we obtain for
$h=2$
the terms
$j_{1}=2$
and
$j_{2}=1$
, and for
$h=3$
we obtain
$j_{1}=j_{2}=j_{3}=1$
). Expanding
$f_{a}$
and
$f_{b}$
using (3.12), the expression for the Rosenbluth potentials in (3.22) and (3.23) and the identity in (3.24), the collision operator in (2.2) can be rewritten in terms of products of
$\boldsymbol{M}_{a}^{lk}$
and
$\boldsymbol{M}_{b}^{lk}$
as
$$\begin{eqnarray}C(f_{a},f_{b})=f_{aM}\mathop{\sum }_{l,k,n,q=0}^{\infty }\mathop{\sum }_{m=0}^{k}\mathop{\sum }_{r=0}^{q}\frac{L_{km}^{l}}{\unicode[STIX]{x1D70E}_{k}^{l}}\frac{L_{qr}^{n}}{\unicode[STIX]{x1D70E}_{q}^{n}}c_{ab}^{lkmnqr},\end{eqnarray}$$
with
$$\begin{eqnarray}c_{ab}^{lkmnqr}=\mathop{\sum }_{u=0}^{\min (2,l,n)}\unicode[STIX]{x1D708}_{\ast abu}^{lm,nr}(v^{2})\mathop{\sum }_{i=0}^{\min (l,n)-u}d_{i}^{l-u,n-u}\unicode[STIX]{x1D654}^{l+n-2(i+u)}(\hat{v})\boldsymbol{\cdot }(\boldsymbol{M}_{a}^{lk}\boldsymbol{\cdot }^{i+u}\boldsymbol{M}_{b}^{nq})_{TS}.\end{eqnarray}$$
The quantity
$\unicode[STIX]{x1D708}_{\ast abu}^{lm,nr}$
consists of a linear combination of
$I_{+}^{l}$
and
$I_{-}^{l}$
integrals and its derivatives, which can be written as linear combinations of the error function and its derivatives. Their expressions are reported in Ji & Held (Reference Ji and Held2009).
4 Gyrokinetic Coulomb collision operator
In § 3, the Coulomb collision operator is cast in terms of velocity moments of the multipole expansion of the particle distribution functions
$f_{a}$
and
$f_{b}$
. We now express it in terms of the gyrokinetic distribution functions
$\langle \overline{F_{a}}\rangle$
and
$\langle \overline{F_{b}}\rangle$
. As a first step, the gyroangle dependence of the basis functions
$\unicode[STIX]{x1D654}^{lk}$
is found explicitly by using a coordinate transformation from the particle phase-space coordinates
$(\boldsymbol{x},\boldsymbol{v})$
to the guiding-centre coordinate system
$\boldsymbol{Z}$
. This allows us to decouple the fast gyromotion time associated with the gyroangle
$\unicode[STIX]{x1D703}$
from the typical turbulence time scales. The multipole moments
$\boldsymbol{M}_{a}^{lk}$
and
$\boldsymbol{M}_{b}^{lk}$
can then be written in terms of moments of the guiding-centre distribution function
$\langle F_{a}\rangle$
and
$\langle F_{b}\rangle$
for arbitrary values of
$k_{\bot }\unicode[STIX]{x1D70C}_{s}$
. As a second step, the gyrocentre coordinate system
$\overline{\boldsymbol{Z}}$
is introduced by using the coordinate transformation
$T$
in (2.13). As shown in § 2, for a gyrokinetic equation up to second-order accurate in
$\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D6FF}}$
, only the lowest-order collision operator
$C_{0}$
needs to be retained. This allows us to straightforwardly obtain the gyrokinetic collision operator from the guiding-centre one by a simple coordinate relabelling.
We first derive the polar and azimuthal angle (gyroangle) dependence of the
$\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})$
tensor in terms of scalar spherical harmonics. This is useful to analytically perform the gyroaverage of the collision operator in the Boltzmann equation, equation (2.30). For this purpose, as a first step, we show that the Laplacian of
$\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})$
vanishes, i.e. that
$\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})$
are harmonic tensors. By applying the operator
$\unicode[STIX]{x1D6FB}_{\boldsymbol{v}}^{2}$
to
$\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})$
, and recalling that
$\unicode[STIX]{x1D6FB}_{\boldsymbol{v}}^{2}(1/v)=0$
for
$v\not =0$
, we obtain
$$\begin{eqnarray}\unicode[STIX]{x1D6FB}_{\boldsymbol{v}}^{2}\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})=\frac{2(-1)^{l}(2l+1)v^{2l+1}}{(2l-1)!!}\left[(l+1)\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{v}}\right)^{l}\frac{1}{v}+\boldsymbol{v}\boldsymbol{\cdot }\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{v}}\right)^{l+1}\frac{1}{v}\right]=0,\end{eqnarray}$$
since
$$\begin{eqnarray}\boldsymbol{v}\boldsymbol{\cdot }\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{v}}\right)^{l+1}\frac{1}{v}=-(l+1)\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{v}}\right)^{l}\frac{1}{v},\end{eqnarray}$$
as can be proved by induction (Weinert Reference Weinert1980). The angular dependence of
$\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})$
can now be found by expressing the Laplacian of (4.1) in spherical coordinates. Using the fact that
$\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})=v^{l}\unicode[STIX]{x1D654}^{l}(\hat{v})$
, we obtain
$$\begin{eqnarray}\displaystyle 0 & = & \displaystyle \unicode[STIX]{x1D6FB}_{\boldsymbol{v}}^{2}\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})=\unicode[STIX]{x1D6FB}_{\boldsymbol{v}}^{2}[v^{l}\unicode[STIX]{x1D654}^{l}(\hat{v})]\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D654}^{l}(\hat{v})\left(\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}v^{2}}+\frac{2}{v}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}v}\right)v^{l}-v^{l-2}L^{2}\unicode[STIX]{x1D654}^{l}(\hat{v}),\end{eqnarray}$$
where
$L^{2}$
is the angular part of the operator
$\unicode[STIX]{x1D6FB}_{\boldsymbol{v}}^{2}$
multiplied by
$v^{2}$
$$\begin{eqnarray}L^{2}=\frac{1}{\sin \unicode[STIX]{x1D711}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D711}}\left(\sin \unicode[STIX]{x1D711}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D711}}\right)+\frac{1}{\sin \unicode[STIX]{x1D711}^{2}}\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}^{2}},\end{eqnarray}$$
with
$\unicode[STIX]{x1D711}$
and
$\unicode[STIX]{x1D703}$
chosen, respectively, as the pitch angle and the gyroangle variables, both defined in (2.15). Evaluating the
$v$
derivatives in (4.3), the following differential equation for
$\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})$
is obtained
$$\begin{eqnarray}L^{2}\unicode[STIX]{x1D654}^{l}(\hat{v})=l(l+1)\unicode[STIX]{x1D654}^{l}(\hat{v}).\end{eqnarray}$$
We identify (4.5) as the eigenvalue equation for the scalar spherical harmonics
$Y_{lm}(\unicode[STIX]{x1D711},\unicode[STIX]{x1D703})$
(Arfken, Weber & Harris Reference Arfken, Weber and Harris2013), which can be written in terms of associated Legendre polynomials
$P_{l}^{m}(\cos \unicode[STIX]{x1D711})$
as (Abramowitz et al.
Reference Abramowitz, Stegun and Miller1965)
$$\begin{eqnarray}Y_{lm}(\unicode[STIX]{x1D711},\unicode[STIX]{x1D703})=(-1)^{m}\sqrt{\frac{(2l+1)}{4\unicode[STIX]{x03C0}}\frac{(l-m)!}{(l+m)!}}P_{l}^{m}(\cos \unicode[STIX]{x1D711})\text{e}^{\text{i}m\unicode[STIX]{x1D703}},\end{eqnarray}$$
with
$$\begin{eqnarray}P_{l}^{m}(x)=(1-x^{2})^{m/2}\frac{\text{d}^{m}}{\text{d}x^{m}}[P_{l}(x)],\end{eqnarray}$$
and
$P_{l}(x)=(\text{d}^{l}/\text{d}x^{l})[(x^{2}-1)^{l}]/(2^{l}l!)$
a Legendre polynomial. Therefore, using (4.5), and denoting
$\unicode[STIX]{x1D65A}^{lm}$
the basis elements of
$\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})$
(an elementary derivation of the basis tensors
$\unicode[STIX]{x1D65A}^{lm}$
is shown in appendix A), we write
$\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})$
as
$$\begin{eqnarray}\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})=v^{l}\sqrt{\frac{2\unicode[STIX]{x03C0}^{3/2}l!}{2^{l}(l+1/2)!}}\mathop{\sum }_{m=-l}^{l}Y_{lm}(\unicode[STIX]{x1D711},\unicode[STIX]{x1D703})\unicode[STIX]{x1D65A}^{lm}.\end{eqnarray}$$
Having derived the gyroangle dependence of the
$\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})$
tensors, we now compute the fluid moments
$\boldsymbol{M}_{a}^{lk}$
in terms of
$v_{\Vert }$
and
$\unicode[STIX]{x1D707}$
moments of the guiding-centre distribution function
$\langle F_{a}\rangle$
. In order to perform the velocity integration in the definition of the moments
$\boldsymbol{M}^{lk}$
in (3.14) at arbitrary
$k_{\bot }\unicode[STIX]{x1D70C}$
in guiding-centre phase-space coordinates, we use the identity
$f(\boldsymbol{x},\boldsymbol{v})=\int f(\boldsymbol{x}^{\prime },\boldsymbol{v})\unicode[STIX]{x1D6FF}(\boldsymbol{x}-\boldsymbol{x}^{\prime })\,\text{d}\boldsymbol{x}^{\prime }$
. By imposing
$\boldsymbol{x}^{\prime }=\boldsymbol{R}+\unicode[STIX]{x1D746}$
, writing the volume element in phase space as
$\text{d}\boldsymbol{x}^{\prime }\,\text{d}\boldsymbol{v}=(B_{\Vert }^{\ast }/m)\,\text{d}\boldsymbol{R}\,\text{d}v_{\Vert }\,\text{d}\unicode[STIX]{x1D707}\,\text{d}\unicode[STIX]{x1D703}$
and using (2.14), we obtain
$$\begin{eqnarray}n_{a}\boldsymbol{M}_{a}^{lk}(\boldsymbol{x})=\int F_{a}(\boldsymbol{R},v_{\Vert },\unicode[STIX]{x1D707},\unicode[STIX]{x1D703})\unicode[STIX]{x1D654}^{lk}\left(\frac{\boldsymbol{v}}{v_{tha}}\right)\unicode[STIX]{x1D6FF}(\boldsymbol{x}-\boldsymbol{R}-\unicode[STIX]{x1D746}_{a})\frac{B_{\Vert }^{\ast }}{m}\,\text{d}\boldsymbol{R}\,\text{d}v_{\Vert }\,\text{d}\unicode[STIX]{x1D707}\,\text{d}\unicode[STIX]{x1D703},\end{eqnarray}$$
from (3.14). Expressing
$\boldsymbol{v}=\boldsymbol{v}(\boldsymbol{Z})$
, as shown by (2.15), and performing the integral over
$\boldsymbol{R}$
in (4.9), it follows that
$$\begin{eqnarray}n_{a}\boldsymbol{M}_{a}^{lk}(\boldsymbol{x})=\int F_{a}(\boldsymbol{x}-\unicode[STIX]{x1D746}_{a},v_{\Vert },\unicode[STIX]{x1D707},\unicode[STIX]{x1D703})\unicode[STIX]{x1D654}^{lk}\left[\frac{\boldsymbol{v}(\boldsymbol{x}-\unicode[STIX]{x1D746}_{a},v_{\Vert },\unicode[STIX]{x1D707},\unicode[STIX]{x1D703})}{v_{tha}}\right]\frac{B_{\Vert }^{\ast }}{m}\,\text{d}v_{\Vert }\,\text{d}\unicode[STIX]{x1D707}\,\text{d}\unicode[STIX]{x1D703}.\end{eqnarray}$$
The orderings in (2.26) and (2.27) for the guiding-centre distribution function
$F_{a}$
allow us to approximate
$F_{a}\simeq \langle F_{a}\rangle _{\boldsymbol{R}}$
(Jorge et al.
Reference Jorge, Ricci and Loureiro2017), effectively neglecting
$\unicode[STIX]{x1D716}^{2}$
effects in
$\boldsymbol{M}_{a}^{lk}$
, hence in the collision operator
$C(f_{a},f_{b})$
. To make further analytical progress, and in line with previous gyrokinetic literature (Li & Ernst Reference Li and Ernst2011; Pan & Ernst Reference Pan and Ernst2019), we represent
$F_{a}(\boldsymbol{R},v_{\Vert },\unicode[STIX]{x1D707},\unicode[STIX]{x1D703})$
by its Fourier transform
$F_{a}(\boldsymbol{k},v_{\Vert },\unicode[STIX]{x1D707},\unicode[STIX]{x1D703})=\int F_{a}(\boldsymbol{R},v_{\Vert },\unicode[STIX]{x1D707},\unicode[STIX]{x1D703})\text{e}^{-\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{R}}\,\text{d}\boldsymbol{R}$
, and write
$$\begin{eqnarray}n_{a}\boldsymbol{M}_{a}^{lk}(\boldsymbol{x})=\int \langle F_{a}(\boldsymbol{k},v_{\Vert },\unicode[STIX]{x1D707},\unicode[STIX]{x1D703})\rangle _{\boldsymbol{R}}\unicode[STIX]{x1D654}^{lk}\left[\frac{\boldsymbol{v}(\boldsymbol{x}-\unicode[STIX]{x1D746}_{a},v_{\Vert },\unicode[STIX]{x1D707},\unicode[STIX]{x1D703})}{v_{tha}}\right]\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}}\text{e}^{-\text{i}\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D70C}_{a}}\frac{B_{\Vert }^{\ast }}{m}\,\text{d}\boldsymbol{k}\,\text{d}v_{\Vert }\,\text{d}\unicode[STIX]{x1D707}\,\text{d}\unicode[STIX]{x1D703}.\end{eqnarray}$$
By aligning the
$\boldsymbol{k}$
coordinate system in the integral of (4.11) with the axes
$(\boldsymbol{b},\boldsymbol{e}_{1},\boldsymbol{e}_{2})$
, i.e.
$\boldsymbol{k}=k_{\Vert }\boldsymbol{b}+k_{\bot }(\cos \unicode[STIX]{x1D703}\boldsymbol{e}_{1}+\sin \unicode[STIX]{x1D703}\boldsymbol{e}_{2})$
, we write
$\exp (-\text{i}\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D746})=\exp (-\text{i}k_{\bot }\unicode[STIX]{x1D70C}\cos \unicode[STIX]{x1D703})$
. We then use the Jacobi–Anger expansion (Andrews Reference Andrews1992)
$$\begin{eqnarray}\text{e}^{-\text{i}k_{\bot }\unicode[STIX]{x1D70C}\cos \unicode[STIX]{x1D703}}=\mathop{\sum }_{p=-\infty }^{\infty }(-i)^{p}J_{p}(k_{\bot }\unicode[STIX]{x1D70C})\text{e}^{-\text{i}p\unicode[STIX]{x1D703}},\end{eqnarray}$$
with
$J_{p}$
the Bessel function of order
$p$
, and rewrite (4.11) as
$$\begin{eqnarray}\displaystyle n_{a}\boldsymbol{M}_{a}^{lk}(\boldsymbol{x}) & = & \displaystyle \mathop{\sum }_{p=-\infty }^{\infty }(-1)^{p}\int J_{p}(k_{\bot }\unicode[STIX]{x1D70C})\langle F_{a}(\boldsymbol{k},v_{\Vert },\unicode[STIX]{x1D707},\unicode[STIX]{x1D703})\rangle _{\boldsymbol{R}}\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}}\nonumber\\ \displaystyle & & \displaystyle \times \,\unicode[STIX]{x1D654}^{lk}\left[\frac{\boldsymbol{v}(\boldsymbol{x}-\unicode[STIX]{x1D746}_{a},v_{\Vert },\unicode[STIX]{x1D707},\unicode[STIX]{x1D703})}{v_{tha}}\right]\text{e}^{-\text{i}p\unicode[STIX]{x1D703}}\frac{B_{\Vert }^{\ast }}{m}\,\text{d}\boldsymbol{k}\,\text{d}v_{\Vert }\,\text{d}\unicode[STIX]{x1D707}\,\text{d}\unicode[STIX]{x1D703}.\end{eqnarray}$$
The velocity
$\boldsymbol{v}$
in the argument of
$\unicode[STIX]{x1D654}^{lk}$
in (4.13) is then expanded as
$$\begin{eqnarray}\boldsymbol{v}(\boldsymbol{x}-\unicode[STIX]{x1D746}_{a},v_{\Vert },\unicode[STIX]{x1D707},\unicode[STIX]{x1D703})=\boldsymbol{v}(\boldsymbol{x},v_{\Vert },\unicode[STIX]{x1D707},\unicode[STIX]{x1D703})+O(\unicode[STIX]{x1D716}_{B}).\end{eqnarray}$$
The second term in (4.14) introduces
$\unicode[STIX]{x1D716}_{B}\ll \unicode[STIX]{x1D716}$
terms in the collision operator and is therefore neglected. An example of a numerical implementation using a similar Fourier representation can be found in Pan & Ernst (Reference Pan and Ernst2019).
Using (4.8) to express
$\unicode[STIX]{x1D654}^{lk}(\boldsymbol{v})$
in terms of spherical harmonics, we perform the gyroangle integration in (4.13). By rewriting the spherical harmonics
$Y_{lm}(\unicode[STIX]{x1D711},0)$
in terms of associated Legendre polynomials
$P_{l}^{m}(\cos \unicode[STIX]{x1D711})$
using (4.6), the gyroaverage of the product
$Y^{lk}(\boldsymbol{v}/v_{tha})\text{e}^{-\text{i}p\unicode[STIX]{x1D703}}$
can be performed, yielding
$$\begin{eqnarray}\displaystyle n_{a}\boldsymbol{M}_{a}^{lk}(\boldsymbol{x}) & = & \displaystyle \mathop{\sum }_{p=-\infty }^{\infty }(-1)^{p}\int J_{p}(k_{\bot }\unicode[STIX]{x1D70C}_{a})\langle F_{a}(\boldsymbol{k},v_{\Vert },\unicode[STIX]{x1D707},\unicode[STIX]{x1D703})\rangle _{\boldsymbol{R}}\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}}\nonumber\\ \displaystyle & & \displaystyle \times \,\left\langle \unicode[STIX]{x1D654}^{lk}\left(\frac{\boldsymbol{v}}{v_{tha}}\right)\text{e}^{-\text{i}p\unicode[STIX]{x1D703}}\right\rangle \frac{B_{\Vert }^{\ast }}{m}\,\text{d}\boldsymbol{k}\,\text{d}v_{\Vert }\,\text{d}\unicode[STIX]{x1D707}2\unicode[STIX]{x03C0},\end{eqnarray}$$
with
$$\begin{eqnarray}\displaystyle \left\langle \unicode[STIX]{x1D654}^{lk}\left(\frac{\boldsymbol{v}}{v_{tha}}\right)\text{e}^{-\text{i}p\unicode[STIX]{x1D703}}\right\rangle & = & \displaystyle L_{k}^{l+1/2}\left(\frac{v^{2}}{v_{tha}^{2}}\right)\left(\frac{v}{v_{tha}}\right)^{l}\sqrt{\frac{\unicode[STIX]{x03C0}^{1/2}l!}{2^{l}(l-1/2)!}}\nonumber\\ \displaystyle & & \displaystyle \times \,\mathop{\sum }_{m=-l}^{l}(-1)^{m}\sqrt{\frac{(l-m)!}{(l+m)!}}P_{l}^{m}(\cos \unicode[STIX]{x1D711})\unicode[STIX]{x1D65A}^{lm}\unicode[STIX]{x1D6FF}_{m,p}.\end{eqnarray}$$
We note that the
$p=0$
case of (4.16) corresponds to the gyroaveraged formulas in Ji, Held & Sovinec (Reference Ji, Held and Sovinec2009), Ji, Held & Jhang (Reference Ji, Held and Jhang2013), Ji & Held (Reference Ji and Held2014) used to derive closures for fluid models at zeroth order in
$\unicode[STIX]{x1D716}$
. Finally, by defining the Bessel–Fourier operator
$$\begin{eqnarray}j_{m}[F_{a}]\equiv \int J_{m}(k_{\bot }\unicode[STIX]{x1D70C}_{a})\langle F_{a}(\boldsymbol{k},v_{\Vert },\unicode[STIX]{x1D707},\unicode[STIX]{x1D703})\rangle _{\boldsymbol{R}}\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}}\,\text{d}\boldsymbol{k},\end{eqnarray}$$
the expression for the fluid moments
$\boldsymbol{M}_{a}^{lk}$
in terms of coupled
$v_{\Vert }$
and
$\unicode[STIX]{x1D707}$
moments of the guiding-centre distribution function
$\langle F_{a}\rangle _{\boldsymbol{R}}$
is obtained
$$\begin{eqnarray}n_{a}\boldsymbol{M}_{a}^{lk}(\boldsymbol{x})=\sqrt{\frac{2\unicode[STIX]{x03C0}^{5/2}l!}{2^{l}(l+1/2)!}}\mathop{\sum }_{m=-l}^{l}\boldsymbol{e}^{lm}(-1)^{m}{\mathcal{M}}_{am}^{lk}(\boldsymbol{x}),\end{eqnarray}$$
with
$$\begin{eqnarray}{\mathcal{M}}_{am}^{lk}(\boldsymbol{x})=\int j_{m}[F_{a}]v^{l}L_{k}^{l+1/2}(x_{a}^{2})Y_{lm}(\unicode[STIX]{x1D711},0)\frac{B_{\Vert }^{\ast }}{m}\,\text{d}v_{\Vert }\,\text{d}\unicode[STIX]{x1D707}.\end{eqnarray}$$
Equation (4.18) can now be used to express the collision operator
$C(f_{a},f_{b})$
in terms of
$v_{\Vert }$
and
$\unicode[STIX]{x1D707}$
integrals of
$\langle F_{a}\rangle$
. Using (3.27), (3.28) and (4.8), and defining
$$\begin{eqnarray}E_{jv}^{lsnt}=\boldsymbol{e}^{l+n-2j\,v}\boldsymbol{\cdot }(\boldsymbol{e}^{ls}\boldsymbol{\cdot }^{j}\boldsymbol{e}^{nt})_{TS},\end{eqnarray}$$
we can write the collision operator in (3.27) and (3.28) as a function of the
${\mathcal{M}}_{am}^{lk}$
moments, i.e.
$$\begin{eqnarray}\displaystyle c_{ab}^{lkmnqr} & = & \displaystyle \mathop{\sum }_{u=0}^{\min (2,l,n)}\mathop{\sum }_{j=0}^{\min (l,n)-u}d_{j}^{l-u,n-u}a_{j+u}^{ln}\mathop{\sum }_{s=-l}^{l}\mathop{\sum }_{t=-n}^{n}\mathop{\sum }_{v=-(l+n-2j-2u)}^{l+n-2j-2u}E_{j+u~v}^{lsnt}\nonumber\\ \displaystyle & & \displaystyle \times \,Y_{l+n-2j-2u~v}(\unicode[STIX]{x1D711},\unicode[STIX]{x1D703})\frac{\unicode[STIX]{x1D708}_{\ast abu}^{lm,nr}(v^{2})}{n_{a}n_{b}}{\mathcal{M}}_{as}^{lk}(\boldsymbol{x}){\mathcal{M}}_{bt}^{nq}(\boldsymbol{x}),\end{eqnarray}$$
with
$$\begin{eqnarray}a_{j}^{ln}=\frac{1}{2^{l+n-j}}\sqrt{\frac{8\unicode[STIX]{x03C0}^{13/2}l!n!(l+n-2j)!}{(l+1/2)!(n+1/2)!(l+n-2j+1/2)!}}.\end{eqnarray}$$
We now focus on the gyroaverage of the collision operator in (4.21) with the gyroaverage operation performed at constant
$\boldsymbol{R}$
. We first note that the gyroangle
$\unicode[STIX]{x1D703}$
dependence in
$c_{ab}^{lkmnqr}$
is present through the spherical harmonic
$Y_{l+n-2j-2u\,v}(\unicode[STIX]{x1D711},\unicode[STIX]{x1D703})$
and through the fluid moments
${\mathcal{M}}_{as}^{lk}$
and
${\mathcal{M}}_{bt}^{nq}$
as the latter are functions of
$\boldsymbol{x}=\boldsymbol{R}+\unicode[STIX]{x1D746}_{a}$
. To make the gyroangle dependence explicit, we write both
${\mathcal{M}}_{as}^{lk}$
and
${\mathcal{M}}_{bt}^{nq}$
in Fourier space as
$$\begin{eqnarray}{\mathcal{M}}_{as}^{lk}(\boldsymbol{x}){\mathcal{M}}_{bt}^{nq}(\boldsymbol{x})=\int \text{d}\boldsymbol{k}\,\text{d}\boldsymbol{k}^{\prime }\text{e}^{\text{i}(\boldsymbol{k}+\boldsymbol{k}^{\prime })\boldsymbol{\cdot }\boldsymbol{R}}{\mathcal{M}}_{as}^{lk}(\boldsymbol{k}){\mathcal{M}}_{bt}^{nq}(\boldsymbol{k}^{\prime })\text{e}^{\text{i}(\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D746}_{a}+\boldsymbol{k}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D746}_{a})}.\end{eqnarray}$$
Using the Jacobi–Anger expansion of (4.12), we find that
$$\begin{eqnarray}\displaystyle & & \displaystyle \langle Y_{lm}(\unicode[STIX]{x1D711},\unicode[STIX]{x1D703}){\mathcal{M}}_{as}^{lk}(\boldsymbol{x}){\mathcal{M}}_{bt}^{nq}(\boldsymbol{x})\rangle _{\boldsymbol{R}}=\int \text{d}\boldsymbol{k}\,\text{d}\boldsymbol{k}^{\prime }\text{e}^{\text{i}(\boldsymbol{k}+\boldsymbol{k}^{\prime })\boldsymbol{\cdot }\boldsymbol{R}}{\mathcal{M}}_{as}^{lk}(\boldsymbol{k}){\mathcal{M}}_{bt}^{nq}(\boldsymbol{k}^{\prime })\nonumber\\ \displaystyle & & \displaystyle \qquad \times \,\sqrt{\frac{2l+1}{4\unicode[STIX]{x03C0}}\frac{(l-m)!}{(l+m)!}}P_{l}^{m}(\cos \unicode[STIX]{x1D711})\mathop{\sum }_{p=-\infty }^{\infty }(-1)^{p+m}\text{e}^{\text{i}(p+m)\unicode[STIX]{x1D6FC}}J_{p}(k_{\bot }\unicode[STIX]{x1D70C}_{a})J_{p+m}(k_{\bot }^{\prime }\unicode[STIX]{x1D70C}_{a}),\end{eqnarray}$$
with
$\unicode[STIX]{x1D6FC}$
the azimuthal angle of the
$\boldsymbol{k}^{\prime }$
vector, i.e. the angle between
$\boldsymbol{k}_{\bot }^{\prime }$
and
$\boldsymbol{k}_{\bot }$
. The gyroaveraged collision operator at arbitrary
$k_{\bot }\unicode[STIX]{x1D70C}$
is therefore given by
$$\begin{eqnarray}\langle C(F_{a},F_{b})\rangle _{\boldsymbol{R}}=f_{aM}\mathop{\sum }_{l,k,n,q=0}^{\infty }\mathop{\sum }_{m=0}^{k}\mathop{\sum }_{r=0}^{q}L_{km}^{l}L_{qr}^{n}\langle c_{ab}^{lkmnqr}\rangle _{\boldsymbol{ R}},\end{eqnarray}$$
with
$$\begin{eqnarray}\displaystyle \langle c_{ab}^{lkmnqr}\rangle _{\boldsymbol{R}} & = & \displaystyle \mathop{\sum }_{u=0}^{\min (2,l,n)}\mathop{\sum }_{j=0}^{\min (l,n)-u}d_{j}^{l-u,n-u}a_{j+u}^{ln}\mathop{\sum }_{s=-l}^{l}\mathop{\sum }_{t=-n}^{n}\mathop{\sum }_{v=-(l+n-2j-2u)}^{l+n-2j-2u}E_{j+u~v}^{lsnt}\nonumber\\ \displaystyle & & \displaystyle \times \,b_{j+u}^{l+nv}P_{l+n-2j-2u}^{v}(\cos \unicode[STIX]{x1D711})\unicode[STIX]{x1D708}_{\ast abu}^{lm,nr}(v^{2})\int {\mathcal{M}}_{as}^{lk}(\boldsymbol{k}){\mathcal{M}}_{bt}^{nq}(\boldsymbol{k}^{\prime })\text{e}^{\text{i}(\boldsymbol{k}+\boldsymbol{k}^{\prime })\boldsymbol{\cdot }\boldsymbol{R}}\nonumber\\ \displaystyle & & \displaystyle \times \mathop{\sum }_{p=-\infty }^{\infty }(-1)^{p+v}\text{e}^{\text{i}(p+v)\unicode[STIX]{x1D6FC}}J_{p}(k_{\bot }\unicode[STIX]{x1D70C}_{a})J_{p+v}(k_{\bot }^{\prime }\unicode[STIX]{x1D70C}_{a})\,\text{d}\boldsymbol{k}\,\text{d}\boldsymbol{k}^{\prime },\end{eqnarray}$$
and
$$\begin{eqnarray}b_{j}^{lv}=i^{v}\sqrt{\frac{2l-4j}{4\unicode[STIX]{x03C0}}\frac{(l-2j-v)!}{(l-2j+v)!}}.\end{eqnarray}$$
We note that, if only first-order
$k_{\bot }\unicode[STIX]{x1D70C}$
terms are kept in the Fourier–Bessel operator of (4.17), the collision operator in (4.25) reduces to the drift-kinetic collision operator found in Jorge et al. (Reference Jorge, Ricci and Loureiro2017).
In (4.25), the gyroaveraged collision operator is cast in terms of
$v_{\Vert }$
and
$\unicode[STIX]{x1D707}$
moments of the guiding-centre distribution function
$\langle F_{a}\rangle$
for arbitrary values of
$k_{\bot }\unicode[STIX]{x1D70C}$
. We now apply the transformation
$T$
, introduced in (2.13), to (4.25) in order to write the gyroaveraged collision operator in terms of
$\overline{v}_{\Vert }$
and
$\overline{\unicode[STIX]{x1D707}}$
moments of the gyrocentre distribution function
$\langle \overline{F_{a}}\rangle _{\overline{\boldsymbol{R}}}$
. As shown in § 2, only the zeroth-order terms in the
$\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D6FF}}$
expansion of
$\langle C(\overline{F_{a}},\overline{F_{b}})\rangle _{\overline{\boldsymbol{R}}}$
are needed in order to adequately describe collisional processes in the gyrokinetic framework. Therefore, using equation (2.13), we apply the zeroth-order transformations
$\boldsymbol{Z}\simeq \overline{\boldsymbol{Z}}$
and
$F_{a}(\boldsymbol{Z})=T\overline{F_{a}}(\boldsymbol{Z})\simeq \overline{F_{a}}(\boldsymbol{Z})$
to the collision operator
$\langle C(F_{a},F_{b})\rangle$
in (4.25), yielding
$$\begin{eqnarray}\langle C(\overline{F_{a}},\overline{F_{b}})\rangle _{\overline{\boldsymbol{R}}}\simeq f_{aM}\mathop{\sum }_{l,k,n,q=0}^{\infty }\mathop{\sum }_{m=0}^{k}\mathop{\sum }_{r=0}^{q}L_{km}^{l}L_{qr}^{n}\langle \overline{c}_{ab}^{lkmnqr}\rangle _{\overline{\boldsymbol{R}}},\end{eqnarray}$$
with
$$\begin{eqnarray}\displaystyle \langle \overline{c}_{ab}^{lkmnqr}\rangle _{\overline{\boldsymbol{R}}} & = & \displaystyle \mathop{\sum }_{u=0}^{\min (2,l,n)}\mathop{\sum }_{j=0}^{\min (l,n)-u}d_{j}^{l-u,n-u}a_{j+u}^{ln}\mathop{\sum }_{s=-l}^{l}\mathop{\sum }_{t=-n}^{n}\mathop{\sum }_{v=-(l+n-2j-2u)}^{l+n-2j-2u}E_{j+u~v}^{lsnt}\nonumber\\ \displaystyle & & \displaystyle \times \,b_{j+u}^{l+nv}P_{l+n-2j-2u}^{v}(\overline{v}_{\Vert }/\overline{v})\unicode[STIX]{x1D708}_{\ast abu}^{lm,nr}(\overline{v}^{2})\int \{_{as}^{lk}(\boldsymbol{k})\{_{bt}^{nq}(\boldsymbol{k}^{\prime })\text{e}^{\text{i}(\boldsymbol{k}+\boldsymbol{k}^{\prime })\boldsymbol{\cdot }\overline{\boldsymbol{R}}}\nonumber\\ \displaystyle & & \displaystyle \times \,\mathop{\sum }_{p=-\infty }^{\infty }(-1)^{p+v}\text{e}^{\text{i}(p+v)\unicode[STIX]{x1D6FC}}J_{p}(k_{\bot }\overline{\unicode[STIX]{x1D70C}}_{a})J_{p+v}(k_{\bot }^{\prime }\overline{\unicode[STIX]{x1D70C}}_{a})\,\text{d}\boldsymbol{k}\,\text{d}\boldsymbol{k}^{\prime },\end{eqnarray}$$
where
$\overline{\unicode[STIX]{x1D70C}_{a}}=\unicode[STIX]{x1D70C}_{a}(\overline{\boldsymbol{R}},\overline{\unicode[STIX]{x1D707}},\overline{\unicode[STIX]{x1D703}})$
,
$\overline{v}^{2}=\overline{v}_{\Vert }^{2}+2B\overline{\unicode[STIX]{x1D707}}/m$
and
$$\begin{eqnarray}\overline{{\mathcal{M}}}_{am}^{lk}=\int j_{m}[\overline{F_{a}}]\overline{v}^{l}L_{k}^{l+1/2}(\overline{v}^{2})Y_{lm}(\overline{\unicode[STIX]{x1D711}},0)\frac{B_{\Vert }^{\ast }}{m}\,\text{d}\overline{v}_{\Vert }\,\text{d}\overline{\unicode[STIX]{x1D707}}.\end{eqnarray}$$
The collision operator in (4.28) represents the gyrokinetic full Coulomb collision operator that can be used in gyrokinetic models that are up to
$O(\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D6FF}}^{2})$
accurate. In (4.28), the integral-differential character of the
$C(f_{a},f_{b})$
operator is replaced by two-dimensional integrals of the gyrocentre distribution function over the velocity coordinates
$\overline{v}_{\Vert }$
and
$\overline{\unicode[STIX]{x1D707}}$
. We note that, in practice, a truncation of the series present in (4.28) and (4.29) requires a numerical study in order to assess their convergence.
5 Hermite–Laguerre expansion of the Coulomb operator
In this section, we expand the distribution function into an orthogonal Hermite–Laguerre polynomial basis and compute the Hermite–Laguerre moments of the Coulomb collision operator in (4.28). An expansion of the drift-kinetic (Jorge et al.
Reference Jorge, Ricci and Loureiro2017) and gyrokinetic (Mandell et al.
Reference Mandell, Dorland and Landreman2018; Frei et al.
Reference Frei, Jorge and Ricci2019) equations in Hermite–Laguerre polynomials has been recently introduced, showing that this is an advantageous approach to the study of plasma waves and instabilities (Jorge, Ricci & Loureiro Reference Jorge, Ricci and Loureiro2018; Jorge et al.
Reference Jorge, Ricci, Brunner, Gamba, Konovets, Loureiro, Perrone and Teixeira2019). A key reason for using a basis of Hermite–Laguerre polynomials in gyrokinetics is that these polynomials are orthogonal with respect to a Maxwellian, and can be directly related to the Bessel functions used in evaluating gyroaverage operators such as the ones present in (4.29) and (4.30). We therefore expand
$\langle \overline{F_{a}}\rangle _{\overline{\boldsymbol{R}}}$
as
$$\begin{eqnarray}\langle \overline{F_{a}}\rangle _{\overline{\boldsymbol{R}}}=f_{Ma}\mathop{\sum }_{p,j}\frac{\overline{N}_{a}^{pj}}{\sqrt{2^{p}p!}}H_{p}(\overline{s}_{\Vert a})L_{j}(\overline{s}_{\bot a}^{2}),\end{eqnarray}$$
where
$H_{p}$
are physicists’ Hermite polynomials of order
$p$
defined by the Rodrigues formula
$$\begin{eqnarray}H_{p}(x)=(-1)^{p}\text{e}^{x^{2}}\frac{\text{d}^{p}}{\text{d}x^{p}}\text{e}^{-x^{2}},\end{eqnarray}$$
and normalized via
$$\begin{eqnarray}\int _{-\infty }^{\infty }\,\text{d}xH_{p}(x)H_{p^{\prime }}(x)\text{e}^{-x^{2}}=2^{p}p!\sqrt{\unicode[STIX]{x03C0}}\unicode[STIX]{x1D6FF}_{pp^{\prime }},\end{eqnarray}$$
and
$L_{j}$
the Laguerre polynomials of order
$j$
defined by the Rodrigues formula
$$\begin{eqnarray}L_{j}(x)=\frac{\text{e}^{x}}{j!}\frac{\text{d}^{j}}{\text{d}x^{j}}(\text{e}^{-x}x^{j}),\end{eqnarray}$$
and orthonormal with respect to the weight
$\text{e}^{-x}$
$$\begin{eqnarray}\int _{0}^{\infty }\,\text{d}xL_{j}(x)L_{j^{\prime }}(x)\text{e}^{-x}=\unicode[STIX]{x1D6FF}_{jj^{\prime }}.\end{eqnarray}$$
In (5.1), we introduce the normalized parallel velocity
$$\begin{eqnarray}\overline{s}_{\Vert a}=\frac{\overline{v}_{\Vert }}{v_{tha}},\end{eqnarray}$$
and the perpendicular velocity coordinate
$$\begin{eqnarray}\overline{s}_{\bot a}^{2}=\frac{\overline{\unicode[STIX]{x1D707}}B}{T_{a}}.\end{eqnarray}$$
Due to the orthogonality of the Hermite–Laguerre polynomial basis, the coefficients
$N_{a}^{pj}$
of the expansion in (5.1) can be computed as
$$\begin{eqnarray}\overline{N}_{a}^{pj}=\int \frac{H_{p}(\overline{s}_{\Vert a})L_{j}(\overline{s}_{\bot a}^{2})\langle \overline{F_{a}}\rangle _{\overline{\boldsymbol{R}}}}{\sqrt{2^{p}p!}}\frac{B}{m_{a}}\,\text{d}\overline{v}_{\Vert }\,\text{d}\overline{\unicode[STIX]{x1D707}}\,\text{d}\overline{\unicode[STIX]{x1D703}}.\end{eqnarray}$$
We note that the integrand of
$\overline{N}_{a}^{pj}$
in (5.8) contains the multiplicative factor
$B/m_{a}$
, as opposed to the Jacobian containing the factor
$B_{\Vert }^{\ast }/m_{a}$
. In the following, we also use the Hermite–Laguerre moments of
$\langle \overline{F_{a}}\rangle _{\overline{\boldsymbol{R}}}$
with
$B_{\Vert }^{\ast }$
in the integrand instead of
$B$
. These are denoted as
$\overline{N}_{a}^{\ast pj}$
, i.e.
$$\begin{eqnarray}\displaystyle \overline{N}_{a}^{\ast pj} & = & \displaystyle \int \frac{H_{p}(\overline{s}_{\Vert a})L_{j}(\overline{s}_{\bot a}^{2})\langle \overline{F_{a}}\rangle _{\overline{\boldsymbol{R}}}}{\sqrt{2^{p}p!}}\frac{B_{\Vert }^{\ast }}{m_{a}}\,\text{d}\overline{v}_{\Vert }\,\text{d}\overline{\unicode[STIX]{x1D707}}\,\text{d}\overline{\unicode[STIX]{x1D703}}\nonumber\\ \displaystyle & = & \displaystyle \overline{N}_{a}^{pj}\left(1+\frac{\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \boldsymbol{v}_{\boldsymbol{E}}}{\unicode[STIX]{x1D6FA}_{a}}\right)+\frac{v_{tha}}{\sqrt{2}\unicode[STIX]{x1D6FA}_{a}}\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \boldsymbol{b}(\sqrt{p+1}\overline{N}_{a}^{p+1\,j}+\sqrt{p}\overline{N}_{a}^{p-1\,j}),\qquad\end{eqnarray}$$
while in the weak-flow regime the term proportional to
$\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \boldsymbol{v}_{\boldsymbol{E}}/\unicode[STIX]{x1D6FA}_{a}$
in (5.9) is set to zero.
In order to express the collision operator in terms of the moments
$\overline{N}_{a}^{pj}$
given in (5.8) and to evaluate its Hermite–Laguerre moments, we first consider the gyrokinetic moments
$\overline{{\mathcal{M}}}_{am}^{lk}$
and write the integral that defines them in (4.30) as a function of the gyrocentre moments
$\overline{N}_{a}^{pj}$
of (5.8). As a first step, we project both the Fourier–Bessel operator
$j_{m}[\overline{F_{a}}]$
and the spherical harmonics
$Y_{lm}$
on the Hermite–Laguerre basis. We remark that the
$\overline{\unicode[STIX]{x1D707}}$
and
$k_{\bot }$
dependence in the Fourier–Bessel operator
$j_{m}$
, (4.17), can be decomposed by introducing
$\unicode[STIX]{x1D70C}_{tha}=v_{tha}/\unicode[STIX]{x1D6FA}_{a}$
and noting that
$|\unicode[STIX]{x1D746}_{a}|=\sqrt{\overline{\unicode[STIX]{x1D707}}B/T_{a}}\unicode[STIX]{x1D70C}_{tha}=\overline{s}_{\bot a}\unicode[STIX]{x1D70C}_{tha}$
. This allows the use of the following identity between Bessel and Legendre functions (Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik2007)
$$\begin{eqnarray}J_{m}(2b_{a}\overline{s}_{\bot a})=\unicode[STIX]{x1D70E}_{m}b_{a}^{|m|}\overline{s}_{\bot a}^{|m|}\text{e}^{-b_{a}^{2}}\mathop{\sum }_{r=0}^{\infty }\frac{L_{r}^{|m|}(\overline{s}_{\bot a}^{2})}{(|m|+r)!}b_{a}^{2r},\end{eqnarray}$$
with
$b_{a}=k_{\bot }\unicode[STIX]{x1D70C}_{tha}/2$
,
$\unicode[STIX]{x1D70E}_{0}=1$
and
$\unicode[STIX]{x1D70E}_{m}=\text{sgn}(m)^{m}$
for
$m\not =0$
. The Fourier–Bessel operator in (4.17), with the identity in (5.10) and the Hermite–Laguerre expansion of (5.1), can then be written as
$$\begin{eqnarray}j_{m}[\overline{F_{a}}]=f_{Ma}\mathop{\sum }_{p=0}^{\infty }\mathop{\sum }_{j=0}^{\infty }\mathop{\sum }_{r=0}^{\infty }\frac{H_{p}(\overline{s}_{\Vert a})L_{j}(\overline{s}_{\bot a}^{2})}{\sqrt{2^{p}p!}}\frac{L_{r}^{m}(\overline{s}_{\bot a}^{2})\overline{s}_{\bot a}^{m}}{(m+r)!}\int \overline{N}_{a}^{pj}(\boldsymbol{k})b_{a}^{m+2r}\text{e}^{-b_{a}^{2}}\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}}\,\text{d}\boldsymbol{k}.\end{eqnarray}$$
As a second step, we consider
$$\begin{eqnarray}Y_{lm}(\unicode[STIX]{x1D711},0)=(-1)^{m}\sqrt{\frac{2l+1}{4\unicode[STIX]{x03C0}}\frac{(l-m)!}{(l+m)!}}P_{l}^{m}(\cos \unicode[STIX]{x1D711}),\end{eqnarray}$$
which are used in the definition in (4.6). In order to expand the associated Legendre polynomials
$P_{l}^{m}(\cos \unicode[STIX]{x1D711})$
appearing in (5.12) in a Hermite–Laguerre basis, we generalize the basis transformation between a Legendre–associated Laguerre and a Hermite–Laguerre basis presented in Jorge et al. (Reference Jorge, Ricci and Loureiro2017) to a transformation between an associated Legendre–associated Laguerre and a Hermite–Laguerre basis, that is
$$\begin{eqnarray}\frac{\overline{v}^{l}}{v_{tha}^{l}}P_{l}^{m}\left(\frac{\overline{v}_{\Vert }}{\overline{v}}\right)L_{k}^{l+1/2}\left(\frac{\overline{v}^{2}}{v_{tha}^{2}}\right)=\mathop{\sum }_{p=0}^{l+2k}\mathop{\sum }_{j=0}^{k+\lfloor l/2\rfloor }T_{lkm}^{pj}H_{p}\left(\frac{\overline{v}_{\Vert a}}{v_{tha}}\right)L_{j}\left(\frac{\overline{\unicode[STIX]{x1D707}}B}{T_{a}}\right)\left(\frac{\overline{\unicode[STIX]{x1D707}}B}{T_{a}}\right)^{m/2}.\end{eqnarray}$$
For the derivation and expression of the
$T_{lkm}^{pj}$
coefficients, see appendix B. The inverse transformation coefficients
$\left(T^{-1}\right)_{pj}^{lkm}$
are defined as
$$\begin{eqnarray}H_{p}\left(\frac{\overline{v}_{\Vert a}}{v_{tha}}\right)L_{j}\left(\frac{\overline{\unicode[STIX]{x1D707}}B}{T_{a}}\right)\left(\frac{\overline{\unicode[STIX]{x1D707}}B}{T_{a}}\right)^{m/2}=\mathop{\sum }_{l=0}^{p+2j}\mathop{\sum }_{k=0}^{j+\lfloor p/2\rfloor }(T^{-1})_{pj}^{lkm}\frac{\overline{v}^{l}}{v_{tha}^{l}}P_{l}^{m}\left(\frac{\overline{v}_{\Vert }}{\overline{v}}\right)L_{k}^{l+1/2}\left(\frac{\overline{v}^{2}}{v_{tha}^{2}}\right).\end{eqnarray}$$
The gyrocentre moments
$\overline{{\mathcal{M}}}_{am}^{lk}$
in (4.30) can then be rewritten using the identities in (5.11) and (5.13) and
$$\begin{eqnarray}L_{r}^{m}(x)L_{j}(x)x^{m}=\mathop{\sum }_{s=0}^{m+r+j}d_{rjs}^{m}L_{s}(x),\end{eqnarray}$$
with the
$d_{mjs}^{r}$
coefficients given by
$$\begin{eqnarray}d_{rjs}^{m}=\mathop{\sum }_{r_{1}=0}^{r}\mathop{\sum }_{j_{1}=0}^{j}\mathop{\sum }_{s_{1}=0}^{s}L_{rr_{1}}^{-1/2}L_{jj_{1}}^{m-1/2}L_{ss_{1}}^{-1/2}(r_{1}+j_{1}+s_{1}+m)!,\end{eqnarray}$$
yielding the following expression
$$\begin{eqnarray}\displaystyle \overline{{\mathcal{M}}}_{am}^{lk}(\boldsymbol{k})=\mathop{\sum }_{g=0}^{\infty }\mathop{\sum }_{h=0}^{l+2k}\mathop{\sum }_{u=0}^{k+\lfloor l/2\rfloor }\mathop{\sum }_{s=0}^{m+r+u}M_{lkmg}^{hus}\overline{N}_{a}^{\ast hs}(\boldsymbol{k})b_{a}^{2g+m}\text{e}^{-b_{a}^{2}}, & & \displaystyle\end{eqnarray}$$
where we defined
$$\begin{eqnarray}M_{lkmg}^{hus}=(-1)^{m}\frac{T_{lkm}^{hu}d_{gus}^{m}\sqrt{2^{p}p!}}{(m+g)!}\sqrt{\frac{2l+1}{4\unicode[STIX]{x03C0}}\frac{(l-m)!}{(l+m)!}}.\end{eqnarray}$$
Using the form for
$\overline{{\mathcal{M}}}_{am}^{lk}$
in (5.17), the collision operator in (4.28) can therefore be expressed in terms of Hermite–Laguerre moments
$N^{pj}$
of the distribution function. We note that the moments
$\overline{{\mathcal{M}}}_{am}^{lk}$
in (5.17) reduce to the ones in Jorge et al. (Reference Jorge, Ricci and Loureiro2017) in the lowest-order drift-kinetic limit
$k_{\bot }\unicode[STIX]{x1D70C}_{tha}=0$
.
We now take Hermite–Laguerre moments of the collision operator
$\langle C(\overline{F_{a}},\overline{F_{b}})\rangle$
, i.e. we evaluate
$$\begin{eqnarray}C_{ab}^{pj}(\overline{\boldsymbol{R}})=\int \langle C(\overline{F_{a}},\overline{F_{b}})\rangle _{\overline{\boldsymbol{R}}}\frac{H_{p}(\overline{s}_{\Vert a})L_{j}(\overline{s}_{\bot a^{2}})}{\sqrt{2^{p}p!}}\frac{B_{\Vert }^{\ast }}{m_{a}}\,\text{d}\overline{v}_{\Vert }\,\text{d}\overline{\unicode[STIX]{x1D707}}\,\text{d}\overline{\unicode[STIX]{x1D703}}.\end{eqnarray}$$
Writing the gyroaveraged collision operator
$\langle C(\overline{F_{a}},\overline{F_{b}})\rangle$
in (4.28) using (4.29) and (5.17), and expanding the Bessel functions
$J_{p}(k_{\bot }\unicode[STIX]{x1D70C}_{a})$
and
$J_{p+m}(k_{\bot }^{\prime }\unicode[STIX]{x1D70C}_{a})$
using (5.10), the following form for the
$\langle \overline{c}_{ab}^{lkmnqr}\rangle _{\overline{\boldsymbol{R}}}$
term appearing in
$\langle C(\overline{F_{a}},\overline{F_{b}})\rangle _{\overline{\boldsymbol{R}}}$
is obtained
$$\begin{eqnarray}\displaystyle \langle \overline{c}_{ab}^{lkmnqr}\rangle _{\overline{\boldsymbol{R}}} & = & \displaystyle \mathop{\sum }_{u=0}^{\min (2,l,n)}\mathop{\sum }_{i=0}^{\min (l,n)-u}\mathop{\sum }_{d=-l-n+2i+2u}^{l+n-2i-2u}\mathop{\sum }_{z=0}^{\infty }\mathop{\sum }_{p,p^{\prime }=0}^{\infty }\int D_{abuidzpp^{\prime }}^{lkmnqr}(\boldsymbol{k},\boldsymbol{k}^{\prime })\nonumber\\ \displaystyle & & \displaystyle \times \,P_{l+n-2i-2u}^{d}\left(\frac{\overline{v}_{\Vert }}{\overline{v}}\right)\overline{s}_{\bot a}^{d+2z}L_{p}^{z}(\overline{s}_{\bot a}^{2})L_{p^{\prime }}^{z+d}(\overline{s}_{\bot a}^{2})\unicode[STIX]{x1D708}_{\ast abu}^{lm,nr}(\overline{v}^{2})\text{e}^{\text{i}(\boldsymbol{k}+\boldsymbol{k}^{\prime })\boldsymbol{R}}\,\text{d}\boldsymbol{k}\,\text{d}\boldsymbol{k}^{\prime }.\end{eqnarray}$$
In (5.20), we defined the
$D_{abuidzpp^{\prime }}^{lkmnqr}$
term
$$\begin{eqnarray}D_{abuidzpp^{\prime }}^{lkmnqr}(\boldsymbol{k},\boldsymbol{k}^{\prime })=\mathop{\sum }_{s=-l}^{l}\mathop{\sum }_{t=-n}^{n}E_{i+ud}^{lsnt}B_{a}^{zdp^{\prime }}\frac{d_{i}^{l-u,n-u}a_{i+u}^{ln}\text{e}^{-b_{a}^{2}-b_{a}^{\prime 2}}}{(p+z)!(z+d+p^{\prime })!}{\mathcal{N}}_{abuidz}^{lkmnqr}(\boldsymbol{k},\boldsymbol{k}^{\prime }),\end{eqnarray}$$
with
$B_{a}^{pvz^{\prime }\!}=b_{a}^{p+2z}b_{a}^{p+v+2z^{\prime }\!}$
and
$b_{a}^{\prime }=k_{\bot }^{\prime }\unicode[STIX]{x1D70C}_{tha}/2$
, while the convolution operator
${\mathcal{N}}_{abuidz}^{lkmnqr}(\boldsymbol{k},\boldsymbol{k}^{\prime })$
is given by
$$\begin{eqnarray}\displaystyle {\mathcal{N}}_{abuidz}^{lkmnqr}(\boldsymbol{k},\boldsymbol{k}^{\prime })=(-1)^{z+d}\text{e}^{\text{i}(z+d)\unicode[STIX]{x1D6FC}}b_{i+u}^{l+nd}\overline{{\mathcal{M}}}_{as}^{lk}(\boldsymbol{k})\overline{{\mathcal{M}}}_{bt}^{nq}(\boldsymbol{k}^{\prime }), & & \displaystyle\end{eqnarray}$$
with
$\overline{{\mathcal{M}}}_{as}^{lk}$
the moments of the distribution function defined in (5.17).
Finally, the result in (5.20) is used in (5.19) in order to find the Hermite–Laguerre moments
$C_{ab}^{pj}$
of the full Coulomb collision operator expressed in (4.28). This yields
$$\begin{eqnarray}C_{ab}^{pj}=\mathop{\sum }_{l,k,n,q=0}^{\infty }\mathop{\sum }_{m=0}^{k}\mathop{\sum }_{r=0}^{q}\frac{L_{km}^{l}L_{qr}^{n}}{\sqrt{2^{p}p!}}C_{ab,lkm}^{pj,nqr},\end{eqnarray}$$
with
$$\begin{eqnarray}\displaystyle C_{ab,lkm}^{pj,nqr}(\boldsymbol{k},\boldsymbol{k}^{\prime })=\mathop{\sum }_{u=0}^{\min (2,l,n)}\mathop{\sum }_{i=0}^{\min (l,n)-u}\mathop{\sum }_{d=-l-n+2i+2u}^{l+n-2i-2u}\mathop{\sum }_{z,p,p^{\prime }=0}^{\infty }D_{abuidzpp^{\prime }}^{lkmnqr}(\boldsymbol{k},\boldsymbol{k}^{\prime })I, & & \displaystyle\end{eqnarray}$$
and
$$\begin{eqnarray}I=\int f_{aM}P_{l+n-2i-2u}^{d}(\overline{v}_{\Vert }/\overline{v})\unicode[STIX]{x1D708}_{\ast abu}^{lm,nr}(\overline{v}^{2})H_{p}(\overline{s}_{\Vert a})L_{j}(\overline{s}_{\bot a}^{2})\overline{s}_{\bot a}^{d+2z}L_{p}^{z}(\overline{s}_{\bot a}^{2})L_{p^{\prime }}^{z+d}(\overline{s}_{\bot a}^{2})\frac{B_{\Vert }^{\ast }}{m_{a}}\,\text{d}\overline{v}_{\Vert }\,\text{d}\overline{\unicode[STIX]{x1D707}}.\end{eqnarray}$$
The integral factor
$I$
can be performed analytically by first rewriting the product of two Laguerre polynomials as a single one using
$$\begin{eqnarray}L_{r}^{m}(x)L_{j}(x)=\mathop{\sum }_{s=0}^{r+j}\overline{d}_{rjs}^{m}L_{s}(x),\end{eqnarray}$$
with
$$\begin{eqnarray}\overline{d}_{rjs}^{m}=\mathop{\sum }_{r_{1}=0}^{r}\mathop{\sum }_{j_{1}=0}^{j}\mathop{\sum }_{s_{1}=0}^{s}L_{rr_{1}}^{-1/2}L_{jj_{1}}^{m-1/2}L_{ss_{1}}^{-1/2}(r_{1}+j_{1}+s_{1})!,\end{eqnarray}$$
expressing the resulting Hermite–Laguerre basis in terms of Legendre–associated Laguerre using (5.14), and writing the phase-space volume
$(B_{\Vert }^{\ast }/m)\,\text{d}\overline{v}_{\Vert }\,\text{d}\overline{\unicode[STIX]{x1D707}}$
as
$\overline{v}^{2}\,\text{d}\overline{v}\,\text{d}\overline{\unicode[STIX]{x1D709}}$
with
$\overline{\unicode[STIX]{x1D709}}=\overline{v}_{\Vert }/\overline{v}$
. This yields
$$\begin{eqnarray}I=\mathop{\sum }_{h=0}^{p+z+j}\mathop{\sum }_{g=0}^{g+p^{\prime }}\mathop{\sum }_{s=0}^{p+2g}\mathop{\sum }_{t=0}^{g+\lfloor p/2\rfloor }d_{pjh}^{std}\overline{d}_{p^{\prime }hg}^{z+d}(T^{-1})_{pg}^{std}C_{\ast abu}^{st,lm,nr}\frac{(s+d)!}{(s-d)!}\frac{\unicode[STIX]{x1D6FF}_{l+n-2i-2u,s}}{4\unicode[STIX]{x03C0}(s+1/2)}.\end{eqnarray}$$
Ji & Held (Reference Ji and Held2009) present an analytical closed expression ready to be numerically implemented of the factor
$C_{\ast abu}^{st,lm,nr}=\int f_{Ma}\unicode[STIX]{x1D708}_{\ast abu}^{lm,nr}(v^{2})L_{t}^{s+1/2}(v^{2})v^{s}\,\text{d}\boldsymbol{v}$
. We note that the long-wavelength limit can be found by setting
$d=z=0$
in the collision operator (5.24). This yields the Hermite–Laguerre moments of the collision operator moments found in Jorge et al. (Reference Jorge, Ricci and Loureiro2017).
6 Small mass-ratio approximation
In this section, we simplify the electron–ion and the ion–electron collision operator in (2.2) by taking advantage of the small electron-to-ion mass ratio
$m_{e}/m_{i}$
, and derive their expressions in the gyrokinetic regime. We first consider the electron–ion collision operator.
In
$(\boldsymbol{x},\boldsymbol{v})$
phase-space coordinates, the electron–ion Coulomb collision operator can be greatly simplified by taking advantage of the fact that the ion thermal speed, is small in comparison to the electron thermal speed, for
$T_{e}\sim T_{i}$
. To first order in
$m_{e}/m_{i}$
, the electron–ion Coulomb collision operator, also called the Lorentz pitch-angle scattering operator, can be written as (Helander & Sigmar Reference Helander and Sigmar2005)
$$\begin{eqnarray}C_{ei}=\frac{n_{i}L_{ei}}{v_{the}^{3}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{c}_{e}}\boldsymbol{\cdot }\left[\frac{1}{c_{e}}\frac{\unicode[STIX]{x2202}f_{e}}{\unicode[STIX]{x2202}\boldsymbol{c}_{e}}-\frac{\boldsymbol{c}_{e}}{c_{e}^{3}}\left(\boldsymbol{c}_{e}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}f_{e}}{\unicode[STIX]{x2202}\boldsymbol{c}_{e}}\right)\right],\end{eqnarray}$$
with
$\boldsymbol{c}_{e}=\boldsymbol{v}/v_{the}$
. We expand
$f_{e}$
according to (3.12). We note that the expansion in (3.12) is an eigenbasis of the pitch-angle scattering operator
$C_{ei}$
with eigenvalue
$l(l+1)$
(Ji & Held Reference Ji and Held2008). Therefore, we write
$$\begin{eqnarray}C_{ei}=-f_{eM}\mathop{\sum }_{l,k}\frac{n_{i}L_{ei}}{v_{the}^{3}c_{e}^{3}}\frac{l(l+1)}{\sqrt{\unicode[STIX]{x1D70E}_{k}^{l}}}L_{k}^{l+1/2}(c_{e}^{2})\unicode[STIX]{x1D654}^{l}(\boldsymbol{c}_{e})\boldsymbol{\cdot }\boldsymbol{M}_{e}^{lk}(\boldsymbol{x}).\end{eqnarray}$$
We now Fourier transform the moments
$\boldsymbol{M}_{e}^{lk}$
in (6.2) as
$\boldsymbol{M}_{e}^{lk}(\overline{\boldsymbol{R}})=\int \boldsymbol{M}_{e}^{lk}(\boldsymbol{k})\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\overline{\boldsymbol{R}}}\,\text{d}\boldsymbol{k}$
and write the gyroaveraged collision operator
$C_{ei}$
as
$$\begin{eqnarray}\langle C_{ei}\rangle _{\overline{\boldsymbol{R}}}=-\int \text{d}\boldsymbol{k}\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\overline{\boldsymbol{R}}}f_{eM}\mathop{\sum }_{l,k}\frac{n_{i}L_{ei}}{v_{the}^{3}c_{e}^{3}}\frac{l(l+1)}{\sqrt{\unicode[STIX]{x1D70E}_{k}^{l}}}L_{k}^{l+1/2}(c_{e}^{2})\langle \unicode[STIX]{x1D654}^{l}(\boldsymbol{c}_{e})\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D746}_{e}}\rangle _{\overline{\boldsymbol{R}}}\boldsymbol{\cdot }\boldsymbol{M}_{e}^{lk}(\boldsymbol{k}).\end{eqnarray}$$
Using the Jacobi–Anger expansion of (4.12), (5.10), the inverse basis transformation (5.14) and the identities
$J_{-p}(x)=(-1)^{p}J_{p}(x)$
and
$$\begin{eqnarray}L_{r}^{m}(x)=\mathop{\sum }_{j=0}^{r}\binom{m+r-j-1}{r-j}L_{j}(x),\end{eqnarray}$$
we obtain
$$\begin{eqnarray}\displaystyle \langle \unicode[STIX]{x1D654}^{l}(\boldsymbol{v})\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D746}_{e}}\rangle _{\boldsymbol{R}} & = & \displaystyle \mathop{\sum }_{m=-l}^{l}\mathop{\sum }_{r=0}^{\infty }\mathop{\sum }_{j=0}^{r}\mathop{\sum }_{s=0}^{2j}\mathop{\sum }_{t=0}^{j}\sqrt{\frac{\unicode[STIX]{x03C0}^{1/2}l!}{2^{l}(l-1/2)!}\frac{(l-m)!}{(l+m)!}}\frac{i^{m}\unicode[STIX]{x1D65A}^{lm}}{(m+r)!}\frac{(m+r-j-1)!}{(r-j)!(m-1)!}\nonumber\\ \displaystyle & & \displaystyle \times \,v_{the}^{l}(T^{-1})_{0j}^{stm}b_{e}^{2r+m}e^{-b_{e}^{2}}c_{e}^{l+s}P_{l}^{m}(\cos \unicode[STIX]{x1D711})P_{s}^{m}(\cos \unicode[STIX]{x1D711})L_{t}^{s+1/2}(c_{e}^{2}),\end{eqnarray}$$
with
$b_{e}=k_{\bot }\unicode[STIX]{x1D70C}_{the}/2$
. The collision operator in (6.3) represents the gyrokinetic electron–ion collision operator.
Equation (6.5) provides an expression of the pitch-angle scattering operator
$\langle C_{ei}\rangle$
in (6.3) suitable for projection onto a Hermite–Laguerre basis, i.e.
$$\begin{eqnarray}\displaystyle C_{ei}^{pj}=\int \langle C_{ei}\rangle \frac{\displaystyle H_{p}\left(\frac{\overline{v}_{\Vert }}{v_{tha}}\right)L_{j}\left(\frac{\overline{\unicode[STIX]{x1D707}}B}{T_{a}}\right)}{\sqrt{2^{p}p!}}\frac{B_{\Vert }^{\ast }}{m_{a}}\,\text{d}\overline{v}_{\Vert }\,\text{d}\overline{\unicode[STIX]{x1D707}}\,\text{d}\overline{\unicode[STIX]{x1D703}}=2\unicode[STIX]{x03C0}\mathop{\sum }_{l=0}^{p+2j}\mathop{\sum }_{k=0}^{j+\lfloor p/2\rfloor }\frac{(T^{-1})_{pj}^{lk0}v_{the}^{3}}{\sqrt{2^{p}p!}}I_{ei}^{lk}, & & \displaystyle\end{eqnarray}$$
where we define
$$\begin{eqnarray}I_{ei}^{lk}=\int \langle C_{ei}\rangle c_{e}^{l}P_{l}(\cos \unicode[STIX]{x1D711})L_{k}^{l+1/2}(c_{e}^{2})c_{e}^{2}\,\text{d}c_{e}\,\text{d}\cos \unicode[STIX]{x1D711}.\end{eqnarray}$$
An analytical form for the integral factor
$I_{ei}^{lk}$
can be derived using the expression for
$\langle C_{ei}\rangle$
, equation (6.3) and (6.5), yielding
$$\begin{eqnarray}\displaystyle I_{ei}^{lk}(\boldsymbol{k}) & = & \displaystyle -\mathop{\sum }_{u,v}\frac{n_{e}n_{i}L_{ei}}{v_{the}^{6-u}}\frac{u(u+1)}{\unicode[STIX]{x03C0}\sqrt{\unicode[STIX]{x1D70E}_{v}^{u}}}\mathop{\sum }_{m=-u}^{u}\boldsymbol{M}_{e}^{lk}(\boldsymbol{k})\boldsymbol{\cdot }\boldsymbol{e}^{um}\mathop{\sum }_{r=0}^{\infty }\mathop{\sum }_{i=0}^{r}\mathop{\sum }_{s=0}^{2i}\mathop{\sum }_{t=0}^{i}(T^{-1})_{0i}^{stm}e^{-b_{e}^{2}}\nonumber\\ \displaystyle & & \displaystyle \times \,\sqrt{\frac{u!}{2^{u}(u-1/2)!}\frac{(u-m)!}{(u+m)!}}\frac{i^{m}b_{e}^{2r+m}}{(m+r)!}\frac{(m+r-i-1)!}{(r-i)!(m-1)!}I_{Lkt}^{lsuv}I_{Pm}^{lus},\end{eqnarray}$$
with
$I_{Lkt}^{lsuv}$
and
$I_{Pm}^{lus}$
defined by
$$\begin{eqnarray}I_{Lkt}^{lsuv}=\int L_{k}^{l+1/2}(x)L_{t}^{s+1/2}(x)\text{e}^{-x}x^{(l+u+v)/2-1}\,\text{d}x,\end{eqnarray}$$
and
$$\begin{eqnarray}I_{Pm}^{lus}=\int _{-1}^{1}P_{l}(x)P_{u}^{m}(x)P_{s}^{m}(x)\frac{\text{d}x}{2},\end{eqnarray}$$
respectively. The electron fluid moments
$\boldsymbol{M}_{e}^{lk}$
can be cast in terms of Hermite–Laguerre moments
$\overline{N}_{e}^{lk}$
using the expressions in (4.18), (5.9) and (5.17). The factor
$I_{Lkt}^{lsuv}$
can be analytically evaluated by expanding the associated Laguerre polynomials using (3.8), which leads to
$$\begin{eqnarray}I_{Lkt}^{lsuv}=\mathop{\sum }_{m_{1}=0}^{k}\mathop{\sum }_{m_{2}=0}^{t}L_{km_{1}}^{l}L_{tm_{2}}^{s}(m_{1}+m_{2}+(l+u+v)/2-1)!\end{eqnarray}$$
Similarly, the factor integral
$I_{Pm}^{lus}$
can be calculated using an extended version of Gaunt’s formula (Gaunt Reference Gaunt1929), yielding (Mavromatis & Alassar Reference Mavromatis and Alassar1999)
$$\begin{eqnarray}\displaystyle I_{Pm}^{lus}=(-1)^{m}\left(\begin{array}{@{}ccc@{}}l & u & s\\ 0 & 0 & 0\end{array}\right)\left(\begin{array}{@{}ccc@{}}l & u & s\\ 0 & m & -m\end{array}\right)\sqrt{\frac{(s+m)!(u+m)!}{(s-m)!(u-m)!}}. & & \displaystyle\end{eqnarray}$$
We note that, in (6.12), the Wigner 3-j symbol
$\left(\begin{array}{@{}ccc@{}}j_{1} & j_{2} & j_{3}\\ m_{1} & m_{2} & m_{3}\end{array}\right)$
is related to the Clebsch–Gordan coefficients
$\langle j_{1}m_{1}j_{2}m_{2}|j_{3}m_{3}\rangle$
via (Olver et al.
Reference Olver, Lozier, Boisvert and Clark2010)
$$\begin{eqnarray}\left(\begin{array}{@{}ccc@{}}j_{1} & j_{2} & j_{3}\\ m_{1} & m_{2} & m_{3}\end{array}\right)=\frac{(-1)^{j_{1}-j_{2}-m_{3}}}{\sqrt{2j_{3}+1}}\langle j_{1}m_{1}j_{2}m_{2}|j_{3}(-m_{3})\rangle ,\end{eqnarray}$$
with the Clebsch–Gordan coefficients given by
$$\begin{eqnarray}\displaystyle & & \displaystyle \langle j_{1}m_{1}j_{2}m_{2}|j_{3}m_{3}\rangle =\unicode[STIX]{x1D6FF}_{m_{3},m_{1}+m_{2}}\sqrt{\frac{(2j_{3}+1)(j_{3}+j_{1}-j_{2})!(j_{3}-j_{1}+j_{2})!(j_{1}+j_{2}-j_{3})!}{(j_{1}+j_{2}+j_{3}+1)!}}\nonumber\\ \displaystyle & & \displaystyle \qquad \times \,\sqrt{(j_{3}+m_{3})!(j_{3}-m3)!(j_{1}-m_{1})!(j_{1}+m_{1})!(j_{2}-m_{2})!(j_{3}-m_{3})!}\nonumber\\ \displaystyle & & \displaystyle \qquad \times \,\mathop{\sum }_{k}\frac{(-1)^{k}}{k!(j_{1}+j_{2}-j_{3}-k)!(j_{1}-m_{1}-k)!(j_{2}+m_{2}-k)!}\nonumber\\ \displaystyle & & \displaystyle \qquad \times \,\frac{1}{(j_{3}-j_{2}+m_{1}+k)!(j_{3}-j_{1}-m_{2}+k)!},\end{eqnarray}$$
where the summation in (6.14) is extended over all integers
$k$
that make every factorial in the sum non-negative (Bohm & Loewe Reference Bohm and Loewe1993).
We now turn to the ion–electron collision operator
$C_{ie}$
. To first order in
$m_{e}/m_{i}$
, this is given by (Helander & Sigmar Reference Helander and Sigmar2005)
$$\begin{eqnarray}C_{ie}=\unicode[STIX]{x1D708}_{ei}\frac{m_{e}}{m_{i}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{v}}\boldsymbol{\cdot }\left(\boldsymbol{v}f_{i}+\frac{T_{e}}{m_{i}}\frac{\unicode[STIX]{x2202}f_{i}}{\unicode[STIX]{x2202}\boldsymbol{v}}\right),\end{eqnarray}$$
where the electron–ion friction force is neglected for simplicity. We simplify (6.15) by using (2.27), therefore approximating the distribution function
$f_{i}$
by its gyroaveraged component
$f_{i}\simeq \langle \overline{F}_{i}\rangle _{\overline{\boldsymbol{R}}}$
, and retaining the lowest-order terms in the
$\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D6FF}}$
expansion. This allows us to convert the
$C_{ie}$
operator in (6.15) to the gyrocentre variables
$\overline{\boldsymbol{Z}}$
using the chain rule at lowest order in
$\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D6FF}}$
, i.e. to express (6.15) in
$\boldsymbol{Z}$
coordinates using the guiding-centre transformation in (2.15)–(2.17) and approximate
$\boldsymbol{Z}\simeq \overline{\boldsymbol{Z}}$
. The velocity derivatives can be written as a function of
$\overline{\boldsymbol{Z}}$
using the chain rule, yielding
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\langle \overline{F}_{i}\rangle _{\overline{\boldsymbol{R}}}}{\unicode[STIX]{x2202}\boldsymbol{v}}=\boldsymbol{b}\frac{\unicode[STIX]{x2202}\langle \overline{F}_{i}\rangle _{\overline{\boldsymbol{R}}}}{\unicode[STIX]{x2202}\overline{v}_{\Vert }}+\boldsymbol{c}\left(\sqrt{\frac{2m_{a}\overline{\unicode[STIX]{x1D707}}}{B}}\frac{\unicode[STIX]{x2202}\langle \overline{F}_{i}\rangle _{\overline{\boldsymbol{R}}}}{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D707}}}-\frac{1}{\unicode[STIX]{x1D6FA}_{i}}\boldsymbol{a}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\overline{\boldsymbol{R}}}\langle \overline{F}_{i}\rangle _{\overline{\boldsymbol{R}}}\right),\end{eqnarray}$$
where we define
$\boldsymbol{c}=(\cos \overline{\unicode[STIX]{x1D703}}\boldsymbol{e}_{1}+\sin \overline{\unicode[STIX]{x1D703}}\boldsymbol{e}_{2})$
and
$\boldsymbol{a}=\boldsymbol{c}\times \boldsymbol{b}=(-\sin \overline{\unicode[STIX]{x1D703}}\boldsymbol{e}_{1}+\cos \overline{\unicode[STIX]{x1D703}}\boldsymbol{e}_{2})$
. The ion–electron collision operator can therefore be written as
$$\begin{eqnarray}\displaystyle C_{ie} & = & \displaystyle \unicode[STIX]{x1D708}_{ei}\frac{m_{e}}{m_{i}}\left[3\langle \overline{F}_{i}\rangle _{\overline{\boldsymbol{R}}}+v_{\Vert }\frac{\unicode[STIX]{x2202}\langle \overline{F}_{i}\rangle _{\overline{\boldsymbol{R}}}}{\unicode[STIX]{x2202}\overline{v}_{\Vert }}+2\overline{\unicode[STIX]{x1D707}}\frac{\unicode[STIX]{x2202}\langle \overline{F}_{i}\rangle _{\overline{\boldsymbol{R}}}}{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D707}}}-\sqrt{\frac{2B\overline{\unicode[STIX]{x1D707}}}{m_{i}}}\frac{\boldsymbol{a}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\overline{\boldsymbol{R}}}}{\unicode[STIX]{x1D6FA}_{i}}\langle \overline{F}_{i}\rangle _{\overline{\boldsymbol{R}}}\right.\nonumber\\ \displaystyle & & \displaystyle +\,\frac{T_{e}}{m_{i}}\left(\frac{\unicode[STIX]{x2202}^{2}\langle \overline{F}_{i}\rangle _{\overline{\boldsymbol{R}}}}{\unicode[STIX]{x2202}\overline{v}_{\Vert }^{2}}+\frac{2m_{i}\overline{\unicode[STIX]{x1D707}}}{B}\frac{\unicode[STIX]{x2202}^{2}\langle \overline{F}_{i}\rangle _{\overline{\boldsymbol{R}}}}{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D707}}^{2}}+\frac{\boldsymbol{a}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\overline{\boldsymbol{R}}}\boldsymbol{a}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\overline{\boldsymbol{R}}}}{\unicode[STIX]{x1D6FA}_{i}^{2}}\langle \overline{F}_{i}\rangle _{\overline{\boldsymbol{R}}}\right.\nonumber\\ \displaystyle & & \displaystyle -\left.\left.2\sqrt{\frac{2m_{i}\overline{\unicode[STIX]{x1D707}}}{B}}\boldsymbol{a}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\overline{\boldsymbol{R}}}\frac{\unicode[STIX]{x2202}\langle \overline{F}_{i}\rangle _{\overline{\boldsymbol{R}}}}{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D707}}}+\frac{m_{i}}{B}\frac{\unicode[STIX]{x2202}\langle \overline{F}_{i}\rangle _{\overline{\boldsymbol{R}}}}{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D707}}}\right)\right].\end{eqnarray}$$
We now Fourier transform both
$T_{e}$
and
$\langle \overline{F}_{i}\rangle _{\overline{\boldsymbol{R}}}$
and gyroaverage
$C_{ie}$
, yielding
$$\begin{eqnarray}\displaystyle \langle C_{ie}\rangle _{\overline{\boldsymbol{R}}} & = & \displaystyle \unicode[STIX]{x1D708}_{ei}\frac{m_{e}}{m_{i}}\int \text{e}^{\text{i}(\boldsymbol{k}+\boldsymbol{k}^{\prime })\boldsymbol{\cdot }\boldsymbol{R}}\left[\langle \overline{F}_{i}\rangle _{\overline{\boldsymbol{R}}}+\frac{\unicode[STIX]{x2202}\langle \overline{F}_{i}\rangle _{\overline{\boldsymbol{R}}}}{\unicode[STIX]{x2202}\overline{v}_{\Vert }}+2\overline{\unicode[STIX]{x1D707}}\frac{\unicode[STIX]{x2202}\langle \overline{F}_{i}\rangle _{\overline{\boldsymbol{R}}}}{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D707}}}\right.\nonumber\\ \displaystyle & & \displaystyle +\,J_{0}(k_{\bot }^{\prime }\overline{\unicode[STIX]{x1D70C}_{i}})\frac{T_{e}(\boldsymbol{k}^{\prime })}{m_{i}}\left(\frac{\unicode[STIX]{x2202}^{2}\langle \overline{F}_{i}\rangle _{\overline{\boldsymbol{R}}}}{\unicode[STIX]{x2202}\overline{v}_{\Vert }^{2}}+\frac{2m_{i}\overline{\unicode[STIX]{x1D707}}}{B}\frac{\unicode[STIX]{x2202}^{2}\langle \overline{F}_{i}\rangle _{\overline{\boldsymbol{R}}}}{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D707}}^{2}}+\frac{m_{i}}{B}\frac{\unicode[STIX]{x2202}\langle \overline{F}_{i}\rangle _{\overline{\boldsymbol{R}}}}{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D707}}}\right).\nonumber\\ \displaystyle & & \displaystyle +\,\frac{T_{e}(\boldsymbol{k}^{\prime })}{m_{i}}i\frac{\langle \overline{F}_{i}\rangle _{\overline{\boldsymbol{R}}}}{2\unicode[STIX]{x1D6FA}_{i}^{2}}[J_{0}(k_{\bot }^{\prime }\overline{\unicode[STIX]{x1D70C}_{i}})k_{\bot }^{2}+J_{2}(k_{\bot }^{\prime }\overline{\unicode[STIX]{x1D70C}_{i}})\boldsymbol{k}\boldsymbol{k}\boldsymbol{ : }(\boldsymbol{e}_{1}\boldsymbol{e}_{1}-\boldsymbol{e}_{2}\boldsymbol{e}_{2})]\nonumber\\ \displaystyle & & \displaystyle -\left.\frac{T_{e}(\boldsymbol{k}^{\prime })}{m_{i}}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{e}_{2}iJ_{1}(k_{\bot }^{\prime }\overline{\unicode[STIX]{x1D70C}_{i}})\frac{2m_{i}\overline{v}_{\bot }}{B\unicode[STIX]{x1D6FA}_{i}}\frac{\unicode[STIX]{x2202}\langle \overline{F}_{i}\rangle _{\overline{\boldsymbol{R}}}}{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D707}}}\right],\end{eqnarray}$$
where we have used the identities
$\langle \boldsymbol{a}\text{e}^{\text{i}\boldsymbol{k}^{\prime }\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D746}_{i}}}\rangle _{\overline{\boldsymbol{R}}}=iJ_{1}(k_{\bot }^{\prime }\overline{\unicode[STIX]{x1D70C}_{i}})\boldsymbol{e}_{2}$
and
$\langle \boldsymbol{a}\boldsymbol{a}\text{e}^{\text{i}\boldsymbol{k}^{\prime }\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D746}_{i}}}\rangle _{\overline{\boldsymbol{R}}}=(1/2)[J_{0}(k_{\bot }^{\prime }\overline{\unicode[STIX]{x1D70C}_{i}})(\boldsymbol{e}_{1}\boldsymbol{e}_{1}+\boldsymbol{e}_{2}\boldsymbol{e}_{2})+J_{2}(k_{\bot }^{\prime }\overline{\unicode[STIX]{x1D70C}_{i}})(\boldsymbol{e}_{1}\boldsymbol{e}_{1}-\boldsymbol{e}_{2}\boldsymbol{e}_{2})]$
. Finally, we take Hermite–Laguerre moments of the gyroaveraged ion–electron collision operator
$\langle C_{ie}\rangle _{\overline{\boldsymbol{R}}}$
in (6.18), using the expansion of
$\langle \overline{F}_{i}\rangle _{\overline{\boldsymbol{R}}}$
in (5.1), yielding
$$\begin{eqnarray}\displaystyle C_{ie}^{pj} & = & \displaystyle \unicode[STIX]{x1D708}_{ei}\frac{m_{e}}{m_{i}}\int \text{e}^{\text{i}(\boldsymbol{k}+\boldsymbol{k}^{\prime })\boldsymbol{\cdot }\boldsymbol{R}}\mathop{\sum }_{l,k}\left[\unicode[STIX]{x1D608}_{lk}^{pj}+\text{e}^{-b_{i}^{2}}\frac{T_{e}(k_{\bot }^{\prime })}{T_{i}}\mathop{\sum }_{r=0}^{\infty }\frac{b_{i}^{2r}}{r!}\left(\mathop{\sum }_{s=0}^{r+j}d_{rjs}^{0}B_{lkrs}^{pj}\right.\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\left.\mathop{\sum }_{v=0}^{r}\mathop{\sum }_{s=0}^{v+j+1}\frac{d_{vjs}^{1}i\unicode[STIX]{x1D70C}_{thi}^{2}b_{i}^{2}\unicode[STIX]{x1D6FF}_{lp}\unicode[STIX]{x1D6FF}_{ks}}{4(r+1)(r+2)}\boldsymbol{k}\boldsymbol{k}:(\boldsymbol{e}_{1}\boldsymbol{e}_{1}-\boldsymbol{e}_{2}\boldsymbol{e}_{2})\right)\right]\overline{N}_{i}^{lk}(\boldsymbol{k}),\end{eqnarray}$$
with
$\unicode[STIX]{x1D608}_{lk}^{pj}$
given by
$$\begin{eqnarray}\unicode[STIX]{x1D608}_{lk}^{pj}=2j\unicode[STIX]{x1D6FF}_{lp}\unicode[STIX]{x1D6FF}_{kj-1}-(p+2j)\unicode[STIX]{x1D6FF}_{lp}\unicode[STIX]{x1D6FF}_{kj}-\sqrt{p(p-1)}\unicode[STIX]{x1D6FF}_{lp-2}\unicode[STIX]{x1D6FF}_{kj},\end{eqnarray}$$
and
$B_{lkrs}^{pj}$
by
$$\begin{eqnarray}\displaystyle B_{lkrs}^{pj} & = & \displaystyle \sqrt{p(p-1)}\unicode[STIX]{x1D6FF}_{lp-2}\unicode[STIX]{x1D6FF}_{ks}+\frac{T_{i}}{m_{i}}\frac{\text{i}k_{\bot }^{2}}{2\unicode[STIX]{x1D6FA}_{i}^{2}}-\mathop{\sum }_{i=0}^{s-1}(3+2s)\unicode[STIX]{x1D6FF}_{lp}\unicode[STIX]{x1D6FF}_{ki}+\mathop{\sum }_{i=0}^{s-2}2s\unicode[STIX]{x1D6FF}_{lp}\unicode[STIX]{x1D6FF}_{ki}\nonumber\\ \displaystyle & & \displaystyle +\,\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{e}_{2}\frac{2v_{thi}}{\unicode[STIX]{x1D6FA}_{i}}\frac{b_{i}}{r+1}[(1+s)\unicode[STIX]{x1D6FF}_{lp}\unicode[STIX]{x1D6FF}_{ks}-s\unicode[STIX]{x1D6FF}_{lp}\unicode[STIX]{x1D6FF}_{ks-1}].\end{eqnarray}$$
7 Conclusion
In this work, a formulation of the nonlinear gyrokinetic Coulomb collision operator is derived, providing an extension of a previously derived nonlinear Coulomb drift-kinetic collision operator to the gyrokinetic regime. This constitutes a key element necessary to perform quantitative studies of turbulence, flows and, in general, of the plasma dynamics in the periphery of magnetized fusion devices. The gyroaveraged collision operator is cast in terms of parallel and perpendicular velocity integrals of the gyroaveraged distribution function at arbitrary
$k_{\bot }\unicode[STIX]{x1D70C}_{s}$
, yielding the formula in (4.28). In order to provide an analytical formulation of the Coulomb collision operator ready to be used in pseudospectral formulations of the gyrokinetic equation for distribution functions arbitrarily far from equilibrium and for an arbitrary collisionality regime, the Hermite–Laguerre moments of the gyroaveraged collision operator are evaluated, yielding (5.23). Furthermore, the electron-to-ion mass ratio is used to simplify the form of the electron–ion and ion–electron collision operators, yielding (6.6) and (6.19), respectively.
We conclude by noting that the present collision operator is derived by porting the Coulomb operator to the gyrocentre phase-space by using a framework valid up to second order in the expansion parameter
$\unicode[STIX]{x1D716}$
, yielding second-order accurate momentum and energy conservation laws. The use of the techniques developed here to analytically gyroaverage the Coulomb operator and obtain its projection onto an orthogonal polynomial basis should, in principle, be applicable to collision operators of the Fokker–Planck type that add the necessary correction terms in order to ensure exact conservation laws (Brizard Reference Brizard2004; Burby et al.
Reference Burby, Brizard and Qin2015; Sugama et al.
Reference Sugama, Watanabe and Nunami2015).
Acknowledgements
We thank L. M. Perrone for the helpful insight on the relation between spherical basis tensors and spherical harmonics and the anonymous referees for their careful review of our article. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 and 2019–2020 under grant agreement no. 633053, from Portuguese FCT (Fundação para a Ciência e Tecnologia) under grant PD/BD/105979/2014, carried out as part of the training in the framework of the Advanced Program in Plasma Science and Engineering (APPLAuSE,) sponsored by FCT under grant no. PD/00505/2012 at Instituto Superior Técnico, and from the Swiss National Science Foundation. The views and opinions expressed herein do not necessarily reflect those of the European Commission.
Appendix A. Basis tensors
In this appendix, we derive the form of the basis tensors
$\unicode[STIX]{x1D65A}^{lm}$
used in the definition of
$\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})$
in (4.8). We start with the
$l=1$
case, for which (4.8) yields
$$\begin{eqnarray}\unicode[STIX]{x1D654}^{1}(\boldsymbol{v})=\boldsymbol{v}=\sqrt{\frac{4\unicode[STIX]{x03C0}}{3}}v\mathop{\sum }_{m=-1}^{1}Y_{1m}(\unicode[STIX]{x1D719},\unicode[STIX]{x1D703})\boldsymbol{e}^{1m}.\end{eqnarray}$$
The spherical basis vectors
$\boldsymbol{e}^{1m}$
can then be derived from (A 1) by expressing the vector
$\boldsymbol{v}$
in spherical coordinates as
$$\begin{eqnarray}\boldsymbol{v}=v(\sin \unicode[STIX]{x1D719}\cos \unicode[STIX]{x1D703}\boldsymbol{e}_{x}+\sin \unicode[STIX]{x1D719}\sin \unicode[STIX]{x1D703}\boldsymbol{e}_{y}+\cos \unicode[STIX]{x1D719}\boldsymbol{e}_{z}),\end{eqnarray}$$
and using the identities for the spherical harmonics
$$\begin{eqnarray}Y_{1m}(\unicode[STIX]{x1D719},\unicode[STIX]{x1D703})=\left\{\begin{array}{@{}ll@{}}\displaystyle \sqrt{\frac{3}{8\unicode[STIX]{x03C0}}}\sin \unicode[STIX]{x1D719}\text{e}^{-\text{i}\unicode[STIX]{x1D703}},\quad & m=-1,\\ \displaystyle \sqrt{\frac{3}{4\unicode[STIX]{x03C0}}}\cos \unicode[STIX]{x1D719},\quad & m=0,\\ \displaystyle -\sqrt{\frac{3}{8\unicode[STIX]{x03C0}}}\sin \unicode[STIX]{x1D719}\text{e}^{\text{i}\unicode[STIX]{x1D703}},\quad & m=1,\end{array}\right.\end{eqnarray}$$
therefore obtaining
$$\begin{eqnarray}\boldsymbol{e}^{1m}=\left\{\begin{array}{@{}ll@{}}\displaystyle \frac{\boldsymbol{e}_{x}-\text{i}\boldsymbol{e}_{y}}{\sqrt{2}},\quad & m=-1,\\ \boldsymbol{e}_{z},\quad & m=0,\\ \displaystyle -\frac{\boldsymbol{e}_{x}+\text{i}\boldsymbol{e}_{y}}{\sqrt{2}},\quad & m=1.\\ \quad \end{array}\right.\end{eqnarray}$$
We now construct spherical basis tensors
$\unicode[STIX]{x1D65A}^{lm}$
from the spherical basis vectors
$\boldsymbol{e}^{1m}$
leveraging the techniques developed for the angular momentum formalism in quantum mechanics (Zettili & Zahed Reference Zettili and Zahed2009; Snider Reference Snider2017). As a first step, we note that the basis vectors
$\boldsymbol{e}^{1m}$
are eigenvectors of the angular momentum matrix
$\unicode[STIX]{x1D642}_{z}$
$$\begin{eqnarray}\unicode[STIX]{x1D642}_{z}=\text{i}\left(\begin{array}{@{}ccc@{}}0 & -1 & 0\\ 1 & 0 & 0\\ 0 & 0 & 0\end{array}\right),\end{eqnarray}$$
with eigenvalue
$m$
, that is
$$\begin{eqnarray}\unicode[STIX]{x1D642}_{z}\boldsymbol{\cdot }\boldsymbol{e}^{1m}=m\boldsymbol{e}^{1m}.\end{eqnarray}$$
As a second step, we note that the relationship between the basis vectors
$\boldsymbol{e}_{\unicode[STIX]{x1D6FC}}$
for
$\unicode[STIX]{x1D6FC}=(x,y,z)$
and the angular momentum matrices
$\unicode[STIX]{x1D642}_{\unicode[STIX]{x1D6FC}}$
is given by
$$\begin{eqnarray}\unicode[STIX]{x1D642}_{\unicode[STIX]{x1D6FC}}=-\text{i}\boldsymbol{e}_{\unicode[STIX]{x1D6FC}}\boldsymbol{\cdot }\unicode[STIX]{x1D716},\end{eqnarray}$$
with
$\unicode[STIX]{x1D716}$
the standard Levi-Civita tensor. In index notation, (A 7) can be written as
$$\begin{eqnarray}(\unicode[STIX]{x1D642}_{\unicode[STIX]{x1D6FC}})_{kl}=-\text{i}\mathop{\sum }_{j=1}^{3}(e_{\unicode[STIX]{x1D6FC}})_{j}\unicode[STIX]{x1D716}_{jkl}.\end{eqnarray}$$
The raising
$G_{+}$
and lowering
$G_{-}$
operators (corresponding to the ladder operators in quantum mechanics), defined by
$$\begin{eqnarray}G_{\pm }=G_{x}\pm \text{i}G_{y}\end{eqnarray}$$
allows us to obtain the basis vectors
$\boldsymbol{e}^{1\pm 1}$
from
$\boldsymbol{e}^{10}$
using
$$\begin{eqnarray}G_{\pm }\boldsymbol{e}^{10}=\boldsymbol{e}^{1\pm 1}.\end{eqnarray}$$
In addition, we have that
$$\begin{eqnarray}\boldsymbol{e}^{1-1}=(G_{-})^{2}\boldsymbol{e}^{11}.\end{eqnarray}$$
We can now define the spherical tensor basis
$\unicode[STIX]{x1D65A}^{lm}$
that define the irreducible tensors
$\unicode[STIX]{x1D654}^{l}$
. We start with the spherical basis tensor
$$\begin{eqnarray}\unicode[STIX]{x1D65A}^{ll}=\boldsymbol{e}^{11}\boldsymbol{e}^{11}\ldots \boldsymbol{e}^{11},\end{eqnarray}$$
formed by the product of
$l$
basis vectors
$\boldsymbol{e}^{11}$
. Similarly to
$\unicode[STIX]{x1D654}^{l}(\boldsymbol{v})$
, this tensor is of rank
$l$
, symmetric, and totally traceless, as
$\boldsymbol{e}^{11}\boldsymbol{\cdot }\boldsymbol{e}^{11}=0$
. Furthermore, we note that
$\unicode[STIX]{x1D65A}^{ll}$
is an eigenvector with eigenvalue
$l$
of the angular momentum tensor
$\unicode[STIX]{x1D642}_{z}^{l}$
, with
$G_{n}^{l}$
a tensor of rank
$2l$
defined by
$$\begin{eqnarray}\displaystyle \left[\unicode[STIX]{x1D642}_{\unicode[STIX]{x1D6FC}}^{l}\right]_{a_{1}a_{2}\ldots a_{l}b_{1}b_{2}\ldots b_{l}} & = & \displaystyle \mathop{\sum }_{j^{\prime }k^{\prime }\ldots l^{\prime }}\left\{[\unicode[STIX]{x1D642}_{\unicode[STIX]{x1D6FC}}]_{a_{1}b_{1}}\unicode[STIX]{x1D6FF}_{a_{2}b_{2}}\cdots \unicode[STIX]{x1D6FF}_{a_{l}b_{l}}+\unicode[STIX]{x1D6FF}_{a_{1}b_{1}}[\unicode[STIX]{x1D642}_{\unicode[STIX]{x1D6FC}}]_{a_{2}b_{2}}\cdots \unicode[STIX]{x1D6FF}_{a_{l}b_{l}}\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\cdots +\unicode[STIX]{x1D6FF}_{a_{1}b_{1}}\unicode[STIX]{x1D6FF}_{a_{2}b_{2}}\cdots [\unicode[STIX]{x1D642}_{\unicode[STIX]{x1D6FC}}]_{a_{l}b_{l}}\right\}.\end{eqnarray}$$
The remaining basis tensor elements
$\unicode[STIX]{x1D65A}^{lm}$
can be obtained by applying the tensorial lowering operator
$G_{-}^{l}=G_{x}^{l}-iG_{y}^{l}$
to
$\unicode[STIX]{x1D65A}^{ll}$
, namely
$$\begin{eqnarray}\unicode[STIX]{x1D65A}^{lm}=\sqrt{\frac{(l+m)!}{(2l)!(l-m)!}}(G_{-}^{l}\cdot ^{l})^{l-m}\unicode[STIX]{x1D65A}^{ll},\end{eqnarray}$$
with
$m=-l,-l+1,\ldots ,-1,0,1,\ldots ,l$
and
$(G_{-}^{l}\cdot ^{l})^{l-m}~\unicode[STIX]{x1D65A}^{ll}$
a tensor of order
$l$
built by the application of the
$G_{-}^{l}\cdot ^{l}$
operator to
$\unicode[STIX]{x1D65A}^{ll}$
$l-m$
times. The normalization factor in (A 14) is obtained by requiring that the contravariant
$\unicode[STIX]{x1D65A}^{lm}$
and the covariant
$\boldsymbol{e}_{m}^{l}$
basis tensors form an orthonormal basis, i.e.
$$\begin{eqnarray}\unicode[STIX]{x1D65A}^{lm}\boldsymbol{\cdot }\boldsymbol{e}_{m^{\prime }}^{l}=\unicode[STIX]{x1D6FF}_{m,m^{\prime }}.\end{eqnarray}$$
In order to find a covariant basis
$e_{m}^{l}$
, we start with the case
$l=1$
and note that the set of vectors
$\boldsymbol{e}_{m}^{1}=(\boldsymbol{e}^{1m})^{\ast }=(-1)^{m}\boldsymbol{e}^{1-m}$
, with
$(\boldsymbol{e}_{m}^{1})^{\ast }$
the complex conjugate of
$\boldsymbol{e}_{m}^{1}$
satisfies (A 15). We therefore define
$\boldsymbol{e}_{m}^{l}=(\unicode[STIX]{x1D65A}^{lm})^{\ast }$
, and use (A 15) to normalize
$\unicode[STIX]{x1D65A}^{lm}$
. For computational purposes, we note that the tensor
$\unicode[STIX]{x1D65A}^{lm}$
can also be written as a function of the basis vectors
$\boldsymbol{e}^{1m}$
as (Snider Reference Snider2017)
$$\begin{eqnarray}\unicode[STIX]{x1D65A}^{lm}=N_{lm}\mathop{\sum }_{n=0}^{\lfloor (l+m)/2\rfloor }a_{n}^{lm}{\{(\boldsymbol{e}^{11})^{m+n}(\boldsymbol{e}^{1-1})^{n}(\boldsymbol{e}^{10})^{l-m-2n}\}}_{TS},\end{eqnarray}$$
where
$N_{lm}=\sqrt{(l+m)!(l-m)!2^{l-m}/(2l)!}$
and
$a_{n}^{lm}=l!/[2^{n}n!(m+n)!(l-m-2n)!]$
.
Appendix B. Basis transformation
In this section, we derive a closed-form expression for the
$T_{lkm}^{pj}$
and
$(T^{-1})_{pj}^{lkm}$
coefficients defined in (5.13) and (5.14). By multiplying (5.13) by a Hermite and a Laguerre polynomial and by an exponential of the form
$\text{e}^{-\overline{v}^{2}}$
, and integrating over the whole
$\overline{v}_{\Vert }$
and
$\overline{\unicode[STIX]{x1D707}}$
space, we obtain the following integral expression for
$T_{lkm}^{pj}$
$$\begin{eqnarray}T_{lkm}^{pj}=\frac{v_{tha}^{m-l}}{2^{p}p!\sqrt{\unicode[STIX]{x03C0}}}\int \frac{\overline{v}^{l}}{\overline{v}_{\bot }^{m}}P_{l}^{m}\left(\frac{\overline{v}_{\Vert }}{\overline{v}}\right)L_{k}^{l+1/2}\left(\frac{\overline{v}^{2}}{v_{tha}^{2}}\right)H_{p}\left(\frac{\overline{v}_{\Vert a}}{v_{tha}}\right)L_{j}\left(\frac{\overline{v}_{\bot }^{2}}{v_{tha}^{2}}\right)\text{e}^{-v^{2}/v_{tha}^{2}}\frac{\text{d}\boldsymbol{v}}{2\unicode[STIX]{x03C0}}.\end{eqnarray}$$
We first write the integrand in (B 1) in terms of
$\overline{\unicode[STIX]{x1D709}}=\overline{v}_{\Vert }/\overline{v}$
and
$\overline{v}$
coordinates using the basis transformation in (5.14), yielding
$$\begin{eqnarray}\displaystyle T_{lkm}^{pj} & = & \displaystyle \mathop{\sum }_{l^{\prime }=0}^{p+2j}\mathop{\sum }_{k^{\prime }=0}^{j+\lfloor p/2\rfloor }\frac{(l+1/2)k!}{(l+k+1/2)!}T_{l^{\prime }k^{\prime }}^{pj}\nonumber\\ \displaystyle & & \displaystyle \times \int _{-1}^{1}\frac{P_{l}^{m}(\overline{\unicode[STIX]{x1D709}})P_{l^{\prime }}(\overline{\unicode[STIX]{x1D709}})}{(1-\overline{\unicode[STIX]{x1D709}})^{2}}\,\text{d}\overline{\unicode[STIX]{x1D709}}\int _{0}^{\infty }x_{a}^{(l+l^{\prime }-m+1)/2}L_{k}^{l+1/2}(x_{a})L_{k^{\prime }}^{l^{\prime }+1/2}(x_{a})\,\text{d}x_{a},\end{eqnarray}$$
where we used the fact that
$(T^{-1})_{lk}^{pj}=T_{lk}^{pj}\sqrt{\unicode[STIX]{x03C0}}2^{p}p!k!(l+1/2)/(k+l+1/2)!$
(Jorge et al.
Reference Jorge, Ricci and Loureiro2017). The
$\unicode[STIX]{x1D709}$
integral in (B 2) is performed by expanding
$P_{l}$
as a finite sum of the form
$$\begin{eqnarray}P_{l}(x)=\mathop{\sum }_{s=0}^{l}c_{s}^{l}x^{s},\end{eqnarray}$$
with the coefficients
$c_{s}^{l}=2^{l}[(l+s-1)/2]!/[s!(l-s)!((s-l-1)/2)!]$
, and using the relation between associated Legendre functions
$P_{l}^{m}(x)$
and Legendre polynomials
$P_{l}(x)$
$$\begin{eqnarray}P_{l}^{m}(x)=(-1)^{m}(1-x^{2})^{m/2}\frac{\text{d}^{m}P_{l}(x)}{\text{d}x^{m}}.\end{eqnarray}$$
The
$x$
integral in (B 2) is performed by using the expansion of the associated Laguerre polynomials in (3.8). The
$T_{lkm}^{pj}$
coefficient can then be written as
$$\begin{eqnarray}\displaystyle T_{lkm}^{pj} & = & \displaystyle \mathop{\sum }_{l^{\prime }=0}^{p+2j}\mathop{\sum }_{k^{\prime }=0}^{j+\lfloor p/2\rfloor }T_{l^{\prime }k^{\prime }}^{pj}\frac{(l^{\prime }+1/2)k^{\prime }!}{(l^{\prime }+k^{\prime }+1/2)!}\mathop{\sum }_{m_{1}=0}^{k}\mathop{\sum }_{m_{2}=0}^{k^{\prime }}\mathop{\sum }_{s_{1}=0}^{l}\mathop{\sum }_{s_{2}=0}^{l^{\prime }}L_{km_{1}}^{l}L_{k^{\prime }m_{2}}^{l^{\prime }}\nonumber\\ \displaystyle & & \displaystyle \times \,\frac{c_{s_{1}}^{l}c_{s_{2}}^{l^{\prime }}}{2}\frac{s_{1}!}{(s_{1}-m)!}\frac{[1+(-1)^{s_{1}+s_{2}-m}]}{s_{1}+s_{2}+1-m}\left(m_{1}+m_{2}+\frac{l+l^{\prime }-m+1}{2}\right)!\end{eqnarray}$$
The inverse transformation coefficients
$(T^{-1})_{pj}^{lkm}$
defined by (5.14) can be found similarly, yielding
$$\begin{eqnarray}(T^{-1})_{pj}^{lkm}=\frac{2^{p}p!\sqrt{\unicode[STIX]{x03C0}}k!(l+1/2)(l-m)!}{(k+l+1/2)!(l+m)!}T_{lkm}^{pj}.\end{eqnarray}$$
