Nonlinear Gyrokinetic Coulomb Collision Operator

TL;DR

This paper derives a nonlinear gyrokinetic Coulomb collision operator suitable for plasma edge regions in magnetic confinement fusion, enabling accurate modeling of far-from-equilibrium collisions at small spatial scales.

Contribution

It introduces a novel gyrokinetic Coulomb collision operator expressed in gyrocenter coordinates, with a numerical implementation framework and expansions in velocity space.

Findings

Operator accurately describes far-from-equilibrium collisions

Provides a gyroaverage formula in gyrocenter phase space

Includes Hermite-Laguerre velocity space expansion

Abstract

A gyrokinetic Coulomb collision operator is derived, which is particularly useful to describe the plasma dynamics at the periphery region of magnetic confinement fusion devices. The derived operator is able to describe collisions occurring in distribution functions arbitrarily far from equilibrium with variations on spatial scales at and below the particle Larmor radius. A multipole expansion of the Rosenbluth potentials is used in order to derive the dependence of the full Coulomb collision operator on the particle gyroangle. The full Coulomb collision operator is then expressed in gyrocentre phase-space coordinates, and a closed formula for its gyroaverage in terms of the moments of the gyrocenter distribution function in a form ready to be numerically implemented is provided. Furthermore, the collision operator is projected onto a Hermite-Laguerre velocity space polynomial basis and expansions in the small electron-to-ion mass ratio are provided.

Equations

$$\frac{\partial f_a}{\partial t}+\dot{\mathbf{x}}\cdot \frac{\partial f_a}{\partial \mathbf{x}}+\dot{\mathbf{v}}\cdot \frac{\partial f_a}{\partial \mathbf{v}}=\sum _{b}C(f_a,f_b),$$

$$\begin{array}{rl}C(f_a,f_b)& =L_{ab}\sum _{j,k}\frac{\partial }{\partial v_k}\left[\frac{\partial }{\partial v_j}\left(\frac{{\partial }^2{G}_b}{\partial v_k\partial v_j}{f}_a\right)-2\left(1+\frac{m_a}{m_b}\right)\frac{\partial H_b}{\partial v_k}{f}_a\right],\end{array}$$

$$\begin{array}{rl}{H}_b& =\int \frac{f_b({\mathbf{v}}^{\prime })}{|\mathbf{v}-{\mathbf{v}}^{\prime }|}d{\mathbf{v}}^{\prime },\end{array}$$

$$\begin{array}{r}{G}_b=\int f_b({\mathbf{v}}^{\prime })|\mathbf{v}-{\mathbf{v}}^{\prime }|d{\mathbf{v}}^{\prime },\end{array}$$

$$ϵ\sim \frac{|{\mathbf{v}}_{\mathbf{E}}|}{c_s}\sim \frac{k_{\parallel }}{k_{\perp }}\ll 1,$$

$$\frac{e{\varphi }_0}{T_e}\sim 1,$$

$$\frac{{\varphi }_1}{{\varphi }_0}\sim {ϵ}_{\delta }\ll 1.$$

$$\mathbf{E}\sim {\nabla }_{\perp }{\varphi }_0\sim {\nabla }_{\perp }{\varphi }_1.$$

$${\rho }_s\left|\frac{{\nabla }_{\perp }{\varphi }_1}{{\varphi }_1}\right|\sim 1,$$

$${\rho }_s\left|\frac{{\nabla }_{\perp }{\varphi }_0}{{\varphi }_0}\right|\sim {ϵ}_{\delta }.$$

$$k_{\perp }{\rho }_s\frac{e{\varphi }_0}{T_e}\sim \frac{e{\varphi }_1}{T_e}\sim ϵ.$$

$$\frac{{\nu }_i}{{\Omega }_i}\sim {ϵ}_{\nu }\sim {ϵ}^2,$$

$$\overline{\mathbf{Z}}=T\mathbf{Z}=\mathbf{Z}+{ϵ}_{\delta }{\mathbf{Z}}_1+...,$$

$$F_a(\mathbf{Z})=f_a(\mathbf{x}(\mathbf{Z}),\mathbf{v}(\mathbf{Z})).$$

$$\mathbf{v}=v_{\parallel }\mathbf{b}+{\mathbf{v}}_{E0}+v_{\perp }{}^{\prime }(\cos \theta {\mathbf{e}}_1+\sin \theta {\mathbf{e}}_2)=v\left[\cos \phi \mathbf{b}+\sin \phi (\cos \theta {\mathbf{e}}_1+\sin \theta {\mathbf{e}}_2)\right],$$

$$\mu =\frac{m_a{v}_{\perp }^{{}^{\prime }2}}{2B},$$

$$\mathbf{x}=\mathbf{R}+{\rho }_a,$$

$${\rho }_a={\rho }_a(\mathbf{R},\mu ,\theta )=\sqrt{\frac{2m_a\mu }{q_a^{2}B}}(-\sin \theta {\mathbf{e}}_1+\cos \theta {\mathbf{e}}_2)$$

$$\overline{F_a}(\overline{\mathbf{Z}})=F_a(\mathbf{Z})=T{\overline{F}}_a(\mathbf{Z}),$$

$$\frac{\partial \overline{F_a}}{\partial t}+\dot{\overline{\mathbf{Z}}}\cdot \frac{\partial \overline{F_a}}{\partial \overline{\mathbf{Z}}}=\sum _{b}C(\overline{F_a},\overline{F_b}).$$

$${⟨\chi ⟩}_{\overline{\mathbf{R}}}=\frac{1}{2\pi }{\int }_0^{2\pi }\chi (\overline{\mathbf{Z}})d\overline{\theta },$$

$$\frac{\partial }{\partial t}{⟨\overline{F_a}⟩}_{\overline{\mathbf{R}}}+{⟨\dot{\overline{\mathbf{Z}}}\cdot \frac{\partial \overline{F_a}}{\partial \overline{\mathbf{Z}}}⟩}_{\overline{\mathbf{R}}}=\sum _b{⟨C(\overline{F_a},\overline{F_b})⟩}_{\overline{\mathbf{R}}}.$$

$$\frac{\partial }{\partial t}{⟨\overline{F_a}⟩}_{\overline{\mathbf{R}}}+\dot{\overline{\mathbf{R}}}\cdot \frac{\partial {⟨\overline{F_a}⟩}_{\overline{\mathbf{R}}}}{\partial \overline{\mathbf{R}}}+\dot{\overline{v_{\parallel }}}\frac{\partial {⟨\overline{F_a}⟩}_{\overline{\mathbf{R}}}}{\partial \overline{v_{\parallel }}}=\sum _b{⟨C(\overline{F_a},\overline{F_b})⟩}_{\overline{\mathbf{R}}}.$$

$$\frac{\partial \widetilde{\overline{F_a}}}{\partial t}+\dot{\overline{\mathbf{R}}}\cdot \frac{\partial \widetilde{\overline{F_a}}}{\partial \overline{\mathbf{R}}}+\dot{\overline{v_{\parallel }}}\frac{\partial \widetilde{\overline{F_a}}}{\partial \overline{v_{\parallel }}}+\dot{\overline{\theta }}\frac{\partial \widetilde{\overline{F_a}}}{\partial \overline{\theta }}=\sum _{b}C(\overline{F_a},\overline{F_b})-{⟨C(\overline{F_a},\overline{F_b})⟩}_{\overline{\mathbf{R}}}.$$

$$\widetilde{\overline{F_a}}\simeq \frac{1}{{\Omega }_a}\sum _b{\int }_0^{\overline{\theta }}\left[C\left({⟨\overline{F_a}⟩}_{\overline{\mathbf{R}}},{⟨\overline{F_b}⟩}_{\overline{\mathbf{R}}}\right)-{⟨C\left({⟨\overline{F_a}⟩}_{\overline{\mathbf{R}}},{⟨\overline{F_b}⟩}_{\overline{\mathbf{R}}}\right)⟩}_{\overline{\mathbf{R}}}\right]d{\overline{\theta }}^{\prime }.$$

$$\frac{\widetilde{\overline{F_e}}}{{⟨\overline{F_e}⟩}_{\overline{\mathbf{R}}}}\sim {\left(\frac{T_i}{T_e}\right)}^{3/2}\sqrt{\frac{m_e}{m_i}}{ϵ}_{\nu }\sim {\left(\frac{T_i}{T_e}\right)}^{3/2}\sqrt{\frac{m_e}{m_i}}{ϵ}^2,$$

