3.1 Pseudo-spectral in phase space

Our approach is perhaps best understood by analogy with pseudo-spectral Fourier methods. Consider a Fourier decomposition for the spatial coordinates in the gyrokinetic equation. By using field-line-following coordinates to avoid irregular domains (Beer, Cowley & Hammett Reference Beer, Cowley and Hammett1995), Fourier techniques can be readily applied; periodic (and otherwise appropriate) domains are guaranteed by the asymptotics of the problem (Abel et al. Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013). Such a Fourier decomposition has the advantage of allowing spectrally accurate evaluation of the derivatives because the spatial derivative has Fourier harmonics as its eigenfunctions [ $\text{d}/\text{d}x\exp (\text{i}kx)=\text{i}k\exp (\text{i}kx)$ ]. Conversely, achieving spectral accuracy in evaluating derivatives with a finite-difference, real-space representation on an $n$ -point grid requires the use of $n$ -point stencils, which is expensive. There are costs, of course, that reduce the advantages of a Fourier representation. Nonlinearity, in particular, necessitates the evaluation of convolutions in a strict Fourier basis: quadratic nonlinearities introduce spectral convolutions with $O(n^{2})$ terms. Pseudo-spectral algorithms reduce this operator count to $O(n\log n)$ by using fast discrete transforms to evaluate the nonlinear terms in real space, without the loss of spectral accuracy.

Our new approach takes this a step further by using pseudo-spectral methods not only in configuration space but also in velocity space. Instead of working with $f(\boldsymbol{R},\boldsymbol{v})$ , we Fourier transform in the spatial coordinate $\boldsymbol{R}$ and Laguerre–Hermite transform in the velocity coordinate $\boldsymbol{v}$ to develop spectral equations governing the evolution of $f(\boldsymbol{k},\boldsymbol{p})$ . Here, $\boldsymbol{k}$ is the usual Fourier wavenumber, and $\boldsymbol{p}=(\ell ,m)$ is the ‘wavenumber’ in the Laguerre–Hermite dual to velocity space. The advantages of the Fourier decomposition discussed in the previous paragraph apply to velocity space as well. Just as a Fourier decomposition in space has the advantage of spectrally accurate evaluation of spatial derivatives, our Laguerre–Hermite spectral velocity decomposition allows spectrally accurate evaluation of the velocity derivatives that appear in (2.1) and in the collision operator. Nonlinearities may be evaluated pseudo-spectrally relatively inexpensively in $(\boldsymbol{R},\boldsymbol{v})$ space with the use of efficient transforms and appropriately chosen grids. We expect this fully pseudo-spectral formulation of gyrokinetics to be a useful addition to the community toolbox.

There are disadvantages to our approach, too. In addition to derivatives in velocity space, there are multiplications by velocity in the gyrokinetic equation. The canonical example is the free streaming of particles along the magnetic field $({\sim}v_{\Vert }\,\text{d}/\text{d}z)$ , which introduces a term like $\text{i}k_{z}v_{\Vert }$ to the fluctuation equations which describe a general fluctuation $\unicode[STIX]{x1D6FF}f\,(\boldsymbol{k},\boldsymbol{v})$ . In $(\boldsymbol{R},v)$ space, a term like this can be evaluated with a finite-difference scheme (Kotschenreuther et al. Reference Kotschenreuther, Rewoldt and Tang1995; Jenko et al. Reference Jenko, Dorland, Kotschenreuther and Rogers2000; Candy & Waltz Reference Candy and Waltz2003). In $(\boldsymbol{k},\boldsymbol{v})$ space, this can be evaluated with an integrating factor. In our $(\boldsymbol{k},\boldsymbol{p})$ space formulation, this term couples the equations describing the evolution of the various Hermite polynomial moments of $\unicode[STIX]{x1D6FF}f\,(v_{\Vert })$ . This is the mathematical manifestation of the physics of phase mixing, and it introduces the need to close the finite set of Hermite moments used. At high enough resolution, collisions (even if weak) provide physical closure, but we anticipate the key benefit of using a variable spectral velocity representation will be in the low resolution limit. The model can be regarded as a generalized gyrofluid system, with the flexibility to increase the number of moments to achieve any desired level of accuracy. In the limit of a large number of moments, the model is equivalent to traditional grid-based gyrokinetic algorithms.

3.2 Laguerre and Hermite functions

We begin by defining expansion and projection functions in our basis. The Laguerre expansion and projection functions are given (respectively) by

(3.1a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D713}^{\ell }(\unicode[STIX]{x1D707}B)=(-1)^{\ell }\text{e}^{-\unicode[STIX]{x1D707}B}\text{L}_{\ell }(\unicode[STIX]{x1D707}B),\quad \unicode[STIX]{x1D713}_{\ell }(\unicode[STIX]{x1D707}B)=(-1)^{\ell }\text{L}_{\ell }(\unicode[STIX]{x1D707}B),\end{eqnarray}$$

where

(3.2) $$\begin{eqnarray}\text{L}_{\ell }(x)=\frac{\text{e}^{x}}{\ell !}\frac{\text{d}^{\ell }}{\text{d}x^{\ell }}x^{\ell }\text{e}^{-x}\end{eqnarray}$$

are the Laguerre polynomials. The Hermite expansion and projection functions are given (respectively) by

(3.3a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D719}^{m}(v_{\Vert })=\frac{\text{e}^{-v_{\Vert }^{2}/2}\text{He}_{m}(v_{\Vert })}{\sqrt{(2\unicode[STIX]{x03C0})^{3}m!}},\quad \unicode[STIX]{x1D719}_{m}(v_{\Vert })=\frac{\text{He}_{m}(v_{\Vert })}{\sqrt{m!}},\end{eqnarray}$$

where

(3.4) $$\begin{eqnarray}\text{He}_{m}(x)=(-1)^{m}\text{e}^{x^{2}/2}\frac{\text{d}^{m}}{\text{d}x^{m}}\text{e}^{-x^{2}/2}\end{eqnarray}$$

are the probabilists’ Hermite polynomials.Footnote ⁷

These definitions yield orthonormality conditions,

(3.5a,b ) $$\begin{eqnarray}\int _{0}^{\infty }\text{d}\unicode[STIX]{x1D707}\,B\unicode[STIX]{x1D713}^{k}(\unicode[STIX]{x1D707}B)\unicode[STIX]{x1D713}_{\ell }(\unicode[STIX]{x1D707}B)=\unicode[STIX]{x1D6FF}_{k\ell }\quad 2\unicode[STIX]{x03C0}\int _{-\infty }^{\infty }\text{d}v_{\Vert }\,\unicode[STIX]{x1D719}^{m}(v_{\Vert })\unicode[STIX]{x1D719}_{n}(v_{\Vert })=\unicode[STIX]{x1D6FF}_{mn},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}_{ij}$ is the Kronecker delta. The expansion functions satisfy recurrence relations given by

(3.6) $$\begin{eqnarray}\displaystyle & \displaystyle (\ell +1)\unicode[STIX]{x1D713}^{\ell +1}(\unicode[STIX]{x1D707}B)=(\unicode[STIX]{x1D707}B-2\ell -1)\unicode[STIX]{x1D713}^{\ell }(\unicode[STIX]{x1D707}B)-\ell \unicode[STIX]{x1D713}^{\ell -1}(\unicode[STIX]{x1D707}B), & \displaystyle\end{eqnarray}$$

