Electromagnetic full-$f$ gyrokinetics in the tokamak edge with discontinuous Galerkin methods

Abstract

We present an energy-conserving discontinuous Galerkin scheme for the full-$f$ electromagnetic gyrokinetic system in the long-wavelength limit. We use the symplectic formulation and solve directly for $\unicode[STIX]{x2202}A_{\Vert }/\unicode[STIX]{x2202}t$, the inductive component of the parallel electric field, using a generalized Ohm’s law derived directly from the gyrokinetic equation. Linear benchmarks are performed to verify the implementation and show that the scheme avoids the Ampère cancellation problem. We perform a nonlinear electromagnetic simulation in a helical open-field-line system as a rough model of the tokamak scrape-off layer using parameters from the National Spherical Torus Experiment (NSTX). This is the first published nonlinear electromagnetic gyrokinetic simulation on open field lines. Comparisons are made to a corresponding electrostatic simulation.

1 Introduction

Understanding turbulent transport physics in the tokamak edge and scrape-off layer (SOL) is critical to developing a successful fusion reactor. The dynamics in these regions plays a key role in determining the L–H transition, the pedestal height and the heat load to the vessel walls. While the edge is often modelled by Braginskii-type fluid models that have provided valuable results and insights (Xu et al. 2008; Tamain et al. 2010; Ricci et al. 2012; Francisquez, Zhu & Rogers 2017; Zhu, Francisquez & Rogers 2017), a kinetic treatment will inevitably be necessary for reliable quantitative predictions in some cases (Jenko & Dorland 2001; Cohen & Xu 2008). Gyrokinetic theory and direct numerical simulation have become important tools for studying turbulence and transport in fusion plasmas, especially in the core region (Parker, Lee & Santoro 1993a; Kotschenreuther, Rewoldt & Tang 1995; Dimits et al. 2000; Dorland et al. 2000; Jenko 2000; Lin et al. 2000; Jost et al. 2001; Candy & Waltz 2003; Idomura, Tokuda & Kishimoto 2003; Watanabe & Sugama 2005; Jolliet et al. 2007; Idomura et al. 2008; Peeters et al. 2009; Lanti et al. 2019). In the edge and SOL, gyrokinetic simulations are particularly challenging because the large, intermittent fluctuations in the SOL make assumptions of scale separation between equilibrium and fluctuations not strongly valid. This necessitates a full-$f$ approach that self-consistently evolves the full distribution function, $f$ (as opposed to the $\unicode[STIX]{x1D6FF}f$ approach commonly used in the core, where one assumes $f=F_{0}+\unicode[STIX]{x1D6FF}f$ with a fixed background $F_{0}$ so that only $\unicode[STIX]{x1D6FF}f$ perturbations must be evolved, and the parallel electric field nonlinearity is frequently neglected). Steady progress in gyrokinetic edge/SOL modelling has been made with both particle-in-cell (PIC) (Ku, Chang & Diamond 2009; Korpilo et al. 2016; Ku et al. 2016) and continuum (Shi et al. 2019; Dorf et al. 2016; Pan et al. 2018) methods. Another challenge is the magnetic geometry of the edge/SOL region, which requires treatment of open and closed magnetic field-line regions and the resulting plasma interactions with material walls on open field lines. The X-point in a diverted geometry is an additional complication which makes the use of field-aligned coordinates challenging. Currently, only the XGC1 hybrid-Lagrangian PIC code (Ku et al. 2016) can simulate gyrokinetic turbulence in a three-dimensional diverted geometry with an X-point.

The edge/SOL region also features steep pressure gradients, especially in the H-mode transport barrier and SOL regions, which contributes to the importance of electromagnetic effects. In this regime, the parallel electron dynamics is no longer fast relative to the drift turbulence, so electrons can no longer be treated adiabatically (Scott 1997). This leads to coupling of the perpendicular vortex motions and kinetic shear Alfvén waves, which results in field-line bending (Xu et al. 2010). Including electromagnetic effects in gyrokinetic simulations has proved numerically and computationally challenging, both in the core and in the edge. The so-called Ampère cancellation problem is one of the main numerical issues that has troubled primarily PIC codes (Reynders 1993; Cummings 1994). Various $\unicode[STIX]{x1D6FF}f$ PIC schemes to address the cancellation problem have been developed and there are interesting recent advances in this area (Chen & Parker 2003; Mishchenko, Hatzky & Könies 2004; Hatzky, Könies & Mishchenko 2007; Mishchenko et al. 2014; Startsev & Lee 2014; Bao, Lin & Lu 2018). Meanwhile, some continuum $\unicode[STIX]{x1D6FF}f$ core codes avoided the cancellation problem completely (Rewoldt, Tang & Hastie 1987; Kotschenreuther et al. 1995), while others had to address somewhat minor issues resulting from it (Jenko 2000; Candy & Waltz 2003). With respect to the cancellation problem, one possible reason for the differences might be that, in continuum codes, the fields and particles are discretized on the same grid, whereas in PIC codes the particle positions do not coincide with the field grid. Because particle positions are randomly located relative to the field grid, one might need to be more careful in some way when treating the interaction of the particles and electromagnetic fields.

