2.1 Basic equations

We solve the full-$f$ electromagnetic gyrokinetic (EMGK) equation in the symplectic formulation (Brizard & Hahm Reference Brizard and Hahm2007), which describes the evolution of the gyrocentre distribution function $f_{s}(\boldsymbol{Z},t)=f_{s}(\boldsymbol{R},v_{\Vert },\unicode[STIX]{x1D707},t)$ for each species $s$, where $\boldsymbol{Z}$ is a phase-space coordinate composed of the guiding centre position $\boldsymbol{R}=(x,y,z)$, the parallel velocity $v_{\Vert }$ and the magnetic moment $\unicode[STIX]{x1D707}=m_{s}v_{\bot }^{2}/(2B)$. In terms of the gyrocentre Hamiltonian and the Poisson bracket in gyrocentre coordinates, and also including collisions $C[f_{s}]$ and sources $S_{s}$ (which do not derive from the bracket), the gyrokinetic equation is given byFootnote ¹

(2.1)$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}f_{s}}{\unicode[STIX]{x2202}t}+\{\,f_{s},H_{s}\}-\frac{q_{s}}{m_{s}}\frac{\unicode[STIX]{x2202}A_{\Vert }}{\unicode[STIX]{x2202}t}\frac{\unicode[STIX]{x2202}f_{s}}{\unicode[STIX]{x2202}v_{\Vert }}=C[f_{s}]+S_{s},\end{eqnarray}$$

or equivalently,

(2.2)$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}f_{s}}{\unicode[STIX]{x2202}t}+\dot{\boldsymbol{R}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f_{s}+\dot{v}_{\Vert }^{H}\frac{\unicode[STIX]{x2202}f_{s}}{\unicode[STIX]{x2202}v_{\Vert }}-\frac{q_{s}}{m_{s}}\frac{\unicode[STIX]{x2202}A_{\Vert }}{\unicode[STIX]{x2202}t}\frac{\unicode[STIX]{x2202}f_{s}}{\unicode[STIX]{x2202}v_{\Vert }}=C[f_{s}]+S_{s},\end{eqnarray}$$

where the gyrokinetic Poisson bracket is given by

(2.3)$$\begin{eqnarray}\{F,G\}=\frac{\boldsymbol{B}^{\ast }}{mB_{\Vert }^{\ast }}\boldsymbol{\cdot }\left(\unicode[STIX]{x1D735}F\frac{\unicode[STIX]{x2202}G}{\unicode[STIX]{x2202}v_{\Vert }}-\frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}v_{\Vert }}\unicode[STIX]{x1D735}G\right)-\frac{\hat{\boldsymbol{b}}}{qB_{\Vert }^{\ast }}\times \unicode[STIX]{x1D735}F\boldsymbol{\cdot }\unicode[STIX]{x1D735}G,\end{eqnarray}$$

and we take the gyrocentre Hamiltonian to be

(2.4)$$\begin{eqnarray}H_{s}={\textstyle \frac{1}{2}}m_{s}v_{\Vert }^{2}+\unicode[STIX]{x1D707}B+q_{s}\unicode[STIX]{x1D719}.\end{eqnarray}$$

Here, we have taken the long-wavelength (drift-kinetic) limit to neglect gyroaveraging of the electrostatic potential $\unicode[STIX]{x1D719}$, and we have also dropped higher-order terms in the Hamiltonian that appear in e.g. Brizard & Hahm (Reference Brizard and Hahm2007); extensions to include gyroaveraging will be included in later work, but these additions will not change the overall scheme presented here. The nonlinear phase-space characteristics are given by

(2.5)$$\begin{eqnarray}\displaystyle & \displaystyle \dot{\boldsymbol{R}}=\{\boldsymbol{R},H_{s}\}=\frac{\boldsymbol{B}^{\ast }}{B_{\Vert }^{\ast }}v_{\Vert }+\frac{\hat{\boldsymbol{b}}}{q_{s}B_{\Vert }^{\ast }}\times (\unicode[STIX]{x1D707}\unicode[STIX]{x1D735}B+q_{s}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}), & \displaystyle\end{eqnarray}$$

(2.6)$$\begin{eqnarray}\displaystyle & \displaystyle \dot{v}_{\Vert }=\dot{v}_{\Vert }^{H}-\frac{q_{s}}{m_{s}}\frac{\unicode[STIX]{x2202}A_{\Vert }}{\unicode[STIX]{x2202}t}=\{v_{\Vert },H_{s}\}-\frac{q_{s}}{m_{s}}\frac{\unicode[STIX]{x2202}A_{\Vert }}{\unicode[STIX]{x2202}t}=-\frac{\boldsymbol{B}^{\ast }}{m_{s}B_{\Vert }^{\ast }}\boldsymbol{\cdot }(\unicode[STIX]{x1D707}\unicode[STIX]{x1D735}B+q_{s}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719})-\frac{q_{s}}{m_{s}}\frac{\unicode[STIX]{x2202}A_{\Vert }}{\unicode[STIX]{x2202}t}.\qquad & \displaystyle\end{eqnarray}$$

Here, $B_{\Vert }^{\ast }=\hat{\boldsymbol{b}}\boldsymbol{\cdot }\boldsymbol{B}^{\ast }$ is the parallel component of the effective magnetic field $\boldsymbol{B}^{\ast }=\boldsymbol{B}+(m_{s}v_{\Vert }/q_{s})\unicode[STIX]{x1D735}\times \hat{\boldsymbol{b}}+\unicode[STIX]{x1D6FF}\boldsymbol{B}$, where $\boldsymbol{B}=B\hat{\boldsymbol{b}}$ is the equilibrium magnetic field and $\unicode[STIX]{x1D6FF}\boldsymbol{B}=\unicode[STIX]{x1D735}\times (A_{\Vert }\hat{\boldsymbol{b}})\approx \unicode[STIX]{x1D735}A_{\Vert }\times \hat{\boldsymbol{b}}$ is the perturbed magnetic field (assuming that the equilibrium magnetic field varies on spatial scales longer than perturbations so that $A_{\Vert }\unicode[STIX]{x1D735}\times \hat{\boldsymbol{b}}$ can be neglected). We neglect higher-order parallel compressional fluctuations of the magnetic field, so that $\unicode[STIX]{x1D6FF}\boldsymbol{B}=\unicode[STIX]{x1D6FF}\boldsymbol{B}_{\bot }$. The species charge and mass are $q_{s}$ and $m_{s}$, respectively. In (2.6), note that we have separated $\dot{v}_{\Vert }$ into a term that comes from the Hamiltonian, $\dot{v}_{\Vert }^{H}=\{v_{\Vert },H_{s}\}$, and another term proportional to the inductive component of the parallel electric field, $(q/m)\unicode[STIX]{x2202}A_{\Vert }/\unicode[STIX]{x2202}t$. We use this notation for convenience, and so that the time derivative of the parallel vector potential $A_{\Vert }$ appears explicitly. Further, we will assume a field-aligned coordinate system (e.g. Beer, Cowley & Hammett Reference Beer, Cowley and Hammett1995), and we will take the perpendicular directions to be $x$ and $y$, and the parallel direction to be $z$.

In the absence of collisions $C[f_{s}]$ and sources $S_{s}$, equation (2.1) can be recognized as a Liouville equation, which shows that the distribution function is conserved along the nonlinear characteristics. Liouville’s theorem also shows that phase-space volume is conserved,

(2.7)$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}{\mathcal{J}}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }({\mathcal{J}}\dot{\boldsymbol{R}})+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}v_{\Vert }}({\mathcal{J}}\dot{v}_{\Vert }^{H})-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}v_{\Vert }}\left({\mathcal{J}}\frac{q_{s}}{m_{s}}\frac{\unicode[STIX]{x2202}A_{\Vert }}{\unicode[STIX]{x2202}t}\right)=0,\end{eqnarray}$$

where ${\mathcal{J}}=B_{\Vert }^{\ast }$ is the Jacobian of the gyrocentre coordinates, and we will make the approximation $\hat{\boldsymbol{b}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \hat{\boldsymbol{b}}\approx 0$ so that $B_{\Vert }^{\ast }\approx B$.

We can now write the gyrokinetic equation in conservative form,

(2.8)$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}({\mathcal{J}}f_{s})}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }({\mathcal{J}}\dot{\boldsymbol{R}}f_{s})+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}v_{\Vert }}({\mathcal{J}}\dot{v}_{\Vert }^{H}f_{s})-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}v_{\Vert }}\left({\mathcal{J}}\frac{q_{s}}{m_{s}}\frac{\unicode[STIX]{x2202}A_{\Vert }}{\unicode[STIX]{x2202}t}f_{s}\right)={\mathcal{J}}C[f_{s}]+{\mathcal{J}}S_{s}.\end{eqnarray}$$

Here, we have used the symplectic formulation of electromagnetic gyrokinetics, where the parallel velocity is used as an independent variable (as opposed to the Hamiltonian formulation which uses the parallel canonical momentum $p_{\Vert }$ as an independent variable) (Hahm, Lee & Brizard Reference Hahm, Lee and Brizard1988; Brizard & Hahm Reference Brizard and Hahm2007). Notably, in the symplectic formulation, the time derivative of $A_{\Vert }$ appears explicitly in the gyrokinetic equation, equation (2.8), and $A_{\Vert }$ appears in $\boldsymbol{B}^{\ast }$ but not in the Hamiltonian.

The electrostatic potential is determined by the quasi-neutrality condition in the long-wavelength limit, given by

(2.9)$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{g}+\unicode[STIX]{x1D70E}_{pol}=\unicode[STIX]{x1D70E}_{g}-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{P}=0,\end{eqnarray}$$