(3.7) $$\begin{eqnarray}\displaystyle & \displaystyle \sqrt{m+1}\unicode[STIX]{x1D719}^{m+1}(v_{\Vert })=v_{\Vert }\unicode[STIX]{x1D719}^{m}(v_{\Vert })-\sqrt{m}\unicode[STIX]{x1D719}^{m-1}(v_{\Vert }). & \displaystyle\end{eqnarray}$$

We will also need the following derivative relations:

(3.8a,b ) $$\begin{eqnarray}B\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{\ell }}{\unicode[STIX]{x2202}B}=-(\ell +1)(\unicode[STIX]{x1D713}^{\ell +1}+\unicode[STIX]{x1D713}^{\ell }),\quad \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}^{m}}{\unicode[STIX]{x2202}v_{\Vert }}=-\sqrt{m+1}\,\unicode[STIX]{x1D719}^{m+1}.\end{eqnarray}$$

Recall that all velocities here are normalized to the species thermal velocity. In the following sections, we suppress the polynomial arguments for concision.

3.3 Laguerre–Hermite expansion

The Laguerre and Hermite expansion functions defined above form an orthogonal basis for gyrotropic $L^{2}(\mathbb{R}^{3},\text{e}^{-v_{\Vert }^{2}/2-\unicode[STIX]{x1D707}B}\,\text{d}^{3}v)$ , which is the Hilbert space of functions $f$ satisfying

(3.9) $$\begin{eqnarray}\int \text{d}^{3}v\,|\,f(v_{\Vert },\unicode[STIX]{x1D707})|^{2}\text{e}^{-v_{\Vert }^{2}/2-\unicode[STIX]{x1D707}B}<\infty .\end{eqnarray}$$

Thus as long as the perturbed distribution functions $g$ and $h$ satisfy this requirement (which should hold under the $\unicode[STIX]{x1D6FF}f$ ordering), we can expand them as

(3.10a,b ) $$\begin{eqnarray}g=\mathop{\sum }_{\ell =0}^{\infty }\mathop{\sum }_{m=0}^{\infty }\unicode[STIX]{x1D713}^{\ell }\unicode[STIX]{x1D719}^{m}G_{\ell ,m},\quad h=\mathop{\sum }_{\ell =0}^{\infty }\mathop{\sum }_{m=0}^{\infty }\unicode[STIX]{x1D713}^{\ell }\unicode[STIX]{x1D719}^{m}H_{\ell ,m}.\end{eqnarray}$$

The spectral amplitudes are defined by the projection of $g$ and $h$ onto the Laguerre–Hermite basis, given by

(3.11) $$\begin{eqnarray}\displaystyle & \displaystyle G_{\ell ,m}=2\unicode[STIX]{x03C0}\int _{-\infty }^{\infty }\text{d}v_{\Vert }\int _{0}^{\infty }\text{d}\unicode[STIX]{x1D707}\,B\unicode[STIX]{x1D713}_{\ell }\unicode[STIX]{x1D719}_{m}g, & \displaystyle\end{eqnarray}$$

(3.12) $$\begin{eqnarray}\displaystyle & \displaystyle H_{\ell ,m}=2\unicode[STIX]{x03C0}\int _{-\infty }^{\infty }\text{d}v_{\Vert }\int _{0}^{\infty }\text{d}\unicode[STIX]{x1D707}\,B\unicode[STIX]{x1D713}_{\ell }\unicode[STIX]{x1D719}_{m}h=G_{\ell ,m}+\frac{Z_{s}}{\unicode[STIX]{x1D70F}_{s}}{\mathcal{J}}_{\ell }\unicode[STIX]{x1D6F7}\unicode[STIX]{x1D6FF}_{m0}, & \displaystyle\end{eqnarray}$$

with ${\mathcal{J}}_{\ell }$ defined in (3.18) below. Note that $G_{\ell ,m}=H_{\ell ,m}$ for $m\neq 0$ . Physically, the $G_{\ell ,m}$ (and $H_{\ell ,m}$ ) are orthonormal fluid moments of $g$ (and $h$ ). Each $G_{\ell ,m}=G_{\ell ,m}(\boldsymbol{R},t)$ is a guiding centre function of space and time. We will frequently work with the perpendicular Fourier components of $G_{\ell ,m}(\boldsymbol{k}_{\bot },z,t)$ without changing notation when the context is clear.

We also note that there is a direct correspondence with the gyrofluid moments defined by Beer & Hammett (Reference Beer and Hammett1996),

(3.13) $$\begin{eqnarray}(G_{0,0},G_{0,1},\sqrt{2}G_{0,2},\sqrt{6}G_{0,3},G_{1,0},G_{1,1})=(n,u_{\Vert },T_{\Vert },q_{\Vert },T_{\bot },q_{\bot }),\end{eqnarray}$$

where each moment on the right appears exactly as defined by Beer, with the understanding that the normalizations have all been applied. Thus our choice of the Laguerre–Hermite basis allows us to extend the Beer model to arbitrary number of moments.

Working from (3.10), we can expand the following derivatives that appear in the gyrokinetic equation:

(3.14) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{\boldsymbol{b}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}h=\mathop{\sum }_{\ell ,m=0}^{\infty }\unicode[STIX]{x1D713}^{\ell }\unicode[STIX]{x1D719}^{m}\hat{\boldsymbol{b}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}H_{\ell ,m}-(\ell +1)\left(\unicode[STIX]{x1D713}^{\ell +1}+\unicode[STIX]{x1D713}^{\ell }\right)\unicode[STIX]{x1D719}^{m}H_{\ell ,m}\,\hat{\boldsymbol{b}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\ln B, & \displaystyle\end{eqnarray}$$

(3.15) $$\begin{eqnarray}\displaystyle & \displaystyle \hat{\boldsymbol{b}}\times \unicode[STIX]{x1D735}h=\mathop{\sum }_{\ell ,m=0}^{\infty }\unicode[STIX]{x1D713}^{\ell }\unicode[STIX]{x1D719}^{m}\,\hat{\boldsymbol{b}}\times \unicode[STIX]{x1D735}\,H_{\ell ,m}, & \displaystyle\end{eqnarray}$$

(3.16) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}v_{\Vert }}=-\mathop{\sum }_{\ell ,m=0}^{\infty }\sqrt{m+1}\,\unicode[STIX]{x1D713}^{\ell }\unicode[STIX]{x1D719}^{m+1}\,H_{\ell ,m}. & \displaystyle\end{eqnarray}$$

Electric fields are everywhere gyroaveraged in the gyrokinetic equation, with the gyroaveraging operation expressed in Fourier space as multiplication by the Bessel function $\text{J}_{0}$ . We can expand $\text{J}_{0}$ in terms of our Laguerre functions as

(3.17) $$\begin{eqnarray}\text{J}_{0}(\sqrt{2\unicode[STIX]{x1D707}Bb})=\mathop{\sum }_{\ell =0}^{\infty }\unicode[STIX]{x1D713}_{\ell }{\mathcal{J}}_{\ell }(b)\equiv \mathop{\sum }_{\ell =0}^{\infty }\unicode[STIX]{x1D713}_{\ell }\frac{b^{\ell }}{\ell !}\frac{\unicode[STIX]{x2202}^{\ell }}{\unicode[STIX]{x2202}b^{\ell }}\langle \text{J}_{0}\rangle ,\end{eqnarray}$$

where $b=b_{s}=k_{\bot }^{2}\unicode[STIX]{x1D70C}_{s}^{2}$ and ${\mathcal{J}}_{\ell }$ can be interpreted as the amplitude of the gyroaveraging operator in Fourier–Laguerre space, defined by

