Laguerre–Hermite pseudo-spectral velocity formulation of gyrokinetics

Abstract

First-principles simulations of tokamak turbulence have proven to be of great value in recent decades. We develop a pseudo-spectral velocity formulation of the turbulence equations that smoothly interpolates between the highly efficient but lower resolution three-dimensional (3-D) gyrofluid representation and the conventional but more expensive 5-D gyrokinetic representation. Our formulation is a projection of the nonlinear gyrokinetic equation onto a Laguerre–Hermite velocity-space basis. We discuss issues related to collisions, closures and entropy. While any collision operator can be used in the formulation, we highlight a model operator that has a particularly sparse Laguerre–Hermite representation, while satisfying conservation laws and the H theorem. Free streaming, magnetic drifts and nonlinear phase mixing each give rise to closure problems, which we discuss in relation to the instabilities of interest and to free energy conservation. We show that the model is capable of reproducing gyrokinetic results for linear instabilities and zonal flow dynamics. Thus the final model is appropriate for the study of instabilities, turbulence and transport in a wide range of geometries, including tokamaks and stellarators.

1 Introduction

Tokamaks and stellarators are magnetic confinement fusion (MCF) concepts which employ toroidal magnetic fields to confine thermonuclear plasmas. Steep gradients of density, momentum and temperature are intrinsic to most MCF concepts. For the most part, the plasma column itself is macroscopically stable, but these gradients make it prone to microscopic instabilities and non-thermal fluctuations. These micro-instabilities are well described by gyrokinetics (Catto 1978; Antonsen & Lane 1980; Frieman & Chen 1982). Many gyrokinetic models, algorithms, and codes have been developed and deployed over a period of decades. The first simulation algorithms used for gyrokinetics were particle-in-cell (PIC) algorithms (Lee 1983; Kotschenreuther 1990; Parker & Lee 1992; Dimits & Lee 1993; Denton & Kotschenreuther 1995; Dimits et al. 1996). Later, approaches that involved fluid moments of the gyrokinetic equation were developed. These gyrofluid (or Landau fluid) models used sophisticated (but linear) closures to model kinetic effects such as Landau damping and finite Larmor radius (FLR) effects (Hammett & Perkins 1990; Hammett, Dorland & Perkins 1992; Dorland & Hammett 1993; Hammett et al. 1993; Beer & Hammett 1996; Snyder & Hammett 2001); however, these reduced models struggled to accurately predict thermal transport levels, in part due to inaccurate zonal flow dynamics (Rosenbluth & Hinton 1998; Dimits et al. 2000). The gyrokinetic codes that became the most widely used emerged around this time (Kotschenreuther, Rewoldt & Tang 1995; Dorland et al. 2000; Jenko & Dorland 2001; Candy & Waltz 2003). These Eulerian codes (GS2, GYRO and GENE) solve the nonlinear gyrokinetic equations on a fixed, five-dimensional mesh in phase space $(\boldsymbol{x},\boldsymbol{v})$ . Since the late 1990s, these codes have been used and cited in hundreds of papers to analyse small-scale instabilities and turbulence in tokamaks, stellarators and laboratory dipole experiments; gyrokinetic codes have also recently been repurposed for studying turbulence in astrophysical plasmas. ¹ There are also many newer approaches, including semi-Lagrangian (Grandgirard et al. 2007) and modern PIC approaches (Jolliet et al. 2007). Each approach has distinct advantages, and most codes have been benchmarked extensively with one another.

A shared weakness of these approaches is the inability to run reliably with low velocity resolution. This is despite the fact that early $\unicode[STIX]{x1D6FF}f$ PIC simulations produced accurate estimates of ion heat flux for as few as 2–4 particles per cell (Dimits et al. 2000), indicating that one does not always need to resolve the fine details of the perturbed distribution function to calculate key fluctuation properties accurately. (Here, a $\unicode[STIX]{x1D6FF}f$ simulation refers to a simulation in which the full distribution function $F$ is split into a time independent background $F_{0}$ and a perturbed distribution function $\unicode[STIX]{x1D6FF}f$ that is evolved.) However, using this few particles in $\unicode[STIX]{x1D6FF}f$ PIC is problematic due to Monte Carlo shot noise (Nevins et al. 2005; Lin et al. 2007) and inability to resolve highly anisotropic spatial structures (Idomura, Wakatani & Tokuda 2000; Nevins et al. 2005). Further, the basic conservation laws are not built in, and there is no inexpensive and rigorous way to measure or control the growth of noise as the simulation progresses. Full- $f$ PIC simulations for tokamaks and stellarators, on the other hand, lack the flexibility to run with relatively constant precision at low resolution (i.e. with few particles). Conventional Eulerian gyrokinetic algorithms also lack the flexibility to run with very low velocity-space resolution while controlling discretization artefacts associated with velocity-space filamentation, although some codes such as AstroGK (Numata et al. 2010) use pseudo-spectral grids to increase the accuracy of integrals.