To this point, all published nonlinear electromagnetic gyrokinetic results have focused on the core region, mostly within the $\unicode[STIX]{x1D6FF}f$ formulation neglecting the $E_{\Vert }$ nonlinearity, although the ORB5 PIC code includes the $E_{\Vert }$ nonlinearity and is effectively full-$f$ (Lanti et al. 2019). The XGC1 code is also full-$f$ and is focused on both the core and the edge/SOL; it has an option for a gyrokinetic ion/drift-fluid massless electron hybrid model (Hager et al. 2017), with a fully kinetic implicit electromagnetic scheme based on Chen & Chacon (2015) recently implemented and under further development (Ku et al. 2018b). Other gyrokinetic codes working on the SOL are not yet electromagnetic. Thus, to our knowledge, the results presented here are the first nonlinear electromagnetic full-$f$ gyrokinetic turbulence simulations on open field lines.

In this paper we present a numerical scheme for simulating the full-$f$ electromagnetic gyrokinetic system using a continuum approach. We use an energy-conserving discontinuous Galerkin (DG) scheme for the discretization of the gyrokinetic system in phase space, building on the work of Liu & Shu (2000), Shi (2017), Shi et al. (2017) and Hakim et al. (2019). DG methods are attractive because they are highly local (enabling fairly straightforward parallelization schemes), allow high-order accuracy and enforce local conservation laws (Durran 2010). The present target of the scheme is simulating the edge and SOL of tokamaks, although the scheme could in principle be used for whole-device modelling, including the core. Our scheme has been implemented as part of the gyrokinetics solver (Shi et al. 2017, 2019; Bernard et al. 2019) of the Gkeyll computational plasma framework, which also includes solvers for the Vlasov–Maxwell system (Cagas et al. 2017; Juno et al. 2018) and multi-moment fluid equations (Wang et al. 2015).

The paper is organized as follows. In § 2, we describe the electromagnetic gyrokinetic system and some of its conservation properties. Section 3 describes the discontinuous Galerkin phase-space discretization of the system, and also presents proofs that the scheme preserves particle and energy conservation. The time-discretization scheme is handled in § 4. In § 5 we present some linear electromagnetic benchmarks that validate the scheme and also demonstrate the avoidance of the cancellation problem. We present nonlinear results showing the first electromagnetic gyrokinetic turbulence simulation on open field lines in § 6, along with comparisons to a corresponding electrostatic simulation. We summarize and address future work in § 7.

2 The electromagnetic gyrokinetic system

2.1 Basic equations

We solve the full-$f$ electromagnetic gyrokinetic (EMGK) equation in the symplectic formulation (Brizard & Hahm 2007), 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 by¹

(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 (2007); 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 1995), 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 1988; Brizard & Hahm 2007). 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. 2018a; Shi et al. 2019). 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 1993; Cummings 1994; Chen & Parker 2001)

(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 1993; Chen & Parker 2001), 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 1958; Dougherty 1964),

(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. 2020). 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. 2008).

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 The discrete EMGK system

In this section we describe the phase-space discretization of the electromagnetic gyrokinetic system used in Gkeyll.

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 (2000) (originally for the two-dimensional incompressible Euler and Navier–Stokes equations) to arbitrary Hamiltonian systems (Hakim et al. 2019; Shi 2017; Shi et al. 2017). However, unlike the nodal approach used in Shi (2017) and Shi et al. (2017), 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 2011). 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 2000; Shi 2017; Shi et al. 2017; Hakim et al. 2019). 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.²

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. 1995), 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 (

3.5

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

3.5

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

Proposition 1. The discrete system conserves total number of particles (the $L_{1}$ norm).

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. (2019). 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 3. The discrete system exactly conserves the $L_{2}$ norm of the distribution function when using a central flux, while the distribution function $L_{2}$ norm monotonically decays when using an upwind flux.

Proof. The proof is given as Proposition 3.3 in Hakim et al. (2019).◻

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

Proof. The proof is given as Proposition 3.4 in Hakim et al. (2019).◻

4 Time-discretization scheme

So far we have considered only the discretization of the phase space for the system, and we have considered the conservation properties of the scheme in the continuous-time limit. Indeed, in the discrete-time system the conservation properties are no longer exact due to truncation error in the non-reversible time-stepping methods that we consider. However, the errors will be independent of the phase-space discretization, and errors can be reduced by taking a smaller time step or by using a high-order time-stepping scheme to improve convergence. Following the approach of the Runge–Kutta discontinuous Galerkin method (Cockburn & Shu 1998, 2001; Shu 2009), we have implemented several explicit multi-stage strong stability-preserving Runge–Kutta high-order schemes (Gottlieb, Shu & Tadmor 2001; Shu 2002); the results in this paper use a three-stage, third-order scheme (SSP-RK3). These schemes have the property that a high-order scheme can be composed of several forward-Euler stages. Thus we will detail our time-stepping scheme for a single forward-Euler stage, which can then be combined into a multi-stage high-order scheme. Note that, although we present the time-discretization scheme in this section in terms of our DG phase-space discretization, the scheme could be generalized to any spatial discretization.

Given $f_{h}^{n}=f_{h}(t=t_{n})$ and $A_{\Vert h}^{n}=A_{\Vert h}(t=t_{n})$ at time $t_{n}$, the steps of the forward-Euler scheme to advance to time $t_{n+1}=t_{n}+\unicode[STIX]{x0394}t$ are as follows:

1. (i) Calculate  using equation (3.5), and then  using equation (3.8).

   (4.1)

   
   (4.2)

   

2. (ii) Calculate the partial EMGK update $(\unicode[STIX]{x2202}({\mathcal{J}}f_{h})^{\star }/\unicode[STIX]{x2202}t)^{n}$ using equation (3.14).

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

3. (iii) Calculate $(\unicode[STIX]{x2202}A_{\Vert h}/\unicode[STIX]{x2202}t)^{n}$ from (3.13) [for $p_{v}=1$, this is only a provisional value, which we will denote as $(\widetilde{\unicode[STIX]{x2202}A_{\Vert h}/\unicode[STIX]{x2202}t})^{n}$].

   (4.4)$$\begin{eqnarray}\displaystyle & & \displaystyle \int _{{\mathcal{K}}_{j}^{R}}\text{d}\boldsymbol{R}\unicode[STIX]{x1D735}_{\bot }\left(\frac{\unicode[STIX]{x2202}A_{\Vert h}}{\unicode[STIX]{x2202}t}\right)^{n}\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 }\left(\frac{\unicode[STIX]{x2202}A_{\Vert h}}{\unicode[STIX]{x2202}t}\right)^{n}\unicode[STIX]{x1D711}^{(j)}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\int _{{\mathcal{K}}_{j}^{R}}\text{d}\boldsymbol{R}\,\unicode[STIX]{x1D711}^{(j)}\left(\frac{\unicode[STIX]{x2202}A_{\Vert h}}{\unicode[STIX]{x2202}t}\right)^{n}\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}^{n}\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 }\left(\frac{\unicode[STIX]{x2202}({\mathcal{J}}f_{s\,h})}{\unicode[STIX]{x2202}t}^{\star }\right)^{n}.\end{eqnarray}$$