$$\frac{\widetilde{{\overline{F}}_i}}{{⟨{\overline{F}}_i⟩}_{\overline{\mathbf{R}}}}\sim \frac{{\nu }_i}{{\Omega }_i}\sim {ϵ}_{\nu }\sim {ϵ}^2.$$

$$C({\overline{F}}_e,\overline{F_b})\sim {\nu }_e{\overline{F}}_e\sim \sqrt{\frac{m_i}{m_e}}{ϵ}_{\nu }{\Omega }_i{\overline{F}}_e+O\left({ϵ}_{\nu }{ϵ}_{\delta }{\Omega }_i{\overline{F}}_e\right)\sim \sqrt{\frac{m_i}{m_e}}{ϵ}^2{\Omega }_i{\overline{F}}_e+O\left({ϵ}^3{\Omega }_i{\overline{F}}_e\right).$$

$$C({\overline{F}}_i,\overline{F_b})\sim {\nu }_i{\overline{F}}_i\sim {ϵ}_{\nu }{\Omega }_i{\overline{F}}_i+O\left({ϵ}_{\nu }{ϵ}_{\delta }{\Omega }_i{\overline{F}}_i\right)\sim {ϵ}^2{\Omega }_i{\overline{F}}_i+O\left({ϵ}^3{\Omega }_i{\overline{F}}_i\right).$$

$$\frac{\partial }{\partial t}{⟨\overline{F_a}⟩}_{\overline{\mathbf{R}}}+\dot{\overline{\mathbf{R}}}\cdot \frac{\partial }{\partial \overline{\mathbf{R}}}{⟨\overline{F_a}⟩}_{\overline{\mathbf{R}}}+\dot{\overline{v_{\parallel }}}\frac{\partial }{\partial \overline{v_{\parallel }}}{⟨\overline{F_a}⟩}_{\overline{\mathbf{R}}}=\sum _b{⟨C_0({⟨\overline{F_a}⟩}_{\overline{\mathbf{R}}},{⟨\overline{F_b}⟩}_{\overline{\mathbf{R}}})⟩}_{\overline{\mathbf{R}}},$$

Full text

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 refueling and the removal of fusion ashes (Ricci, 2015). Understanding the plasma dynamics in the periphery is therefore crucial for the success of the whole fusion program (Connor et al., 1998).

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, $\omega\ll\Omega_{i}$, and in the limits $k_{\parallel}\lambda_{mfpe}\ll 1$ and $k_{\perp}\rho_{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., 2009; Tamain et al., 2009; Ricci et al., 2012; Halpern et al., 2016; Stegmeir et al., 2016; Zhu et al., 2018; Paruta et al., 2018), 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_{\perp}\rho_{i}\sim 1$ fluctuations are important (Hahm et al., 2009). This is especially relevant in the high-temperature tokamak H-mode regime (Zweben et al., 2007), the regime of operation relevant for ITER and future devices. Despite recent progress (Chang et al., 2017; Shi et al., 2017; Pan et al., 2018), 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_{\perp}\rho_{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., 2009).

Since the first formulations of the gyrokinetic theory, there have been significant research efforts to take collisions into account (Catto & Tsang, 1977; McCoy et al., 1981; Start, 2002; Brizard, 2004; Abel et al., 2008b; Barnes et al., 2009; Li & Ernst, 2011; Dorf et al., 2012; Estève et al., 2015; Hakim et al., 2019; Pan & Ernst, 2019). The first effort devoted to a gyrokinetic collision operator can be traced back to the work of Catto & Tsang (1977), later improved by Abel et al. (2008a) 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 (2011) and Madsen (2013). However, as turbulence in the tokamak periphery is essentially nonlinear, the relative level of fluctuations in this region being of order unity (Scott, 2002), and the level of collisions is not sufficient for a local thermalization, the distribution function may significantly deviate from a local Maxwellian distribution (Tskhakaya, 2012). 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, 2004; Sugama et al., 2015; Burby et al., 2015). 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 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 et al., 2018; Frei et al., 2019), 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, 2011; Jorge et al., 2017)] without requiring ad-hoc truncation of infinite series.

This paper is organized as follows. Section 2 derives the gyrokinetic equation and Section 3 presents the multipole expansion of the Coulomb collision operator. In Section 4, the Coulomb operator is ported to a gyrocenter coordinate system, while Section 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 Section 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}(\mathbf{x},\mathbf{v})$ is given by the Boltzmann equation

[TABLE]

with $C(f_{a},f_{b})$ the Coulomb (also known as Landau) collision operator (Landau, 1936; Rosenbluth et al., 1957). 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

[TABLE]

with

[TABLE]

and

[TABLE]

the Rosenbluth potentials satisfying $\nabla^{2}_{v}G_{b}=2H_{b}$ and $\nabla^{2}_{v}H_{b}=-4\pi f_{b}$ (Rosenbluth et al., 1957). In Eq. 2, $L_{ab}=q_{a}^{2}q_{b}^{2}\lambda_{ab}/(8\pi\epsilon_{0}^{2}m_{a}^{2})=\nu_{ab}v_{tha}^{3}/n_{b}$ is introduced, where $\lambda_{ab}$ and $\nu_{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, 2009; Frei et al., 2019). More precisely, denoting typical turbulent frequencies as $\omega\sim|\partial_{t}\log n|\sim|\partial_{t}\log T_{e}|$ with $n_{e}$ and $T_{e}$ the electron density and temperature respectively, and typical wavenumbers $\mathbf{k}=k_{\parallel}\mathbf{b}+\mathbf{k}_{\perp}$, being $\mathbf{k}\sim|\nabla\log n_{e}|\sim|\nabla\log T_{e}|$ and $\mathbf{b}=\mathbf{B}/B$ the magnetic field unit vector, we assume

[TABLE]

where $c_{s}=\sqrt{T_{e}/m_{i}}$ is the sound speed, $\rho_{s}=c_{s}/\Omega_{i}$ the sound Larmor radius, $\Omega_{i}=eB/m_{i}$ the ion gyrofrequency, and $\mathbf{v}_{\mathbf{E}}=\mathbf{E}\times\mathbf{B}/B^{2}$ the $\mathbf{E}\times\mathbf{B}$ drift velocity with $\mathbf{E}=-\nabla\phi-\partial_{t}\mathbf{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 $\phi=\phi_{0}+\phi_{1}$ (Dimits et al., 1992; Qin et al., 2006; Frei et al., 2019), i.e. into a possibly large-scale drift-kinetic component, $\phi_{0}$, satisfying

[TABLE]

and its small-amplitude gyrokinetic, $\phi_{1}$, component

[TABLE]

A similar decomposition into large and small scale fluctuations is applied to the magnetic vector potential $\mathbf{A}=\mathbf{A}_{0}+\mathbf{A}_{1}$, with $|\mathbf{A}_{1}|/|\mathbf{A}_{0}|\sim\epsilon_{\delta}$. Both $\phi_{0}$ and $\phi_{1}$ are assumed to yield a similar contribution to the total electric field

[TABLE]

Therefore, by ordering typical gradient lengths of $\phi_{1}$ to be comparable to $\rho_{s}$, we set

[TABLE]

which, using Eqs. 7 and 8, constraints typical gradient lengths of $\phi_{0}$ to be much larger than $\rho_{s}$, as

[TABLE]

In the following, we set $\epsilon_{\delta}\sim\epsilon$, which, using Eq. 5, yields

[TABLE]

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 $\epsilon_{B}\sim{\rho_{s}}/{L_{B}}$. We note that the collision operator developed here is valid for both $\epsilon_{B}\sim\epsilon^{2}$ and $\epsilon_{B}\sim\epsilon^{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., 2009). Finally, the collision frequency is ordered as

[TABLE]

with $\nu_{i}=\nu_{ii}$ the ion-ion collision frequency. For $T_{e}\sim T_{i}$, the ordering in Eq. 12 implies that $k_{\parallel}\lambda_{mfpe}\sim k_{\parallel}\lambda_{mfpi}\sim k_{\perp}\rho_{s}/\epsilon$, with $\lambda_{mfpa}=v_{tha}/\nu_{a}$.

By taking advantage of the ordering in Eqs. (5-7) and Eq. 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. ($\mathbf{x},\mathbf{v}$), to five. While linear and nonlinear gyrokinetic equations of motion were originally derived using recursive techniques (Taylor & Hastie, 1968; Rutherford & Frieman, 1968; Catto, 1978), more recent derivations of the gyrokinetic equation based on Hamiltonian Lie perturbation theory (Cary, 1981) ensure the existence of phase-space volume and magnetic moment conservation laws (Hahm, 1988; Brizard & Hahm, 2007; Hahm et al., 2009; Frei et al., 2019). 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 ($\phi_{1}$ and $A_{\parallel 1}$) are neglected (Cary & Brizard, 2009). Within this approximation, the coordinate transformation from particle phase-space coordinates ($\mathbf{x},\mathbf{v})$ to guiding-center coordinates $\mathbf{Z}=(\mathbf{R},v_{\parallel},\mu,\theta)$ is derived, where $\mathbf{R}$ is the guiding-center, $v_{\parallel}$ the parallel velocity, $\mu$ the magnetic moment, and $\theta$ the gyroangle. The second step introduces small-scale and small-amplitude electromagnetic fluctuations, $\phi_{1}$ and $\mathbf{A}_{1}$. For this purpose, a gyrocenter coordinate system $\overline{\mathbf{Z}}=(\overline{\mathbf{R}},\overline{v}_{\parallel},\overline{\mu},\overline{\theta})$ is constructed perturbatively from the guiding-center coordinates $\mathbf{Z}$ via a transformation $T$ of the form

[TABLE]

where $\mathbf{Z}_{1}$ contains terms proportional to $\phi_{1}$ and $\mathbf{A}_{1}$, such that $\overline{\mu}=T\mu=\mu+\epsilon_{\delta}\mu_{1}+...$, remains an adiabatic invariant [see, e.g., Brizard & Hahm (2007)]. 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}(\mathbf{x},\mathbf{v})$ is first expressed in terms of the guiding-center coordinates $\mathbf{Z}$ by defining the guiding centre distribution function $F_{a}(\mathbf{Z})$ as