In this work we present a pseudo-spectral velocity formulation of gyrokinetics. Unlike the Eulerian algorithms described above that discretize velocity space on a mesh, our approach projects the gyrokinetic equation onto a velocity basis composed of Laguerre and Hermite polynomials. ² Projecting the distribution function onto this basis produces fluid-like quantities, which correspond to density, parallel momentum, etc. The result is a formulation that is equivalent to a generalized gyrofluid system with arbitrarily many moments, which at high velocity-space resolution corresponds directly to conventional gyrokinetic Eulerian models. A key advantage of our approach is the flexibility to use very low velocity-space resolution within the same framework, where the system corresponds precisely to the well-established gyrofluid models described above. Further, since our approach produces a system of fluid-quantity conservation laws, one can maintain fundamental properties like number, momentum and energy conservation, even at very low resolution. Any inaccuracy from using low resolution can be isolated to the shortcomings of the closure assumptions used. For example, truncating the system appropriately reproduces the reduced double-adiabatic equations, which conserve number, momentum and the first and second adiabatic invariants. Importantly, since the worst case scenario is the level of accuracy of the original gyrofluid models (which used only 6 moments), our model should be able to obtain fairly precise estimates of most quantities of interest at resolutions only modestly higher than that of the original gyrofluid models. We also present a model collision operator which has a particularly sparse representation in the Laguerre–Hermite basis. This collision operator maintains conservation laws and the H theorem.

The particular application we target is the modelling of MCF turbulence in the context of whole-device modelling (WDM). The WDM idea is to simulate the full operation of a tokamak or stellarator, from the currents in the external magnets to heating systems, to turbulence, etc. with the run-time flexibility to use a range of models for any given process in any given calculation, depending on the fidelity required for that process for that simulation run. Our model can run at relatively low velocity resolution as an inexpensive WDM module, producing a significant performance advantage over other approaches, or at high resolution at costs comparable to the state of the art. Thus, a Laguerre–Hermite pseudo-spectral gyrokinetic module for WDM applications promises good fidelity at modest cost, together with straightforward validation at high resolution.

This paper is organized as follows: in § 2, we introduce the gyrokinetic equation, which is the starting point of our approach. In § 3, we project this equation onto the Laguerre–Hermite basis, discussing details of our model collision operator in § 3.5. Section 4 discusses free energy and its conservation in the context of this formulation. We present linear results of our model in § 5, and discuss conclusions and future work in § 6.

2 Gyrokinetic description of dynamics

The literature of gyrokinetics and gyrokinetic turbulence is vast and mature. Our notation and philosophy follow four particular references: Antonsen & Lane (1980), Frieman & Chen (1982), Barnes et al. (2010) and Abel et al. (2013). Detailed discussions that appear in these papers are not repeated here unless necessary. In the context of the latter, this report describes a potentially useful representation of the nonlinear $\unicode[STIX]{x1D6FF}f$ gyrokinetic equation [equation (108) of Abel et al. (2013)] in the absence of strong equilibrium flows and in the absence of electromagnetic fluctuations. Our immediate goal is to derive a model appropriate for radially local, flux-tube simulations, such as would be used to describe the turbulent fluxes through a given flux surface of a global WDM simulation. Generalization of our final equations for radially global applications would require the derivation of physically plausible, radially non-periodic, gyrokinetic boundary conditions and a more nuanced treatment of short wavelength finite Larmor radius effects. Inclusion of electromagnetic fluctuations challenges the kinds of closure models we favour, even in the radially local, flux-tube limit (Snyder & Hammett 2001). Generalizing the model to include strong flows, on the other hand, is straightforward.

We begin with a non-dimensional, normalized form of the electrostatic gyrokinetic equation:

(2.1) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}g_{s}}{\unicode[STIX]{x2202}t}+\left[v_{ts}v_{\Vert }\hat{\boldsymbol{b}}+\langle \boldsymbol{v}_{E}\rangle _{\boldsymbol{R}}+\frac{\unicode[STIX]{x1D70F}_{s}}{Z_{s}}\boldsymbol{v}_{d}\right]\boldsymbol{\cdot }\unicode[STIX]{x1D735}h_{s}+\langle \boldsymbol{v}_{E}\rangle _{\boldsymbol{R}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}F_{Ms}-v_{ts}\unicode[STIX]{x1D707}(\hat{\boldsymbol{b}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B)\frac{\unicode[STIX]{x2202}h_{s}}{\unicode[STIX]{x2202}v_{\Vert }}=C(h_{s}).\end{eqnarray}$$

This equation describes the evolution of the fluctuating gyrokinetic distribution function $g_{s}=g_{s}(\boldsymbol{R},v_{\Vert },\unicode[STIX]{x1D707},t)$ for a particular species $s$ . All quantities in the above equation (e.g. $h,g,\boldsymbol{v}_{E},v_{\Vert },\unicode[STIX]{x1D707},t$ , etc.) have been non-dimensionalized, with the specific normalizations explained in appendix A. The total distribution function, $F_{s}=F_{0s}+\unicode[STIX]{x1D6FF}f_{s}=F_{Ms}(1-Z_{s}\unicode[STIX]{x1D6F7}/\unicode[STIX]{x1D70F}_{s})+h_{s}$ is a Maxwellian with a small Boltzmann component [ $\propto \unicode[STIX]{x1D6F7}(\boldsymbol{r},t)$ , the electrostatic potential] and a general perturbation, $h_{s}(\boldsymbol{R},v_{\Vert },\unicode[STIX]{x1D707},t)$ . The total gyroaveraged distribution function, $\langle F_{s}\rangle _{\boldsymbol{R}}=F_{0s}+g_{s}$ , with $g_{s}=h_{s}-(Z_{s}/\unicode[STIX]{x1D70F}_{s})\langle \unicode[STIX]{x1D6F7}\rangle _{\boldsymbol{R}}F_{Ms}$ , describes the probability of finding a particle of species $s$ with guiding centre (or gyrocentre) position $\boldsymbol{R}$ , velocity parallel to the magnetic field $v_{\Vert }$ , and magnetic moment $\unicode[STIX]{x1D707}=v_{\bot }^{2}/2B$ , where $v_{\bot }$ is the speed in the plane perpendicular to the magnetic field. The electrostatic potential $\unicode[STIX]{x1D6F7}=\unicode[STIX]{x1D6F7}(\boldsymbol{r},t)$ is a function of particle position $\boldsymbol{r}$ , where $\boldsymbol{r}=\boldsymbol{R}+\unicode[STIX]{x1D746}$ so that $\unicode[STIX]{x1D746}$ is the gyroradius vector that rotates at the gyrofrequency $\unicode[STIX]{x1D6FA}$ and points from the gyrocentre (at $\boldsymbol{R}$ ) to the particle (at $\boldsymbol{r}$ ). The equilibrium magnetic field $\boldsymbol{B}$ has magnitude $B$ and direction $\hat{\boldsymbol{b}}=\boldsymbol{B}/B$ .

The equilibrium distribution function is a Maxwellian (which we take to have no flows),

(2.2) $$\begin{eqnarray}F_{0s}=F_{Ms}=\frac{n_{s}}{(2\unicode[STIX]{x03C0}v_{ts}^{2})^{3/2}}\text{e}^{-v_{\Vert }^{2}/2-\unicode[STIX]{x1D707}B},\end{eqnarray}$$

and the fluctuating gyroaveraged distribution function satisfies $g\ll F_{M}$ . We also have the following dimensionless species parameters: equilibrium density $n_{s}(\unicode[STIX]{x1D6F9})$ , equilibrium temperature $\unicode[STIX]{x1D70F}_{s}(\unicode[STIX]{x1D6F9})$ , mass $m_{s}$ , charge $Z_{s}$ and thermal velocity $v_{ts}=\sqrt{\unicode[STIX]{x1D70F}_{s}/m_{s}}$ . Here $\unicode[STIX]{x1D6F9}=\langle \unicode[STIX]{x1D6F9}\rangle _{\boldsymbol{R}}=\langle \unicode[STIX]{x1D6F9}\rangle _{\boldsymbol{r}}$ is the poloidal flux, which labels flux surfaces and effectively serves as a radial coordinate. In flux-tube geometry, all equilibrium quantities such as $n_{s}$ , $\unicode[STIX]{x1D70F}_{s}$ , $B$ (and therefore $\unicode[STIX]{x1D70C}$ ) are assumed to be constant in each radially local flux-tube domain.

The notation $\langle \cdots \rangle _{\boldsymbol{R}}$ denotes a gyroaverage taken at constant $\boldsymbol{R}$ . ³ The gyroaveraged $\boldsymbol{E}\times \boldsymbol{B}$ velocity is $\langle \boldsymbol{v}_{E}\rangle _{\boldsymbol{R}}=\hat{\boldsymbol{b}}\times \unicode[STIX]{x1D735}\langle \unicode[STIX]{x1D6F7}\rangle _{\boldsymbol{R}}$ , the magnetic drift velocity is $\boldsymbol{v}_{d}=(v_{\Vert }^{2}+\unicode[STIX]{x1D707}B)\hat{\boldsymbol{b}}\times \unicode[STIX]{x1D735}B/B^{2}$ , and we take $\unicode[STIX]{x1D735}\equiv \unicode[STIX]{x2202}/\unicode[STIX]{x2202}\boldsymbol{R}$ throughout. The collision operator, which we define in § 3.5, is denoted by $C(h)$ . ⁴ Finally, the potential is determined by the quasineutrality equation, which we also write in terms of $h$ :

(2.3) $$\begin{eqnarray}\mathop{\sum }_{s}Z_{s}n_{s}\int \text{d}^{3}\boldsymbol{v}\,\langle h_{s}\rangle _{\boldsymbol{r}}=\mathop{\sum }_{s}\frac{n_{s}Z_{s}^{2}}{\unicode[STIX]{x1D70F}_{s}}\unicode[STIX]{x1D6F7}.\end{eqnarray}$$

In the limit of a single hydrogenic ion species with Boltzmann electrons, this reduces to

(2.4) $$\begin{eqnarray}\int \text{d}^{3}\boldsymbol{v}\,\langle h\rangle _{\boldsymbol{ r}}=\unicode[STIX]{x1D70F}_{e}^{-1}[\unicode[STIX]{x1D6F7}-\langle \langle \unicode[STIX]{x1D6F7}\rangle \rangle ]+\unicode[STIX]{x1D6F7},\end{eqnarray}$$

where $\unicode[STIX]{x1D70F}_{e}=T_{e}/T_{i}$ in this limit of $T_{i}=T_{\text{ref}}$ . (Note that several authors choose the reciprocal definition for $\unicode[STIX]{x1D70F}_{e}$ .) The flux-surface average is denoted by $\langle \langle \cdots \rangle \rangle$ . This term comes from the standard correction for Boltzmann electrons (Dorland & Hammett 1993), which ensures that the flux-surface-averaged electron density perturbation vanishes.

Equations (2.1) and (2.4) describe self-consistent electrostatic gyrokinetic dynamics in field-line-following coordinates, including magnetically trapped particles and as yet unspecified collisional physics. These equations can be solved in general geometry (tokamak and stellarator) to find the instabilities, fluctuation spectra and turbulent fluxes of particles and energy. Landau damping, trapped particle effects, and other ‘kinetic’ phenomena are described by this comprehensive approach.

3 Laguerre–Hermite pseudo-spectral formulation

Our basic approach is to project the velocity dependence of the fluctuating distribution functions $h$ and $g$ onto orthonormal polynomials. The result is a spectral representation ⁵ within which the amplitudes of the polynomial basis functions, which can be interpreted as fluid moments, ⁶ become the dynamical variables. Thus instead of keeping track of $g(v_{\Vert },\unicode[STIX]{x1D707}B)$ with a finite-difference algorithm on a discretized $(v_{\Vert },\unicode[STIX]{x1D707}B)$ grid, we evolve coupled fluid-like quantities, such as density, parallel momentum, etc., which emerge from the projection. Specifically, we project onto Hermite polynomials in the $v_{\Vert }$ coordinate and Laguerre polynomials in the $\unicode[STIX]{x1D707}B$ coordinate. We choose this basis in part because the resulting system is identical at low resolution to the gyrofluid models of Dorland & Hammett (1993), and Beer & Hammett (1996). The Laguerre and Hermite polynomials are also orthogonal with respect to a Maxwellian. Further, the Laguerre and Hermite polynomials are eigenfunctions of our model collision operator, described in § 3.5.

At fixed accuracy, pseudo-spectral algorithms are very efficient and thus produce a small footprint in device memory (Boyd 2001). This makes it possible to move a well-resolved gyrokinetic simulation of ion temperature gradient-driven tokamak turbulence from hundreds or thousands of cores on a massively parallel supercomputer to a single graphics processor (GPU). Additionally, the regular memory access patterns of a pseudo-spectral formulation are well matched to the advantages of the GPU architecture (or any comparably high-performance vector hardware).

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 1995), Fourier techniques can be readily applied; periodic (and otherwise appropriate) domains are guaranteed by the asymptotics of the problem (Abel et al. 2013). 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. 1995; Jenko et al. 2000; Candy & Waltz 2003). 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. ⁷

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 (1996),