with the guiding centre charge density (neglecting gyroaveraging in the long-wavelength limit)

(2.10)$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{g}=\mathop{\sum }_{s}q_{s}\int \text{d}\boldsymbol{w}{\mathcal{J}}f_{s}.\end{eqnarray}$$

Here we have defined $\text{d}\boldsymbol{w}=2\unicode[STIX]{x03C0}m_{s}^{-1}\,\text{d}v_{\Vert }\,\text{d}\unicode[STIX]{x1D707}=m_{s}^{-1}\,\text{d}v_{\Vert }\,\text{d}\unicode[STIX]{x1D707}\int \text{d}\unicode[STIX]{x1D6FC}$ as the gyrocentre velocity-space volume element $(\text{d}\boldsymbol{v}=m_{s}^{-1}\,\text{d}v_{\Vert }\,\text{d}\unicode[STIX]{x1D707}\,\text{d}\unicode[STIX]{x1D6FC}{\mathcal{J}})$ with the gyroangle $\unicode[STIX]{x1D6FC}$ integrated away and the Jacobian factored out. The polarization vector is

(2.11)$$\begin{eqnarray}\boldsymbol{P}=-\mathop{\sum }_{s}\int \text{d}\boldsymbol{w}\frac{m_{s}}{B^{2}}{\mathcal{J}}f_{s}\unicode[STIX]{x1D735}_{\bot }\unicode[STIX]{x1D719}\approx -\mathop{\sum }_{s}\frac{m_{s}n_{0s}}{B^{2}}\unicode[STIX]{x1D735}_{\bot }\unicode[STIX]{x1D719}\equiv -\unicode[STIX]{x1D716}_{\bot }\unicode[STIX]{x1D735}_{\bot }\unicode[STIX]{x1D719},\end{eqnarray}$$

where $\unicode[STIX]{x1D735}_{\bot }=\unicode[STIX]{x1D735}-\hat{\boldsymbol{b}}(\hat{\boldsymbol{b}}\boldsymbol{\cdot }\unicode[STIX]{x1D735})$ is the gradient perpendicular to the background magnetic field. We use a linearized polarization density $n_{0}$ that we take to be a constant in time, which is consistent with neglecting a second-order $E\times B$ energy term in the Hamiltonian. While the validity of this approximation in the SOL can be questioned due to large density fluctuations, a linearized polarization density is commonly used for computational efficiency (Ku et al. Reference Ku, Chang, Hager, Churchill, Tynan, Cziegler, Greenwald, Hughes, Parker and Adams2018a; Shi et al. Reference Shi, Hammett, Stoltzfus-Dueck and Hakim2019). Future work will include the nonlinear polarization density along with the second-order $E\times B$ energy term in the Hamiltonian. The quasi-neutrality condition can then be rewritten as the long-wavelength gyrokinetic Poisson equation,

(2.12)$$\begin{eqnarray}-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\mathop{\sum }_{s}\frac{m_{s}n_{0s}}{B^{2}}\unicode[STIX]{x1D735}_{\bot }\unicode[STIX]{x1D719}=\mathop{\sum }_{s}q_{s}\int \text{d}\boldsymbol{w}{\mathcal{J}}f_{s}.\end{eqnarray}$$

Even in the long-wavelength limit with no gyroaveraging, the first-order polarization charge density on the left-hand side of (2.12) incorporates some finite Larmor radius (FLR) effects.

The parallel vector potential $A_{\Vert }$ is determined by the parallel Ampère equation,

(2.13)$$\begin{eqnarray}-\unicode[STIX]{x1D6FB}_{\bot }^{2}A_{\Vert }=\unicode[STIX]{x1D707}_{0}J_{\Vert }=\unicode[STIX]{x1D707}_{0}\mathop{\sum }_{s}q_{s}\int \text{d}\boldsymbol{w}{\mathcal{J}}v_{\Vert }\,f_{s}.\end{eqnarray}$$

Note that we can also take the time derivative of this equation to get a generalized Ohm’s law which can be solved directly for $\unicode[STIX]{x2202}A_{\Vert }/\unicode[STIX]{x2202}t$, the inductive component of the parallel electric field $E_{\Vert }$ (Reynders Reference Reynders1993; Cummings Reference Cummings1994; Chen & Parker Reference Chen and Parker2001)

(2.14)$$\begin{eqnarray}-\unicode[STIX]{x1D6FB}_{\bot }^{2}\frac{\unicode[STIX]{x2202}A_{\Vert }}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D707}_{0}\mathop{\sum }_{s}q_{s}\int \text{d}\boldsymbol{w}v_{\Vert }\frac{\unicode[STIX]{x2202}({\mathcal{J}}f_{s})}{\unicode[STIX]{x2202}t}.\end{eqnarray}$$

Writing the gyrokinetic equation as

(2.15)$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}({\mathcal{J}}f_{s})}{\unicode[STIX]{x2202}t}=\frac{\unicode[STIX]{x2202}({\mathcal{J}}f_{s})}{\unicode[STIX]{x2202}t}^{\star }+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}v_{\Vert }}\left({\mathcal{J}}\frac{q_{s}}{m_{s}}\frac{\unicode[STIX]{x2202}A_{\Vert }}{\unicode[STIX]{x2202}t}f_{s}\right),\end{eqnarray}$$

where $\unicode[STIX]{x2202}({\mathcal{J}}f_{s})^{\star }/\unicode[STIX]{x2202}t$ denotes all the terms in the gyrokinetic equation (including sources and collisions) except the $\unicode[STIX]{x2202}A_{\Vert }/\unicode[STIX]{x2202}t$ term, Ohm’s law can be rewritten (after an integration by parts) as

(2.16)$$\begin{eqnarray}\left(-\unicode[STIX]{x1D735}_{\bot }^{2}+\mathop{\sum }_{s}\frac{\unicode[STIX]{x1D707}_{0}q_{s}^{2}}{m_{s}}\int \text{d}\boldsymbol{w}{\mathcal{J}}f_{s}\right)\frac{\unicode[STIX]{x2202}A_{\Vert }}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D707}_{0}\mathop{\sum }_{s}q_{s}\int \text{d}\boldsymbol{w}v_{\Vert }\frac{\unicode[STIX]{x2202}({\mathcal{J}}f_{s})}{\unicode[STIX]{x2202}t}^{\star }.\end{eqnarray}$$

As we will show in § 4, this form allows for the use of an explicit time-stepping scheme in which one can first compute $\unicode[STIX]{x2202}({\mathcal{J}}f_{s})^{\star }/\unicode[STIX]{x2202}t$ (which does not involve $\unicode[STIX]{x2202}A_{\Vert }/\unicode[STIX]{x2202}t$), then compute $\unicode[STIX]{x2202}A_{\Vert }/\unicode[STIX]{x2202}t$ and finally compute $\unicode[STIX]{x2202}({\mathcal{J}}f_{s})/\unicode[STIX]{x2202}t$. Note, however, that in some PIC approaches (Reynders Reference Reynders1993; Chen & Parker Reference Chen and Parker2001), one must expand the right-hand side of (2.16) by inserting the gyrokinetic equation so that the right-hand side involves only moments of $f_{s}$ without time derivatives. In our continuum scheme we can compute $\unicode[STIX]{x2202}({\mathcal{J}}f_{s})^{\star }/\unicode[STIX]{x2202}t$ directly and then perform the integration. Further, note that although we are using the symplectic ($v_{\Vert }$) formulation of EMGK, our Ohm’s law from (2.16) contains two integral terms which must cancel exactly. This is the root of the cancellation problem that appears in Ampère’s law in the Hamiltonian ($p_{\Vert }$) formulation, and in appendix A we show that the same cancellation problem could arise from (2.16) if the integrals are not treated consistently.

To model the effect of collisions we use a conservative Lenard–Bernstein (or Dougherty) collision operator (Lenard & Bernstein Reference Lenard and Bernstein1958; Dougherty Reference Dougherty1964),

(2.17)$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{J}}C[f]=\unicode[STIX]{x1D708}\left\{\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}v_{\Vert }}\left[(v_{\Vert }-u_{\Vert }){\mathcal{J}}f+v_{t}^{2}\frac{\unicode[STIX]{x2202}({\mathcal{J}}f)}{\unicode[STIX]{x2202}v_{\Vert }}\right]+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\left[2\unicode[STIX]{x1D707}{\mathcal{J}}f+2\unicode[STIX]{x1D707}\frac{m}{B}v_{t}^{2}\frac{\unicode[STIX]{x2202}({\mathcal{J}}f)}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\right]\right\}, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where

(2.18a,b)$$\begin{eqnarray}nu_{\Vert }=\int \text{d}\boldsymbol{w}{\mathcal{J}}v_{\Vert }\,f,\quad nu_{\Vert }^{2}+3nv_{t}^{2}=\int \text{d}\boldsymbol{w}{\mathcal{J}}(v_{\Vert }^{2}+2\unicode[STIX]{x1D707}B/m)f,\end{eqnarray}$$

with $n=\int \text{d}\boldsymbol{w}{\mathcal{J}}f$. This collision operator contains the effect of drag and pitch-angle scattering, and it conserves number, momentum and energy density. Consistent with our present long-wavelength treatment of the gyrokinetic system, finite Larmor radius effects are ignored. For simplicity we restrict ourselves to the case in which the collision frequency $\unicode[STIX]{x1D708}$ is velocity independent, i.e. $\unicode[STIX]{x1D708}\neq \unicode[STIX]{x1D708}(v)$. Further details about this collision operator, including its conservation properties and its discretization, are left to a separate paper (Francisquez et al. Reference Francisquez, Bernard, Mandell, Hammett and Hakim2020). In this work, we include only the effects of like-species collisions, which neglects electron–ion collisions and the resulting resistivity. A conservative scheme for cross-species collisions has also been implemented and will be included in later work. Extensions to a more complete collision operator are in progress.