[TABLE]

The coordinate transformation $\mathbf{v}(\mathbf{Z})$ is given by

[TABLE]

with $(\mathbf{b},\mathbf{e}_{1},\mathbf{e}_{2})$ a fixed right-handed coordinate set, $v_{E0}=-\nabla\phi_{0}$ the drift-kinetic $\mathbf{E}\times\mathbf{B}$ drift, $\cos\varphi=v_{\parallel}/v$ the cosine of the pitch angle and $\theta$ the gyroangle. The magnetic moment $\mu$ is defined as

[TABLE]

whereas the particle position $\mathbf{x}(\mathbf{Z})$ is written as

[TABLE]

with

[TABLE]

the Larmor radius and $\mathbf{R}$ the guiding-center of the particle. The Jacobian of the guiding-center transformation of Eqs. (15-17) is given by $B_{\parallel}^{*}/m_{a}=(B/m_{a})(1{+\mathbf{b}\cdot\nabla\times\mathbf{v}_{E}/\Omega_{a}}+v_{\parallel}\mathbf{b}\cdot\nabla\times\mathbf{b}/\Omega_{a})$ (Cary & Brizard, 2009).

We note that, in a weak-flow regime where $\mathbf{v}_{E}$ is absent from Eq. 15, the calculation that follows remains valid, except for the Jacobian of the guiding-center transformation $B_{\parallel}^{*}$, which should be replaced by $B_{\parallel}^{*}=(B/m_{a})(1+v_{\parallel}\mathbf{b}\cdot\nabla\times\mathbf{b}/\Omega_{a})$.

To account for the small scale fluctuations and magnetic inhomogeneity, the gyrokinetic distribution function $\overline{F}_{a}(\overline{\mathbf{Z}})$ is then defined as (Brizard & Hahm, 2007)

[TABLE]

with the coordinate transformation between $\overline{\mathbf{Z}}$ and ${\mathbf{Z}}$ given perturbatively by Eq. 13. Indeed, using the chain rule to rewrite the Boltzmann equation, Eq. 1, in terms of gyrocenter $\overline{\mathbf{Z}}$ coordinates, we obtain

[TABLE]

We now introduce the gyroaverage operator $\left<...\right>_{\overline{\mathbf{R}}}$ defined by

[TABLE]

where all gyrocenter coordinates $\overline{\mathbf{Z}}$ but $\overline{\theta}$ are kept fixed during the integration. By applying the gyroaverage operator to Eq. 20, the gyrokinetic equation is obtained

[TABLE]

Equation (22) can be further simplified by noting that, in the gyrokinetic framework, the transformation in Eq. 13 is constructed in such a way that the gyrocenter equations of motion, i.e. the equations that determine $\dot{\overline{\mathbf{Z}}}=(\dot{\overline{\mathbf{R}}},\dot{\overline{v_{\parallel}}},\dot{\overline{\mu}},\dot{\overline{\theta}})$, are gyroangle independent and that $\overline{\mu}$ is an adiabatic invariant satisfying $\dot{\overline{\mu}}=0$ (Brizard & Hahm, 2007). Therefore, the gyrokinetic equation in Eq. 22 can be written as

[TABLE]

In order to further simplify Eq. 23, we estimate the order of magnitude of the gyrophase dependent part of the distribution function $\widetilde{\overline{F_{a}}}=\overline{F_{a}}-\left<\overline{F_{a}}\right>_{\overline{\mathbf{R}}}$, where $\overline{F_{a}}$ obeys Eq. 20 and $\left<\overline{F_{a}}\right>_{\overline{\mathbf{R}}}$ obeys Eq. 23. For this purpose, we note that the equation for the evolution of $\widetilde{\overline{F_{a}}}=\overline{F_{a}}-\left<\overline{F_{a}}\right>_{\overline{\mathbf{R}}}$ can be obtained by subtracting Eq. 23 from Eq. 20, yielding

[TABLE]

To lowest order, $\dot{\overline{\theta}}\partial_{\theta}\widetilde{\overline{F_{a}}}\sim\Omega_{a}\widetilde{\overline{F_{a}}}$ and $\partial_{t}\sim\dot{\overline{\mathbf{R}}}\cdot\nabla_{\overline{\mathbf{R}}}\sim\dot{\overline{v_{\parallel}}}\partial_{{\overline{v_{\parallel}}}}\sim\epsilon\Omega_{i}$. Therefore, the leading order estimate of Eq. 24 gives

[TABLE]

Using the fact that $C(\overline{F_{a}},\overline{F_{b}})\sim\nu_{a}\overline{F}_{a}$, together with Eq. 12, and expanding $\overline{F_{a}}$ as $\overline{F_{a}}=\left<\overline{F_{a}}\right>_{\overline{\mathbf{R}}}+\epsilon_{\nu}\overline{F_{a1}}+...$, we find that

[TABLE]

and

[TABLE]

showing that, up to second order in $\epsilon$, the gyroangle dependence of the distribution function can be neglected in Eq. 23. We remark that a similar estimate for the gyrophase dependent part of the guiding-center distribution function was found in Jorge et al. (2017).

We now evaluate the magnitude of the collisional term in Eq. 24. Using the expansion $C(\overline{F_{a}},\overline{F_{b}})=C_{0}(\overline{F_{a}},\overline{F_{b}})+\epsilon_{\delta}C_{1}(\overline{F_{a}},\overline{F_{b}})+...$ with $C_{0}(\overline{F_{a}},\overline{F_{b}})\sim\nu_{a}\overline{F_{a}}$, and noting that the first order gyrocenter transformation $\mathbf{Z}_{1}=\mathbf{Z}-\overline{\mathbf{Z}}+O(\epsilon_{\delta}^{2})$ in Eq. 13 is mass dependent, i.e. $\overline{\mathbf{Z}_{1e}}\sim\sqrt{m_{e}/m_{i}}\overline{\mathbf{Z}_{1i}}$ [see, e.g., Brizard & Hahm (2007)], the magnitude of the Coulomb collision operator for electrons can be estimated as

[TABLE]

Thus, the third order term in the expansion in Eq. 28 does not contain a $\sqrt{m_{i}/m_{e}}$ factor, in contrast with its lower order counterpart. A similar argument holds for the ions, yielding

[TABLE]