(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 (1992). 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 (1993). 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 (1993), 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 1993). 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 (1938), 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. ⁸

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 (2005). 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 1964), which is itself a generalization of the Lenard–Bernstein collision operator (Lenard & Bernstein 1958). 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 (2007). 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 1977). 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. (2008) 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 calculated ⁹ 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. (2008), 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. 2016).

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 (1990), which was used and extended in the gyrofluid models of Dorland & Hammett (1993), Beer & Hammett (1996) and Snyder & Hammett (2001). Further, Smith (1997) 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. ¹⁰

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 Free energy

The importance of conserved quantities in turbulence modelling is well appreciated (Sugama et al. 1996). In multi-scale gyrokinetics, energy flows between fluctuations and mean fields, and is conserved on the transport scale. The relevant fluctuating conserved quantity in the absence of driving and damping is the free energy (Sugama et al. 1996; Abel et al. 2013). In the electrostatic gyrokinetic formalism, the (normalized) free energy defined by

(4.1) $$\begin{eqnarray}\displaystyle W & = & \displaystyle \mathop{\sum }_{s}\int \text{d}^{3}\boldsymbol{r}\int \text{d}^{3}\boldsymbol{v}\,\frac{n_{s}\unicode[STIX]{x1D70F}_{s}\langle \unicode[STIX]{x1D6FF}f_{s}^{2}\rangle _{\boldsymbol{ r}}}{2F_{Ms}}=\mathop{\sum }_{s}\int \text{d}^{3}\boldsymbol{r}\int \text{d}^{3}\boldsymbol{v}\,\frac{n_{s}\unicode[STIX]{x1D70F}_{s}\langle h_{s}^{2}\rangle _{\boldsymbol{ r}}}{2F_{Ms}}-\mathop{\sum }_{s}\int \text{d}^{3}\boldsymbol{r}\frac{n_{s}Z_{s}^{2}}{2\unicode[STIX]{x1D70F}_{s}}\unicode[STIX]{x1D6F7}^{2}\nonumber\\ \displaystyle & = & \displaystyle \mathop{\sum }_{s}\int \text{d}^{3}\boldsymbol{R}\int \text{d}^{3}\boldsymbol{v}\,\frac{n_{s}\unicode[STIX]{x1D70F}_{s}~h_{s}^{2}}{2F_{Ms}}-\mathop{\sum }_{s}\int \text{d}^{3}\boldsymbol{r}\frac{n_{s}Z_{s}^{2}}{2\unicode[STIX]{x1D70F}_{s}}\unicode[STIX]{x1D6F7}^{2}\end{eqnarray}$$

is conserved in the absence of drive and damping (Krommes & Hu 1993; Howes et al. 2006). The free energy evolves in time according to