2.2 Conservation properties

In the absence of collisions and sources, the Hamiltonian structure of the gyrokinetic system guarantees conservation of arbitrary functions of $f$ along the characteristics,

(2.19)$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}G(f)}{\unicode[STIX]{x2202}t}+\{G(f),H\}-\frac{q}{m}\frac{\unicode[STIX]{x2202}A_{\Vert }}{\unicode[STIX]{x2202}t}\frac{\unicode[STIX]{x2202}G(f)}{\unicode[STIX]{x2202}v_{\Vert }}=0,\end{eqnarray}$$

along with corresponding Casimir invariants $\int \text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}{\mathcal{J}}G(f)$, where $\text{d}\boldsymbol{R}=\text{d}x\,\text{d}y\,\text{d}z$. Thus, the system has an infinite number of conserved quantities such as the total particle number (or $L_{1}$ norm) $N=\int \text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}{\mathcal{J}}f$, the $L_{2}$ norm $M=\int \text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}{\mathcal{J}}f^{2}$ and the kinetic entropy $S=-\int \text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}{\mathcal{J}}f\ln f$ (Idomura et al. Reference Idomura, Ida, Kano, Aiba and Tokuda2008).

The system also conserves total energy, $W=W_{K}+W_{E}+W_{B}=W_{H}-W_{E}+W_{B}$, where the kinetic particle energy (neglecting the kinetic energy of the $E\times B$ flow) is

(2.20)$$\begin{eqnarray}W_{K}=\mathop{\sum }_{s}\int \text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}{\mathcal{J}}({\textstyle \frac{1}{2}}m_{s}v_{\Vert }^{2}+\unicode[STIX]{x1D707}B)f_{s},\end{eqnarray}$$

the (non-vacuum) electrostatic field energy (equivalent to the kinetic energy associated with the $E\times B$ flow of particles) is

(2.21)$$\begin{eqnarray}W_{E}=\mathop{\sum }_{s}\int \text{d}\boldsymbol{R}\frac{1}{2}\frac{m_{s}n_{0s}}{B^{2}}|\unicode[STIX]{x1D735}_{\bot }\unicode[STIX]{x1D719}|^{2}=\int \text{d}\boldsymbol{R}\frac{\unicode[STIX]{x1D716}_{\bot }}{2}|\unicode[STIX]{x1D735}_{\bot }\unicode[STIX]{x1D719}|^{2},\end{eqnarray}$$

the (perturbed) electromagnetic field energy is

(2.22)$$\begin{eqnarray}W_{B}=\int \text{d}\boldsymbol{R}\frac{1}{2\unicode[STIX]{x1D707}_{0}}|\unicode[STIX]{x1D735}_{\bot }A_{\Vert }|^{2},\end{eqnarray}$$

and

(2.23)$$\begin{eqnarray}W_{H}=\mathop{\sum }_{s}\int \text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}{\mathcal{J}}H_{s}f_{s}.\end{eqnarray}$$

(Note that $W_{H}$ is the sum of the particle kinetic energy and twice the potential energy, because every pair of particle interactions is double counted in the raw integral of $q_{s}\unicode[STIX]{x1D719}f_{s}$.)

Assuming the boundary conditions are periodic or that the distribution function vanishes at the boundary so that surface terms vanish, the evolution of these quantities can be calculated as

(2.24)$$\begin{eqnarray}\displaystyle \frac{\text{d}W_{H}}{\text{d}t} & = & \displaystyle \mathop{\sum }_{s}\int \text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}H_{s}\frac{\unicode[STIX]{x2202}({\mathcal{J}}f_{s})}{\unicode[STIX]{x2202}t}+\mathop{\sum }_{s}\int \text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}{\mathcal{J}}f_{s}\frac{\unicode[STIX]{x2202}H_{s}}{\unicode[STIX]{x2202}t}\nonumber\\ \displaystyle & = & \displaystyle -\int \text{d}\boldsymbol{R}J_{\Vert }\frac{\unicode[STIX]{x2202}A_{\Vert }}{\unicode[STIX]{x2202}t}+\int \text{d}\boldsymbol{R}\unicode[STIX]{x1D70E}_{g}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t},\end{eqnarray}$$

(2.25)$$\begin{eqnarray}\displaystyle \frac{\text{d}W_{E}}{\text{d}t} & = & \displaystyle \mathop{\sum }_{s}\int \text{d}\boldsymbol{R}\frac{m_{s}n_{0s}}{B^{2}}\unicode[STIX]{x1D735}_{\bot }\unicode[STIX]{x1D719}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\bot }\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t}=\int \text{d}\boldsymbol{R}\unicode[STIX]{x1D70E}_{g}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t},\end{eqnarray}$$

(2.26)$$\begin{eqnarray}\frac{\text{d}W_{B}}{\text{d}t}=\int \text{d}\boldsymbol{R}\frac{1}{\unicode[STIX]{x1D707}_{0}}\unicode[STIX]{x1D735}_{\bot }A_{\Vert }\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\bot }\frac{\unicode[STIX]{x2202}A_{\Vert }}{\unicode[STIX]{x2202}t}=\int \text{d}\boldsymbol{R}J_{\Vert }\frac{\unicode[STIX]{x2202}A_{\Vert }}{\unicode[STIX]{x2202}t},\end{eqnarray}$$

so that total energy is indeed conserved:

(2.27)$$\begin{eqnarray}\frac{\text{d}W}{\text{d}t}=\frac{\text{d}W_{H}}{\text{d}t}-\frac{\text{d}W_{E}}{\text{d}t}+\frac{\text{d}W_{B}}{\text{d}t}=0.\end{eqnarray}$$

3.1 Discrete equations

We use an energy-conserving discontinuous Galerkin scheme to discretize the gyrokinetic system in phase space. The scheme generalizes the algorithm of Liu & Shu (Reference Liu and Shu2000) (originally for the two-dimensional incompressible Euler and Navier–Stokes equations) to arbitrary Hamiltonian systems (Hakim et al. Reference Hakim, Hammett, Shi and Mandell2019; Shi Reference Shi2017; Shi et al. Reference Shi, Hammett, Stoltzfus-Dueck and Hakim2017). However, unlike the nodal approach used in Shi (Reference Shi2017) and Shi et al. (Reference Shi, Hammett, Stoltzfus-Dueck and Hakim2017), we use a modal DG scheme.

We start by decomposing the global phase-space domain $\unicode[STIX]{x1D6FA}$ into a structured phase-space mesh ${\mathcal{T}}$ with cells ${\mathcal{K}}_{j}\in {\mathcal{T}},j=1,\ldots ,N$. We then introduce a piecewise-polynomial approximation space for the distribution function $f(\boldsymbol{R},v_{\Vert },\unicode[STIX]{x1D707})$,

(3.1)$$\begin{eqnarray}{\mathcal{V}}_{h}^{p}=\{v:v|_{{\mathcal{K}}_{j}}\in \boldsymbol{P}^{p},\forall {\mathcal{K}}_{j}\in {\mathcal{T}}\},\end{eqnarray}$$

where $\boldsymbol{P}^{p}$ is some space of polynomials with maximum degree $p$ (by some measure). That is, $v(z)$ are polynomial functions of $\boldsymbol{z}$ in each cell, and $\boldsymbol{P}^{p}$ is the space of the linear combination of some set of multi-variate polynomials. In this work, we choose $\boldsymbol{P}^{p}$ to be an orthonormalized serendipity polynomial element space (Arnold & Awanou Reference Arnold and Awanou2011). The serendipity basis set has the advantage of using fewer basis functions while giving the same formal convergence order (although it is less accurate) as the Lagrange tensor basis, although note that, for $p=1$, the serendipity basis is equivalent to the Lagrange tensor basis. We can then obtain the discrete weak form of the gyrokinetic equation by multiplying equation (2.8) by any test function $\unicode[STIX]{x1D713}\in {\mathcal{V}}_{h}^{p}$ and integrating (by parts) in each cell

(3.2)$$\begin{eqnarray}\displaystyle & & \displaystyle \int _{{\mathcal{K}}_{j}}\text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}\unicode[STIX]{x1D713}\frac{\unicode[STIX]{x2202}({\mathcal{J}}f_{h})}{\unicode[STIX]{x2202}t}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\oint _{\unicode[STIX]{x2202}{\mathcal{K}}_{j}}\text{d}\boldsymbol{w}\,\text{d}\boldsymbol{s}_{R}\boldsymbol{\cdot }\dot{\boldsymbol{R}}_{h}\unicode[STIX]{x1D713}^{-}\widehat{{\mathcal{J}}f_{h}}+\oint _{\unicode[STIX]{x2202}{\mathcal{K}}_{j}}\text{d}\boldsymbol{R}\,\text{d}s_{w}\left(\dot{v}_{\Vert h}^{H}-\frac{q}{m}\frac{\unicode[STIX]{x2202}A_{\Vert h}}{\unicode[STIX]{x2202}t}\right)\unicode[STIX]{x1D713}^{-}\widehat{{\mathcal{J}}f_{h}}\nonumber\\ \displaystyle & & \displaystyle \qquad -\,\int _{{\mathcal{K}}_{j}}\text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}{\mathcal{J}}f_{h}\dot{\boldsymbol{R}}_{h}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}-\int _{{\mathcal{K}}_{j}}\text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}{\mathcal{J}}f_{h}\left(\dot{v}_{\Vert h}^{H}-\frac{q}{m}\frac{\unicode[STIX]{x2202}A_{\Vert h}}{\unicode[STIX]{x2202}t}\right)\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}v_{\Vert }}\nonumber\\ \displaystyle & & \displaystyle \quad =\int _{{\mathcal{K}}_{j}}\text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}\unicode[STIX]{x1D713}({\mathcal{J}}C[f_{h}]+{\mathcal{J}}S).\end{eqnarray}$$