4. (iv) ($p_{v}=1$ only) Use the provisional $(\widetilde{\unicode[STIX]{x2202}A_{\Vert h}/\unicode[STIX]{x2202}t})^{n}$ from step 3 to calculate the upwinding direction in the surface terms in (3.12), and then calculate $(\unicode[STIX]{x2202}A_{\Vert h}/\unicode[STIX]{x2202}t)^{n}$.

   (4.5)$$\begin{eqnarray}\displaystyle & & \displaystyle \int _{{\mathcal{K}}_{j}^{R}}\text{d}\boldsymbol{R}\unicode[STIX]{x1D735}_{\bot }\left(\frac{\unicode[STIX]{x2202}A_{\Vert h}}{\unicode[STIX]{x2202}t}\right)^{n}\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 }\left(\frac{\unicode[STIX]{x2202}A_{\Vert h}}{\unicode[STIX]{x2202}t}\right)^{n}\unicode[STIX]{x1D711}^{(j)}\nonumber\\ \displaystyle & & \displaystyle \qquad -\,\int _{{\mathcal{K}}_{j}^{R}}\text{d}\boldsymbol{R}\,\unicode[STIX]{x1D711}^{(j)}\left(\frac{\unicode[STIX]{x2202}A_{\Vert h}}{\unicode[STIX]{x2202}t}\right)^{n}\mathop{\sum }_{s}\frac{\unicode[STIX]{x1D707}_{0}q_{s}^{2}}{m_{s}}\mathop{\sum }_{i}\oint _{\unicode[STIX]{x2202}{\mathcal{K}}_{i}^{v}}\text{d}s_{w}\bar{v}_{\Vert }^{-}\widehat{{\mathcal{J}}f_{s\,h}}^{n}\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 }\left(\frac{\unicode[STIX]{x2202}({\mathcal{J}}f_{s\,h})}{\unicode[STIX]{x2202}t}^{\star }\right)^{n}-\mathop{\sum }_{i}\oint _{\unicode[STIX]{x2202}{\mathcal{K}}_{i}^{v}}\text{d}s_{w}\bar{v}_{\Vert }^{-}\dot{v}_{\Vert h}^{H\,n}\widehat{{\mathcal{J}}f_{s\,h}}^{n}\right].\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

5. (v) Calculate the full GK update, $((\unicode[STIX]{x2202}({\mathcal{J}}f_{h}))/\unicode[STIX]{x2202}t)^{n}$, using equation (3.15). For $p_{v}=1$, the provisional $(\widetilde{(\unicode[STIX]{x2202}A_{\Vert h})/\unicode[STIX]{x2202}t})^{n}$ from step 3 should again be used to calculate the upwinding direction in the surface terms for consistency.

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

6. (vi) Advance $f_{h}$ and $A_{\Vert h}$ to time $t_{n+1}$.

   (4.7)$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{J}}f_{h}^{n+1}={\mathcal{J}}f_{h}^{n}+\unicode[STIX]{x0394}t\left(\frac{\unicode[STIX]{x2202}({\mathcal{J}}f_{h})}{\unicode[STIX]{x2202}t}\right)^{n}, & \displaystyle\end{eqnarray}$$
   (4.8)$$\begin{eqnarray}\displaystyle & \displaystyle A_{\Vert h}^{n+1}=A_{\Vert h}^{n}+\unicode[STIX]{x0394}t\left(\frac{\unicode[STIX]{x2202}A_{\Vert h}}{\unicode[STIX]{x2202}t}\right)^{n}. & \displaystyle\end{eqnarray}$$

Note that the parallel Ampère equation, equation (3.10), is only used to solve for the initial condition of $A_{\Vert h}(t=0)$. For all other times, equation (4.8) is used to advance $A_{\Vert h}$. This prevents the system from being over-determined and ensures consistency between $A_{\Vert h}$ and $\unicode[STIX]{x2202}A_{\Vert h}/\unicode[STIX]{x2202}t$.