(3.18) $$\begin{eqnarray}{\mathcal{J}}_{\ell }\equiv \frac{b^{\ell }}{\ell !}\frac{\unicode[STIX]{x2202}^{\ell }}{\unicode[STIX]{x2202}b^{\ell }}\langle \text{J}_{0}\rangle .\end{eqnarray}$$

Here, $\langle \text{J}_{0}\rangle$ denotes the velocity integral of $\text{J}_{0}$ . In the infinite series limit, we can use the exact expression, $\langle \text{J}_{0}\rangle =\int \text{d}\unicode[STIX]{x1D707}\,B\text{J}_{0}(\sqrt{2\unicode[STIX]{x1D707}Bb})\text{e}^{-\unicode[STIX]{x1D707}B}=\text{e}^{-b/2}$ . This expression gives FLR terms that are consistent with the FLR approach of Brizard (Reference Brizard1992). In the finite series limit, however, with the sum truncated at some ${\mathcal{L}}$ , equation (3.17) is only $O(b^{{\mathcal{L}}})$ accurate. For applications with the maximum required value of $b\lesssim {\mathcal{L}}$ , this remains a very good approach. In appendix B we show that this results in accurate FLR terms in the dispersion relation for approximately $b\lesssim {\mathcal{L}}$ . We also consider an alternative choice in the finite series limit, $\langle \text{J}_{0}\rangle =\unicode[STIX]{x1D6E4}_{0}^{1/2}$ , which is consistent with the FLR approach of Dorland & Hammett (Reference Dorland and Hammett1993). We show in appendix B that this approach gives better asymptotic behaviour at large $b$ for the dispersion relation, for a range of values of ${\mathcal{L}}$ . Extending the ideas of Dorland & Hammett (Reference Dorland and Hammett1993), one can probably choose a Padé approximation for $\langle \text{J}_{0}\rangle$ accurate to order $b\sim {\mathcal{L}}$ , using the additional freedom to satisfy additional constraints. We have not explored this possibility exhaustively.

Finally, we note that the procedure used to derive the fluid equations that follow is agnostic to the specific definition of $\langle \text{J}_{0}\rangle$ , and only relies on the fact that we use (3.17) to expand $\text{J}_{0}$ in the Laguerre basis. This is also true of the results pertaining to free energy conservation in § 4.

3.4 Laguerre–Hermite fluid equations

We now derive the infinite set of coupled three-dimensional (3-D) fluid equations governing the evolution of the Laguerre–Hermite spectral amplitudes $G_{\ell ,m}$ . The procedure is straightforward: we work term by term in the gyrokinetic equation, first using the identities above to expand each term as an infinite Laguerre–Hermite sum, and then using the orthogonality relations to project the terms onto the Laguerre–Hermite basis. To illustrate this procedure, we show as an example how we project the magnetic drift term in the gyrokinetic equation in appendix C.

This effectively reduces (without approximation) the 5-D gyrokinetic equation to an infinite collection of coupled 3D fluid equations [(3.21), below]. Coupling arises from finite Larmor radius (FLR) effects, the $\boldsymbol{E}\times \boldsymbol{B}$ and toroidal drifts, collisions and flows along field lines, including flows associated with magnetic trapping.

The projection of the nonlinear term in (2.1) involves the convective derivative,

(3.19) $$\begin{eqnarray}\displaystyle \frac{\text{d}G_{\ell ,m}}{\text{d}t} & \equiv & \displaystyle \frac{\unicode[STIX]{x2202}G_{\ell ,m}}{\unicode[STIX]{x2202}t}+\mathop{\sum }_{k=0}^{\infty }\mathop{\sum }_{n=|k-\ell |}^{k+\ell }\unicode[STIX]{x1D6FC}_{k\ell n}\,({\mathcal{J}}_{n}\boldsymbol{v}_{E})\boldsymbol{\cdot }\unicode[STIX]{x1D735}G_{k,m}\nonumber\\ \displaystyle & = & \displaystyle \frac{\unicode[STIX]{x2202}G_{\ell ,m}}{\unicode[STIX]{x2202}t}+\mathop{\sum }_{k=0}^{\infty }~\mathop{\sum }_{n=|k-\ell |}^{k+\ell }\unicode[STIX]{x1D6FC}_{k\ell n}\,({\mathcal{J}}_{n}\boldsymbol{v}_{E})\boldsymbol{\cdot }\unicode[STIX]{x1D735}H_{k,m},\end{eqnarray}$$

where the equivalence of the expressions with $G$ and $H$ follows from the fact that $\langle \boldsymbol{v}_{E}\rangle \boldsymbol{\cdot }\unicode[STIX]{x1D735}h=\langle \boldsymbol{v}_{E}\rangle \boldsymbol{\cdot }\unicode[STIX]{x1D735}g$ since $\langle \boldsymbol{v}_{E}\rangle \boldsymbol{\cdot }\unicode[STIX]{x1D735}\langle \unicode[STIX]{x1D6F7}\rangle =0$ . The convolution in (3.19) arises from FLR-induced coupling, and in particular is related to for nonlinear, FLR-induced phase mixing (Dorland & Hammett Reference Dorland and Hammett1993). The convolution coefficients are given by

(3.20) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}_{k\ell n}=\int _{0}^{\infty }\text{d}\unicode[STIX]{x1D707}\,B\unicode[STIX]{x1D713}_{k}\unicode[STIX]{x1D713}^{\ell }\unicode[STIX]{x1D713}_{n}=\mathop{\sum }_{j}\frac{(k+\ell -j)!2^{2j-k-\ell +n}}{(k-j)!(\ell -j)!(2j-k-\ell +n)!(k+\ell -n-j)!},\end{eqnarray}$$

as first calculated by Watson (Reference Watson1938), where the summation limits are set by the requirement that all factorials have non-negative arguments. This requirement is consistent with the bandwidth limits in the sum over $n$ in (3.19).

Note that instead of evaluating this term as a convolution, one could use a pseudo-spectral approach by transforming to $(v_{\Vert },\unicode[STIX]{x1D707}B)$ coordinates and evaluating the nonlinearity as a pointwise multiplication. This pseudo-spectral alternative, which is described in appendix D, requires efficient implementations of the transforms described by (3.10) and (3.11). These transforms can be evaluated via matrix multiplication, which is efficient unless one uses quite a large number of Laguerre moments.Footnote ⁸

In terms of Laguerre–Hermite basis functions, equation (2.1) is

(3.21) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\text{d}G_{\ell ,m}}{\text{d}t}+v_{ts}\unicode[STIX]{x1D735}_{\Vert }\left(\sqrt{m+1}H_{\ell ,m+1}+\sqrt{m}H_{\ell ,m-1}\right)\nonumber\\ \displaystyle & & \displaystyle \qquad \quad +\,v_{ts} [-(\ell +1)\sqrt{m+1}H_{\ell ,m+1}-\ell \sqrt{m+1}H_{\ell -1,m+1}\nonumber\\ \displaystyle & & \displaystyle \qquad \quad +\,\ell \,\sqrt{m}\,H_{\ell ,m-1}+(\ell +1)\sqrt{m}\,H_{\ell +1,m-1} ]\unicode[STIX]{x1D735}_{\Vert }\ln B\nonumber\\ \displaystyle & & \displaystyle \qquad \quad +\,\text{i}\unicode[STIX]{x1D714}_{ds} [\sqrt{(m+1)(m+2)}\,H_{\ell ,m+2}+(\ell +1)\,H_{\ell +1,m}\nonumber\\ \displaystyle & & \displaystyle \qquad \quad +\,2\,(\ell +m+1)\,H_{\ell ,m}+\sqrt{m(m-1)}\,H_{\ell ,m-2}+\ell \,H_{\ell -1,m} ]\nonumber\\ \displaystyle & & \displaystyle \qquad =D_{\ell ,m}+C_{\ell ,m}.\end{eqnarray}$$