Solving this equation for all test functions $\unicode[STIX]{x1D713}\in {\mathcal{V}}_{h}^{p}$ in all cells ${\mathcal{K}}_{j}\in {\mathcal{T}}$ yields the discretized distribution function $f_{h}\in {\mathcal{V}}_{h}^{p}$, where the subscript $h$ denotes a discrete quantity in ${\mathcal{V}}_{h}^{p}$. In the surface terms, $\text{d}\boldsymbol{s}_{R}$ is the differential element on a configuration-space surface (pointing outward normal to the surface), $\text{d}s_{w}=2\unicode[STIX]{x03C0}m_{s}^{-1}\,\text{d}\unicode[STIX]{x1D707}\,\boldsymbol{n}\boldsymbol{\cdot }(\unicode[STIX]{x2202}\boldsymbol{Z}/\unicode[STIX]{x2202}v_{\Vert })$ is the differential element on a $v_{\Vert }$ surface and the notation $\unicode[STIX]{x1D713}^{-}$ ($\unicode[STIX]{x1D713}^{+}$) indicates that the function $\unicode[STIX]{x1D713}$ is evaluated just inside (outside) the surface $\unicode[STIX]{x2202}{\mathcal{K}}_{j}$. The notation $\widehat{f}=\widehat{f}(f^{+},f^{-})$ indicates a ‘numerical flux’, which takes a single value at the cell surface and in general can depend on the solution on both sides of the surface since the solution is discontinuous at the surface. Here, we choose to use standard upwind fluxes, which depend on the local value of the phase-space characteristic flow normal to the surface evaluated at each Gaussian quadrature point on the surface. Denoting the flow as $\unicode[STIX]{x1D736}_{h}$, the upwind flux can be expressed as

(3.3)$$\begin{eqnarray}\widehat{f_{h}}={\textstyle \frac{1}{2}}(f_{h}^{+}+f_{h}^{-})-{\textstyle \frac{1}{2}}\text{sgn}(\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D736}_{h})(f_{h}^{+}-f_{h}^{-}),\end{eqnarray}$$

where $\boldsymbol{n}=\text{d}\boldsymbol{s}/|\text{d}\boldsymbol{s}|$ is the unit normal pointing out of the $\unicode[STIX]{x2202}{\mathcal{K}}_{j}$ surface.

We will introduce a subset of ${\mathcal{V}}_{h}^{p}$ where the piecewise polynomials are continuous across cell interfaces, denoted by $\overline{{\mathcal{V}}}_{h}^{p}$. As we will show later, in order to preserve energy conservation in our discrete scheme, we will require that the discrete Hamiltonian be continuous across cell interfaces, i.e. $H_{h}\in \overline{{\mathcal{V}}}_{h}^{p}$ (Liu & Shu Reference Liu and Shu2000; Shi Reference Shi2017; Shi et al. Reference Shi, Hammett, Stoltzfus-Dueck and Hakim2017; Hakim et al. Reference Hakim, Hammett, Shi and Mandell2019). Note that one can show that this ensures that the discrete phase-space characteristics, $\dot{\boldsymbol{R}}_{h}=\{\boldsymbol{R},H_{h}\}$ and $\dot{v}_{\Vert h}^{H}-(q_{s}/m_{s})\unicode[STIX]{x2202}A_{\Vert h}/\unicode[STIX]{x2202}t=\{v_{\Vert },H_{h}\}-(q_{s}/m_{s})\unicode[STIX]{x2202}A_{\Vert h}/\unicode[STIX]{x2202}t$, are also continuous across cell interfaces.Footnote ²

We must also discretize the field equations. We introduce the restriction of the phase-space mesh to configuration space, ${\mathcal{T}}^{R}$, and we denote the configuration-space cells by ${\mathcal{K}}_{j}^{R}\in {\mathcal{T}}^{R}$ for $j=1,\ldots ,N_{R}$, where $N_{R}$ is the number of configuration-space cells. We also restrict ${\mathcal{V}}_{h}^{p}$ to configuration space as

(3.4)$$\begin{eqnarray}{\mathcal{X}}_{h}^{p}={\mathcal{V}}_{h}^{p}\setminus {\mathcal{T}}^{R}.\end{eqnarray}$$

Further, we introduce the subset of polynomials that are piecewise continuous across configuration-space cell interfaces $\overline{{\mathcal{X}}}_{h}^{p}\subset {\mathcal{X}}_{h}^{p}$, along with an additional subset  where continuity is required in the directions perpendicular to the magnetic field, but not in the direction parallel to the field. Assuming a field-aligned coordinate system (e.g. Beer et al. Reference Beer, Cowley and Hammett1995), we will take the perpendicular directions to be $x$ and $y$, and the parallel direction to be $z$.

Since we require $H_{h}$ to be continuous across all cell interfaces, this means that we require $\unicode[STIX]{x1D719}_{h}$ to be continuous, i.e. $\unicode[STIX]{x1D719}_{h}\in \overline{{\mathcal{X}}}_{h}^{p}$. Thus, to solve the Poisson equation we use the (continuous) finite-element method (FEM). While one could ensure $\unicode[STIX]{x1D719}_{h}$ is continuous in all directions by using a three-dimensional FEM solve, we instead use a two-dimensional FEM solve in the $x$ and $y$ directions, followed by a one-dimensional smoothing operation in the $z$ direction. That is, we first solve for  using a two-dimensional FEM solve, and then we use a smoothing/projection operation to ensure continuity in the $z$ direction. We will denote this operation as  and define it below. We can make this splitting because $\unicode[STIX]{x1D735}_{\bot }$ only produces coupling in the $x$ and $y$ (perpendicular) directions.

For the two-dimensional solve, we solve for  by multiplying equation (2.12) by a test function  and integrating (by parts) in each configuration-space cell ${\mathcal{K}}_{j}^{R}$ to obtain the discrete local weak form

(3.5)

where $\unicode[STIX]{x1D709}^{(j)}$ denotes the restriction of $\unicode[STIX]{x1D709}$ to cell $j$ and

(3.6)$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{g\,h}=\mathop{\sum }_{s}q_{s}\int _{{\mathcal{T}}^{v}}\text{d}\boldsymbol{w}{\mathcal{J}}f_{s\,h},\end{eqnarray}$$

with ${\mathcal{T}}^{v}$ the restriction of ${\mathcal{T}}$ to velocity space. The global weak form is then obtained by summing equation (

) over cells in $x$ and $y$ (but not in $z$), which results in cancellation of the surface terms at cell interfaces and leaves only a global $\unicode[STIX]{x2202}{\mathcal{T}}^{R}$ boundary term. Note that, in order to maintain energetic consistency (as we will see below), the introduction of ${\mathcal{P}}_{z}$ necessitates the modification of the right-hand side of (

) with ${\mathcal{P}}_{z}^{\ast }$, the adjoint of ${\mathcal{P}}_{z}$, defined as

(3.7)$$\begin{eqnarray}\int _{{\mathcal{T}}^{R}}\text{d}\boldsymbol{R}\,f{\mathcal{P}}_{z}[g]=\int _{{\mathcal{T}}^{R}}\text{d}\boldsymbol{R}\,{\mathcal{P}}_{z}^{\ast }[f]g.\end{eqnarray}$$

For the smoothing operation , we use a one-dimensional FEM solve in the $z$ direction. This can be written as the solution $\unicode[STIX]{x1D719}_{h}$ of the global (in $z$) weak equality

(3.8)

where $\unicode[STIX]{x1D712}\in \widehat{{\mathcal{X}}}_{h}^{p}\subset {\mathcal{X}}_{h}^{p}$, with $\widehat{{\mathcal{X}}}_{h}^{p}$ a subset of the configuration-space basis where continuity is required only in the $z$ direction. Here, ${\mathcal{T}}_{j}^{z}$ denotes a restriction of the domain that is global in $z$ but cell-wise local in $x$ and $y$. We remark that using an FEM solve for this operation makes ${\mathcal{P}}_{z}$ self-adjoint, so that ${\mathcal{P}}_{z}^{\ast }={\mathcal{P}}_{z}$. Note, however, that one could instead use a different, local smoothing operation that is not self-adjoint, so we will keep the distinction between ${\mathcal{P}}_{z}$ and ${\mathcal{P}}_{z}^{\ast }$. Also note that ${\mathcal{P}}_{z}$ is a projection operator, in that .

The continuous discrete Hamiltonian $H_{h}\in \overline{{\mathcal{V}}}_{h}^{p}$ is then given by

(3.9)

where $v_{\Vert h}^{2}$ is the projection of $v_{\Vert }^{2}$ onto $\overline{{\mathcal{V}}}_{h}^{p}$. Note that this is only necessary when $v_{\Vert }^{2}$ is not in the basis, i.e. when $p_{v}<2$, where $p_{v}$ is the maximum degree of the $v_{\Vert }$ monomials in the basis set.

For the parallel Ampère equation we will take  so that $A_{\Vert h}$ is continuous in $x$ and $y$ but discontinuous in $z$. Multiplying equation (2.13) by a test function  and integrating, we can obtain the discrete weak form of this equation. The local weak form in cell $j$ is