(4.2) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}W}{\unicode[STIX]{x2202}t}=\mathop{\sum }_{s}\int \text{d}^{3}\boldsymbol{r}\int \text{d}^{3}\boldsymbol{v}\,\frac{n_{s}\unicode[STIX]{x1D70F}_{s}}{F_{Ms}}\,\left\langle h_{s}\frac{\unicode[STIX]{x2202}g_{s}}{\unicode[STIX]{x2202}t}\right\rangle _{\boldsymbol{r}}=\mathop{\sum }_{s}\int \text{d}^{3}\boldsymbol{R}\int \text{d}^{3}\boldsymbol{v}\,\frac{n_{s}\unicode[STIX]{x1D70F}_{s}}{F_{Ms}}\,h_{s}\frac{\unicode[STIX]{x2202}g_{s}}{\unicode[STIX]{x2202}t}.\end{eqnarray}$$

Thus, given (2.1) for each species, one calculates $\unicode[STIX]{x2202}W/\unicode[STIX]{x2202}t$ by multiplying each species’ gyrokinetic equation by $n_{s}\unicode[STIX]{x1D70F}_{s}h_{s}/F_{0s}$ , integrating the result over the entire phase space, and summing over species.

We will now explore the conservation of $W$ in the Laguerre–Hermite basis. Expanding and projecting $h_{s}$ in (4.1), we have

(4.3) $$\begin{eqnarray}W=\mathop{\sum }_{s}n_{s}\unicode[STIX]{x1D70F}_{s}\int \text{d}^{3}\boldsymbol{R}\,\mathop{\sum }_{\ell =0}^{{\mathcal{L}}-1}\mathop{\sum }_{m=0}^{{\mathcal{M}}-1}\frac{1}{2}H_{\ell ,m}^{2}-\mathop{\sum }_{s}\frac{n_{s}Z_{s}^{2}}{2\unicode[STIX]{x1D70F}_{s}}\int \text{d}^{3}\boldsymbol{r}\,\unicode[STIX]{x1D6F7}^{2},\end{eqnarray}$$

where for now we have truncated the Laguerre–Hermite expansion at some ${\mathcal{L}}$ and ${\mathcal{M}}$ . We will primarily address free energy conservation under the assumption of closure by truncation, as described in § 3.6. We also briefly discuss subtleties of free energy conservation in the Beer gyrofluid model at the end of this section. Further, we can always recover the fully kinetic limit by taking ${\mathcal{L}},{\mathcal{M}}\rightarrow \infty$ .

The free energy for each species $W_{s}$ evolves according to ¹¹

(4.4) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}W_{s}}{\unicode[STIX]{x2202}t}=n_{s}\unicode[STIX]{x1D70F}_{s}\int \text{d}^{3}\boldsymbol{R}\left[\mathop{\sum }_{\ell =0}^{{\mathcal{L}}-1}\mathop{\sum }_{m=0}^{{\mathcal{M}}-1}H_{\ell ,m}\frac{\unicode[STIX]{x2202}G_{\ell ,m}}{\unicode[STIX]{x2202}t}\right].\end{eqnarray}$$

Most of the terms in this large sum cancel, leaving only contributions that involve perturbations with either $\ell \geqslant {\mathcal{L}}$ or $m\geqslant {\mathcal{M}}$ . These remainder contributions are precisely where closures are required.

Since there are several different types of terms in the $G_{\ell ,m}$ evolution equations, it is convenient to consider the free energy evolution of each separately. The parallel convective derivative terms (including the $\unicode[STIX]{x1D735}_{\Vert }\ln B$ magnetic trapping terms) give

(4.5) $$\begin{eqnarray}\displaystyle \frac{1}{n_{s}\unicode[STIX]{x1D70F}_{s}v_{ts}}\frac{\unicode[STIX]{x2202}W_{\Vert s}}{\unicode[STIX]{x2202}t} & = & \displaystyle -\int \text{d}^{3}\boldsymbol{R}\left[\mathop{\sum }_{\ell =0}^{{\mathcal{L}}-1}\sqrt{{\mathcal{M}}}H_{\ell ,{\mathcal{M}}-1}\unicode[STIX]{x1D735}_{\Vert }H_{\ell ,{\mathcal{M}}}\right]\nonumber\\ \displaystyle & & \displaystyle +\,\left[\mathop{\sum }_{\ell =0}^{{\mathcal{L}}-1}[-(\ell +1)\sqrt{{\mathcal{M}}}H_{\ell ,{\mathcal{M}}-1}H_{\ell ,{\mathcal{M}}}-\ell \sqrt{{\mathcal{M}}}H_{\ell ,{\mathcal{M}}-1}H_{\ell -1,{\mathcal{M}}}]\unicode[STIX]{x1D735}_{\Vert }\ln B\right]\nonumber\\ \displaystyle & & \displaystyle +\,\left[\mathop{\sum }_{m=0}^{{\mathcal{M}}-1}[{\mathcal{L}}\sqrt{m}H_{{\mathcal{L}}-1,m}H_{{\mathcal{L}},m-1}]\unicode[STIX]{x1D735}_{\Vert }\ln B\right],\end{eqnarray}$$

where we have used the identity