Parallel convection, including bounce motion induced by magnetic trapping in the equilibrium magnetic field, is described by the terms proportional to $\unicode[STIX]{x1D735}_{\Vert }\equiv \hat{\boldsymbol{b}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}$ . Toroidicity gives rise to the terms proportional to $\text{i}\unicode[STIX]{x1D714}_{ds}=\text{i}\unicode[STIX]{x1D714}_{d}(\unicode[STIX]{x1D70F}_{s}/Z_{s})$ , where $\text{i}\unicode[STIX]{x1D714}_{d}\equiv (1/B^{2})\hat{\boldsymbol{b}}\times \unicode[STIX]{x1D735}B\boldsymbol{\cdot }\unicode[STIX]{x1D735}$ . Drive terms from equilibrium gradients, denoted by $D_{\ell ,m}$ , are given by

(3.22) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}l@{}}D_{\ell ,m=0}=\text{i}\unicode[STIX]{x1D714}_{\ast }\left[{\displaystyle \frac{1}{L_{ns}}}{\mathcal{J}}_{\ell }\unicode[STIX]{x1D6F7}+{\displaystyle \frac{1}{L_{Ts}}}[\ell {\mathcal{J}}_{\ell -1}\unicode[STIX]{x1D6F7}+2\ell {\mathcal{J}}_{\ell }\unicode[STIX]{x1D6F7}+(\ell +1){\mathcal{J}}_{\ell +1}\unicode[STIX]{x1D6F7}]\right]\\ D_{\ell ,m=2}={\displaystyle \frac{1}{\sqrt{2}}}\text{i}\unicode[STIX]{x1D714}_{\ast }\,{\displaystyle \frac{1}{L_{Ts}}}\,{\mathcal{J}}_{\ell }\unicode[STIX]{x1D6F7}\\ D_{\ell ,m}=0\quad \text{otherwise},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where $\text{i}\unicode[STIX]{x1D714}_{\ast }\equiv -\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F9}\boldsymbol{\cdot }\hat{\boldsymbol{b}}\times \unicode[STIX]{x1D735}$ , and $L_{ns}$ and $L_{Ts}$ are the normalized density and temperature gradient scale lengths, respectively. The collision terms, denoted by $C_{\ell ,m}$ , are presented in § 3.5.

To complete the Laguerre–Hermite equation set, we must use (2.3) to find the potential $\unicode[STIX]{x1D6F7}$ , given $h$ . The left-hand side of (2.3) involves the non-Boltzmann part of the particle space density $\bar{n}=\bar{n}(\boldsymbol{r})$ , which is given by

(3.23) $$\begin{eqnarray}\bar{n}=\int \text{d}^{3}\boldsymbol{v}\langle h\rangle _{\boldsymbol{ r}}=\int \text{d}^{3}\boldsymbol{v}\,\text{J}_{0}h=\mathop{\sum }_{\ell =0}^{\infty }{\mathcal{J}}_{\ell }H_{\ell ,0},\end{eqnarray}$$

where we have used $\langle h\rangle _{\boldsymbol{r}}$ and $\text{J}_{0}h$ interchangeably for the gyroaverage of $h$ , with the understanding that in the latter case we have taken the Fourier transform so that $h=h(\boldsymbol{k}_{\bot })$ . For hydrogenic plasma with Boltzmann electrons, the quasineutrality equation reduces to (2.4), which projects to

(3.24) $$\begin{eqnarray}\mathop{\sum }_{\ell =0}^{\infty }{\mathcal{J}}_{\ell }H_{\ell ,0}=\unicode[STIX]{x1D70F}_{e}^{-1}[\unicode[STIX]{x1D6F7}-\langle \langle \unicode[STIX]{x1D6F7}\rangle \rangle ]+\unicode[STIX]{x1D6F7}.\end{eqnarray}$$

Equivalently, quasineutrality can be expressed in terms of $G_{\ell ,0}$ , as

(3.25) $$\begin{eqnarray}\mathop{\sum }_{\ell =0}^{\infty }{\mathcal{J}}_{\ell }G_{\ell ,0}=\unicode[STIX]{x1D70F}^{-1}[\unicode[STIX]{x1D6F7}-\langle \langle \unicode[STIX]{x1D6F7}\rangle \rangle ]+\left[1-\mathop{\sum }_{\ell =0}^{\infty }{\mathcal{J}}_{\ell }^{2}\right]\unicode[STIX]{x1D6F7}.\end{eqnarray}$$

When we truncate the system and use only ${\mathcal{L}}$ Laguerre moments, we truncate the sum of the left-hand side of (3.24) at ${\mathcal{L}}-1$ . In (3.25), this is equivalent to truncating the sums on both the left- and right-hand side of (3.25) at ${\mathcal{L}}-1$ . This maintains consistency in the FLR treatment of the left and right-hand side, which is most clear when we express quasineutrality in terms of $H$ . Further, note that $\sum {\mathcal{J}}_{\ell }^{2}\approx \unicode[STIX]{x1D6E4}_{0}=\int \text{d}^{3}\boldsymbol{v}\,\text{J}_{0}^{2}F_{M}$ , where $\unicode[STIX]{x1D6E4}_{0}=\unicode[STIX]{x1D6E4}_{0}(b)=\text{I}_{0}(b)\text{e}^{-b}$ , and $\text{I}_{0}(b)=\text{J}_{0}(\text{i}b)$ is the modified Bessel function.

Finally, we stress that the definition of ${\mathcal{J}}_{\ell }$  (3.18) used in (3.25) for a particular species should be the same as is used in the evolution equations (3.21) for that species. This ensures free energy conservation, as noted by Scott (Reference Scott2005). The Beer and Dorland gyrofluid models do not conserve free energy at all orders of $k_{\bot }\unicode[STIX]{x1D70C}$ because the FLR expressions in the quasineutrality equation are only consistent with the FLR terms in the fluid equations at long perpendicular wavelengths.

3.5 Collision operator

To model collisional physics, we generalize the single-species Dougherty collision operator (Dougherty Reference Dougherty1964), which is itself a generalization of the Lenard–Bernstein collision operator (Lenard & Bernstein Reference Lenard and Bernstein1958). Though it is a simplified, model collision operator, the Dougherty operator has several appealing properties. It describes pitch-angle scattering, energy diffusion and slowing down. It vanishes on (and only on) a Maxwellian, and it conserves number, momentum and energy. Thus it satisfies the appropriate H-theorem, and drives the plasma to thermal equilibrium in the long-time limit. These properties were recognized and emphasized by Anderson & O’Neil (Reference Anderson and O’Neil2007). Below, we gyroaverage the Dougherty operator and express the result in normalized $(v_{\Vert },\unicode[STIX]{x1D707})$ coordinates. The result is a gyrokinetic operator that can be compactly expressed in terms of Laguerre–Hermite polynomials, which are eigenfunctions of the differential part of the operator, while preserving all of the above properties.

3.5.1 Definition

The starting point is