(3.10)$$\begin{eqnarray}\int _{{\mathcal{K}}_{j}^{R}}\text{d}\boldsymbol{R}\unicode[STIX]{x1D735}_{\bot }A_{\Vert h}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\bot }\unicode[STIX]{x1D711}^{(j)}-\oint _{\unicode[STIX]{x2202}{\mathcal{K}}_{j}^{R}}\text{d}\boldsymbol{s}_{R}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\bot }A_{\Vert h}\unicode[STIX]{x1D711}^{(j)}=\unicode[STIX]{x1D707}_{0}\int _{{\mathcal{K}}_{j}^{R}}\text{d}\boldsymbol{R}\,\unicode[STIX]{x1D711}^{(j)}J_{\Vert h},\end{eqnarray}$$

where again the surface terms will cancel on summing over cells except at the global $\unicode[STIX]{x2202}{\mathcal{T}}^{R}$ boundary, and

(3.11)$$\begin{eqnarray}J_{\Vert h}=\mathop{\sum }_{s}\frac{q_{s}}{m_{s}}\int _{{\mathcal{T}}^{v}}\text{d}\boldsymbol{w}{\mathcal{J}}\frac{\unicode[STIX]{x2202}H_{s\,h}}{\unicode[STIX]{x2202}v_{\Vert }}f_{s\,h}.\end{eqnarray}$$

Here, note that we have replaced the $v_{\Vert }$ in the $J_{\Vert }$ definition from (2.13) with $(1/m)\unicode[STIX]{x2202}H_{h}/\unicode[STIX]{x2202}v_{\Vert }$; this will be required for energy conservation in the $p_{v}=1$ case, since $\unicode[STIX]{x2202}H_{h}/\unicode[STIX]{x2202}v_{\Vert }\neq mv_{\Vert }$ when $v_{\Vert }^{2}$ is not in the basis. Instead, for $p_{v}=1$, $\unicode[STIX]{x2202}H_{h}/\unicode[STIX]{x2202}v_{\Vert }=m\bar{v}_{\Vert }$, the piecewise-constant projection of $mv_{\Vert }$. As before, we solve equation (3.10) using a two-dimensional FEM solve in the $x$ and $y$ directions. Note, however, that we do not require the smoothing operation in $z$ here because $A_{\Vert h}$ is allowed to be discontinuous in the $z$ direction, since it does not appear in the Hamiltonian in the symplectic formulation of EMGK.

The discrete weak form of Ohm’s law can be obtained by taking the time derivative of (3.10); after some manipulation, which we leave to appendix B, the local weak form becomes

(3.12)$$\begin{eqnarray}\displaystyle & & \displaystyle \!\!\!\!\int _{{\mathcal{K}}_{j}^{R}}\text{d}\boldsymbol{R}\unicode[STIX]{x1D735}_{\bot }\frac{\unicode[STIX]{x2202}A_{\Vert h}}{\unicode[STIX]{x2202}t}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\bot }\unicode[STIX]{x1D711}^{(j)}\nonumber\\ \displaystyle & & \displaystyle \qquad -\,\oint _{\unicode[STIX]{x2202}{\mathcal{K}}_{j}^{R}}\text{d}\boldsymbol{s}_{R}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\bot }\frac{\unicode[STIX]{x2202}A_{\Vert h}}{\unicode[STIX]{x2202}t}\unicode[STIX]{x1D711}^{(j)}-\int _{{\mathcal{K}}_{j}^{R}}\text{d}\boldsymbol{R}\,\unicode[STIX]{x1D711}^{(j)}\frac{\unicode[STIX]{x2202}A_{\Vert h}}{\unicode[STIX]{x2202}t}\mathop{\sum }_{s,i}\frac{\unicode[STIX]{x1D707}_{0}q_{s}^{2}}{m_{s}}\oint _{\unicode[STIX]{x2202}{\mathcal{K}}_{i}^{v}}\text{d}s_{w}\bar{v}_{\Vert }^{-}\widehat{{\mathcal{J}}f_{s\,h}}\nonumber\\ \displaystyle & & \displaystyle \quad =\unicode[STIX]{x1D707}_{0}\mathop{\sum }_{s}q_{s}\int _{{\mathcal{K}}_{j}^{R}}\text{d}\boldsymbol{R}\,\unicode[STIX]{x1D711}^{(j)}\left[\int _{{\mathcal{T}}^{v}}\text{d}\boldsymbol{w}\bar{v}_{\Vert }\frac{\unicode[STIX]{x2202}({\mathcal{J}}f_{s\,h})}{\unicode[STIX]{x2202}t}^{\star }\!-\!\mathop{\sum }_{i}\oint _{\unicode[STIX]{x2202}{\mathcal{K}}_{i}^{v}}\text{d}s_{w}\bar{v}_{\Vert }^{-}\dot{v}_{\Vert h}^{H}\widehat{{\mathcal{J}}f_{s\,h}}\right],\quad \!(p_{v}=1),\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(3.13)$$\begin{eqnarray}\displaystyle & & \displaystyle \!\!\!\!\int _{{\mathcal{K}}_{j}^{R}}\text{d}\boldsymbol{R}\unicode[STIX]{x1D735}_{\bot }\frac{\unicode[STIX]{x2202}A_{\Vert h}}{\unicode[STIX]{x2202}t}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\bot }\unicode[STIX]{x1D711}^{(j)}-\!\oint _{\unicode[STIX]{x2202}{\mathcal{K}}_{j}^{R}}\!\text{d}\boldsymbol{s}_{R}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\bot }\frac{\unicode[STIX]{x2202}A_{\Vert h}}{\unicode[STIX]{x2202}t}\unicode[STIX]{x1D711}^{(j)}+\!\int _{{\mathcal{K}}_{j}^{R}}\text{d}\boldsymbol{R}\,\unicode[STIX]{x1D711}^{(j)}\frac{\unicode[STIX]{x2202}A_{\Vert h}}{\unicode[STIX]{x2202}t}\mathop{\sum }_{s}\frac{\unicode[STIX]{x1D707}_{0}q_{s}^{2}}{m_{s}}\!\!\int _{{\mathcal{T}}^{v}}\!\text{d}\boldsymbol{w}{\mathcal{J}}f_{s\,h}\nonumber\\ \displaystyle & & \displaystyle \quad =\unicode[STIX]{x1D707}_{0}\mathop{\sum }_{s}q_{s}\int _{{\mathcal{K}}_{j}^{R}}\text{d}\boldsymbol{R}\,\unicode[STIX]{x1D711}^{(j)}\int _{{\mathcal{T}}^{v}}\text{d}\boldsymbol{w}v_{\Vert }\frac{\unicode[STIX]{x2202}({\mathcal{J}}f_{s\,h})}{\unicode[STIX]{x2202}t}^{\star },\quad (p_{v}>1),\end{eqnarray}$$

where , and

(3.14)$$\begin{eqnarray}\displaystyle & & \displaystyle \int _{{\mathcal{K}}_{j}}\text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}\unicode[STIX]{x1D713}\frac{\unicode[STIX]{x2202}({\mathcal{J}}f_{h})}{\unicode[STIX]{x2202}t}^{\star }=-\oint _{\unicode[STIX]{x2202}{\mathcal{K}}_{j}}\text{d}\boldsymbol{w}\,\text{d}\boldsymbol{s}_{R}\boldsymbol{\cdot }\dot{\boldsymbol{R}}_{h}\unicode[STIX]{x1D713}^{-}\widehat{{\mathcal{J}}f_{h}}+\int _{{\mathcal{K}}_{j}}\text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}{\mathcal{J}}f_{h}\dot{\boldsymbol{R}}_{h}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\int _{{\mathcal{K}}_{j}}\text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}{\mathcal{J}}f_{h}\dot{v}_{\Vert h}^{H}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}v_{\Vert }}+\int _{{\mathcal{K}}_{j}}\text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}\unicode[STIX]{x1D713}({\mathcal{J}}C[f_{h}]+{\mathcal{J}}S),\end{eqnarray}$$

so that the gyrokinetic equation can be written as

(3.15)$$\begin{eqnarray}\displaystyle & & \displaystyle \int _{{\mathcal{K}}_{j}}\text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}\unicode[STIX]{x1D713}\frac{\unicode[STIX]{x2202}({\mathcal{J}}f_{h})}{\unicode[STIX]{x2202}t}\nonumber\\ \displaystyle & & \displaystyle \quad =\int _{{\mathcal{K}}_{j}}\text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}\unicode[STIX]{x1D713}\frac{\unicode[STIX]{x2202}({\mathcal{J}}f_{h})}{\unicode[STIX]{x2202}t}^{\star }-\oint _{\unicode[STIX]{x2202}{\mathcal{K}}_{j}}\text{d}\boldsymbol{R}\,\text{d}s_{w}\left(\dot{v}_{\Vert h}^{H}-\frac{q}{m}\frac{\unicode[STIX]{x2202}A_{\Vert h}}{\unicode[STIX]{x2202}t}\right)\unicode[STIX]{x1D713}^{-}\widehat{{\mathcal{J}}f_{h}}\nonumber\\ \displaystyle & & \displaystyle \qquad -\,\int _{{\mathcal{K}}_{j}}\text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}{\mathcal{J}}f_{h}\frac{q}{m}\frac{\unicode[STIX]{x2202}A_{\Vert h}}{\unicode[STIX]{x2202}t}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}v_{\Vert }}.\end{eqnarray}$$