5 Linear benchmarks

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 1961). 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 1993)

(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. 1993) 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. 2019; Shi et al. 2019), 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. (1993), although the differences in FLR terms ($b=0$) cause our growth rates to be larger than those in Kim et al. (1993). 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.

6 Nonlinear results

We now present preliminary nonlinear electromagnetic results from Gkeyll. We simulate turbulence on helical, open field lines as a rough model of the tokamak scrape-off layer. These simulations are a direct extension of the work of Shi et al. (2019) to include electromagnetic fluctuations. As such, we use the same simulation geometry and similar NSTX-like parameters. In the non-orthogonal field-aligned geometry, $x$ is the radial coordinate, $z$ is the coordinate along the field lines and $y$ is the binormal coordinate which labels field lines at constant $x$ and $z$. These coordinates map to physical cylindrical coordinates $(R,\unicode[STIX]{x1D711},Z)$ via $R=x$, $\unicode[STIX]{x1D711}=(y/\sin \unicode[STIX]{x1D703}+z\cos \unicode[STIX]{x1D703})/R_{c}$, $Z=z\sin \unicode[STIX]{x1D703}$. (Note that this fixes an error in the $\unicode[STIX]{x1D711}(y,z)$ mapping in (Shi et al. 2019).) The field-line pitch $\sin \unicode[STIX]{x1D703}=B_{v}/B$ is taken to be constant, with $B_{v}$ the vertical component of the magnetic field (analogous to the poloidal field in typical tokamak geometry), and $B$ the total magnitude of the background magnetic field. Further, $R_{c}=R_{0}+a$ is the radius of curvature at the centre of the simulation domain, with $R_{0}$ the device major radius and $a$ the minor radius. As in Shi et al. (2019), we neglect all geometrical factors arising from the non-orthogonal coordinate system, except for the assumption that perpendicular gradients are much stronger than parallel gradients so that we can approximate

(6.1)$$\begin{eqnarray}(\unicode[STIX]{x1D735}\times \hat{\boldsymbol{b}})\boldsymbol{\cdot }\unicode[STIX]{x1D735}f(x,y,z)\approx [(\unicode[STIX]{x1D735}\times \hat{\boldsymbol{b}})\boldsymbol{\cdot }\unicode[STIX]{x1D735}y]\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}y}=-\frac{1}{x}\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}y},\end{eqnarray}$$

where we have used $\boldsymbol{B}=B_{axis}(R_{0}/R)\hat{\boldsymbol{e}}_{z}$ in the last step, with $B_{axis}$ the magnetic field strength at the magnetic axis.

We use this geometry to simulate a flux-tube-like domain on the outboard side that wraps helically around the torus and terminates on conducting plates at each end in $z$. The simulation box is centred at $(x,y,z)=(R_{c},0,0)$ with dimensions $L_{x}=50\unicode[STIX]{x1D70C}_{s0}\approx 14.6~\text{cm}$, $L_{y}=100\unicode[STIX]{x1D70C}_{s0}\approx 29.1~\text{cm}$ and $L_{z}=L_{p}/\sin \unicode[STIX]{x1D703}=8~\text{m}$, where $L_{p}=2.4~\text{m}$ and $\unicode[STIX]{x1D70C}_{s0}=c_{s0}/\unicode[STIX]{x1D6FA}_{i}$. Periodic boundary conditions are used in the $y$ direction, and a Dirichlet boundary condition $\unicode[STIX]{x1D719}=0$ is applied in $x$, which effectively prevents flows into the boundaries in $x$. Conducting-sheath boundary conditions are applied in the $z$ direction (Shi et al. 2017, 2019), which partially reflect one species (typically electrons) and fully absorbs the other species depending on the sign of the sheath potential. This involves solving the gyrokinetic Poisson equation to evaluate the potential at the $z$ boundary, corresponding to the sheath entrance, and using the resulting sheath potential to determine a cutoff velocity below which particles are reflected by the sheath. Notably, our sheath boundary condition allows current fluctuations in and out of the sheath, which we will find to important later in this section. This is different from the standard logical sheath boundary condition (Parker et al. 1993b) that imposes that there is no net current to the sheath by assuming that the ion and electron currents at the sheath entrance are equal at all times. The velocity-space grid has extents $4v_{ts}\leqslant v_{\Vert }\leqslant 4v_{ts}$ and $0\leqslant \unicode[STIX]{x1D707}\leqslant 6T_{s0}/B_{0}$, where $v_{ts}=\sqrt{T_{s0}/m_{s}}$ and $B_{0}=B_{axis}R_{0}/R_{c}$. We use piecewise-linear ($p=1$) basis functions, with $(N_{x},N_{y},N_{z},N_{v_{\Vert }},N_{\unicode[STIX]{x1D707}})=(16,32,10,10,5)$. Note that, although the domain that we simulate is a flux tube, the simulations are not performed in the local limit; the simulations include radial variation of the magnetic field and the profiles, and are thus effectively global.