(3.26) $$\begin{eqnarray}C(f)=\unicode[STIX]{x1D708}\,\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{v}}\boldsymbol{\cdot }\left[\unicode[STIX]{x1D63F}[f](\boldsymbol{v})\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}\boldsymbol{v}}+P[f](\boldsymbol{v})f\right],\end{eqnarray}$$

where $\unicode[STIX]{x1D63F}$ is a velocity-space diffusion tensor that is a functional of $f$ , and $P$ is a field–particle operator, also a functional of $f$ in the general case. The collision frequency $\unicode[STIX]{x1D708}$ is constant, i.e. independent of velocity. Dougherty approximated $\unicode[STIX]{x1D63F}$ by $T/M$ , where $M$ is the particle mass and $T$ is the temperature, allowing $T=T_{0}+\unicode[STIX]{x1D6FF}T$ . For $P[f]$ , he used $\boldsymbol{v}-\boldsymbol{u}$ , with the flow $\boldsymbol{u}=\boldsymbol{u}_{0}+\unicode[STIX]{x1D6FF}\boldsymbol{u}$ . Below, we make the further approximation that $\boldsymbol{u}_{0}=0$ , though this could be relaxed in the future if desired.

With these assumptions, we can then linearize to obtain

(3.27) $$\begin{eqnarray}C(\unicode[STIX]{x1D6FF}f)=\unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{v}}\boldsymbol{\cdot }\left[\frac{T_{0}}{M}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FF}f}{\unicode[STIX]{x2202}\boldsymbol{v}}+\frac{\unicode[STIX]{x1D6FF}T}{M}\frac{\unicode[STIX]{x2202}F_{M}}{\unicode[STIX]{x2202}\boldsymbol{v}}+\boldsymbol{v}\unicode[STIX]{x1D6FF}f-\unicode[STIX]{x1D6FF}\boldsymbol{u}F_{M}\right].\end{eqnarray}$$

By construction, $C(F_{M})=0$ . Number is automatically conserved because $C$ has the form of a velocity-space divergence. Momentum and energy are conserved by the field–particle terms, which are $(\propto \unicode[STIX]{x1D6FF}T,\unicode[STIX]{x1D6FF}\boldsymbol{u})$ . The forms of the restoring terms are explicitly constructed to preserve the properties listed above.

In gyrokinetic variables, recall from above that the perturbed distribution function can be written as $\unicode[STIX]{x1D6FF}f=h-Z_{s}/\unicode[STIX]{x1D70F}_{s}\unicode[STIX]{x1D6F7}F_{M}$ . The velocity-space structure of the Boltzmann component $(\propto \unicode[STIX]{x1D6F7})$ is Maxwellian, and thus $C(\unicode[STIX]{x1D6FF}f)=C(h)$ . Note that in general $C(h)\neq C(g)$ due to the velocity dependence of the gyroaveraging operator.

We gyroaverage this collision operator in the standard way (Catto & Tsang Reference Catto and Tsang1977). Collisions occur at fixed position $\boldsymbol{r}$ rather than at fixed gyrocentre, $\boldsymbol{R}$ . Thus, the collision operator must be evaluated at fixed $\boldsymbol{r}$ . Only at the end does one transform back to guiding centre position $\boldsymbol{R}$ , as required by the fact that (2.1) describes the evolution of $g(\boldsymbol{R})$ .

Keeping track of these dependences is slightly awkward, because $h=h(\boldsymbol{R})$ . Following Abel et al. (Reference Abel, Barnes, Cowley, Dorland and Schekochihin2008) we express the distribution function $h$ in $\boldsymbol{k}$ space to facilitate the required manipulations:

(3.28) $$\begin{eqnarray}\displaystyle \langle C_{\boldsymbol{r}}(h)\rangle _{\boldsymbol{R}} & = & \displaystyle \left\langle C\left(\mathop{\sum }_{\boldsymbol{k}}\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{R}}h_{k}\right)\right\rangle _{\boldsymbol{R}}=\mathop{\sum }_{\boldsymbol{k}}\langle \text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{r}}C(\text{e}^{-\text{i}\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D746}}h_{k})\rangle _{\boldsymbol{R}}\nonumber\\ \displaystyle & = & \displaystyle \mathop{\sum }_{\boldsymbol{k}}\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{R}}\langle \text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D746}}C(\text{e}^{-\text{i}\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D746}}h_{k})\rangle _{\boldsymbol{R}}\equiv C(h).\end{eqnarray}$$

In this form, it is less difficult to perform the gyroaverages. The phase factors combine to describe classical diffusion $\propto -\unicode[STIX]{x1D708}(k_{\bot }^{2}\unicode[STIX]{x1D70C}^{2})h$ . (Number, momentum and energy are conserved before and after gyroaveraging, as discussed below and in appendix E.) Setting aside the field–particle (conservation) terms, in cylindrical $(v_{\Vert },\unicode[STIX]{x1D707})$ coordinates the result of gyroaveraging the velocity-space derivative (or test particle diffusion) part of the collision operator is

(3.29) $$\begin{eqnarray}C(h)=\unicode[STIX]{x1D708}\left[\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}v_{\Vert }}\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}v_{\Vert }}+v_{\Vert }\right)+2\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\left(\frac{\unicode[STIX]{x1D707}}{B}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}+\unicode[STIX]{x1D707}\right)-k_{\bot }^{2}\unicode[STIX]{x1D70C}^{2}\right]h.\end{eqnarray}$$

Pitch-angle scattering is a particularly important dynamical process in many experimental scenarios. Our collision operator, even as it is expressed in cylindrical coordinates (3.29), describes pitch-angle scattering. Parenthetically, note that one can make pitch-angle and energy scattering manifest by changing coordinates to $(v,\unicode[STIX]{x1D709})$ , with $\unicode[STIX]{x1D709}=v_{\Vert }/v$ :

(3.30) $$\begin{eqnarray}C(h)=\unicode[STIX]{x1D708}\left\{\frac{1}{v^{2}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}\left[(1-\unicode[STIX]{x1D709}^{2})\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}\right]+\frac{1}{v^{2}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}v}\left[v^{2}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}v}+v^{3}\right]-k_{\bot }^{2}\unicode[STIX]{x1D70C}^{2}\right\}h.\end{eqnarray}$$

The field–particle terms are calculatedFootnote ⁹ by integrals of $h$ at fixed $\boldsymbol{r}$ . The full expression required for the momentum-conserving terms in the collision operator is

(3.31) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D6FF}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{v}}{v_{ts}^{2}}=\bar{\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{v}=[\bar{u}_{\Vert }v_{\Vert }\text{J}_{0}+\bar{u}_{\bot }v_{\bot }J_{1}].\end{eqnarray}$$

Here, the perturbed parallel flow of particles (not guiding centres) is

(3.32) $$\begin{eqnarray}\bar{u}_{\Vert }\equiv \int \text{d}^{3}\boldsymbol{v}\,\text{J}_{0}v_{\Vert }h.\end{eqnarray}$$

The standard factor of $\text{J}_{0}$ in the integrand arises from integrating at fixed $\boldsymbol{r}$ , as dictated by the spatial locality of collisions, which requires a gyroaverage of $h(\boldsymbol{R})$ . The integral for perpendicular momentum conservation is also standard, involving

(3.33) $$\begin{eqnarray}\bar{u}_{\bot }\equiv \int \text{d}^{3}\boldsymbol{v}\,J_{1}v_{\bot }h,\end{eqnarray}$$

where $J_{1}(a)=-\text{d}/\text{d}a\,\text{J}_{0}(a)$ . Without the field–particle terms, the collision operator also fails to conserve energy. The total perturbed energy consists of both parallel and perpendicular energy,

(3.34) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D6FF}T}{mv_{ts}^{2}}=\bar{T}=\frac{1}{3}(\bar{T}_{\Vert }+2\bar{T}_{\bot })=\frac{1}{3}\int \text{d}^{3}\boldsymbol{v}\,[(v_{\Vert }^{2}-1)+2(\unicode[STIX]{x1D707}B-1)]\text{J}_{0}h,\end{eqnarray}$$