Note that some special attention is required to ensure that upwinding of the numerical fluxes is handled consistently in (3.12) and (3.15) in the $p_{v}=1$ case. The upwind flow for the $v_{\Vert }$ surface terms is $\dot{v}_{\Vert h}^{H}-(q/m)\unicode[STIX]{x2202}A_{\Vert h}/\unicode[STIX]{x2202}t$; this is somewhat problematic because we cannot readily solve for $\unicode[STIX]{x2202}A_{\Vert h}/\unicode[STIX]{x2202}t$ from (3.12) without first knowing the upwind direction, which depends on $\unicode[STIX]{x2202}A_{\Vert h}/\unicode[STIX]{x2202}t$. Thus, for $p_{v}=1$ only, we use an approximate $\widetilde{\unicode[STIX]{x2202}A_{\Vert h}/\unicode[STIX]{x2202}t}$, calculated using equation (3.13) (which contains no surface term contributions), to compute the upwind direction for the $v_{\Vert }$ surface terms in (3.12) and (3.15). (One could extend this algorithm by iterating with a new estimate of the upwind direction based on the previous estimate of $\unicode[STIX]{x2202}A_{\Vert \,h}/\unicode[STIX]{x2202}t$, but we leave that for future work. The present algorithm seems to work well for the cases tested so far, and we expect that $\widetilde{\unicode[STIX]{x2202}A_{\Vert h}/\unicode[STIX]{x2202}t}$ results in the correct upwind direction most of the time.)

In our modal DG scheme, integrals in the above weak forms are computed analytically using a quadrature-free scheme that results in exact integrations (of the discrete integrands). (This means there are no aliasing errors, and that integration by parts operations that led to these integrals are treated exactly, for the specified discrete representation of $f_{h}$ and other factors in the integrand.) This is important for ensuring the conservation properties of the scheme, since the conservation laws in the EMGK system are indirect, involving integrals of the gyrokinetic equation. The fact that integrations are exact also has important implications for the cancellation problem. Since integrals in the discrete Ohm’s law are computed exactly, the discretization errors (which are solely embedded in the discrete integrands) cancel exactly, avoiding the cancellation problem. For more details about the modal scheme, the analytical integrations and the avoidance of the cancellation problem, we have included in appendix C a derivation of a semi-discrete Alfvén wave dispersion relation that results from our scheme.

3.2 Discrete conservation properties

Now we would like to show that the discrete system (in the continuous-time limit) preserves various conservation laws of the continuous system. As with the continuous system, we will consider the conservation properties in the absence of collisions, sources and sinks, and we will assume that the boundary conditions are either periodic or that the distribution function vanishes at the boundary.

Proof. Taking $\unicode[STIX]{x1D713}=1$ in the discrete weak form of the gyrokinetic equation, equation (3.2), and summing over all cells, we have

(3.16)$$\begin{eqnarray}\displaystyle & & \displaystyle \mathop{\sum }_{j}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\int _{{\mathcal{K}}_{j}}\text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}{\mathcal{J}}f_{h}+\mathop{\sum }_{j}\oint _{\unicode[STIX]{x2202}{\mathcal{K}}_{j}}\text{d}\boldsymbol{w}\,\text{d}\boldsymbol{s}_{R}\boldsymbol{\cdot }\dot{\boldsymbol{R}}_{h}\widehat{{\mathcal{J}}f_{h}}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\mathop{\sum }_{j}\oint _{\unicode[STIX]{x2202}{\mathcal{K}}_{j}}\text{d}\boldsymbol{R}\,\text{d}s_{w}\left(\dot{v}_{\Vert h}^{H}-\frac{q}{m}\frac{\unicode[STIX]{x2202}A_{\Vert h}}{\unicode[STIX]{x2202}t}\right)\widehat{{\mathcal{J}}f_{h}}=0\nonumber\\ \displaystyle & & \displaystyle \quad \Rightarrow \quad \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\int _{{\mathcal{T}}}\text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}{\mathcal{J}}f_{h}=0,\end{eqnarray}$$

where the surface terms cancel exactly at cell interfaces because the integrands (both the phase-space characteristics and the numerical fluxes) are continuous across the interfaces. ◻

Proposition 2. The discrete system conserves a discrete total energy, $W_{h}=W_{H\,h}-W_{E\,h}+W_{B\,h}$, where

(3.17)

(3.18)

and

(3.19)$$\begin{eqnarray}W_{B\,h}=\int _{{\mathcal{T}}}\,\text{d}\boldsymbol{R}\frac{1}{2\unicode[STIX]{x1D707}_{0}}|\unicode[STIX]{x1D735}_{\bot }A_{\Vert h}|^{2}.\end{eqnarray}$$

Proof. The proof follows from Proposition 3.2 in Hakim et al. (Reference Hakim, Hammett, Shi and Mandell2019). We start by calculating

(3.20)$$\begin{eqnarray}\frac{\text{d}W_{H\,h}}{\text{d}t}=\mathop{\sum }_{s,j}\int _{{\mathcal{K}}_{j}}\,\text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}H_{s\,h}\frac{\unicode[STIX]{x2202}({\mathcal{J}}f_{s\,h})}{\unicode[STIX]{x2202}t}+{\mathcal{J}}f_{s\,h}\frac{\unicode[STIX]{x2202}H_{s\,h}}{\unicode[STIX]{x2202}t}.\end{eqnarray}$$

The first term can be calculated by taking $\unicode[STIX]{x1D713}=H_{h}$ in (3.2) and summing over cells and species, since $\unicode[STIX]{x1D713}\in {\mathcal{V}}_{h}^{p}$ and $H_{h}\in \overline{{\mathcal{V}}}_{h}^{p}\subset {\mathcal{V}}_{h}^{p}$:

(3.21)$$\begin{eqnarray}\displaystyle & & \displaystyle \mathop{\sum }_{s,j}\int _{{\mathcal{K}}_{j}}\text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}H_{s\,h}\frac{\unicode[STIX]{x2202}({\mathcal{J}}f_{s\,h})}{\unicode[STIX]{x2202}t}+\mathop{\sum }_{s,j}\oint _{\unicode[STIX]{x2202}{\mathcal{K}}_{j}}\text{d}\boldsymbol{w}\,\text{d}\boldsymbol{s}_{R}\boldsymbol{\cdot }\dot{\boldsymbol{R}}_{h}H_{s\,h}^{-}\widehat{{\mathcal{J}}f_{s\,h}}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\mathop{\sum }_{s,j}\oint _{\unicode[STIX]{x2202}{\mathcal{K}}_{j}}\text{d}\boldsymbol{R}\,\text{d}s_{w}\left(\dot{v}_{\Vert h}^{H}-\frac{q_{s}}{m_{s}}\frac{\unicode[STIX]{x2202}A_{\Vert h}}{\unicode[STIX]{x2202}t}\right)H_{s\,h}^{-}\widehat{{\mathcal{J}}f_{s\,h}}\nonumber\\ \displaystyle & & \displaystyle \quad -\,\mathop{\sum }_{s,j}\int _{{\mathcal{K}}_{j}}\text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}{\mathcal{J}}f_{s\,h}\left(\dot{\boldsymbol{R}}_{h}\boldsymbol{\cdot }\unicode[STIX]{x1D735}H_{s\,h}+\dot{v}_{\Vert h}^{H}\frac{\unicode[STIX]{x2202}H_{s\,h}}{\unicode[STIX]{x2202}v_{\Vert }}\right)\nonumber\\ \displaystyle & & \displaystyle \quad +\,\mathop{\sum }_{s,j}\int _{{\mathcal{K}}_{j}}\text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}{\mathcal{J}}f_{s\,h}\frac{q_{s}}{m_{s}}\frac{\unicode[STIX]{x2202}A_{\Vert h}}{\unicode[STIX]{x2202}t}\frac{\unicode[STIX]{x2202}H_{s\,h}}{\unicode[STIX]{x2202}v_{\Vert }}=0.\end{eqnarray}$$

Here, we see why we must require $H_{h}$ to be continuous; we want the surface terms to vanish, which means the integrands must be continuous across cell interfaces so that the contributions from either side of the interface cancel exactly when we sum over cells. The numerical flux $\widehat{{\mathcal{J}}f_{h}}$ is by definition continuous across the interface, and we have already noted above that the phase-space characteristics $\dot{\boldsymbol{R}}_{h}$ and $\dot{v}_{\Vert h}^{H}-(q/m)\unicode[STIX]{x2202}A_{\Vert h}/\unicode[STIX]{x2202}t$ are also continuous across cell interfaces. This leaves the Hamiltonian, which we require to be continuous so that the surface terms do indeed vanish. Further, the first volume term vanishes exactly because $\dot{\boldsymbol{R}}_{h}\boldsymbol{\cdot }\unicode[STIX]{x1D735}H_{h}+\dot{v}_{\Vert h}^{H}\unicode[STIX]{x2202}H_{h}/\unicode[STIX]{x2202}v_{\Vert }=\{H_{h},H_{h}\}=0$ by definition of the Poisson bracket. This leaves

(3.22)$$\begin{eqnarray}\mathop{\sum }_{s,j}\int _{{\mathcal{K}}_{j}}\text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}H_{s\,h}\frac{\unicode[STIX]{x2202}({\mathcal{J}}f_{s\,h})}{\unicode[STIX]{x2202}t}=-\mathop{\sum }_{s,j}\int _{{\mathcal{K}}_{j}}\text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}{\mathcal{J}}f_{s\,h}\frac{q_{s}}{m_{s}}\frac{\unicode[STIX]{x2202}A_{\Vert h}}{\unicode[STIX]{x2202}t}\frac{\unicode[STIX]{x2202}H_{s\,h}}{\unicode[STIX]{x2202}v_{\Vert }}=-\int _{{\mathcal{T}}^{R}}\text{d}\boldsymbol{R}\frac{\unicode[STIX]{x2202}A_{\Vert h}}{\unicode[STIX]{x2202}t}J_{\Vert h},\end{eqnarray}$$

where, here, we see why we have defined $J_{\Vert h}$ using the derivative of $H_{h}$ instead of $v_{\Vert }$, as noted after equation (3.11). We now have the desired result for this term. For the second term in (3.20), we have