The simulation parameters are similar to those used in Shi et al. (2019), roughly approximating an H-mode deuterium plasma in the NSTX SOL: $B_{axis}=0.5~\text{T}$, $R_{0}=0.85~\text{m}$, $a=0.5~\text{m}$, $T_{e0}=T_{i0}=40~\text{eV}$. For the particle source, we use the same form as in Shi et al. (2019) but we increase the source particle rate by a factor of 10 to access a higher $\unicode[STIX]{x1D6FD}$ regime where electromagnetic effects will be more important. The source is localized in the region $x<x_{S}+3\unicode[STIX]{x1D706}_{S}$, with $x_{S}=R_{c}-0.05~\text{m}$ and $\unicode[STIX]{x1D706}_{S}=5\times 10^{-3}~\text{m}$. The location $x=x_{S}+3\unicode[STIX]{x1D706}_{S}$, which separates the source region from the SOL region, can be thought of as the separatrix. A floor of one tenth the peak particle source rate is used near the midplane to prevent regions of $n\ll n_{0}$ from developing at large $x$. Unlike in Shi et al. (2019) we do not use numerical heating to keep $f>0$ despite the fact that our DG algorithm does not guarantee positivity. While the simulations appear to be robust to negative $f$ in some isolated regions, lowering the source floor in the SOL region leads to simulation failures due to positivity issues at large $x$. A more sophisticated algorithm for ensuring positivity is left to future work. We also artificially lower the collision frequency to one tenth the physical value to offset the increased particle source rate so that the time-step limit from collisions does not become too restrictive. Further, we model only ion–ion and electron–electron collisions, leaving cross-species collisions to future work.

Figure 4. Snapshots from an electromagnetic simulation on open, helical field lines. From left to right, we show the density, temperature and plasma $\unicode[STIX]{x1D6FD}$ of electrons (a–c) and ions (d–f). The snapshots are taken at the midplane $(z=0)$ at $t=620~\unicode[STIX]{x03BC}\text{s}$. The dashed line indicates the boundary between the source and SOL regions. A blob with mushroom structure is being ejected from the source region and propagating radially outward into the SOL region.

Figure 5. Snapshots (at $z=0$, $t=620~\unicode[STIX]{x03BC}\text{s}$) of the electrostatic potential $\unicode[STIX]{x1D719}$, parallel magnetic vector potential $A_{\Vert }$ and normalized magnetic fluctuation amplitude $|\unicode[STIX]{x1D6FF}B_{\bot }|/B_{0}=|\unicode[STIX]{x1D735}_{\bot }A_{\Vert }|/B_{0}$ (a–c), along with the components of the parallel electric field $E_{\Vert }=-\unicode[STIX]{x1D735}_{\Vert }\unicode[STIX]{x1D719}-\unicode[STIX]{x2202}A_{\Vert }/\unicode[STIX]{x2202}t$ (d–f).

Figure 6. Radial profile of the normalized magnetic fluctuation amplitude, $|\unicode[STIX]{x1D6FF}B_{\bot }|/B_{0}=|\unicode[STIX]{x1D735}_{\bot }A_{\Vert }|/B_{0}$, averaged in $y$, $z$ and time using data near the midplane $(|z|<0.4~\text{m})$ over a period of $400~\unicode[STIX]{x03BC}\text{s}$. On average, we observe magnetic fluctuations of the order of $0.5{-}1\,\%$. The source region is shaded.

Using the novel electromagnetic scheme described in this paper, we ran a simulation in this configuration to $t=1~\text{ms}$, with a quasi-steady state being reached around $t=600~\unicode[STIX]{x03BC}\text{s}$ when the sources balance losses to the end plates. For reference, the ion transit time is $\unicode[STIX]{x1D70F}_{i}=(L_{z}/2)/v_{ti}\approx 50~\unicode[STIX]{x03BC}\text{s}$. In figure 4 we show snapshots of the density, temperature and $\unicode[STIX]{x1D6FD}$ of electrons (a–c) and ions (d–f). The snapshots are taken at the midplane ($z=0$) at $t=620~\unicode[STIX]{x03BC}\text{s}$. We can see a blob with a mushroom structure being ejected from the source region. We also show in figure 5 snapshots of the electromagnetic fields taken at the same time and location. We show the electrostatic potential $\unicode[STIX]{x1D719}$, the parallel magnetic vector potential $A_{\Vert }$ and the normalized magnetic fluctuation amplitude $|\unicode[STIX]{x1D6FF}B_{\bot }|/B_{0}=|\unicode[STIX]{x1D735}_{\bot }A_{\Vert }|/B_{0}$ (a–c), along with the components of the parallel electric field $E_{\Vert }=-\unicode[STIX]{x1D735}_{\Vert }\unicode[STIX]{x1D719}-\unicode[STIX]{x2202}A_{\Vert }/\unicode[STIX]{x2202}t$ (d–f). Note that only $\unicode[STIX]{x1D719}$, $A_{\Vert }$ and $\unicode[STIX]{x2202}A_{\Vert }/\unicode[STIX]{x2202}t$ are evolved quantities in the simulation, with the other quantities derived. We see that $\unicode[STIX]{x2202}A_{\Vert }/\unicode[STIX]{x2202}t$ is of comparable magnitude to $\unicode[STIX]{x1D735}_{\Vert }\unicode[STIX]{x1D719}$, indicating that the dynamics is in the electromagnetic regime. Significant magnetic fluctuations of over $2.5\,\%$ can be seen in $|\unicode[STIX]{x1D6FF}B_{\bot }|/B_{0}$ in this snapshot. We also show in figure 6 the time and spatially averaged profile of magnetic fluctuations vs $x$, which shows that, on average, we observe magnetic fluctuations of the order of $0.5\,\%{-}1\,\%$. This radial profile, and similar ones that will follow, is computed by averaging in $y$ and $z$ using data near the midplane $(|z|<0.4~\text{m})$ over the period of $600~\unicode[STIX]{x03BC}\text{s}{-}1~\text{ms}$.