Equations (28) and (29) show that the lowest order collision operator $C_{0}(\overline{F_{a}},\overline{F_{b}})$ is, in fact, $O(\epsilon^{2})$. Therefore, the gyrokinetic equation valid up to second order in $\epsilon$, considered in a large number of edge gyrokinetic models [see, e.g., Qin et al. (2006, 2007); Hahm et al. (2009); Frei et al. (2019)], can thus be written as

[TABLE]

assuming $\dot{\overline{\mathbf{R}}}$ and $\dot{\overline{v_{\parallel}}}$ to be at least $O(\epsilon^{2})$ accurate. We note that although only the lowest order in $\epsilon_{\delta}$ collision operator $C_{0}(\left<\overline{F_{a}}\right>_{\overline{\mathbf{R}}},\left<\overline{F_{b}}\right>_{\overline{\mathbf{R}}})$ is used in Eq. 30, all orders in $k_{\perp}\rho_{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 Eq. 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 Eq. 3 is analogous to the expression of the electrostatic potential due to a charge distribution, a similarity already noted by Rosenbluth et al. (1957). This fact allows us to make use of known electrostatic expansion techniques (Jackson, 1998) to perform a multipole expansion of the Rosenbluth potentials. We first Taylor expand the factor $1/|\mathbf{v}-\mathbf{v}^{\prime}|$ in Eq. 3 around $\mathbf{v}=0$ if $v\leq v^{\prime}$ or around $\mathbf{v}^{\prime}=0$ if $v>v^{\prime}$, yielding

[TABLE]

where we used the identity $\partial_{\mathbf{v}}(1/|\mathbf{v}-\mathbf{v}^{\prime}|)_{v=0}=-\partial_{\mathbf{v}^{\prime}}(1/v^{\prime})$ and where we denote the inner product between all the $l$ indices of two $l$-rank tensors, $\mathbf{T}_{1}^{l}$ and $\mathbf{T}_{2}^{l}$, as $\mathbf{T}_{1}^{l}\cdot\mathbf{T}_{2}^{l}$. Both $v\leq v^{\prime}$ and $v>v^{\prime}$ cases are included in order to take into account the fact that $f_{b}(\mathbf{v}^{\prime})$ is, in general, finite over the entire velocity space $\mathbf{v}^{\prime}$. Denoting $\mathbf{Y}^{l}(\mathbf{v})$ the spherical harmonic tensor (Weinert, 1980)

[TABLE]

we obtain the following expression for $H_{b}$

[TABLE]

In order to simplify Eq. 33, we note that the tensor $\mathbf{Y}^{l}(\mathbf{v})=Y_{\alpha\beta...\gamma}^{l}(\mathbf{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_{\alpha\beta...\gamma}^{l}(\mathbf{v})$ is interchangeable as the velocity derivatives commute for $v\not=0$. The totally traceless feature (i.e. $\sum_{\alpha}Y_{\alpha\alpha...\gamma}^{l}(\mathbf{v})=0$ and the same for any other pair of indices), stems from the fact that the contraction between any two indices in $Y_{\alpha\beta...\gamma}^{l}(\mathbf{v})$ leads to the multiplicative factor $\nabla_{\mathbf{v}}^{2}=\partial_{\mathbf{v}}\cdot\partial_{\mathbf{v}}(1/v)$, which vanishes for $v\not=0$ (we note that $Y^{l}(\mathbf{v})$ vanishes in the limit $\mathbf{v}=0$). Furthermore, by defining the tensor $(\mathbf{v})^{l}_{\text{TS}}$ as the traceless symmetric counterpart of $(\mathbf{v})^{l}$ [e.g., $(\mathbf{v})^{2}_{\text{TS}}=\mathbf{v}\mathbf{v}-\mathbf{I}v^{2}/3$ with $\mathbf{I}$ the identity matrix], we replace the tensors $(\mathbf{v}^{\prime})^{l}$ and $(\mathbf{v})^{l}$ in Eq. 33 by their traceless symmetric counterpart $(\mathbf{v}^{\prime})^{l}_{\text{TS}}$ and $(\mathbf{v})^{l}_{\text{TS}}$ respectively

[TABLE]

as they differ only by terms proportional to the identity matrix that vanish when summed with $\mathbf{Y}^{l}(\mathbf{v})$ and $\mathbf{Y}^{l}(\mathbf{v}^{\prime})$. In fact, for the $l=2$ case we have $(\mathbf{v}^{2}-(\mathbf{v})^{2}_{\text{TS}})\cdot\mathbf{Y}^{2}(\mathbf{v})=(v^{2}/3)\mathbf{I}\cdot\mathbf{Y}^{2}(\mathbf{v})=(v^{2}/3)\sum_{\alpha}Y^{2}_{\alpha\alpha}=0$, and similarly for $l>2$. In addition, following the convention in Snider (2017), scalars ($l=0$) and vectors ($l=1$) are considered to be traceless symmetric quantities. Finally, we relate the tensors $(\mathbf{v})^{l}_{\text{TS}}$ and $\mathbf{Y}^{l}(\mathbf{v})$. For $l=0$ and $l=1$, we have $\mathbf{Y}^{0}(\mathbf{v}^{\prime})=1=(\mathbf{v}^{\prime})^{0}_{\text{TS}}$ and $\mathbf{Y}^{1}(\mathbf{v}^{\prime})=\mathbf{v}^{\prime}=(\mathbf{v}^{\prime})^{1}_{\text{TS}}$. For $l=2$, Eq. 32 gives

[TABLE]

The results obtained for $l=0,1$, and 2 can be generalized, i.e. $(\mathbf{v}^{\prime})^{l}_{\text{TS}}=\mathbf{Y}^{l}(\mathbf{v}^{\prime})$ as proved by induction (Weinert, 1980). The Rosenbluth potential $H_{b}$ can therefore be written as

[TABLE]

The first term in Eq. 36 can be regarded as the potential due to the charge distribution $f_{b}(\mathbf{v}^{\prime})$ inside a sphere of radius $v$, while the second term is the potential due to a finite charge distribution $f_{b}(\mathbf{v}^{\prime})$ at $v^{\prime}\geq v$.

We now look for an expansion of $f_{b}$ that allows us to perform the integrals in Eq. 36 analytically by writing $H_{b}$ as a sum of velocity moments of $f_{b}$. We consider the basis functions (Hirshman & Sigmar, 1976)

[TABLE]

with $L_{k}^{l+1/2}(v)$ an associated Laguerre polynomial (Abramowitz et al., 1965), i.e.

[TABLE]

where

[TABLE]

The basis $\mathbf{Y}^{lk}(\mathbf{v})$ is orthogonal, being the orthogonality relation given by (Banach & Piekarski, 1989; Snider, 2017)

[TABLE]

with $\mathbf{T}^{lk}$ an arbitrary symmetric and traceless tensor, and $\sigma_{k}^{l}$ the normalization constant

[TABLE]

A proof that $\mathbf{Y}^{lk}(\mathbf{v})$ is a complete basis, i.e. that each $l$ and $k$ element of $\mathbf{Y}^{lk}(\mathbf{v})$ is linearly independent and that a linear combination of its elements spans any smooth function $f(\mathbf{v})$, can be found in Banach & Piekarski (1989), where the equivalence between Grad’s moment expansion in tensorial Hermite polynomials (which forms a complete basis) and $\mathbf{Y}^{lk}(\mathbf{v})$ is shown. We then write $f_{b}$ as

[TABLE]

with $f_{Mb}$ a Maxwellian distribution function

[TABLE]

According to Eq. 40, with ${\mathbf{T}^{lk}}=\mathbf{M}^{lk}_{b}$, the coefficients $\mathbf{M}^{lk}_{b}$ are obtained by taking velocity moments of $f_{b}$ of the form

[TABLE]

Finally, we note that Eq. 42 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 Eq. 42 into Eq. 36, we obtain the following expression

[TABLE]

where we define the normalized velocity $x_{b}=v^{2}/v_{thb}^{2}$, the solid angle $\sigma$ such that $d\mathbf{v}=v^{2}dvd\sigma$, and use the relation $\mathbf{Y}^{l}\mathbf{(}v)=v^{l}\mathbf{Y}^{l}(\hat{v})$ with $\mathbf{v}=v\hat{v}$ (Weinert, 1980). Applying the orthogonality relation of Eq. 40 for $k=0$, and expanding the associated Laguerre polynomials using Eq. 38, we write $H_{b}$ as

[TABLE]

where the identity

[TABLE]

is used to simplify Eq. 46.

We note that the expression of $H_{b}$ in Eq. 46 corresponds to the one in Ji & Held (2006), having replaced the $\mathbf{Y}^{l}(\mathbf{v})$ tensors by the $\mathbf{P}^{l}(\mathbf{v})$ tensors defined by the recursion relation [see Eq. (14) of Ji & Held (2006)]

[TABLE]

with $\mathbf{P}^{0}(\mathbf{v})=1$ and $\mathbf{P}^{1}(\mathbf{v})=\mathbf{v}$. We can indeed prove that $\mathbf{Y}^{l}(\mathbf{v})=\mathbf{P}^{l}(\mathbf{v})$ by deriving the tensor $\mathbf{Y}^{l}(\mathbf{v})$ using Eq. 32. This yields

[TABLE]

which is the same recursion relation present in Eq. 48. Since $\mathbf{Y}^{0}(\mathbf{v})=\mathbf{P}^{0}(\mathbf{v})$ and $\mathbf{Y}^{1}(\mathbf{v})=\mathbf{P}^{1}(\mathbf{v})$, the proof is complete.

The integrals in Eq. 46 can be put in terms of upper

[TABLE]

and lower

[TABLE]

incomplete gamma functions (Abramowitz et al., 1965), yielding

[TABLE]

A procedure similar to the one used to obtain Eq. 52 can be followed for the second Rosenbluth potential $G_{b}$ by expanding the distribution function $f_{b}$ appearing in $G_{b}$ according to Eq. 42, therefore obtaining

[TABLE]

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 Eq. 2 in terms of a single spherical harmonic tensor $\mathbf{Y}^{l}(\mathbf{v})$, we make use of the following identity between symmetric traceless tensors (Ji & Held, 2009)

[TABLE]

where $\cdot^{n}$ is the $n$-fold inner product [e.g., for the matrix $\mathbf{A}=A_{ij}$, $(\mathbf{A}\cdot^{1}\mathbf{A})_{ij}=\sum_{k}A_{ki}A_{kj}$]. The $d_{j}^{l,n}$ coefficient can be written in terms of the $t_{j}^{l,n}$ coefficient

[TABLE]

as

[TABLE]

In Eq. 56, the summation range involves the integer partitions of $j$, i.e., the decomposition of $j$ into different sums of $h$ positive integers, here labeled 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 Eq. 42, the expression for the Rosenbluth potentials in Eqs. 52 and 53, and the identity in Eq. 54, the collision operator in Eq. 2 can be rewritten in terms of products of $\mathbf{M}_{a}^{lk}$ and $\mathbf{M}_{b}^{lk}$ as

[TABLE]

with

[TABLE]

The quantity $\nu_{*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 (2009).

4 Gyrokinetic Coulomb Collision Operator

In Section 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 $\left<\overline{F_{a}}\right>$ and $\left<\overline{F_{b}}\right>$. As a first step, the gyroangle dependence of the basis functions $\mathbf{Y}^{lk}$ is found explicitly by using a coordinate transformation from the particle phase-space coordinates $(\mathbf{x},\mathbf{v})$ to the guiding-center coordinate system $\mathbf{Z}$. This allows us to decouple the fast gyromotion time associated with the gyroangle $\theta$ from the typical turbulence time scales. The multipole moments $\mathbf{M}_{a}^{lk}$ and $\mathbf{M}_{b}^{lk}$ can then be written in terms of moments of the guiding-center distribution function $\left<F_{a}\right>$ and $\left<F_{b}\right>$ for arbitrary values of $k_{\perp}\rho_{s}$. As a second step, the gyrocenter coordinate system $\overline{\mathbf{Z}}$ is introduced by using the coordinate transformation $T$ in Eq. 13. As shown in Section 2, for a gyrokinetic equation up to second order accurate in $\epsilon_{\delta}$, 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-center one by a simple coordinate relabeling.

We first derive the polar and azimuthal angle (gyroangle) dependence of the $\mathbf{Y}^{l}(\mathbf{v})$ tensor in terms of scalar spherical harmonics. This is useful to analytically perform the gyroaverage of the collision operator in the Boltzmann equation, Eq. 30. For this purpose, as a first step, we show that the Laplacian of $\mathbf{Y}^{l}(\mathbf{v})$ vanishes, i.e. that $\mathbf{Y}^{l}(\mathbf{v})$ are harmonic tensors. By applying the operator $\nabla_{\mathbf{v}}^{2}$ to $\mathbf{Y}^{l}(\mathbf{v})$, and recalling that $\nabla_{\mathbf{v}}^{2}(1/v)=0$ for $v\not=0$, we obtain

[TABLE]

since

[TABLE]

as can be proved by induction (Weinert, 1980). The angular dependence of $\mathbf{Y}^{l}(\mathbf{v})$ can now be found by expressing the Laplacian of Eq. 59 in spherical coordinates. Using the fact that $\mathbf{Y}^{l}(\mathbf{v})=v^{l}\mathbf{Y}^{l}(\hat{v})$, we obtain

[TABLE]

where $L^{2}$ is the angular part of the operator $\nabla^{2}_{\mathbf{v}}$ multiplied by $v^{2}$

[TABLE]

with $\varphi$ and $\theta$ chosen, respectively, as the pitch angle and the gyroangle variables, both defined in Eq. 15. Evaluating the $v$ derivatives in Eq. 61, the following differential equation for $\mathbf{Y}^{l}(\mathbf{v})$ is obtained

[TABLE]

We identify Eq. 63 as the eigenvalue equation for the scalar spherical harmonics $Y_{lm}(\varphi,\theta)$ (Arfken et al., 2013), which can be written in terms of associated Legendre polynomials $P_{l}^{m}(\cos\varphi)$ as (Abramowitz et al., 1965)

[TABLE]

with

[TABLE]

and $P_{l}(x)=(d^{l}/dx^{l})[(x^{2}-1)^{l}]/(2^{l}l!)$ a Legendre polynomial. Therefore, using Eq. 63, and denoting $\mathbf{e}^{lm}$ the basis elements of $\mathbf{Y}^{l}(\mathbf{v})$ (an elementary derivation of the basis tensors $\mathbf{e}^{lm}$ is shown in Appendix A), we write $\mathbf{Y}^{l}(\mathbf{v})$ as

[TABLE]

Having derived the gyroangle dependence of the $\mathbf{Y}^{l}(\mathbf{v})$ tensors, we now compute the fluid moments $\mathbf{M}_{a}^{lk}$ in terms of $v_{\parallel}$ and $\mu$ moments of the guiding-center distribution function $\left<F_{a}\right>$. In order to perform the velocity integration in the definition of the moments $\mathbf{M}^{lk}$ in Eq. 44 at arbitrary $k_{\perp}\rho$ in guiding-center phase-space coordinates, we use the identity $f(\mathbf{x},\mathbf{v})=\int f(\mathbf{x}^{\prime},\mathbf{v})\delta(\mathbf{x}-\mathbf{x}^{\prime})d\mathbf{x}^{\prime}$. By imposing $\mathbf{x}^{\prime}=\mathbf{R}+\mathbf{\rho}$, writing the volume element in phase-space as $d\mathbf{x}^{\prime}d\mathbf{v}=(B_{\parallel}^{*}/m)d\mathbf{R}dv_{\parallel}d\mu d\theta$, and using Eq. 14, we obtain

[TABLE]

from Eq. 44. Expressing $\mathbf{v}=\mathbf{v}(\mathbf{Z})$, as shown by Eq. 15, and performing the integral over $\mathbf{R}$ in Eq. 67, it follows that

[TABLE]

The orderings in Eqs. 26 and 27 for the guiding-center distribution function $F_{a}$ allows us to approximate $F_{a}\simeq\left<F_{a}\right>_{\mathbf{R}}$ (Jorge et al., 2017), effectively neglecting $\epsilon^{2}$ effects in $\mathbf{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, 2011; Pan & Ernst, 2019), we represent $F_{a}(\mathbf{R},v_{\parallel},\mu,\theta)$ by its Fourier transform $F_{a}(\mathbf{k},v_{\parallel},\mu,\theta)=\int F_{a}(\mathbf{R},v_{\parallel},\mu,\theta)e^{-i\mathbf{k}\cdot\mathbf{R}}d\mathbf{R}$, and write

[TABLE]

By aligning the $\mathbf{k}$ coordinate system in the integral of Eq. 69 with the axes $(\mathbf{b},\mathbf{e}_{1},\mathbf{e}_{2})$, i.e. $\mathbf{k}=k_{\parallel}\mathbf{b}+k_{\perp}(\cos\theta\mathbf{e}_{1}+\sin\theta\mathbf{e}_{2})$, we write $\exp(-i\mathbf{k}\cdot\mathbf{\rho})=\exp(-ik_{\perp}\rho\cos\theta)$. We then use the Jacobi-Anger expansion (Andrews, 1992)

[TABLE]

with $J_{p}$ the Bessel function of order $p$, and rewrite Eq. 69 as

[TABLE]

The velocity $\mathbf{v}$ in the argument of $\mathbf{Y}^{lk}$ in Eq. 71 is then expanded as

[TABLE]

The second term in Eq. 72 introduces $\epsilon_{B}\ll\epsilon$ 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 (2019).

Using Eq. 66 to express $\mathbf{Y}^{lk}(\mathbf{v})$ in terms of spherical harmonics, we perform the gyroangle integration in Eq. 71. By rewriting the spherical harmonics $Y_{lm}(\varphi,0)$ in terms of associated Legendre polynomials $P_{l}^{m}(\cos\varphi)$ using Eq. 64, the gyroaverage of the product $Y^{lk}(\mathbf{v}/v_{tha})e^{-ip\theta}$ can be performed, yielding

[TABLE]

with

[TABLE]

We note that, the $p=0$ case of Eq. 74 corresponds to the gyroaveraged formulas in Ji et al. (2009, 2013); Ji & Held (2014) used to derive closures for fluid models at zeroth order in $\epsilon$. Finally, by defining the Bessel-Fourier operator

[TABLE]

the expression for the fluid moments $\mathbf{M}_{a}^{lk}$ in terms of coupled $v_{\parallel}$ and $\mu$ moments of the guiding-center distribution function $\left<F_{a}\right>_{\mathbf{R}}$ is obtained

[TABLE]

with

[TABLE]

Equation (76) can now be used to express the collision operator $C(f_{a},f_{b})$ in terms of $v_{\parallel}$ and $\mu$ integrals of $\left<F_{a}\right>$. Using Eqs. 58, 66 and 57, and defining

[TABLE]

we can write the collision operator in Eqs. 57 and 58 as a function of the $\mathcal{M}^{lk}_{am}$ moments, i.e.

[TABLE]

with

[TABLE]

We now focus on the gyroaverage of the collision operator in Eq. 79 with the gyroaverage operation performed at constant $\mathbf{R}$. We first note that the gyroangle $\theta$ dependence in $c^{lkmnqr}_{ab}$ is present through the spherical harmonic $Y_{l+n-2j-2u~{}v}(\varphi,\theta)$ and through the fluid moments $\mathcal{M}^{lk}_{as}$ and $\mathcal{M}^{nq}_{bt}$ as the latter are functions of $\mathbf{x}=\mathbf{R}+{\mathbf{\rho}_{a}}$. To make the gyroangle dependence explicit, we write both $\mathcal{M}^{lk}_{as}$ and $\mathcal{M}^{nq}_{bt}$ in Fourier space as

[TABLE]

Using the Jacobi-Anger expansion of Eq. 70, we find that

[TABLE]

with $\alpha$ the azimuthal angle of the $\mathbf{k}^{\prime}$ vector, i.e., the angle between $\mathbf{k}_{\perp}^{\prime}$ and $\mathbf{k}_{\perp}$. The gyroaveraged collision operator at arbitrary $k_{\perp}\rho$ is therefore given by

[TABLE]

with

[TABLE]

and

[TABLE]

We note that, if only first order $k_{\perp}\rho$ terms are kept in the Fourier-Bessel operator of Eq. 75, the collision operator in Eq. 83 reduces to the drift-kinetic collision operator found in Jorge et al. (2017).

In Eq. 83, the gyroaveraged collision operator is cast in terms of $v_{\parallel}$ and $\mu$ moments of the guiding-center distribution function $\left<F_{a}\right>$ for arbitrary values of $k_{\perp}\rho$. We now apply the transformation $T$, introduced in Eq. 13, to Eq. 83 in order to write the gyroaveraged collision operator in terms of $\overline{v}_{\parallel}$ and $\overline{\mu}$ moments of the gyrocenter distribution function $\left<\overline{F_{a}}\right>_{\overline{\mathbf{R}}}$. As shown in Section 2, only the zeroth order terms in the $\epsilon_{\delta}$ expansion of $\left<C(\overline{F_{a}},\overline{F_{b}})\right>_{\overline{\mathbf{R}}}$ are needed in order to adequately describe collisional processes in the gyrokinetic framework. Therefore, using Eq. 13, we apply the zeroth order transformations $\mathbf{Z}\simeq\overline{\mathbf{Z}}$ and $F_{a}(\mathbf{Z})=T\overline{F_{a}}(\mathbf{Z})\simeq\overline{F_{a}}(\mathbf{Z})$ to the collision operator $\left<C(F_{a},F_{b})\right>$ in Eq. 83, yielding

[TABLE]

with

[TABLE]

where $\overline{\rho_{a}}=\rho_{a}(\overline{\mathbf{R}},\overline{\mu},\overline{\theta})$, $\overline{v}^{2}=\overline{v}_{\parallel}^{2}+2B\overline{\mu}/m$, and

[TABLE]

The collision operator in Eq. 86 represents the gyrokinetic full Coulomb collision operator that can be used in gyrokinetic models that are up to $O(\epsilon_{\delta}^{2})$ accurate. In Eq. 86, the integral-differential character of the $C(f_{a},f_{b})$ operator is replaced by two-dimensional integrals of the gyrocenter distribution function over the velocity coordinates $\overline{v}_{\parallel}$ and $\overline{\mu}$. We note that, in practice, a truncation of the series present in Eqs. 86 and 87 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 Eq. 86. An expansion of the drift-kinetic (Jorge et al., 2017) and gyrokinetic (Mandell et al., 2018; Frei et al., 2019) equation in Hermite-Laguerre polynomials has been recently introduced, showing that this is an advantageous approach to the study of plasma waves and instabilities (Jorge et al., 2018, 2019). 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 Eqs. 87 and 88. We therefore expand $\left<\overline{F_{a}}\right>_{\overline{\mathbf{R}}}$ as

[TABLE]

where $H_{p}$ are physicists’ Hermite polynomials of order $p$ defined by the Rodrigues’ formula

[TABLE]

and normalized via

[TABLE]

and $L_{j}$ the Laguerre polynomials of order $j$ defined by the Rodrigues’ formula

[TABLE]

and orthonormal with respect to the weight $e^{-x}$

[TABLE]

In Eq. 89, we introduce the normalized parallel velocity

[TABLE]

and the perpendicular velocity coordinate

[TABLE]

Due to the orthogonality of the Hermite-Laguerre polynomial basis, the coefficients $N_{a}^{pj}$ of the expansion in Eq. 89 can be computed as

[TABLE]

We note that the integrand of $\overline{N}_{a}^{pj}$ in Eq. 96 contains the multiplicative factor $B/m_{a}$, as opposed to the Jacobian containing the factor $B_{\parallel}^{*}/m_{a}$. In the following, we also use the Hermite-Laguerre moments of $\left<\overline{F_{a}}\right>_{\overline{\mathbf{R}}}$ with $B_{\parallel}^{*}$ in the integrand instead of B. These are denoted as $\overline{N}_{a}^{*pj}$, i.e.

[TABLE]

while in the weak flow regime the term proportional to $\mathbf{b}\cdot\nabla\times\mathbf{v}_{\mathbf{E}}/\Omega_{a}$ in Eq. 97 is set to zero.

In order to express the collision operator in terms of the moments $\overline{N}_{a}^{pj}$ given in Eq. 96 and evaluate its Hermite-Laguerre moments, we first consider the gyrokinetic moments $\overline{\mathcal{M}}_{am}^{lk}$ and write the integral that defines them in Eq. 88 as a function of the gyrocenter moments $\overline{N}_{a}^{pj}$ of Eq. 96. 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{\mu}$ and $k_{\perp}$ dependence in the Fourier-Bessel operator $j_{m}$, Eq. 75, can be decomposed by introducing $\rho_{tha}=v_{tha}/\Omega_{a}$ and noting that $|\mathbf{\rho}_{a}|=\sqrt{\overline{\mu}B/T_{a}}\rho_{tha}=\overline{s}_{\perp a}\rho_{tha}$. This allows the use of the following identity between Bessel and Legendre functions (Gradshteyn & Ryzhik, 2007)

[TABLE]

with $b_{a}=k_{\perp}\rho_{tha}/2$, $\sigma_{0}=1$ and $\sigma_{m}=\text{sgn}(m)^{m}$ for $m\not=0$. The Fourier-Bessel operator in Eq. 75, with the identity in Eq. 98 and the Hermite-Laguerre expansion of Eq. 89, can then be written as

[TABLE]

As a second step, we consider

[TABLE]

which are used in the definition in Eq. 64. In order to expand the associated Legendre polynomials $P_{l}^{m}(\cos\varphi)$ appearing in Eq. 100 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. (2017) to a transformation between an Associated Legendre-Associated Laguerre and a Hermite-Laguerre basis, that is

[TABLE]

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

[TABLE]

The gyrocenter moments $\overline{\mathcal{M}}^{lk}_{am}$ in Eq. 88 can then be rewritten using the identities in Eqs. 99 and 101 and

[TABLE]

with the $d^{r}_{mjs}$ coefficients given by

[TABLE]

yielding the following expression

[TABLE]

where we defined

[TABLE]

Using the form for $\overline{\mathcal{M}}_{am}^{lk}$ in Eq. 105, the collision operator in Eq. 86 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 Eq. 105 reduce to the ones in Jorge et al. (2017) in the lowest order drift-kinetic limit $k_{\perp}\rho_{tha}=0$.

We now take Hermite-Laguerre moments of the collision operator $\left<C(\overline{F_{a}},\overline{F_{b}})\right>$, i.e. we evaluate

[TABLE]

Writing the gyroaveraged collision operator $\left<C(\overline{F_{a}},\overline{F_{b}})\right>$ in Eq. 86 using Eqs. 105 and 87, and expanding the Bessel functions $J_{p}(k_{\perp}\rho_{a})$ and $J_{p+m}(k_{\perp}^{\prime}{\rho_{a}})$ using Eq. 98, the following form for the $\left<\overline{c}_{ab}^{lkmnqr}\right>_{\overline{\mathbf{R}}}$ term appearing in $\left<C(\overline{F_{a}},\overline{F_{b}})\right>_{\overline{\mathbf{R}}}$ is obtained

[TABLE]

In Eq. 108, we defined the $D_{abuidz{pp^{\prime}}}^{lkmnqr}$ term

[TABLE]

with $B_{a}^{pvz^{\prime}}=b_{a}^{p+2z}b_{{a}}^{p+v+2z^{\prime}}$ and $b_{{a}}^{\prime}=k_{\perp}^{\prime}\rho_{th{a}}/2$, while the convolution operator $\mathcal{N}_{abuidz}^{lkmnqr}(\mathbf{k},\mathbf{k}^{\prime})$ is given by

[TABLE]

with $\overline{\mathcal{M}}_{as}^{lk}$ the moments of the distribution function defined in Eq. 105.

Finally, the result in Eq. 108 is used in Eq. 107 in order to find the Hermite-Laguerre moments $C_{ab}^{pj}$ of the full Coulomb collision operator expressed in Eq. 86. This yields

[TABLE]

with

[TABLE]

and

[TABLE]

The integral factor $I$ can be performed analytically by first rewriting the product of two Laguerre polynomials as a single one using

[TABLE]

with

[TABLE]

expressing the resulting Hermite-Laguerre basis in terms of Legendre-Associated Laguerre using Eq. 102, and writing the phase-space volume $(B_{\parallel}^{*}/m)d\overline{v}_{\parallel}d\overline{\mu}$ as $\overline{v}^{2}d\overline{v}d\overline{\xi}$ with $\overline{\xi}=\overline{v}_{\parallel}/\overline{v}$. This yields

[TABLE]

Ji & Held (2009) present an analytical closed expression ready to be numerically implemented of the factor $C_{*abu}^{st,lm,nr}=\int f_{Ma}\nu_{*abu}^{lm,nr}(v^{2})L_{t}^{s+1/2}(v^{2})v^{s}d\mathbf{v}$. We note that the long-wavelength limit can be found by setting $d=z=0$ in the collision operator Eq. 112. This yields the Hermite-Laguerre moments of the collision operator moments found in Jorge et al. (2017).

6 Small Mass-Ratio Approximation

In this section, we simplify the electron-ion and the ion-electron collision operator in Eq. 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 $(\mathbf{x},\mathbf{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 Lorentz pitch-angle scattering operator, can be written as (Helander & Sigmar, 2005)

[TABLE]

with $\mathbf{c}_{e}=\mathbf{v}/v_{the}$. We expand $f_{e}$ according to Eq. 42. We note that the expansion in Eq. 42 is an eigenbasis of the pitch-angle scattering operator $C_{ei}$ with eigenvalue $l(l+1)$ (Ji & Held, 2008). Therefore, we write

[TABLE]

We now Fourier transform the moments $\mathbf{M}_{e}^{lk}$ in Eq. 118 as $\mathbf{M}_{e}^{lk}(\overline{\mathbf{R}})=\int\mathbf{M}_{e}^{lk}(\mathbf{k})e^{i\mathbf{k}\cdot\overline{\mathbf{R}}}d\mathbf{k}$ and write the gyroaveraged collision operator $C_{ei}$ as

[TABLE]

Using the Jacobi-Anger expansion of Eq. 70, Eq. 98, the inverse basis transformation Eq. 102, and the identities $J_{-p}(x)=(-1)^{p}J_{p}(x)$ and

[TABLE]

we obtain

[TABLE]

with $b_{e}=k_{\perp}\rho_{the}/2$. The collision operator in Eq. 119 represents the gyrokinetic electron-ion collision operator.

Equation (121) provides an expression of the pitch-angle scattering operator $\left<C_{ei}\right>$ in Eq. 119 suitable for projection onto a Hermite-Laguerre basis, i.e.

[TABLE]

where we define

[TABLE]

An analytical form for the integral factor $I_{ei}^{lk}$ can be derived using the expression for $\left<C_{ei}\right>$, Eq. 119, and Eq. 121, yielding

[TABLE]

with $I_{Lkt}^{lsuv}$ and $I_{Pm}^{lus}$ defined by

[TABLE]

and

[TABLE]

respectively. The electron fluid moments $\mathbf{M}_{e}^{lk}$ can be cast in terms of Hermite-Laguerre moments $\overline{N}_{e}^{lk}$ using the expressions in Eqs. (76), (97), and (105). The factor $I_{Lkt}^{lsuv}$ can be analytically evaluated by expanding the associated Laguerre polynomials using Eq. 38, which leads to

[TABLE]

Similarly, the factor integral $I_{Pm}^{lus}$ can be calculated using an extended version of Gaunt’s formula (Gaunt, 1929), yielding (Mavromatis & Alassar, 1999)

[TABLE]

We note that, in Eq. 128, the Wigner 3-j symbol $\begin{pmatrix}j_{1}&j_{2}&j_{3}\\ m_{1}&m_{2}&m_{3}\end{pmatrix}$ is related to the Clebsch-Gordan coefficients $\left<j_{1}m_{1}j_{2}m_{2}|j_{3}m_{3}\right>$ via (Olver et al., 2010)

[TABLE]

with the Clebsch-Gordan coefficients given by

[TABLE]

where the summation in Eq. 130 is extended over all integers $k$ that make every factorial in the sum nonnegative (Bohm & Loewe, 1993).

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, 2005)

[TABLE]

where the electron-ion friction force is neglected for simplicity. We simplify Eq. 131 by using Eq. 27, therefore approximating the distribution function $f_{i}$ by its gyroaveraged component $f_{i}\simeq\left<\overline{F}_{i}\right>_{\overline{\mathbf{R}}}$, and retaining the lowest-order terms in the $\epsilon_{\delta}$ expansion. This allows us to convert the $C_{ie}$ operator in Eq. 131 to the gyrocenter variables $\overline{\mathbf{Z}}$ using the chain rule at lowest order in $\epsilon_{\delta}$, i.e. to express Eq. 131 in $\mathbf{Z}$ coordinates using the guiding-center transformation in Eqs. 15, 16 and 17 and approximate $\mathbf{Z}\simeq\overline{\mathbf{Z}}$. The velocity derivatives can be written as a function of $\overline{\mathbf{Z}}$ using the chain rule, yielding

[TABLE]

where we define $\mathbf{c}=(\cos\overline{\theta}\mathbf{e}_{1}+\sin\overline{\theta}\mathbf{e}_{2})$ and $\mathbf{a}=\mathbf{c}\times\mathbf{b}=(-\sin\overline{\theta}\mathbf{e}_{1}+\cos\overline{\theta}\mathbf{e}_{2})$. The ion-electron collision operator can therefore be written as

[TABLE]

We now Fourier transform both $T_{e}$ and $\left<\overline{F}_{i}\right>_{\overline{\mathbf{R}}}$ and gyroaverage $C_{ie}$, yielding

[TABLE]

where we have used the identities $\left<\mathbf{a}e^{i\mathbf{k}^{\prime}\cdot\mathbf{\overline{}}{\rho_{i}}}\right>_{\overline{\mathbf{R}}}=iJ_{1}(k_{\perp}^{\prime}\overline{\rho_{i}})\mathbf{e}_{2}$ and $\left<\mathbf{a}\mathbf{a}e^{i\mathbf{k}^{\prime}\cdot\mathbf{\overline{}}{\rho_{i}}}\right>_{\overline{\mathbf{R}}}=(1/2)[J_{0}(k_{\perp}^{\prime}\overline{\rho_{i}})\left(\mathbf{e}_{1}\mathbf{e}_{1}+\mathbf{e}_{2}\mathbf{e}_{2}\right)+J_{2}(k_{\perp}^{\prime}\overline{\rho_{i}})\left(\mathbf{e}_{1}\mathbf{e}_{1}-\mathbf{e}_{2}\mathbf{e}_{2}\right)]$. Finally, we take Hermite-Laguerre moments of the gyroaveraged ion-electron collision operator $\left<C_{ie}\right>_{\overline{\mathbf{R}}}$ in Eq. 134, using the expansion of $\left<\overline{F}_{i}\right>_{\overline{\mathbf{R}}}$ in Eq. 89, yielding

[TABLE]

with $A_{lk}^{pj}$ given by

[TABLE]

and $B_{lkrs}^{pj}$ by

[TABLE]

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_{\perp}\rho_{s}$, yielding the formula in Eq. 86. 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 Eq. 111. Furthermore, the electron-to-ion mass ratio is used to simplify the form of the electron-ion and ion-electron collision operators, yielding Eq. 122 and Eq. 135, respectively.

We conclude by noting that the present collision operator is derived by porting the Coulomb operator to the gyrocenter phase-space by using a framework valid up to second order in the expansion parameter $\epsilon$, 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, 2004; Sugama et al., 2015; Burby et al., 2015).

8 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 $\mathbf{e}^{lm}$ used in the definition of $\mathbf{Y}^{l}(\mathbf{v})$ in Eq. 66. We start with the $l=1$ case, for which Eq. 66 yields

[TABLE]

The spherical basis vectors $\mathbf{e}^{1m}$ can then be derived from Eq. 138 by expressing the vector $\mathbf{v}$ in spherical coordinates as

[TABLE]

and using the identities for the spherical harmonics

[TABLE]

therefore obtaining

[TABLE]

We now construct spherical basis tensors $\mathbf{e}^{lm}$ from the spherical basis vectors $\mathbf{e}^{1m}$ leveraging the techniques developed for the angular momentum formalism in quantum mechanics (Zettili & Zahed, 2009; Snider, 2017). As a first step, we note that the basis vectors $\mathbf{e}^{1m}$ are eigenvectors of the angular momentum matrix $G_{z}$

[TABLE]

with eigenvalue $m$, that is

[TABLE]

As a second step, we note that the relationship between the basis vectors $\mathbf{e}_{\alpha}$ for $\alpha=(x,y,z)$ and the angular momentum matrices $G_{\alpha}$ is given by

[TABLE]

with $\epsilon$ the standard Levi-Civita tensor. In index notation, Eq. 144 can be written as

[TABLE]

The raising $G_{+}$ and lowering $G_{-}$ operators (corresponding to the ladder operators in quantum mechanics), defined by

[TABLE]

allows us to obtain the basis vectors $\mathbf{e}^{1\pm 1}$ from $\mathbf{e}^{10}$ using

[TABLE]

In addition, we have that

[TABLE]

We can now define the spherical tensor basis $\mathbf{e}^{lm}$ that define the irreducible tensors $\mathbf{Y}^{l}$. We start with the spherical basis tensor

[TABLE]

formed by the product of $l$ basis vectors $\mathbf{e}^{11}$. Similarly to $\mathbf{Y}^{l}(\mathbf{v})$, this tensor is of rank $l$, symmetric, and totally traceless, as $\mathbf{e}^{11}\cdot\mathbf{e}^{11}=0$. Furthermore, we note that $\mathbf{e}^{ll}$ is an eigenvector with eigenvalue $l$ of the angular momentum tensor $G_{z}^{l}$, with $G_{n}^{l}$ a tensor of rank $2l$ defined by

[TABLE]

The remaining basis tensor elements $\mathbf{e}^{lm}$ can be obtained by applying the tensorial lowering operator $G^{l}_{-}=G_{x}^{l}-iG_{y}^{l}$ to $\mathbf{e}^{ll}$, namely

[TABLE]

with $m=-l,-l+1,...,-1,0,1,...,l$ and $\left(G^{l}_{-}\cdot^{l}\right)^{l-m}\mathbf{e}^{ll}$ a tensor of order $l$ built by the application of the $G^{l}_{-}\cdot^{l}$ operator to $\mathbf{e}^{ll}$ $l-m$ times. The normalization factor in Eq. 151 is obtained by requiring that the contravariant $\mathbf{e}^{lm}$ and the covariant $\mathbf{e}^{l}_{m}$ basis tensors form an orthonormal basis, i.e.

[TABLE]

In order to find a covariant basis $e^{l}_{m}$, we start with the case $l=1$ and note that the set of vectors $\mathbf{e}^{1}_{m}=(\mathbf{e}^{1m})^{*}=(-1)^{m}\mathbf{e}^{1-m}$, with $(\mathbf{e}^{1}_{m})^{*}$ the complex conjugate of $\mathbf{e}^{1}_{m}$ satisfies Eq. 152. We therefore define $\mathbf{e}^{l}_{m}=(\mathbf{e}^{lm})^{*}$, and use Eq. 152 to normalize $\mathbf{e}^{lm}$. For computational purposes, we note that the tensor $\mathbf{e}^{lm}$ can also be written as a function of the basis vectors $\mathbf{e}^{1m}$ as (Snider, 2017)

[TABLE]

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})^{lkm}_{pj}$ coefficients defined in Eqs. 101 and 102. By multiplying Eq. 101 by a Hermite and a Laguerre polynomial and by an exponential of the form $e^{-\overline{v}^{2}}$, and integrating over the whole $\overline{v}_{\parallel}$ and $\overline{\mu}$ space, we obtain the following integral expression for $T_{lkm}^{pj}$

[TABLE]

We first write the integrand in Eq. 154 in terms of $\overline{\xi}=\overline{v}_{\parallel}/\overline{v}$ and $\overline{v}$ coordinates using the basis transformation in Eq. 102, yielding

[TABLE]

where we used the fact that $(T^{-1})_{lk}^{pj}=T_{lk}^{pj}\sqrt{\pi}2^{p}p!k!(l+1/2)/(k+l+1/2)!$ (Jorge et al., 2017). The $\xi$ integral in Eq. 155 is performed by expanding $P_{l}$ as a finite sum of the form

[TABLE]

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)$

[TABLE]

The $x$ integral in Eq. 155 is performed by using the expansion of the associated Laguerre polynomials in Eq. 38. The $T_{lkm}^{pj}$ coefficient can then be written as

[TABLE]

The inverse transformation coefficients $(T^{-1})^{lkm}_{pj}$ defined by Eq. 102 can be found similarly, yielding

[TABLE]