(3.23)

Thus we have

(3.24)

which is consistent with equation (2.24).

Next, we calculate

(3.25)

where we have used  in (3.5) to make the second equality, noting that the surface term vanishes upon summing over cells because  is continuous in the perpendicular directions. Here, we see why we modified the right-hand side of (3.5) with ${\mathcal{P}}_{z}^{\ast }$, so that the resulting term in (

3.25

) matches the one in (3.23).

Finally, we calculate

(3.26)$$\begin{eqnarray}\frac{\text{d}W_{B\,h}}{\text{d}t}=\mathop{\sum }_{j}\int _{{\mathcal{K}}_{j}^{R}}\text{d}\boldsymbol{R}\frac{1}{\unicode[STIX]{x1D707}_{0}}\unicode[STIX]{x1D735}_{\bot }A_{\Vert h}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\bot }\frac{\unicode[STIX]{x2202}A_{\Vert h}}{\unicode[STIX]{x2202}t}=\int _{{\mathcal{T}}^{R}}\text{d}\boldsymbol{R}\frac{\unicode[STIX]{x2202}A_{\Vert h}}{\unicode[STIX]{x2202}t}J_{\Vert h},\end{eqnarray}$$

where we have used $\unicode[STIX]{x1D711}^{(j)}=(1/\unicode[STIX]{x1D707}_{0})\unicode[STIX]{x2202}A_{\Vert h}/\unicode[STIX]{x2202}t$ in (3.10) to make the second equality, again noting that the surface term vanishes upon summing over cells because  is continuous in the perpendicular directions.

We now have conservation of discrete total energy

(3.27)$$\begin{eqnarray}\frac{\text{d}W_{h}}{\text{d}t}=\frac{\text{d}W_{H\,h}}{\text{d}t}-\frac{\text{d}W_{E\,h}}{\text{d}t}+\frac{\text{d}W_{B\,h}}{\text{d}t}=0.\end{eqnarray}$$

We note that this proof did not rely on the particular choice of numerical flux function. ◻

Proposition 4. If the discrete distribution function $f_{h}$ remains positive definite, then the discrete scheme grows the discrete entropy monotonically,

(3.28)$$\begin{eqnarray}-\frac{\text{d}}{\text{d}t}\int _{{\mathcal{T}}}\text{d}\boldsymbol{R}\,\text{d}\boldsymbol{w}{\mathcal{J}}f_{h}\ln (f_{h})\geqslant 0.\end{eqnarray}$$

5.1 Kinetic Alfvén wave

As a first benchmark of our electromagnetic scheme, we consider the kinetic Alfvén wave. In a slab (straight background magnetic field) geometry, with stationary ions (assuming $\unicode[STIX]{x1D714}\gg k_{\Vert }v_{ti}$), the gyrokinetic equation for electrons reduces to

(5.1)$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}f_{e}}{\unicode[STIX]{x2202}t}=\{H_{e},f_{e}\}-\frac{e}{m}\frac{\unicode[STIX]{x2202}f_{e}}{\unicode[STIX]{x2202}v_{\Vert }}\frac{\unicode[STIX]{x2202}A_{\Vert }}{\unicode[STIX]{x2202}t}=-v_{\Vert }\frac{\unicode[STIX]{x2202}f_{e}}{\unicode[STIX]{x2202}z}-\frac{e}{m}\frac{\unicode[STIX]{x2202}f_{e}}{\unicode[STIX]{x2202}v_{\Vert }}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}z}+\frac{\unicode[STIX]{x2202}A_{\Vert }}{\unicode[STIX]{x2202}t}\right).\end{eqnarray}$$

Taking a single Fourier mode with perpendicular wavenumber $k_{\bot }$ and parallel wavenumber $k_{\Vert }$, the field equations become

(5.2)$$\begin{eqnarray}\displaystyle & \displaystyle k_{\bot }^{2}\frac{m_{i}n_{0}}{B^{2}}\unicode[STIX]{x1D719}=en_{0}-e\int \text{d}v_{\Vert }\,f_{e}, & \displaystyle\end{eqnarray}$$

(5.3)$$\begin{eqnarray}\displaystyle & \displaystyle k_{\bot }^{2}A_{\Vert }=-\unicode[STIX]{x1D707}_{0}e\int \text{d}v_{\Vert }\,v_{\Vert }f_{e}, & \displaystyle\end{eqnarray}$$

(5.4)$$\begin{eqnarray}\displaystyle & \displaystyle \left(k_{\bot }^{2}+\frac{\unicode[STIX]{x1D707}_{0}e^{2}}{m_{e}}\int \text{d}v_{\Vert }\,f_{e}\right)\frac{\unicode[STIX]{x2202}A_{\Vert }}{\unicode[STIX]{x2202}t}=-\unicode[STIX]{x1D707}_{0}e\int \text{d}v_{\Vert }~v_{\Vert }\{H_{e},f_{e}\}. & \displaystyle\end{eqnarray}$$

After linearizing the gyrokinetic equation by assuming a uniform Maxwellian background with density $n_{0}$ and temperature $m_{e}v_{te}^{2}$, so that $f_{e}=F_{Me}+\unicode[STIX]{x1D6FF}f_{e}$, the dispersion relation becomes

(5.5)$$\begin{eqnarray}\unicode[STIX]{x1D714}^{2}\left[1+\frac{\unicode[STIX]{x1D714}}{\sqrt{2}k_{\Vert }v_{te}}Z\left(\frac{\unicode[STIX]{x1D714}}{\sqrt{2}k_{\Vert }v_{te}}\right)\right]=\frac{k_{\Vert }^{2}v_{te}^{2}}{\hat{\unicode[STIX]{x1D6FD}}}\left[1+k_{\bot }^{2}\unicode[STIX]{x1D70C}_{s}^{2}+\frac{\unicode[STIX]{x1D714}}{\sqrt{2}k_{\Vert }v_{te}}Z\left(\frac{\unicode[STIX]{x1D714}}{\sqrt{2}k_{\Vert }v_{te}}\right)\right],\end{eqnarray}$$

where $\hat{\unicode[STIX]{x1D6FD}}=(\unicode[STIX]{x1D6FD}_{e}/2)m_{i}/m_{e}$, with $\unicode[STIX]{x1D6FD}_{e}=2\unicode[STIX]{x1D707}_{0}n_{0}T_{e}/B^{2}$, $v_{te}=\sqrt{T_{e}/m_{e}}$ is the electron thermal speed, $\unicode[STIX]{x1D70C}_{s}$ is the ion sound gyroradius and $Z(x)$ is the plasma dispersion function (Fried & Conte Reference Fried and Conte1961). In the limit $k_{\bot }\unicode[STIX]{x1D70C}_{s}\ll 1$ the wave becomes the standard shear Alfvén wave from magnetohydrodyamics (MHD), which is an undamped wave with frequency $\unicode[STIX]{x1D714}=k_{\Vert }v_{A}$, where $v_{A}=v_{te}/\hat{\unicode[STIX]{x1D6FD}}^{1/2}$ is the Alfvén velocity. For larger values of $k_{\bot }\unicode[STIX]{x1D70C}_{s}$, the mode is damped by kinetic effects.

Figure 1. Real frequencies (a) and damping rates (b) for the kinetic Alfvén wave vs $k_{\bot }\unicode[STIX]{x1D70C}_{s}$. Solid lines are the exact values from (5.5) for three different values of $\hat{\unicode[STIX]{x1D6FD}}=(\unicode[STIX]{x1D6FD}_{e}/2)m_{i}/m_{e}$, and black dots are the numerical results from Gkeyll.

Figure 2. Values of $\unicode[STIX]{x1D719}_{h}$ (blue) and $\unicode[STIX]{x2202}A_{\Vert h}/\unicode[STIX]{x2202}t$ (yellow) for the case with $\hat{\unicode[STIX]{x1D6FD}}=10$ and $k_{\bot }\unicode[STIX]{x1D70C}_{s}=0.01$. The amplitude of $E_{\Vert h}$ (green) is ${\sim}10^{-9}$.

In figure 1, we show the real frequencies (a) and damping rates (b) obtained by solving equation (5.5) for a few values of $\hat{\unicode[STIX]{x1D6FD}}$. We also show numerical results from Gkeyll, which match the analytic results very well. These results are a good indication that our scheme avoids the Ampère cancellation problem, which can cause large errors for modes with $\hat{\unicode[STIX]{x1D6FD}}/k_{\bot }^{2}\unicode[STIX]{x1D70C}_{s}^{2}\gg 1$ (see appendix A); we see no such errors, even for the case with $\hat{\unicode[STIX]{x1D6FD}}/k_{\bot }^{2}\unicode[STIX]{x1D70C}_{s}^{2}=10^{5}$. Each Gkeyll simulation was run using piecewise-linear basis functions ($p=1$) in a reduced dimensionality mode with one configuration-space dimension and one velocity-space dimension, with $(N_{z},N_{v_{\Vert }})=(32,64)$ the number of cells in each dimension. The perpendicular dimensions ($x$ and $y$), which appear only in the field equations in this simple system, were handled by replacing $\unicode[STIX]{x1D735}_{\bot }$ by $k_{\bot }$, as in (5.2) and (5.3). We use periodic boundary conditions in $z$ and zero-flux boundary conditions in $v_{\Vert }$.