where again, the Bessel function ultimately arises from the locality of collisions. Putting it all together, the collision operator is

(3.35) $$\begin{eqnarray}\displaystyle C(h) & = & \displaystyle \unicode[STIX]{x1D708}\left\{\left[\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}v_{\Vert }}\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}v_{\Vert }}+v_{\Vert }\right)+2\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\left(\frac{\unicode[STIX]{x1D707}}{B}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}+\unicode[STIX]{x1D707}\right)-k_{\bot }^{2}\unicode[STIX]{x1D70C}^{2}\right]h\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,[\bar{T}[(v_{\Vert }^{2}-1)+2(\unicode[STIX]{x1D707}B-1)]\text{J}_{0}+\bar{\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{v}]F_{M}\right\}.\end{eqnarray}$$

The dimensionless collision frequency $\unicode[STIX]{x1D708}$ is defined by (A 3).

This collision operator is a good physical model of like-particle collisions, which are important to gyrokinetic dynamics. It captures the physics of the collision operator presented in Abel et al. (Reference Abel, Barnes, Cowley, Dorland and Schekochihin2008), except that our collision frequency does not have velocity dependence.

3.5.2 Laguerre–Hermite projection

This collision operator has an additional attractive feature: the Laguerre and Hermite polynomials are eigenfunctions of the differential part of the operator. It can therefore be evaluated efficiently. Applying the parallel velocity components of the cylindrical velocity-space Laplacian in (3.35) to the Hermite basis functions gives

(3.36) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}v_{\Vert }}\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}v_{\Vert }}+v_{\Vert }\right)\unicode[STIX]{x1D719}^{m}=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}v_{\Vert }}({\unicode[STIX]{x1D719}^{\prime }}^{m}+\sqrt{m+1}\unicode[STIX]{x1D719}^{m+1}+\sqrt{m}\unicode[STIX]{x1D719}^{m-1})=-m\unicode[STIX]{x1D719}^{m}.\end{eqnarray}$$

This means that the Hermite basis functions are eigenfunctions of the parallel velocity components of the cylindrical velocity-space Laplacian, with eigenvalue $-m$ . Similarly, for the perpendicular components we have

(3.37) $$\begin{eqnarray}2\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\left(\frac{\unicode[STIX]{x1D707}}{B}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}+\unicode[STIX]{x1D707}\right)\unicode[STIX]{x1D713}^{\ell }=\frac{2\ell }{B}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}(\unicode[STIX]{x1D713}^{\ell }+\unicode[STIX]{x1D713}^{\ell -1})=-2\ell \unicode[STIX]{x1D713}^{\ell },\end{eqnarray}$$

so the Laguerre basis functions are eigenfunctions of the perpendicular components of the Laplacian. The result is that the projection of the velocity-space Laplacian operating on $h$ is sparse in our basis. Finally, the conservation terms only project onto the $m=0,1$ , and 2 Hermite moments; we leave the details to appendix E.1.

Thus the Laguerre–Hermite projection of our collision operator is

(3.38) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}l@{}}C_{\ell ,0}=-\unicode[STIX]{x1D708}\,(b+2\ell +0)\,H_{\ell ,0}+\unicode[STIX]{x1D708}\left(\sqrt{b}({\mathcal{J}}_{\ell }+{\mathcal{J}}_{\ell -1})\bar{u}_{\bot }\right.\\ \qquad +\,\left.2[\ell {\mathcal{J}}_{\ell -1}+2\ell {\mathcal{J}}_{\ell }+(\ell +1){\mathcal{J}}_{\ell +1}]\bar{T}\right),\\ C_{\ell ,1}=-\unicode[STIX]{x1D708}\,(b+2\ell +1)H_{\ell ,1}+\unicode[STIX]{x1D708}\,{\mathcal{J}}_{\ell }\bar{u}_{\Vert },\\ C_{\ell ,2}=-\unicode[STIX]{x1D708}\,(b+2\ell +2)\,H_{\ell ,2}+\unicode[STIX]{x1D708}\,\sqrt{2}\,{\mathcal{J}}_{\ell }\bar{T},\\ C_{\ell ,m}=-\unicode[STIX]{x1D708}(b+2\ell +m)H_{\ell ,m},\quad (m>2)\end{array}\right\} & & \displaystyle\end{eqnarray}$$

with the projections of (3.32)–(3.34) given by

(3.39) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{u}_{\Vert }=\int \text{d}^{3}\boldsymbol{v}\,\text{J}_{0}v_{\Vert }h=\mathop{\sum }_{\ell =0}^{\infty }{\mathcal{J}}_{\ell }H_{\ell ,1}, & \displaystyle\end{eqnarray}$$

(3.40) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{u}_{\bot }=\int \text{d}^{3}\boldsymbol{v}\,J_{1}v_{\bot }h=\sqrt{b_{s}}\mathop{\sum }_{\ell =0}^{\infty }({\mathcal{J}}_{\ell }+{\mathcal{J}}_{\ell -1})H_{\ell ,0}, & \displaystyle\end{eqnarray}$$

(3.41) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{T}_{\Vert }=\int \text{d}^{3}\boldsymbol{v}\,\text{J}_{0}(v_{\Vert }^{2}-1)h=\sqrt{2}\mathop{\sum }_{\ell =0}^{\infty }{\mathcal{J}}_{\ell }H_{\ell ,2}, & \displaystyle\end{eqnarray}$$

(3.42) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{T}_{\bot }=\int \text{d}^{3}\boldsymbol{v}\,\text{J}_{0}(\unicode[STIX]{x1D707}B-1)h=\mathop{\sum }_{\ell =0}^{\infty }[\ell {\mathcal{J}}_{\ell -1}+2\ell {\mathcal{J}}_{\ell }+(\ell +1){\mathcal{J}}_{\ell +1}]H_{\ell ,0}. & \displaystyle\end{eqnarray}$$

From (3.41)–(3.42) one can find that $\bar{T}=(\bar{T}_{\Vert }+2\bar{T}_{\bot })/3$ in general. However, if the perpendicular temperature fluctuation is not being evolved (i.e. if ${\mathcal{L}}=1$ , a severe and unrecommended truncation), one should use $\bar{T}=\bar{T}_{\Vert }$ . With this change of definition, our operator reduces to the energy-conserving form of the Kirkwood operator.

3.6 Summary and closure considerations