(4.6) $$\begin{eqnarray}\unicode[STIX]{x1D735}_{\Vert }(FG)=B\unicode[STIX]{x1D735}_{\Vert }\left(\frac{FG}{B}\right)+FG\unicode[STIX]{x1D735}_{\Vert }\ln B\end{eqnarray}$$

repeatedly, and also used the fact that $B\unicode[STIX]{x1D735}_{\Vert }(FG/B)$ has the form of a total divergence, which vanishes upon integration over all space (keeping in mind that $B$ is the Jacobian of our coordinates).

Each of these terms describes the coupling of one resolved Laguerre–Hermite moment $(\ell <{\mathcal{L}},m<{\mathcal{M}})$ to one unresolved moment. In the first sum on the right-hand side, there are ${\mathcal{L}}$ unresolved moments, each with $m={\mathcal{M}}$ ; each of these moments has contributions from $\int \text{d}^{3}\boldsymbol{v}\,v_{\Vert }^{{\mathcal{M}}}v_{\bot }^{2\ell }g$ ; all relate to the largest scale in $g(v_{\Vert })$ that cannot be resolved. There are comparable limitations on any discrete $v$ -space representation of $g(v_{\Vert })$ . Similarly, in the second and third sums on the right-hand side, there are ${\mathcal{M}}$ and ${\mathcal{L}}$ additional unresolved moments that must be closed, respectively. When closing by truncation, we take $H_{\ell ,{\mathcal{M}}}=H_{{\mathcal{L}},m-1}=H_{\ell ,{\mathcal{M}}}=H_{\ell -1,{\mathcal{M}}}=0$ , so all terms in these three sums vanish. Thus closure by truncation gives $\unicode[STIX]{x2202}W_{\Vert s}/\unicode[STIX]{x2202}t=0$ .

Note that in the kinetic limit, $\unicode[STIX]{x2202}W_{\Vert s}/\unicode[STIX]{x2202}t\rightarrow 0$ as long as $G_{{\mathcal{L}},{\mathcal{M}}}\rightarrow 0$ ‘fast enough’ as ${\mathcal{L}},{\mathcal{M}}\rightarrow \infty$ , which is the expected result. This will be the case for all of the parts of the free energy except the drive and damping parts.

The toroidal drift terms give

(4.7) $$\begin{eqnarray}\displaystyle \frac{Z_{s}}{n_{s}\unicode[STIX]{x1D70F}_{s}^{2}}\frac{\unicode[STIX]{x2202}W_{ds}}{\unicode[STIX]{x2202}t} & = & \displaystyle -\int \text{d}^{3}\boldsymbol{R}\left[\mathop{\sum }_{\ell =0}^{{\mathcal{L}}-1}\sqrt{{\mathcal{M}}({\mathcal{M}}+1)}H_{\ell ,{\mathcal{M}}-1}\text{i}\unicode[STIX]{x1D714}_{d}H_{\ell ,{\mathcal{M}}+1}\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,\sqrt{{\mathcal{M}}({\mathcal{M}}-1)}H_{\ell ,{\mathcal{M}}-2}\text{i}\unicode[STIX]{x1D714}_{d}H_{\ell ,{\mathcal{M}}}\right]+\mathop{\sum }_{m=0}^{{\mathcal{M}}-1}{\mathcal{L}}H_{{\mathcal{L}}-1,m}\text{i}\unicode[STIX]{x1D714}_{d}H_{{\mathcal{L}},m},\end{eqnarray}$$

where here $\text{i}\unicode[STIX]{x1D714}_{d}(FG)$ has the form of a total divergence and vanishes upon integration over all space. There are $2{\mathcal{L}}$ unresolved moments in the first sum, and ${\mathcal{M}}$ unresolved moments in the second sum. As before, closure by truncation ensures that all of these terms vanish, giving $\unicode[STIX]{x2202}W_{ds}/\unicode[STIX]{x2202}t=0$ .

The nonlinear terms give

(4.8) $$\begin{eqnarray}\frac{1}{n_{s}\unicode[STIX]{x1D70F}_{s}}\frac{\unicode[STIX]{x2202}W_{NL,s}}{\unicode[STIX]{x2202}t}=-\int \text{d}^{3}\boldsymbol{R}\,\mathop{\sum }_{\ell =0}^{{\mathcal{L}}-1}\mathop{\sum }_{m=0}^{{\mathcal{M}}-1}H_{\ell ,m}\mathop{\sum }_{k={\mathcal{L}}}^{\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}$$

This remainder term describes nonlinear coupling of resolved moments $H_{\ell <{\mathcal{L}},m}$ to unresolved moments $H_{k\geqslant {\mathcal{L}},m}$ via FLR corrections to the $\boldsymbol{E}\times \boldsymbol{B}$ drift. Once again, closure by truncation ensures that all of these terms vanish, giving $\unicode[STIX]{x2202}W_{NL,s}/\unicode[STIX]{x2202}t=0$ .

For appropriate choices of density and temperature gradients, the drive terms can cause the free energy to increase:

(4.9) $$\begin{eqnarray}\displaystyle \frac{1}{n_{s}\unicode[STIX]{x1D70F}_{s}}\frac{\unicode[STIX]{x2202}W_{\ast s}}{\unicode[STIX]{x2202}t} & = & \displaystyle \int \text{d}^{3}\boldsymbol{R}\,\mathop{\sum }_{\ell =0}^{{\mathcal{L}}-1}\left[H_{\ell ,0}\text{i}\unicode[STIX]{x1D714}_{\ast }\!\left(\frac{1}{L_{n}}{\mathcal{J}}_{\ell }\unicode[STIX]{x1D6F7}+\frac{1}{L_{T}}[\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)\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,\frac{1}{\sqrt{2}}\frac{1}{L_{T}}H_{\ell ,2}\text{i}\unicode[STIX]{x1D714}_{\ast }{\mathcal{J}}_{\ell }\unicode[STIX]{x1D6F7}\right]\nonumber\\ \displaystyle & = & \displaystyle \int \text{d}^{3}\boldsymbol{r}\,\left(\frac{\bar{n}}{L_{n}}+\frac{\bar{T}}{L_{T}}\right)\text{i}\unicode[STIX]{x1D714}_{\ast }\unicode[STIX]{x1D6F7}\neq 0.\end{eqnarray}$$

These terms require no closure. Although one could contemplate taking ${\mathcal{J}}_{{\mathcal{L}}}=0$ in these terms proportional to $\unicode[STIX]{x1D714}_{\ast }$ (to be consistent with closure by truncation where $G_{{\mathcal{L}},0}=H_{{\mathcal{L}},0}$ implies ${\mathcal{J}}_{{\mathcal{L}}}=0$ ), we show in appendix B that this can result in spurious instability. Thus in these $\unicode[STIX]{x1D714}_{\ast }$ terms, we revert to (3.18) for the definition of ${\mathcal{J}}_{{\mathcal{L}}}$ .

Because the Laguerre–Hermite moments are eigenfunctions of the velocity-space derivatives of our collision operator, each moment contributes $-\unicode[STIX]{x1D708}(b+2\ell +m)H_{\ell ,m}^{2}$ to $\text{d}W/\text{d}t$ . The field–particle (conservation) terms are more complicated, but as shown in appendix E, never result in a net increase $W$ after summing over all $(\ell ,m)$ . Thus, the collision operator can only decrease the total free energy or leave it unchanged:

(4.10) $$\begin{eqnarray}\displaystyle \frac{1}{n_{s}\unicode[STIX]{x1D70F}_{s}}\frac{\unicode[STIX]{x2202}W_{Cs}}{\unicode[STIX]{x2202}t} & = & \displaystyle \int \text{d}^{3}\boldsymbol{R}\int \text{d}^{3}\boldsymbol{v}\,\frac{h_{s}}{F_{Ms}}\,C(h_{s})\nonumber\\ \displaystyle & = & \displaystyle -\unicode[STIX]{x1D708}\int \text{d}^{3}\boldsymbol{R}\left[\mathop{\sum }_{\ell =0}^{{\mathcal{L}}-1}\mathop{\sum }_{m=0}^{{\mathcal{M}}-1}[(b+2\ell +m)|H_{\ell ,m}|^{2}]-(|\bar{u}_{\Vert }|^{2}+|\bar{u}_{\bot }|^{2}+3|\bar{T}|^{2})\right]\nonumber\\ \displaystyle & {\leqslant} & \displaystyle 0.\end{eqnarray}$$

Upon denoting the total drive by ${\mathcal{D}}\equiv \sum _{s}n_{s}\unicode[STIX]{x1D70F}_{s}\unicode[STIX]{x2202}W_{\ast s}/\unicode[STIX]{x2202}t$ and the total effect of collisions by ${\mathcal{C}}\equiv -\sum _{s}n_{s}\unicode[STIX]{x1D70F}_{s}\unicode[STIX]{x2202}W_{Cs}/\unicode[STIX]{x2202}t$ , we have

(4.11) $$\begin{eqnarray}\displaystyle \frac{\text{d}W}{\text{d}t} & = & \displaystyle {\mathcal{D}}-{\mathcal{C}}+\mathop{\sum }_{s}n_{s}\unicode[STIX]{x1D70F}_{s}\left(v_{ts}\frac{\unicode[STIX]{x2202}W_{\Vert }}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x1D70F}_{s}}{Z_{s}}\frac{\unicode[STIX]{x2202}W_{d}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}W_{NL}}{\unicode[STIX]{x2202}t}\right)\nonumber\\ \displaystyle & = & \displaystyle {\mathcal{D}}-{\mathcal{C}}\quad (G_{\ell ,m}=H_{\ell ,m}=0\text{ for }\ell \geqslant {\mathcal{L}}\text{ or }m\geqslant {\mathcal{M}}),\end{eqnarray}$$

where in the second line we have indicated the desired result: the total free energy $W$ is conserved in the absence of driving and collisional damping when the system is closed by truncation, and the free energy can only decrease from collisions. Further, we obtain the expected result in the kinetic limit,

(4.12) $$\begin{eqnarray}\lim _{{\mathcal{L}},{\mathcal{M}}\rightarrow \infty }\frac{\text{d}W}{\text{d}t}={\mathcal{D}}-{\mathcal{C}}.\end{eqnarray}$$