We also show in figure 2 the fields $\unicode[STIX]{x1D719}_{h}$ and $\unicode[STIX]{x2202}A_{\Vert h}/\unicode[STIX]{x2202}t$ for the case with $\hat{\unicode[STIX]{x1D6FD}}=10$ and $k_{\bot }\unicode[STIX]{x1D70C}_{s}=0.01$, which gives $\hat{\unicode[STIX]{x1D6FD}}/k_{\bot }^{2}\unicode[STIX]{x1D70C}_{s}^{2}=10^{5}$. For these parameters the system is near the MHD limit, which means we should expect $E_{\Vert }=-\unicode[STIX]{x2202}\unicode[STIX]{x1D719}/\unicode[STIX]{x2202}z-\unicode[STIX]{x2202}A_{\Vert }/\unicode[STIX]{x2202}t\approx 0$. While this condition is never enforced, getting the physics correct requires the scheme to allow $\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{h}/\unicode[STIX]{x2202}z\approx -\unicode[STIX]{x2202}A_{\Vert h}/\unicode[STIX]{x2202}t$. The fact that our scheme allows discontinuities in $A_{\Vert }$ in the parallel direction is an advantage in this case. Because $\unicode[STIX]{x1D719}_{h}$ is piecewise linear here, $\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{h}/\unicode[STIX]{x2202}z$ is piecewise constant; this is necessarily discontinuous for non-trivial solutions. Thus the scheme produces a piecewise constant $\unicode[STIX]{x2202}A_{\Vert h}/\unicode[STIX]{x2202}t$ in this MHD-limit case, as shown in figure 2, resulting in $E_{\Vert h}\approx 0$. If our scheme did not allow discontinuities in $A_{\Vert h}$, a continuous $\unicode[STIX]{x2202}A_{\Vert h}/\unicode[STIX]{x2202}t$ would never be able to exactly cancel a discontinuous $\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{h}/\unicode[STIX]{x2202}z$, and the resulting $E_{\Vert h}\neq 0$ would make the solution inaccurate. Notably, this would be the case had we chosen the Hamiltonian ($p_{\Vert }$) formulation of the gyrokinetic system, which uses $p_{\Vert }=mv_{\Vert }+qA_{\Vert }$ as the parallel velocity coordinate. This is because $A_{\Vert }$ is included in the Hamiltonian in the $p_{\Vert }$ formulation, which would require continuity of $A_{\Vert h}$ (and thereby $\unicode[STIX]{x2202}A_{\Vert h}/\unicode[STIX]{x2202}t$) to conserve energy in our discretization scheme.

5.2 Kinetic ballooning mode

We use the kinetic ballooning mode (KBM) instability in the local limit as a second linear benchmark of our electromagnetic scheme. The dispersion relation is given by solving (Kim, Horton & Dong Reference Kim, Horton and Dong1993)

(5.6)$$\begin{eqnarray}\unicode[STIX]{x1D714}[\unicode[STIX]{x1D70F}+k_{\bot }^{2}+\unicode[STIX]{x1D6E4}_{0}(b)-P_{0}]\unicode[STIX]{x1D719}=[\unicode[STIX]{x1D70F}(\unicode[STIX]{x1D714}-\unicode[STIX]{x1D714}_{\ast e})-k_{\Vert }P_{1}]\frac{\unicode[STIX]{x1D714}}{k_{\Vert }}A_{\Vert },\end{eqnarray}$$

(5.7)$$\begin{eqnarray}\displaystyle \frac{2k_{\Vert }^{2}k_{\bot }^{2}}{\unicode[STIX]{x1D6FD}_{i}}A_{\Vert } & = & \displaystyle k_{\Vert }[k_{\Vert }P_{1}-\unicode[STIX]{x1D70F}(\unicode[STIX]{x1D714}-\unicode[STIX]{x1D714}_{\ast e})]\unicode[STIX]{x1D719}\nonumber\\ \displaystyle & & \displaystyle -\,[k_{\Vert }^{2}P_{2}-\unicode[STIX]{x1D70F}(\unicode[STIX]{x1D714}(\unicode[STIX]{x1D714}-\unicode[STIX]{x1D714}_{\ast e})-2\unicode[STIX]{x1D714}_{de}(\unicode[STIX]{x1D714}-\unicode[STIX]{x1D714}_{\ast e}(1+\unicode[STIX]{x1D702}_{e})))]A_{\Vert },\end{eqnarray}$$

where

(5.8)$$\begin{eqnarray}P_{m}=\int _{0}^{\infty }\text{d}v_{\bot }v_{\bot }\int _{-\infty }^{\infty }\text{d}v_{\Vert }\frac{1}{\sqrt{2\unicode[STIX]{x03C0}}}\text{e}^{-(v_{\Vert }^{2}+v_{\bot }^{2})/2}(v_{\Vert })^{m}\frac{\unicode[STIX]{x1D714}-\unicode[STIX]{x1D714}_{\ast i}[1+\unicode[STIX]{x1D702}_{i}(v^{2}/2-3/2)]}{\unicode[STIX]{x1D714}-k_{\Vert }v_{\Vert }-\unicode[STIX]{x1D714}_{di}(v_{\Vert }^{2}+v_{\bot }^{2}/2)}J_{0}^{2}(v_{\bot }\sqrt{b}),\end{eqnarray}$$

with $\unicode[STIX]{x1D70F}=T_{i}/T_{e}$, $\unicode[STIX]{x1D714}_{\ast e}=k_{y}$, $\unicode[STIX]{x1D714}_{\ast i}=-k_{y}$, $\unicode[STIX]{x1D702}_{s}=L_{n}/L_{Ts}$ and $\unicode[STIX]{x1D6E4}_{0}(b)=I_{0}(b)e^{-b}$ with $I_{0}(b)=J_{0}(ib)$ the modified Bessel function. Here, the wavenumbers $k_{y}$ and $k_{\Vert }$ are normalized to $\unicode[STIX]{x1D70C}_{i}$ and $L_{n}$, respectively, and the frequencies $\unicode[STIX]{x1D714}$ and $\unicode[STIX]{x1D714}_{\ast }$ are normalized to $v_{ti}/L_{n}$. In the local limit, $\unicode[STIX]{x1D714}_{ds}=\unicode[STIX]{x1D714}_{\ast s}L_{n}/R$ and $k_{\bot }=k_{y}$ do not vary along the field line. Note that in (5.6) we have modified the FLR terms from (Kim et al. Reference Kim, Horton and Dong1993) so that we can take $b=k_{\bot }^{2}\rightarrow 0$ while keeping $k_{\bot }\neq 0$ to neglect all FLR effects except for the first-order polarization term, which is consistent with our long-wavelength Poisson equation.

The local limit can be achieved by simulating a helical flux tube with no magnetic shear, which gives a system with constant magnetic curvature that corresponds to $\unicode[STIX]{x1D714}_{d}=\text{const}$. This geometry has been previously used for SOL turbulence studies with Gkeyll (Bernard et al. Reference Bernard, Shi, Gentle, Hakim, Hammett, Stoltzfus-Dueck and Taylor2019; Shi et al. Reference Shi, Hammett, Stoltzfus-Dueck and Hakim2019), except in this section we take the boundary condition along the field lines to be periodic. We will provide further details about the helical geometry and the coordinates in the next section.

Figure 3. Growth rates for the KBM instability in the local limit, as a function of $\unicode[STIX]{x1D6FD}_{i}$, with $k_{\bot }\unicode[STIX]{x1D70C}_{i}=0.5,k_{\Vert }L_{n}=0.1,R/L_{n}=5,R/L_{Ti}=12.5,R/L_{Te}=10$ and $\unicode[STIX]{x1D70F}=1$. The black dots are numerical results from Gkeyll, and the coloured lines are the result of numerically solving the analytic dispersion relation given by (5.6)–(5.7).

We show the results of Gkeyll simulations of the KBM instability in the local-limit helical geometry for several values of $\unicode[STIX]{x1D6FD}_{i}$ in figure 3. The results agree well with the analytic result obtained by numerically solving equations (5.6)–(5.7). The parameters $k_{\bot }\unicode[STIX]{x1D70C}_{i}=0.5,k_{\Vert }L_{n}=0.1,R/L_{n}=5,R/L_{Ti}=12.5,R/L_{Te}=10,\unicode[STIX]{x1D70F}=1$ are chosen to match those used in figure 1 of Kim et al. (Reference Kim, Horton and Dong1993), although the differences in FLR terms ($b=0$) cause our growth rates to be larger than those in Kim et al. (Reference Kim, Horton and Dong1993). These are fully five-dimensional simulations with the real deuterium–electron mass ratio using piecewise-linear ($p=1$) basis functions, with $(N_{x},N_{y},N_{z},N_{v_{\Vert }},N_{\unicode[STIX]{x1D707}})=(1,16,16,32,16)$. The boundary conditions are periodic in the three configuration-space dimensions and zero flux in the velocity dimensions. The initial distribution function of each species is composed of a background Maxwellian with gradients in the density and temperature corresponding to the desired $L_{n}$ and $L_{Ts}$, plus a perturbed Maxwellian (for the electrons only) with small sinusoidal variations in the density corresponding to the desired $k_{y}$ and $k_{\Vert }$. Note that, since we are using a full-$f$ representation, the presence of a background gradient in the distribution function means that we must apply the periodic boundary conditions by first subtracting off the initial background distribution function, then applying periodicity to the perturbations only and then adding back the background distribution. Finally, we note that since Gkeyll is designed primarily for nonlinear calculations, the fact that Fourier modes are not eigenfunctions of the DG discretization of the system makes these linear tests somewhat difficult for Gkeyll. This may play a role in the small deviation of the results from the analytical theory. Because of this, Fourier modes other than the one initialized can grow and pollute the results. In particular, we have not included results from the ion temperature gradient (ITG) branch because we find that a mode with $k_{\Vert }=0$ grows and overcomes the initialized finite $k_{\Vert }$ mode before its growth rate has converged.