We now have a full set of Laguerre–Hermite spectral equations, given by (3.21)–(3.24) with collision terms given by (3.38). Solving this gyrokinetic system with many $(\ell ,m)$ moments is rigorously equivalent to solving the gyrokinetic system of (2.1)–(2.4) with high resolution in $(v_{\Vert },\unicode[STIX]{x1D707})$ coordinates. There are natural advantages to using each representation. We expect that pseudo-spectral algorithms, which have access to both representations, will have significant advantages over purely spectral and purely $v$ -space algorithms, and also over Lagrangian (such as particle-in-cell) and semi-Lagrangian formulations, especially with realistic values of collisionality. This is because the Laguerre–Hermite formulation expresses critical conservation laws (density, momentum and energy) with the first few moments, while the higher moments are damped progressively more strongly, since $C\sim -\unicode[STIX]{x1D708}\,(b+2\ell +m)$ . Our collision operator reflects this prioritization, as it expresses the key conservation laws at long wavelength, but sharply attenuates higher moments at short wavelengths, as a result of classical diffusion, pitch-angle scattering, energy diffusion and slowing down. This short wavelength attenuation is qualitatively correct.

Of course, in practice one cannot solve an infinite set of Laguerre–Hermite evolution equations. In the limit of high Laguerre–Hermite resolution, the system can be simply truncated at some large ${\mathcal{L}}$ and ${\mathcal{M}}$ , and closed by setting $G_{\ell ,m}=H_{\ell ,m}=0$ , for all perturbations with either $\ell \geqslant {\mathcal{L}}$ or $m\geqslant {\mathcal{M}}$ that appear in the equations for the resolved moments. This also implies ${\mathcal{J}}_{\ell }=0$ in these terms for $\ell \geqslant {\mathcal{L}}$ when $m=0$ . We will term this closure approach ‘closure by truncation’. In this case, the unresolved moments correspond to comparable limitations on any discrete $v$ -space representation of $g(v)$ , where fine scales in $g(v)$ above the grid resolution cannot be resolved. The collision operator regulates these fine velocity scales by smoothing the distribution function. Our collision operator fulfils this purpose by acting increasingly strongly on higher Laguerre and Hermite moments, limiting their amplitude. Thus for a given collisionality, there is a physical cutoff at some ${\mathcal{L}}$ and ${\mathcal{M}}$ beyond which fine scales in velocity space are completely wiped out by collisions, which justifies truncation of the moment series. This is the simplest high-resolution closure, but not the only option. One can also obtain an asymptotically correct collisional closure by assuming the collision term becomes dominant in the unresolved moment equations (Loureiro et al. Reference Loureiro, Dorland, Fazendeiro, Kanekar, Mallet, Vilelas and Zocco2016).

In the limit of low Laguerre–Hermite resolution, the closure situation is more complicated. Unresolved moments are not expected to be negligible at collisionalities of interest, so closure by truncation at low resolution will generally give poor results. One possible approach is to follow the gyrofluid closure approach pioneered by Hammett & Perkins (Reference Hammett and Perkins1990), which was used and extended in the gyrofluid models of Dorland & Hammett (Reference Dorland and Hammett1993), Beer & Hammett (Reference Beer and Hammett1996) and Snyder & Hammett (Reference Snyder and Hammett2001). Further, Smith (Reference Smith1997) extended the Hammett–Perkins approach to an arbitrary number of moments in the slab limit. Thus in the slab limit, one could use Smith’s scheme to generate closures for the Hermite moments. No linear closure is needed for the Laguerre moments in the slab limit, so in this limit the linear closure problem is solved.Footnote ¹⁰

Extending Smith’s generalized closure approach to toroidal geometry is particularly challenging due to the presence of branch cuts in the kinetic dispersion relation. In particular, closures for the Hermite moments are complicated by the fact that phase mixing arises from both parallel Landau damping and the curvature drift resonance. Thus we leave the task of developing generalized closures in toroidal geometry for later work. We do however discuss closures that correspond to Beer’s toroidal gyrofluid equations in appendix F. This allows us to reproduce Beer’s results using our generalized Laguerre–Hermite model.

4.1 Summary of free energy issues

We have shown that the Laguerre–Hermite projection produces a set of coupled equations of moments, each of which has three spatial dimensions and evolves in time. If one truncates the collisionless system of equations [by simply neglecting all terms in the Laguerre–Hermite series with $\ell \geqslant {\mathcal{L}}$ or $m\geqslant {\mathcal{M}}$ (provided ${\mathcal{L}}>1$ and ${\mathcal{M}}>2$ )], and considers the case for which there are no driving gradients, these equations conserve the truncated gyrokinetic free energy given in (4.3). The orthogonality of the Laguerre and Hermite basis functions gives the structure that is missing from the Beer gyrofluid equations, which did not conserve a free energy functional (see appendix G).

In this truncated system, if one then includes collisions, the truncated gyrokinetic free energy is explicitly damped, even as the truncated expressions for number [(3.23), except with the $\ell$ sum running from $0$ to ${\mathcal{L}}-1$ ], parallel and perpendicular momenta [(3.39) and (3.40), with the same $\ell$ truncation] and energy [(3.34), using the projections in (3.42) and (3.41), with the same $\ell$ truncation] are conserved.

With non-zero, fixed-in-time gradients, the free energy density $W(t)$ is no longer formally bounded, even with finite collisions. This might happen, for example, if one had strong gradients, weak collisions, and only a few Laguerre–Hermite moments. In the the context of the gyrokinetic free energy, we have described the desired properties of Landau fluid closures that might be used to model the transfer of free energy density to unresolved scales in velocity space ( $\ell \geqslant {\mathcal{L}}$ or $m\geqslant {\mathcal{M}}$ ). In such cases, collisions, which are more important at small scales in velocity space, would ultimately provide a physical sink.

5.1 Local linear ITG

We first examine an ion temperature gradient (ITG) instability in the local limit, where $k_{\Vert }$ , $B$ and $\unicode[STIX]{x1D714}_{d}$ are treated as constants. We select a case used to validate the Beer $4+2$ gyrofluid model (Beer & Hammett Reference Beer and Hammett1996), which first appeared in Dong, Horton & Kim (Reference Dong, Horton and Kim1992). The relevant parameters (in our units) are $q=2$ , $R_{0}/a=5$ , $k_{\Vert }a=a/qR_{0}$ , $a/L_{n}=1$ and $a/L_{T}=1.5,2$ , and 3. We also set $\unicode[STIX]{x1D708}=0.01$ . Figure 1 shows linear growth rates over a range of $k_{y}\unicode[STIX]{x1D70C}_{i}$ for two choices of Laguerre–Hermite velocity resolution. For comparison, we also show the results from the Beer gyrofluid model. We see that the growth rates converge as the Laguerre–Hermite resolution is increased.

We calculate a normalized integrated velocity-space spectrum,

(5.1) $$\begin{eqnarray}P(\ell ,m)=\frac{\displaystyle \int \text{d}^{3}\boldsymbol{R}\,|G_{\ell ,m}|^{2}}{\displaystyle \int \text{d}^{3}\boldsymbol{R}\,|G_{0,0}|^{2}},\end{eqnarray}$$

to examine the structure of the distribution function in the Laguerre–Hermite basis. $P$ can also be interpreted as a measure of the free energy in each moment. Figure 2 shows this spectrum for the higher resolution case in the above $a/L_{T}=2$ calculation. The amplitude $P$ is highest at small $\ell$ and $m$ , which is expected since free energy is injected by the gradients at these large scales. The upper right quadrant of the plot, where both $\ell$ and $m$ are large, is more than six orders of magnitude smaller than the large amplitudes at small $\ell$ and $m$ . For best contrast we have truncated the colour scale to seven orders of magnitude, but note that the amplitude in the top right corner, corresponding to the finest resolved scale in the velocity space, is $P(15,31)\sim 10^{-10}$ . This is a good indication that this calculation is well resolved in velocity space. Further, we see that $P(7,15)\sim 10^{-5}$ . This is the finest resolved scale in the lower resolution case, which has a spectrum (not shown) that looks qualitatively similar to the lower left quadrant of figure 2, and which also has $P(7,15)\sim 10^{-5}$ . The fact that the lower resolution case agrees with the higher resolution case in both growth rates and spectra suggests that the low resolution case has sufficient resolution.