In figures 7 and 8 we show projections of the three-dimensional magnetic field-line trajectories. These plots are created by integrating the field-line equations for the total (background plus fluctuation) magnetic field. In figure 7, each field line starts at $z=-4~\text{m}$ and either $x=1.33~\text{m}$ or $x=1.38~\text{m}$ for a range of $y$ values and is traced to $z=4~\text{m}$. The starting points (at $z=-4~\text{m}$) are marked with circles, while the ending points (at $z=4~\text{m}$) are marked with crosses. The trajectories have been projected onto the $x-y$ plane, and we have also plotted the ion density at $z=0~\text{m}$ in the background. From left to right, we show a short time series of snapshots, with $t=230$, 240 and $250~\unicode[STIX]{x03BC}\text{s}$. At $t=230~\unicode[STIX]{x03BC}\text{s}$, a blob is starting to emerge from the source region at $y\approx 0.04~\text{m}$. The field lines that start at $x=1.33~\text{m}$ are beginning to be stretched radially outward as the blob emerges. In the $t=240~\unicode[STIX]{x03BC}\text{s}$ snapshot, we see that the blob is now propagating radially outward into the SOL region and the $x=1.33~\text{m}$ field lines have been stretched further. The field lines that start at $x=1.38~\text{m}$ are now also starting to be stretched near $y\approx 0.2~\text{m}$, and they are stretched even more in the $t=250~\unicode[STIX]{x03BC}\text{s}$ snapshot as the blob continues to propagate. We can also see the remnants of another blob that was ejected near $y=-0.1~\text{m}$ in previous frames. In the $t=230~\unicode[STIX]{x03BC}\text{s}$ snapshot, the field lines have been stretched by this blob, but by $t=250~\unicode[STIX]{x03BC}\text{s}$ the field lines in this region have returned closer to their equilibrium position. This behaviour of blobs bending and stretching the field lines is an inherently full-$f$ phenomenon. The blobs have a higher density and temperature than the background, so they raise the local plasma $\unicode[STIX]{x1D6FD}$ as they propagate. This causes the field lines to move with the plasma, allowing the fields lines to be deformed and stretched by the radially propagating blobs and ultimately leading to larger magnetic fluctuations.