The question of free energy conservation is more complicated when using sophisticated closures in the style of Hammett & Perkins (1990), as in Beer’s gyrofluid model. The calculation of free energy conservation in our Laguerre–Hermite model in this section demonstrates that the Beer gyrofluid model conserves free energy in the absence of driving and collisions when all closure terms are dropped (after making some reinterpretations addressed in appendix G). However, one would like to be able to show that the free energy does not increase when Beer’s closure terms are included. This is complicated by the fact that the Beer closure scheme includes both dissipative and non-dissipative (reactive) terms. While the dissipative terms by nature can only decrease the free energy, the non-dissipative terms are more troubling. Nonetheless, the dissipative and non-dissipative terms collectively serve to mimic physical damping from phase mixing, so one would expect that the aggregate contribution of these terms is to decrease the free energy. It is then enough to show that the eigenvalues of the closed system are all damped in the absence of driving and collisions.

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 Linear results

We now show some preliminary linear calculations, as a proof of concept of our Laguerre–Hermite formulation of gyrokinetics. These calculations have been performed with GX, a new gyrokinetic code that uses our Laguerre–Hermite spectral velocity discretization. Numerical details and further results from GX will be reported separately. All of the results from this section use closure by truncation.

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 1996), which first appeared in Dong, Horton & Kim (1992). 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. 2000). 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. 2000). 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 2000) 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. 1993; Waltz, Kerbel & Milovich 1994). Rosenbluth & Hinton (1998) 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. 2000). The gyrofluid closure was modified to be able to capture some of the zonal flow residual with limited success (Beer & Hammett 1998).

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.

6 Summary and conclusion

In this report, we have outlined a new pseudo-spectral velocity formulation of flux-tube gyrokinetics. This is achieved by projecting the gyrokinetic equation onto a velocity basis composed of Laguerre and Hermite polynomials. A key advantage of the resulting model is the flexibility in choice of velocity-space resolution. At the lowest velocity-space resolution, the model corresponds directly to gyrofluid models. Within the same framework, the model corresponds to conventional gyrokinetic approaches at high velocity resolution. Between these limits, one has the freedom to tailor the resolution to one’s needs: one can smoothly increase resolution to improve accuracy, or one can minimize resolution to produce a performance advantage. There is also an opportunity for a dynamic fidelity refinement approach, where the resolution can be adjusted dynamically during a simulation.

The flexibility to use relatively lower velocity resolution than is used in standard gyrokinetic approaches is due to the fact that our Laguerre–Hermite formulation expresses critical conservation laws (density, momentum, energy and free energy) even at very low resolution. We have presented a model collision operator that reflects this prioritization: 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. Our collision operator is also efficiently expressed in our basis, and it satisfies the H theorem.

The main disadvantage of our approach is the need for closures due to the coupling of the spectral moments, which arises as the mathematical manifestation of the physics of phase mixing. We recognize that the value of our approach will only be fully realized when a generalized closure scheme has been found that can give gyrokinetic fidelity with sub-kinetic velocity resolution. The relative success of gyrofluid closures indicates that this is feasible, and we can use these closures directly in our model at the lowest resolution (see appendix F). We thus leave as important future work the generalization of gyrofluid-like closures to arbitrary spectral velocity resolution. This will bolster the flexibility of the model, as one will be able to capture important kinetic effects with whatever fidelity is required by the problem at hand. This includes modelling linear phase mixing from Landau damping and toroidal drifts as well as nonlinear effects such as nonlinear FLR phase mixing (Dorland & Hammett 1993; Schekochihin et al. 2009; Tatsuno et al. 2009; Chen et al. 2010; Plunk et al. 2010; Howes et al. 2011; Kawamori 2013; Numata & Loureiro 2015), and phase un-mixing from the plasma echo (Kanekar 2014; Parker & Dellar 2015; Schekochihin et al. 2016).

We also leave the extension of this model to include electromagnetic fluctuations for the future. While we could simply include electrons as a second species and find the perpendicular magnetic field fluctuations from Ampere’s law, such an approach is potentially expensive. For many applications, the Snyder & Hammett (2001) model will be the best way forward, for example. Importantly, a model with $\unicode[STIX]{x1D6FF}A_{\Vert }$ and $\unicode[STIX]{x1D6FF}B_{\Vert }$ presents no new challenges to the Laguerre–Hermite projections.

Our 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 grid resolution, even with the simplest possible closure scheme. More advanced closure schemes could lead to similar accuracy at lower cost in cases where there are not sharp features in velocity space. Nonlinear simulations of ion temperature gradient turbulence, which make use of the pseudo-spectral approach described in appendix D, will be reported in a later work. In the nonlinear system, turbulence is an additional source of regularization in velocity space, which should allow the use of lower velocity resolution in the pseudo-spectral system than in a standard grid-based approach. We have also coupled our model to the TRINITY full-torus gyrokinetic transport framework (Barnes et al. 2010). In this setting the flexibility of our model enables a significant performance advantage over standard approaches, allowing efficient and tractable whole-device modelling of tokamaks and stellarators (Highcock et al. 2016).

Finally, we note that our model is inherently parallelizable and has a small memory footprint due to its pseudo-spectral approach. This allows us to move a well-resolved gyrokinetic simulation of ion temperature gradient-driven tokamak turbulence from hundreds or thousands of cores on an massively parallel HPC device to a single graphics processor (GPU). Thus our new Laguerre–Hermite pseudo-spectral formulation has produced a high-fidelity gyrokinetic simulation code that runs at a few teraflops on a desktop computer with an inexpensive GPU.