5.2 Cyclone linear ITG

We also examine an ITG instability in the non-local limit, for which we use the Cyclone base case parameters, a widely benchmarked test case (Dimits et al. Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000). Figure 3 shows the linear growth rates over a range of $k_{y}\unicode[STIX]{x1D70C}_{i}$ for the same two choices of Laguerre–Hermite velocity resolution as above. For comparison, we also show the results from the Beer $4+2$ gyrofluid model, along with results from the gyrokinetic code GS2. GS2 solves the gyrokinetic equation using a polar velocity grid in energy and pitch angle coordinates, $\unicode[STIX]{x1D700}$ and $\unicode[STIX]{x1D706}=\unicode[STIX]{x1D707}/\unicode[STIX]{x1D700}$ . In the limit of large velocity resolution, our Laguerre–Hermite approach should agree with the grid-based approach of GS2. From the figure we see that the two approaches do indeed agree in the higher resolution case. We also see that the convergence in resolution is slower here in the non-local limit than in the local limit above, as the lower resolution case gives growth rates too large for higher $k_{y}\unicode[STIX]{x1D70C}_{i}$ . This suggests that more velocity resolution is needed in the non-local limit than in the local limit.

Figure 1. Growth rates of an ITG instability in the local limit. The results from the Laguerre–Hermite formulation are shown at two choices of velocity resolution (blue and green), as well results from the Beer $4+2$ gyrofluid model (red). The Laguerre–Hermite results show good convergence in velocity resolution, especially at small $k_{y}\unicode[STIX]{x1D70C}_{i}$ .

Figure 2. The normalized integrated velocity-space spectrum, defined by (5.1), resulting from the higher resolution $a/L_{T}=2$ case in figure 1. The amplitudes decrease by several orders of magnitude moving from small to large $\ell$ and $m$ , indicating that the calculation is well resolved in velocity space.

Figure 3. Linear growth rates for the Cyclone base case (Dimits et al. Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000). The results from the Laguerre–Hermite formulation are shown at two choices of velocity resolution (blue and green), as well as results from the gyrokinetic code GS2 (black) and the Beer $4+2$ gyrofluid model (red).

Figure 4. The normalized integrated velocity-space spectrum resulting from the higher resolution case in figure 3. While this spectrum is still well resolved, the amplitudes decrease more gradually than in the local limit spectrum from figure 2.

This is confirmed by the velocity spectrum for this case, shown in figure 4. Whereas in the local limit a significant portion of the velocity space had amplitudes more than six orders of magnitude smaller than the largest amplitudes, in the non-local limit the amplitudes decrease much more gradually. This indicates that the distribution function has more structure in velocity space in this case. This also shows why the lower resolution case failed to produce accurate growth rates: the amplitudes in the lower left quadrant of figure 4 do not decrease by as much as in the local case, with $P(7,15)\sim 10^{-4}$ .

5.3 Rosenbluth–Hinton zonal flow residual

As a final test, we examine zonal flow dynamics in our Laguerre–Hermite formulation. Zonal flows are nonlinearly driven and nonlinearly damped (Rogers, Dorland & Kotschenreuther Reference Rogers, Dorland and Kotschenreuther2000) toroidally and poloidally symmetric sheared $\boldsymbol{E}\times \boldsymbol{B}$ flows that have been shown to play a key role in determining the turbulence saturation level (Hammett et al. Reference Hammett, Beer, Dorland, Cowley and Smith1993; Waltz, Kerbel & Milovich Reference Waltz, Kerbel and Milovich1994). Rosenbluth & Hinton (Reference Rosenbluth and Hinton1998) first showed that zonal flows are not linearly damped by collisionless processes, showing that in a simplified equilibrium model the residual flow is given by

(5.2) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D6F7}(t=\infty )}{\unicode[STIX]{x1D6F7}(t=0)}=\frac{1}{1+1.6q^{2}/\sqrt{\unicode[STIX]{x1D716}}},\end{eqnarray}$$

where $q$ is the safety factor and $\unicode[STIX]{x1D716}=r/R_{0}$ is the inverse aspect ratio. The original Beer gyrofluid model did not accurately capture the residual because of the damping introduced by the closures. This caused discrepancies with gyrokinetic models in nonlinear simulations (Dimits et al. Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000). The gyrofluid closure was modified to be able to capture some of the zonal flow residual with limited success (Beer & Hammett Reference Beer, Hammett, Connor, Sindoni and Vaclavik1998).

Figure 5. Rosenbluth–Hinton residual flow for various values of $\unicode[STIX]{x1D716}=r/R_{0}$ . The residuals calculated by the Laguerre–Hermite model with ${\mathcal{L}}=16$ and ${\mathcal{M}}=32$ agree well with those calculated by GS2. We also show the expected theoretical result, equation (5.2).

Figure 6. The normalized integrated velocity-space spectrum for the $\unicode[STIX]{x1D716}=0.2$ case from figure 5. High amplitudes at large $m$ along the axis are the result of trapped particle dynamics that produces sharp features in $v_{\Vert }$ .

Figure 5 shows residuals for several values of $\unicode[STIX]{x1D716}$ , for the $k_{x}\unicode[STIX]{x1D70C}_{i}=0.01$ mode with $q=1.4$ , as calculated by our Laguerre–Hermite model with ${\mathcal{L}}=16$ and ${\mathcal{M}}=32$ . We also show the residual calculated by GS2 and the expected theoretical value from (5.2). Our result for the average residual agrees well with GS2, though we note that we observe larger oscillations and less damping of the potential in the Laguerre–Hermite model. We attribute this to resolution issues. The velocity spectrum, shown in figure 6 for the $\unicode[STIX]{x1D716}=0.2$ case, shows overall larger amplitudes compared to the spectra from the linear instability calculations. We also see that there are high amplitudes at large $m$ along the axis, indicating that sharp features are being generated in $g(v_{\Vert })$ . This is to be expected, since the zonal flow residual is the result of complicated kinetic effects involving trapped particle dynamics, and resolving the trapped–passing boundary (which the velocity grid in GS2 is explicitly designed for) is important. Naturally, our spectral approach is not as well suited to resolving these sharp features. In the problems we will target, however, collisions and nonlinearity will help to smear away sharp features in the distribution function.

5.4 Summary of results

These preliminary results show that our model is capable of reproducing gyrokinetic results for linear instabilities and zonal flow dynamics at spectral resolution comparable to conventional gyrokinetic velocity resolution, even with the simplest possible closure scheme. In the linear instability comparisons we see that the Beer $4+2$ gyrofluid model is nearly as accurate at much lower resolution. Since our spectral system is effectively a generalized gyrofluid model, implementing gyrofluid-like closures in our system will be straightforward. Thus the Beer results are encouraging, since they indicate that using more advanced closure schemes in our system can lead to gyrokinetic accuracy at sub-kinetic velocity resolution. However, there are cases where we will still need to use higher velocity resolution, such as in the zonal flow results. In this case, the failure of the Beer model to produce the undamped residual shows the shortcomings of using low velocity resolution with gyrofluid-like closures when the correct solution has sharp features in velocity space.