Figure 7. Three-dimensional magnetic field-line trajectories at $t=230$, 240 and $250~\unicode[STIX]{x03BC}\text{s}$, projected onto the $x-y$ plane. The ion density at $z=0~\text{m}$ is plotted in the background. Each field line starts at $z=-4~\text{m}$ and either $x=1.33~\text{m}$ or $x=1.38~\text{m}$ for a range of $y$ values and is traced to $z=4~\text{m}$. The starting points are marked with circles and the ending points are marked with crosses. Focusing on the blob that is being ejected near $y=0~\text{m}$, we see that field lines are stretched and bent by the blob as it propagates into the SOL region. In previous frames (not shown) a blob was also ejected near $y=-0.1~\text{m}$. At $t=230~\unicode[STIX]{x03BC}\text{s}$ the field lines are still stretched from this event, but they return closer to their equilibrium position by $t=250~\unicode[STIX]{x03BC}\text{s}$. (A full animation of this time series is included as supplementary materials online at https://doi.org/10.1017/S0022377820000070.)

Figure 8. Three-dimensional magnetic field-line trajectories at $t=240~\unicode[STIX]{x03BC}\text{s}$, projected onto the $x-y$ plane in (a) and the $x-z$ plane in (b). The ion density is plotted in the background, at $z=0~\text{m}$ in (a) and averaged over $|y|<0.02~\text{m}$ in (b). Each field line starts at $y=0~\text{m}$ and $z=-4~\text{m}$ for a range of $x$ values and is traced to $z=4~\text{m}$. The starting points are marked with circles and the ending points are marked with crosses. Each field line is coloured the same in both (a,b). The field lines in the near SOL are stretched radially outward by a blob near $y=0~\text{m}$.

In figure 8 we show a slightly different view of the field-line trajectories at $t=240~\unicode[STIX]{x03BC}\text{s}$. Field lines are still traced from the bottom ($z=-4~\text{m}$) to the top ($z=4~\text{m}$), but now each field line starts at $y=0~\text{m}$ for a range of $x$. The starting points are again marked with circles and the ending points are marked with crosses. We have projected the three-dimensional trajectories onto the $x-y$ plane in figure 8(a), and onto the $x-z$ plane in figure 8(b). In (a) we again plot the ion density at $z=0~\text{m}$ in the background; in (b) the ion density has been averaged over $|y|<0.02~\text{m}$. As can be seen in figure 8(b), the blob propagating near $y\approx 0~\text{m}$ has stretched several field lines radially outward near the midplane. These bowed-out field lines originate from a range of $x$ values, $1.3~\text{m}\lesssim x\lesssim 1.35~\text{m}$, and have all been dragged along with the blob as it was ejected from the source region and propagated radially outward. We also see some degree of line tying in these plots, with many of the field lines ending at a similar point in $x-y$ space to where they began, despite being stretched near the midplane. The field lines are not perfectly line tied, however; if they were, the crosses would perfectly align with their corresponding circles in the $x-y$ projections. Because our sheath boundary condition allows current fluctuations at the sheath interface, we can model the finite resistance of the sheath, which makes line tying only partial (Kunkel & Guillory 1966). This allows the footpoints of the field lines to slip at the sheath interface (Ryutov 2006). Examining figures 7 and 8, we see evidence of this in the simulation, with most of the end points moving slowly and smoothly in the vicinity of their origin, especially at larger $x$. In the source region, however, there are other field lines whose end points suddenly jump further away from their origin. This suggests that we are seeing line breaking (reconnection) due to electron inertia effects and numerical diffusion, as field lines are pushed close together by large perturbations in the source region.

We have also run a corresponding electrostatic simulation in this configuration for direct comparison. This simulation is identical in configuration to the $L_{z}=8~\text{m}$ case from Shi et al. (2019) except for the increased particle source rate and lack of cross-species collisions. In figure 9 we show a comparison of radial profiles of density, temperature and $\unicode[STIX]{x1D6FD}$ for the electromagnetic and electrostatic cases. The profiles for the electromagnetic case are shallower in the SOL region and steeper in the source region. This suggests that there is less radial particle and heat transport into the SOL region in the electromagnetic case. This is in part confirmed by the profile of the radial particle flux in figure 10, showing approximately 40 % less particle transport in the electromagnetic case. The total particle flux $\unicode[STIX]{x1D6E4}_{n,r}$ includes the $E\times B$ particle flux, $\unicode[STIX]{x1D6E4}_{n,r,E\times B}=\langle {\tilde{n}}_{e}\tilde{v}_{r}\rangle$, with $v_{r}=E_{r}/B=-(1/B)\unicode[STIX]{x2202}\unicode[STIX]{x1D719}/\unicode[STIX]{x2202}y$. In the electromagnetic case, $\unicode[STIX]{x1D6E4}_{n,r}$ also includes the particle flux due to magnetic flutter, $\unicode[STIX]{x1D6E4}_{n,r,\text{flutter}}=\langle \widetilde{n_{e}u_{\Vert e}}b_{r}\rangle$, with $b_{r}=(1/B)\unicode[STIX]{x2202}A_{\Vert }/\unicode[STIX]{x2202}y$. The tilde indicates the fluctuation of a time-varying quantity, defined as $\tilde{A}=A-\bar{A}$ with $\bar{A}$ the time average of $A$. The brackets $\langle A\rangle$ denote an average in $y$, $z$ (near the midplane) and time. We also note that the electromagnetic profiles might be even shallower in the SOL region if not for the floor on the source used to prevent positivity issues in the distribution function at large $x$.

Figure 9. Radial profiles of density (a), temperature (b) and $\unicode[STIX]{x1D6FD}$ (c) for electrons (solid) and ions (dashed). Profiles from the electromagnetic case (EM) are blue, and the electrostatic profiles (ES) are yellow. The profiles are averaged in $y$, $z$ and time using data near the midplane $(|z|<0.4~\text{m})$ over a period of $400~\unicode[STIX]{x03BC}\text{s}$. The electromagnetic case shows shallower profiles in the SOL region, indicating that there is less radial particle and heat transport (as confirmed by figure 10).

Figure 10. Radial profile of the radial electron particle flux $\unicode[STIX]{x1D6E4}_{n,r}$, averaged in $y$, $z$ and time using data near the midplane $(|z|<0.4~\text{m})$ over a period of $400~\unicode[STIX]{x03BC}\text{s}$. The transport in the electromagnetic case (EM, blue) is approximately 40 % lower than in the electrostatic case (ES, yellow). This contributes to the shallower electron density profile in the electromagnetic case, as seen in figure 9. The electromagnetic case contains radial transport from both $E\times B$ drift (dashed) and magnetic flutter (dot-dashed).

Figure 11. Comparison of fluctuation statistics for the electron density (a–c), electrostatic potential (d–f) and radial electron particle flux (g,h) between the electromagnetic case (EM, blue) and a corresponding electrostatic case (ES, yellow). From left to right, we show radial profiles of the normalized root-mean-square fluctuation amplitude, skewness and excess kurtosis. These plots were computed by averaging in $y$, $z$ and time using data near the midplane $(|z|<0.4~\text{m})$ over a period of $400~\unicode[STIX]{x03BC}\text{s}$. The electromagnetic case shows higher electron density fluctuation amplitude, skewness and excess kurtosis. This is an indication that the electromagnetic case has more intermittent density fluctuations. The skewness and kurtosis of the particle flux also indicate that the transport is more intermittent in the electromagnetic case. The statistics for the potential are relatively similar for the electrostatic and electromagnetic cases.

We also compare fluctuation statistics between the electromagnetic and electrostatic cases in figure 11. Statistics of the electron density are shown on the top row, the middle row shows statistics of electrostatic potential fluctuations and the bottom row shows statistics of the radial electron particle flux. Despite the fact that the electromagnetic case shows less particle transport, the root-mean-square relative density fluctuations are larger in the electromagnetic case by up to a factor of two. The electromagnetic case also has higher skewness and excess kurtosis, indicating that the density fluctuations in the electromagnetic case are more intermittent. The skewness and kurtosis of the particle flux also indicate more intermittency in the electromagnetic case. This contributes to the reduced transport shown in the electromagnetic case, since the transport events are rarer even if the fluctuation levels are larger. Meanwhile, the fluctuation statistics for the potential are relatively similar between the electromagnetic and electrostatic cases. The statistics of the density and potential appear more coupled in the electrostatic case, consistent with the electrons behaving adiabatically when electromagnetic effects are neglected. In the electromagnetic case the density fluctuations are more intermittent and higher amplitude than the potential fluctuations.

Finally, we note that in terms of computational cost, the electromagnetic simulation is less than twice as expensive as the corresponding electrostatic simulation on a per-time-step basis. On 128 cores, the time per time step was 0.41 s/step for the electrostatic simulation and 0.68 s/step for the electromagnetic simulation. The increased cost is due to the additional field solves required for Ohm’s law, along with additional terms in the gyrokinetic equation. However, due to time-step restrictions on an electrostatic simulation due to the electrostatic shear Alfvén mode (also known as the $\unicode[STIX]{x1D714}_{H}$ mode) (Lee 1987), the electromagnetic simulation makes up some of the additional cost by taking slightly larger time steps. The total wall-clock time (on 128 cores) for the electrostatic simulation was approximately 65 h, and the electromagnetic simulation took about 82 h. Altogether, the cost of these simulations is relatively modest, and the addition of electromagnetic effects only makes the simulations marginally (${\sim}25$ %) more expensive. We also note that the new version of Gkeyll, which uses a quadrature-free modal DG scheme, is approximately 10 times faster than the previous version of Gkeyll used in Shi et al. (2019), which used nodal DG with Gaussian quadrature. More details about the improvements from the quadrature-free modal scheme will be reported elsewhere.

7 Summary and conclusion

In this paper we have presented an energy-conserving scheme for the full-$f$ electromagnetic gyrokinetic system. We choose the symplectic formulation of EMGK, which uses the parallel velocity as an independent variable. This leads to the time derivative of the parallel vector potential, $\unicode[STIX]{x2202}A_{\Vert }/\unicode[STIX]{x2202}t$, appearing explicitly in the gyrokinetic equation. We handle this term directly by solving an Ohm’s law. We presented the conservation properties of the EMGK system.

We described the discontinuous Galerkin scheme used to discretize the EMGK system in phase space. We proved that the scheme preserves particle conservation, and that the scheme also preserves energy conservation provided that the discrete Hamiltonian is continuous. This is achieved by using the (continuous) finite-element method for the field solves. We also detailed a basic forward-Euler time-stepping scheme to be used in the stages of a multi-stage high-order SSP-RK scheme. The time-stepping scheme updates the gyrokinetic equation in two parts, with the result of the first part [denoted $\unicode[STIX]{x2202}({\mathcal{J}}f_{s})^{\star }/\unicode[STIX]{x2202}t$] being used in Ohm’s law to solve directly for $\unicode[STIX]{x2202}A_{\Vert }/\unicode[STIX]{x2202}t$ so that it can then be used in the second part of the gyrokinetic update.

We have implemented the scheme in the gyrokinetic solver of the Gkeyll computational plasma framework. We provided two linear benchmarks to validate the electromagnetic scheme: a kinetic Alfvén wave calculation and a local kinetic ballooning mode instability calculation. In both cases results from Gkeyll agree well with analytic results. The success of these calculations, especially in cases with $\hat{\unicode[STIX]{x1D6FD}}/k_{\bot }^{2}\unicode[STIX]{x1D70C}_{s}^{2}\gg 1$, indicates that the scheme avoids the Ampère cancellation problem. Further, the discontinuous Galerkin nature of the scheme enables (but does not enforce) exact cancellation of the components of $E_{\Vert }$, allowing the system to capture the MHD limit with $E_{\Vert }=0$.

We presented a nonlinear electromagnetic full-$f$ gyrokinetic simulation of turbulence on helical, open field lines as a rough model of the tokamak scrape-off layer. This simulation is the first nonlinear electromagnetic gyrokinetic simulation on open field lines. We showed data illustrating the interplay between blobs propagating into the SOL and the resulting bending and stretching of magnetic field lines. We also made quantitative comparisons between the electromagnetic simulation and a corresponding electrostatic simulation. Notably, the electromagnetic simulation exhibits less transport and shallower density and temperature profiles in the SOL, as well as larger root-mean-square density fluctuations with more intermittency. Future work will examine the particle and energy balance in these nonlinear simulations to confirm the conservation properties proved in § 3.2 after accounting for sources and losses to the walls.

A number of extensions have been left to future work. The capability to simulate more realistic magnetic geometries using a non-orthogonal field-aligned coordinate system is currently in progress, so that effects of magnetic shear and non-constant curvature can be included. The inclusion of both open and closed-field-line regions, including the X-point in diverted geometries, is an additional complication that must be addressed. Gyroaveraging is another important effect that must be implemented to improve fidelity. Further, a robust solution to the issue of maintaining positivity of the distribution function has been implemented for Hamiltonian terms and is in progress for the collision operator. This could, for example, alleviate the need to use a source floor in the nonlinear simulations presented in § 6, which could further enhance the differences between the electromagnetic and electrostatic profiles.

The modest cost of the nonlinear full-$f$ gyrokinetic simulations that we have presented make the prospect of using Gkeyll for whole-device modelling a feasible goal. The inclusion of electromagnetic effects will be crucial to the fidelity of these efforts. Such a tool would be invaluable to future studies of turbulent transport in fusion devices, both from a theoretical perspective and also as a model for predicting and analysing the performance of current and future experiments.