$\dot{\mbox{\boldmath${R}$}}=\dot{x}\,\mbox{\boldmath${e}$}_{x}+\dot{y}\,\mbox{\boldmath${e}$}_{y}+\dot{z}\,\mbox{\boldmath${e}$}_{z}$,
with components

$\displaystyle\dot{x}$
$\displaystyle=\\{x,H\\}=\dot{\mbox{\boldmath${R}$}}\boldsymbol{\cdot}\nabla
x=\frac{B^{*x}}{B_{\parallel}^{*}}v_{\parallel}+\frac{\mathbf{\hat{b}}}{qB_{\parallel}^{*}}\times\left(\mu\nabla
B+q\nabla\Phi\right)\boldsymbol{\cdot}\nabla x$
$\displaystyle=\frac{1}{JB_{\parallel}^{*}}\left[-\frac{mv_{\parallel}^{2}}{qB}\frac{\partial(B\,b_{y})}{\partial
z}+\frac{mv_{\parallel}^{2}+\mu B}{qB}b_{y}\frac{\partial B}{\partial
z}-b_{z}\frac{\partial\Phi}{\partial y}+b_{y}\frac{\partial\Phi}{\partial
z}+v_{\parallel}\left(\frac{\partial(A_{\parallel}b_{z})}{\partial
y}-\frac{\partial(A_{\parallel}b_{y})}{\partial z}\right)\right]$ (5.34)
$\displaystyle\dot{y}$
$\displaystyle=\\{y,H\\}=\dot{\mbox{\boldmath${R}$}}\boldsymbol{\cdot}\nabla
y=\frac{B^{*y}}{B_{\parallel}^{*}}v_{\parallel}+\frac{\mathbf{\hat{b}}}{qB_{\parallel}^{*}}\times\left(\mu\nabla
B+q\nabla\Phi\right)\boldsymbol{\cdot}\nabla y$
$\displaystyle=\frac{1}{JB_{\parallel}^{*}}\left[-\frac{mv_{\parallel}^{2}}{qB}\left(\frac{\partial(B\,b_{z})}{\partial
x}-\frac{\partial(B\,b_{x})}{\partial z}\right)+\frac{mv_{\parallel}^{2}+\mu
B}{qB}\left(b_{z}\frac{\partial B}{\partial x}-b_{x}\frac{\partial B}{\partial
z}\right)+b_{z}\frac{\partial\Phi}{\partial
x}-b_{x}\frac{\partial\Phi}{\partial z}\right.$
$\displaystyle\qquad\qquad\quad\left.+v_{\parallel}\left(\frac{\partial(A_{\parallel}b_{x})}{\partial
z}-\frac{\partial(A_{\parallel}b_{z})}{\partial x}\right)\right]$ (5.35)
$\displaystyle\dot{z}$
$\displaystyle=\\{z,H\\}=\dot{\mbox{\boldmath${R}$}}\boldsymbol{\cdot}\nabla
z=\frac{B^{*z}}{B_{\parallel}^{*}}v_{\parallel}+\frac{\mathbf{\hat{b}}}{qB_{\parallel}^{*}}\times\left(\mu\nabla
B+q\nabla\Phi\right)\boldsymbol{\cdot}\nabla z$
$\displaystyle=\frac{1}{JB_{\parallel}^{*}}\left[\mathcal{C}v_{\parallel}+\frac{mv_{\parallel}^{2}}{qB}\frac{\partial(B\,b_{y})}{\partial
x}-\frac{mv_{\parallel}^{2}+\mu B}{qB}b_{y}\frac{\partial B}{\partial
x}+b_{x}\frac{\partial\Phi}{\partial y}-b_{y}\frac{\partial\Phi}{\partial
x}+v_{\parallel}\left(\frac{\partial(A_{\parallel}b_{y})}{\partial
x}-\frac{\partial(A_{\parallel}b_{x})}{\partial y}\right)\right]$ (5.36)

The parallel acceleration is given by

$\displaystyle\dot{v}_{\parallel}$
$\displaystyle=-\frac{1}{JB_{\parallel}^{*}}\left[-\frac{mv_{\parallel}}{q}\frac{\partial
b_{y}}{\partial z}+\left(\frac{\partial(A_{\parallel}b_{z})}{\partial
y}-\frac{\partial(A_{\parallel}b_{y})}{\partial
z}\right)\right]\left(\frac{\mu}{m}\frac{\partial B}{\partial
x}+\frac{q}{m}\frac{\partial\Phi}{\partial x}\right)$
$\displaystyle-\frac{1}{JB_{\parallel}^{*}}\left[\frac{mv_{\parallel}}{q}\left(\frac{\partial
b_{x}}{\partial z}-\frac{\partial b_{z}}{\partial
x}\right)+\left(\frac{\partial(A_{\parallel}b_{x})}{\partial
z}-\frac{\partial(A_{\parallel}b_{z})}{\partial
x}\right)\right]\frac{q}{m}\frac{\partial\Phi}{\partial y}$
$\displaystyle-\frac{1}{JB_{\parallel}^{*}}\left[\mathcal{C}+\frac{mv_{\parallel}}{q}\frac{\partial
b_{y}}{\partial x}+\left(\frac{\partial(A_{\parallel}b_{y})}{\partial
x}-\frac{\partial(A_{\parallel}b_{x})}{\partial
y}\right)\right]\left(\frac{\mu}{m}\frac{\partial B}{\partial
z}+\frac{q}{m}\frac{\partial\Phi}{\partial z}\right)-\frac{q}{m}\frac{\partial
A_{\parallel}}{\partial t}.$ (5.37)

Here we have not neglected any terms due to smallness of parallel derivatives
compared to perpendicular derivatives. While this is a common approximation
made in local gyrokinetic codes, we note that dropping such terms can break
Liouville’s theorem, Eq. 2.82, so that the gyrokinetic equation can no longer
be written in conservative form (phase-space volume is no longer conserved
exactly). In particular, Liouville’s theorem in this case requires

$\nabla\boldsymbol{\cdot}(\mathbf{\hat{b}}\times\nabla
H)=(\nabla\times\mathbf{\hat{b}})\boldsymbol{\cdot}\nabla H$ (5.38)

so that corresponding terms cancel exactly. If parallel derivatives were
dropped, this could become

$\nabla\boldsymbol{\cdot}(\mathbf{\hat{b}}\times\nabla
H)\neq(\nabla\times\mathbf{\hat{b}})_{\perp}\boldsymbol{\cdot}\nabla H,$
(5.39)

which breaks Liouville’s theorem. Since our energy-conserving discontinuous
Galerkin scheme relies on the conservative form of the gyrokinetic equation,
it is important for Liouville’s theorem to be preserved.

####  5.1.5 Field equations in field-aligned coordinates

For the gyrokinetic Poisson equation, Eq. 2.120, we must calculate

$\displaystyle\nabla\boldsymbol{\cdot}(\epsilon_{\perp}\nabla_{\perp}\Phi)$
$\displaystyle=\frac{1}{J}\left[\frac{\partial}{\partial
x}\left(J\epsilon_{\perp}\nabla_{\perp}\Phi\boldsymbol{\cdot}\nabla
x\right)+\frac{\partial}{\partial
y}\left(J\epsilon_{\perp}\nabla_{\perp}\Phi\boldsymbol{\cdot}\nabla
y\right)+\frac{\partial}{\partial
z}\left(J\epsilon_{\perp}\nabla_{\perp}\Phi\boldsymbol{\cdot}\nabla
z\right)\right]$
$\displaystyle\approx\frac{1}{J}\left[\frac{\partial}{\partial
x}\left(J\epsilon_{\perp}\left(\frac{\partial\Phi}{\partial
x}g^{xx}+\frac{\partial\Phi}{\partial
y}g^{xy}\right)\right)+\frac{\partial}{\partial
y}\left(J\epsilon_{\perp}\left(\frac{\partial\Phi}{\partial
x}g^{xy}+\frac{\partial\Phi}{\partial y}g^{yy}\right)\right)\right].$ (5.40)

Here, unlike above, we _do_ neglect the $\partial/\partial z$ terms compared
to the perpendicular derivative terms, so that the required Poisson solve
remains two-dimensional. For energetic consistency, a similar treatment would
need to be made in the corresponding electrostatic field energy term (or in
the second-order $E\times B$ energy term if it is kept in the Hamiltonian).
Similarly, in Ampère’s law, Eq. 2.122, we have

$\nabla_{\perp}^{2}A_{\parallel}=\frac{1}{J}\left[\frac{\partial}{\partial
x}\left(J\left(\frac{\partial A_{\parallel}}{\partial x}g^{xx}+\frac{\partial
A_{\parallel}}{\partial y}g^{xy}\right)\right)+\frac{\partial}{\partial
y}\left(J\left(\frac{\partial A_{\parallel}}{\partial x}g^{xy}+\frac{\partial
A_{\parallel}}{\partial y}g^{yy}\right)\right)\right].$ (5.41)

####  5.1.6 Summary of geometric quantities

The geometry quantities of interest, which appear in either the equations of
motion or the field equations, are

$\displaystyle\texttt{bmag}=B$ (5.42)
$\displaystyle\texttt{cmag}={\mathcal{C}}=\frac{JB}{\sqrt{g_{zz}}}$ (5.43)
$\displaystyle\texttt{b\\_x}=b_{x}=\frac{g_{xz}}{\sqrt{g_{zz}}}$ (5.44)
$\displaystyle\texttt{b\\_y}=b_{y}=\frac{g_{yz}}{\sqrt{g_{zz}}}$ (5.45)
$\displaystyle\texttt{b\\_z}=b_{z}=\sqrt{g_{zz}}$ (5.46)
$\displaystyle\texttt{gxx}=g^{xx}$ (5.47) $\displaystyle\texttt{gxy}=g^{xy}$
(5.48) $\displaystyle\texttt{gyy}=g^{yy}$ (5.49)
$\displaystyle\texttt{jacobPhase}=B_{\parallel}^{*}\approx B$ (5.50)
$\displaystyle\texttt{jacobGeo}=J=\sqrt{\det g_{ij}}.$ (5.51)

Here we have also included the variable names for these quantities in Gkeyll.

###  5.2 Helical SOL configuration including magnetic shear

Thus far, our treatment of field-aligned geometry has been completely general,
except for the assumption of axisymmetry. Now we will examine how a particular
magnetic field geometry affects the choice of field-aligned coordinates and
the resulting metric quantities that appear in the equations.

The helical field in an SMT is given in cylindrical $(R,\varphi,Z)$
coordinates (where $\varphi$ is counter-clockwise viewed from above) by

$\mbox{\boldmath${B}$}=B_{\varphi}\boldsymbol{\hat{\boldsymbol{\varphi}}}+B_{v}\mathbf{\hat{Z}}=\frac{B_{0}R_{0}}{R}\boldsymbol{\hat{\boldsymbol{\varphi}}}+B_{v}\mathbf{\hat{Z}}.$
(5.52)

Note we can also write this as

$\mbox{\boldmath${B}$}=\nabla\Psi\times\nabla\varphi+B_{0}R_{0}\nabla\varphi,$
(5.53)

with $\Psi=R^{2}B_{v}/2$ the vertical magnetic flux function (analogous to the
poloidal flux in a tokamak). The field line pitch varies with radius, and can
be expressed via (Perez et al., 2006)

$q(R)=\frac{HB_{\varphi}}{2\pi RB_{v}}=\frac{B_{0}R_{0}H}{2\pi R^{2}B_{v}},$
(5.54)

where $q$ is analogous to the safety factor in a tokamak, and $H$ is the
vertical height between the bottom and top end-plates where the field lines
terminate. Similarly, we can define the magnetic shear as

$\hat{s}=\frac{R}{q}\frac{\textnormal{d}q}{\textnormal{d}R}.$ (5.55)

Note that if we allow the vertical field to have some simple radial dependence
via $B_{v}=B_{v}(R)=B_{v0}(R/x_{0})^{n}$, the shear is

$\hat{s}=-2-\frac{R}{B_{v}}\frac{\textnormal{d}B_{v}}{\textnormal{d}R}=-2-n.$
(5.56)

When $B_{v}=$ const (as is the case for standard SMTs like the Helimak), we
have $\hat{s}=-2$. The connection length is given by $L_{c}=HB/B_{v}$, which
in general varies with radius. In Fig. 5.1$(a)$ we plot the connection length
$L_{c}$ as a function of radius for several values of $\hat{s}$ with NSTX-like
parameters.

Note that here all quantities, including the field line pitch and the magnetic
shear, are still constant _along_ the field lines. This is in contrast to a
real tokamak SOL, where the pitch and shear can vary significantly along the
field lines, especially near the X-point due to flux expansion. The field line
pitch at the midplane is also nearly constant with radius in a real tokamak
SOL, while here we have the pitch varying with radius. Nonetheless, this
simple geometry will still allow us to study some of the effects of magnetic
shear. Comparing to the actual connection length in the NSTX experiment shown
in Fig. 5.1$(b)$ (with the figure adapted from Fig. 2 in Boedo et al. (2014)),
we see that the variation in $L_{c}$ with radius as $|\hat{s}|$ increases is
approaching the variation of the connection length in the experiment. However,
in the experiment the connection length varies even more than in the
$\hat{s}=-10$ case over a shorter radial length, changing by a factor of
almost four over a few centimeters.

Figure 5.1: $(a)$ Connection length, $L_{c}/2=HB/B_{v}/2$, as a function of
radius for several values of magnetic shear, $\hat{s}$, using NSTX-like
parameters: $R_{0}=0.85$ m, $x_{0}=R_{0}+a=1.35$ m, $B_{0}=0.5$ T, $H=2.4$ m.
Here we choose $B_{v0}$ so that $L_{c}(x_{0})=8$ m, which results in
$B_{v0}\approx 0.1$ T. $(b)$ Connection length from the NSTX experiment,
$L_{\parallel}$, defined as parallel length from midplane to lower divertor
plate. Figure adapted from Boedo et al. (2014).

Now we need field-line-following coordinates $(x,y,z)$, such that

$\mbox{\boldmath${B}$}=\mathcal{C}\nabla x\times\nabla y.$ (5.57)

This can be achieved by choosing the coordinates to be222Alternatively, we
could have chosen the $z$ coordinate to measure distance along the field line,
with $z=Z/\sin\vartheta=ZB/B_{v}$, as we show in Appendix 5.A. In this
approach, the $z$ domain extent is given by the connection length; however,
the connection length can vary with radius, so the $z$ domain extent would
also need to vary with radius. Choosing the vertical height $Z$ for the $z$
coordinate does not have this issue (assuming the vertical height between the
top and bottom end plates does not change with radius), so this is the
approach that we take in the bulk of this chapter.

$\displaystyle x=R,\qquad z=Z,\qquad y=x_{0}\left(\varphi-\frac{2\pi
qZ}{H}\right)=x_{0}\left(\varphi-\frac{B_{\varphi}Z}{B_{v}R}\right).$ (5.58)

The resulting gradient basis vectors are

$\displaystyle\nabla x=\mathbf{\hat{R}},\qquad\nabla y=-\frac{2\pi
x_{0}z}{H}\frac{\textnormal{d}q}{\textnormal{d}x}\mathbf{\hat{R}}+\frac{x_{0}}{x}\boldsymbol{\hat{\boldsymbol{\varphi}}}-\frac{2\pi
x_{0}q}{H}\mathbf{\hat{Z}},\qquad\nabla z=\mathbf{\hat{Z}},$ (5.59)

so that

$\nabla x\times\nabla y=\frac{x_{0}}{x}\left(\frac{2\pi
qx}{H}\boldsymbol{\hat{\boldsymbol{\varphi}}}+\mathbf{\hat{Z}}\right)=\frac{x_{0}}{x}\left(\frac{B_{\varphi}}{B_{v}}\boldsymbol{\hat{\boldsymbol{\varphi}}}+\mathbf{\hat{Z}}\right).$
(5.60)

Now taking $\mathcal{C}=B_{v}x/x_{0}$, we obtain the correct form of the
field,

$\mbox{\boldmath${B}$}=\frac{B_{v}x}{x_{0}}\nabla x\times\nabla
y=B_{\varphi}\boldsymbol{\hat{\boldsymbol{\varphi}}}+B_{v}\mathbf{\hat{Z}}.$
(5.61)

Figure 5.2: Diagram of field-aligned coordinate basis vectors in the
$(\varphi,Z)$ plane for a helical magnetic geometry. Note that vector
magnitudes are not drawn to scale.

We can also calculate the tangent basis vectors

$\mbox{\boldmath${e}$}_{\alpha}=\frac{\partial\mbox{\boldmath${X}$}}{\partial\alpha}=\frac{\partial
R}{\partial\alpha}\mathbf{\hat{R}}+R\frac{\partial\varphi}{\partial\alpha}\boldsymbol{\hat{\boldsymbol{\varphi}}}+\frac{\partial
Z}{\partial\alpha}\mathbf{\hat{Z}},\qquad\alpha=(x,y,z),$ (5.62)

which gives

$\displaystyle\mbox{\boldmath${e}$}_{x}=\mathbf{\hat{R}}+\frac{2\pi
xz}{H}\frac{\textnormal{d}q}{\textnormal{d}x}\boldsymbol{\hat{\boldsymbol{\varphi}}}=\mathbf{\hat{R}}+\frac{B_{\varphi}z}{B_{v}x}\hat{s}\boldsymbol{\hat{\boldsymbol{\varphi}}},\qquad\mbox{\boldmath${e}$}_{y}=\frac{x}{x_{0}}\boldsymbol{\hat{\boldsymbol{\varphi}}},\qquad\mbox{\boldmath${e}$}_{z}=\frac{2\pi
qx}{H}\boldsymbol{\hat{\boldsymbol{\varphi}}}+\mathbf{\hat{Z}}=\frac{B_{\varphi}}{B_{v}}\boldsymbol{\hat{\boldsymbol{\varphi}}}+\mathbf{\hat{Z}}.$
(5.63)

A diagram of the field-aligned basis vectors in the $(\varphi,Z)$ plane is
shown in Fig. 5.2.

Finally, the mapping from the field-aligned coordinates $(x,y,z)$ to physical
Cartesian coordinates $(X,Y,Z)$ is given by

$\displaystyle X=R\cos\varphi=x\cos\left(\frac{y}{x_{0}}+\frac{2\pi
qz}{H}\right)=x\cos\left(\frac{y}{x_{0}}+\frac{B_{\varphi}z}{B_{v}x}\right)$
(5.64) $\displaystyle Y=R\sin\varphi=x\sin\left(\frac{y}{x_{0}}+\frac{2\pi
qz}{H}\right)=x\sin\left(\frac{y}{x_{0}}+\frac{B_{\varphi}z}{B_{v}x}\right)$
(5.65) $\displaystyle Z=z.$ (5.66)

The geometry quantities of interest are

$\displaystyle\texttt{bmag}=B=B_{v}\sqrt{1+\frac{B_{\varphi}^{2}}{B_{v}^{2}}}$
(5.67)
$\displaystyle\texttt{cmag}={\mathcal{C}}=\frac{JB}{\sqrt{g_{zz}}}=\frac{B_{v}x}{x_{0}}$
(5.68)
$\displaystyle\texttt{b\\_x}=b_{x}=\frac{g_{xz}}{\sqrt{g_{zz}}}=\frac{\hat{s}zB_{0}^{2}R_{0}^{2}}{B_{v}Bx^{3}}$
(5.69)
$\displaystyle\texttt{b\\_y}=b_{y}=\frac{g_{yz}}{\sqrt{g_{zz}}}=\frac{B_{0}R_{0}}{Bx_{0}}$
(5.70) $\displaystyle\texttt{b\\_z}=b_{z}=\sqrt{g_{zz}}=\frac{B}{B_{v}}$
(5.71) $\displaystyle\texttt{gxx}=g^{xx}=1$ (5.72)
$\displaystyle\texttt{gxy}=g^{xy}=-\frac{\hat{s}zB_{0}R_{0}x_{0}}{B_{v}x^{3}}$
(5.73)
$\displaystyle\texttt{gyy}=g^{yy}=\frac{B_{0}^{2}R_{0}^{2}x_{0}^{2}}{B_{v}^{2}x^{6}}\left(x^{2}+\hat{s}^{2}z^{2}\right)+\frac{x_{0}^{2}}{x^{2}}$
(5.74) $\displaystyle\texttt{jacobPhase}=B_{\parallel}^{*}\approx B$ (5.75)
$\displaystyle\texttt{jacobGeo}=J=\sqrt{\det g_{ij}}=\frac{x}{x_{0}}.$ (5.76)

Note that in Gkeyll, none of these expressions for the metric quantities are
explicitly implemented; instead the mapping from computational to physical
coordinates, Eqs. 5.64, 5.65 and 5.66, is supplied as an input and then the
metric quantities are computed via automatic differentiation operations. This
enables the flexibility to use more complicated mappings where analytical
expressions for the metric quantities may not be available (see _e.g._ Section
5.3).

We can now compute how various terms in the equations of motion are affected
by the geometry. For example, the components of the magnetic (curvature and
$\nabla B$) drift are

$\displaystyle\mbox{\boldmath${v}$}_{d}\boldsymbol{\cdot}\nabla x=0$ (5.77)
$\displaystyle\mbox{\boldmath${v}$}_{d}\boldsymbol{\cdot}\nabla
y=\frac{B}{JB_{v}}\left[-\frac{mv_{\parallel}^{2}+\mu
B}{qB}\left(\frac{1}{x}+\frac{B_{v}^{2}}{B^{2}}(1+\hat{s})\right)+\frac{mv_{\parallel}^{2}}{qB}\frac{B_{v}^{2}}{B^{2}x}(2+\hat{s})\right]$
(5.78) $\displaystyle\mbox{\boldmath${v}$}_{d}\boldsymbol{\cdot}\nabla
z=\frac{B_{\varphi}}{B}\frac{mv_{\parallel}^{2}+\mu
B}{qB}\left(\frac{1}{x}+\frac{B_{v}^{2}}{B^{2}}(1+\hat{s})\right).$ (5.79)

Comparing these terms to Eq. 4.3 that we used for the simplified helical
geometry in Chapter 4, we see that the first term in
$\mbox{\boldmath${v}$}_{d}\boldsymbol{\cdot}\nabla y$ in Eq. 5.78 is the same
except for a factor of $B/(JB_{v})$. We also have some new terms in
$\mbox{\boldmath${v}$}_{d}\boldsymbol{\cdot}\nabla y$, which are small
corrections when $B_{v}\ll B_{\varphi}\sim B$ as assumed in Chapter 4. We also
have a finite $\mbox{\boldmath${v}$}_{d}\boldsymbol{\cdot}\nabla z$, unlike in
the simplified geometry treatment, though this term is small compared to
$\mbox{\boldmath${v}$}_{d}\boldsymbol{\cdot}\nabla y$ when $B_{v}\ll
B_{\varphi}\sim B$.

####  5.2.1 Simulation results: dependence on magnetic shear in helical
configuration

In this section we present preliminary electrostatic gyrokinetic simulations
in the helical configuration with magnetic shear described above. We perform a
scan of the magnetic shear parameter, taking $\hat{s}=\\{-2,-5,-10\\}$, with
the geometry becoming more sheared as $|\hat{s}|$ increases. We use NSTX-like
geometry parameters: $R_{0}=0.85$ m, $x_{0}=R_{0}+a=1.35$ m, $B_{0}=0.5$ T,
$H=2.4$ m, and we choose $B_{v0}=B_{v}(x_{0})$ so that $L_{c}(x_{0})=8$ m,
resulting in $B_{v0}\approx 0.1$ T. The resulting connection length as a
function of radius is shown in Fig. 5.1. The domain extents in the radial,
binormal, and parallel directions, respectively, are $1.26\leq x\leq 1.42$ m
($L_{x}\approx 56\rho_{\mathrm{s}0}$), $-0.485\leq y\leq 0.485$ m
($L_{y}\approx 100\rho_{\mathrm{s}0}H/L_{c}(x_{0})$), and $-H/2\leq z\leq
H/2$. All other parameters are the same as used in the base case ($\hat{n}=1$)
from Section 4.3, including the source power
$P_{\mathrm{src}}=P_{\mathrm{SOL}}L_{y}/(2\pi R_{c})=0.62$ MW and the source
profile with a Gaussian peak at $x=1.3$ m.

Figure 5.3: Time-averaged midplane density and temperature profiles for
several values of magnetic shear, $\hat{s}$. The profiles steepen as
$|\hat{s}|$ increases and the geometry becomes more sheared.

Figure 5.4: Snapshots of electron density at the midplane for
$\hat{s}=-2,-,5,-10$.

In Fig. 5.4 we show time-averaged density and temperature radial profiles for
each case. The profiles for the $\hat{s}=-2$ case are similar to the profiles
from the base case in Section 4.3, which used a simplified geometry neglecting
magnetic shear and assumed a constant connection length. This is somewhat
expected since Fig. 5.1 shows that the connection length in the $\hat{s}=-2$
case varies little over the domain. As we move to more sheared geometries, the
density profiles steepen, with the peak midplane density more than doubling
between the $\hat{s}=-2$ and $\hat{s}=-10$ cases.

Snapshots of the electron density at the midplane for each case are shown in
Fig. 5.4. While the $\hat{s}=-2$ case looks similar to the cases from Chapter
4 with blobs moving radially outwards, the $\hat{s}=-10$ case shows little
evidence of radial transport. We measure the radial $E\times B$ particle flux
near the midplane in Fig. 5.6 and confirm that indeed particle transport is
reduced as $|\hat{s}|$ increases. Recalling the blob dynamics discussion from
Section 1.3.1, magnetic shear can short-circuit the blob polarization by
allowing currents to close through the thin sheared part of the blob,
resulting in slower blobs. This seems to be consistent with the picture here,
with blob transport getting weaker in more sheared cases.

Figure 5.5: Time-averaged radial electron $E\times B$ particle flux near the
midplane. The cross-field transport decreases as the geometry becomes more
sheared at larger $|\hat{s}|$, resulting in the steeper profiles seen in Fig.
5.4.

Figure 5.6: Time-averaged electron particle and heat fluxes to the bottom
(solid) and top (dashed) end plates. The peak fluxes increase as $|\hat{s}|$
increases, consistent with less cross-field transport upstream.

Figure 5.7: Time-averaged electron radial $E\times B$ particle flux profiles
along the field line, taken just outside the source region at $x\sim 1.32$ m.
There is a slight asymmetry in the profiles, with slightly more radial
transport at $z>0$ than $z<0$ for the $\hat{s}=-5$ and $\hat{s}=-10$ cases.
This is consistent with the differences in the particle fluxes to the bottom
and top end plates shown in Fig. 5.6.

Figure 5.8: Time-averaged electron density profiles along the field line,
taken just outside the source region at $x\sim 1.32$ m. There is a slight
asymmetric shift in the profiles, with slightly higher density near the top of
the domain ($z=H/2=1.4$ m) especially for the more sheared cases.

As a result of weaker cross-field transport, the peak particle and heat fluxes
to the end plate increase in the more sheared cases, as shown in Fig. 5.6.
Note that we have separately shown the fluxes to the top and bottom end
plates. There are some differences in the particle flux profiles between the
ends, especially for the $\hat{s}=-10$ case, where there is a noticeable shift
in the peak to higher $x$ between the bottom and top. When we examine how the
radial particle flux (taken just outside the source region near $x=1.32$ m)
varies along the field line in Fig. 5.8, we see asymmetry as well. There
slightly more radial transport at $z>0$ than $z<0$ for the $\hat{s}=-5$ and
$\hat{s}=-10$ cases. This is consistent with the radial shift in the end plate
particle fluxes to higher $x$ from bottom to top.

The average electron density also shows asymmetry along the field line, as
shown in Fig. 5.8. Here we have again evaluated the profiles just outside the
source region near $x=1.32$ m. There is a slight shift in the profiles to
higher $z$, with the shift larger in the more sheared cases. One possible
reason for the asymmetry is the presence of a vertical component of the
$E\times B$ drift,

$\mbox{\boldmath${v}$}_{E}\boldsymbol{\cdot}\nabla
z=\frac{1}{JB_{\parallel}^{*}}\left(b_{x}\frac{\partial\Phi}{\partial
y}-b_{y}\frac{\partial\Phi}{\partial
x}\right)=\frac{x_{0}}{xB}\left(\frac{\hat{s}zB_{0}^{2}R_{0}^{2}}{B_{v}Bx^{3}}\frac{\partial\Phi}{\partial
y}-\frac{B_{0}R_{0}}{Bx_{0}}\frac{\partial\Phi}{\partial x}\right).$ (5.80)

This term was not present in the earlier simplified geometry simulations from
Section 4.2. It has been suggested that this term is responsible for asymmetry
between top and bottom profiles in the Helimak (Bernard et al., 2020).

###  5.3 Solov’ev model analytical equilibria in the SOL

The Grad-Shafranov equation in cylindrical $(R,Z,\phi)$ coordinates,

$-\mu_{0}RJ_{\phi}=-\nabla^{*}\Psi=-R^{2}\nabla\boldsymbol{\cdot}\frac{1}{R^{2}}\nabla\Psi=\mu_{0}R^{2}p^{\prime}(\Psi)+II^{\prime}(\Psi),$
(5.81)

relates the equilibrium, defined by the poloidal flux function $\Psi$, to the
pressure and current profiles, $p(\Psi)$ and $I(\Psi)$ respectively. Given
$\Psi(R,Z)$, one particular choice for the field-aligned coordinates are
$(x,y,z)=(\Psi,-\alpha,\theta)$, where $x=\Psi$ is the poloidal flux,
$z=\theta$ is a generalized poloidal angle, and $y=-\alpha$ is a field-line-
labeling coordinate defined so that the Clebsch representation of the magnetic
field is given by

$\mbox{\boldmath${B}$}=\nabla\alpha\times\nabla\Psi=\nabla x\times\nabla y,$
(5.82)

with $\mathcal{C}=1$ for this choice of coordinates. Note that for an
axisymmetric system, the background magnetic field can also be expressed as

$\mbox{\boldmath${B}$}=I(\Psi)\nabla\phi+\nabla\Psi\times\nabla\phi,$ (5.83)

where $\phi$ is the toroidal angle and $I(\Psi)=RB_{\phi}$. Thus in practice,
the functions $\Psi(R,Z)$ and $B_{\phi}(R,Z)$ are sufficient to determine the
magnetic geometry. In this section we will take an analytical Solov’ev
solution of the Grad-Shafranov equation for $\Psi(R,Z)$, and show how to
compute the remaining $\theta$ and $\alpha$ coordinates. In principle, the
procedure outlined could be used for an arbitrary $\Psi(R,Z)$ profile,
including one from a numerical equilibrium file generated by _e.g._ EFIT (Lao
et al., 1985, 1990).

Figure 5.9: Flux surfaces (surfaces of constant $\Psi$) for the analytical
Solov’ev equilibrium given by Eq. 5.84, with NSTX-like parameters:
$R_{0}=0.85$, $B_{0}=0.55$, $\kappa_{0}=2$, and $\bar{q}=2$. For the open
field lines, we will take the end-plates to be at $Z=\pm Z_{\text{end}}=\pm
1.5$.

Taking a vacuum field (which is a good approximation in the SOL) so that $I=$
const and approximating $p\propto\Psi$, we obtain an analytical Solov’ev
solution (Chance et al., 1978; Jardin, 2010),

$\Psi(R,Z)=\frac{B_{0}}{2R_{0}^{2}\kappa_{0}\bar{q}}\left(R^{2}Z^{2}+\frac{\kappa_{0}^{2}}{4}(R^{2}-R_{0}^{2})^{2}\right),$
(5.84)

where $R_{0}$ is the major radius of the magnetic axis, $B_{0}$ is the
toroidal field strength at the magnetic axis, $\kappa_{0}$ is the ellipticity
at the axis, and $\bar{q}$ is the safety factor at the axis. The resulting
flux surfaces for NSTX-like parameters are shown in Fig. 5.9. The flux
surfaces are up-down symmetric, and there are closed and open surfaces,
separated by a separatrix given by the surface
$\Psi=\Psi_{\text{sep}}=B_{0}\kappa_{0}R_{0}^{2}/(8\bar{q})$. In this
(unphysical) equilibrium, there are X-points where the separatrix intersects
the $R=0$ axis. However, we will focus only on the open flux surfaces
(sufficiently far away from the X-points), that is, those with
$\Psi>\Psi_{\text{sep}}$. We will assume that the field lines on these open
surfaces terminate on end-plates at $Z=\pm Z_{\text{end}}$. For a given
surface with $\Psi>\Psi_{\text{sep}}$, we can then parametrize the surface
with

$R(\Psi,Z)=\sqrt{R_{0}^{2}+\frac{2Z^{2}}{\kappa_{0}^{2}}+\frac{2}{\kappa_{0}^{2}}\sqrt{\frac{2\kappa_{0}^{3}\bar{q}R_{0}^{2}}{B_{0}}\Psi-\kappa_{0}^{2}R_{0}^{2}Z^{2}+Z^{4}}}.$
(5.85)

We will now calculate the generalized poloidal angle, $\theta$. For this,
first consider a magnetic surface coordinate system $(\Psi,\theta,\phi)$. The
line element in these coordinates is given by

$\textnormal{d}\mbox{\boldmath${\ell}$}=\frac{\partial\mbox{\boldmath${R}$}}{\partial\Psi}\textnormal{d}\Psi+\frac{\partial\mbox{\boldmath${R}$}}{\partial\theta}\textnormal{d}\theta+\frac{\partial\mbox{\boldmath${R}$}}{\partial\phi}\textnormal{d}\phi.$
(5.86)

On a magnetic surface we have $\textnormal{d}\Psi=0$ and in a poloidal plane
we have $\textnormal{d}\phi=0$, so the line element on the surface in the
poloidal plane is simply

$\textnormal{d}\mbox{\boldmath${\ell}$}_{p}=\frac{\partial\mbox{\boldmath${R}$}}{\partial\theta}\textnormal{d}\theta={J}\nabla\phi\times\nabla\Psi\,\textnormal{d}\theta,$
(5.87)

with magnitude

$\textnormal{d}\ell_{p}=|\textnormal{d}\mbox{\boldmath${\ell}$}_{p}|=|{J}\nabla\phi\times\nabla\Psi|\textnormal{d}\theta=\frac{|{J}\nabla\Psi|}{R}\textnormal{d}\theta.$
(5.88)

Now we can find $\theta$ by integrating

$\textnormal{d}\theta=\frac{R}{|{J}\nabla\Psi|}\textnormal{d}\ell_{p}$ (5.89)

along contours of $\Psi$,

$\theta=\int_{\theta_{0}}^{\theta}\textnormal{d}\theta^{\prime}=\int_{\theta_{0}}^{\theta}\frac{R}{|{J}\nabla\Psi|}\textnormal{d}\ell_{p}.$
(5.90)

Using the flux surface parametrization from Eq. 5.85, the differential length
along contours of $\Psi$ is given by

$\textnormal{d}\ell_{p}=\sqrt{1+\left(\frac{\partial R(\Psi,Z)}{\partial
Z}\right)^{2}}\textnormal{d}Z.$ (5.91)

Instead of prescribing the form of the generalized poloidal angle and
subsequently computing the Jacobian
${J}=(\nabla\Psi\times\nabla\theta\boldsymbol{\cdot}\nabla\phi)^{-1}$, we can
instead prescribe the Jacobian to give desired properties of the generalized
poloidal angle (Jardin, 2010). Here, we will choose

${J}=s(\Psi)\frac{R}{|\nabla\Psi|},$ (5.92)

which gives an equal-arc-length poloidal angle in $(-\pi,\pi]$, with

$s(\Psi)=\frac{1}{2\pi}\oint\textnormal{d}\ell_{p}=\frac{1}{\pi}\int_{-Z_{\text{end}}}^{Z_{\text{end}}}\sqrt{1+\left(\frac{\partial
R(\Psi,Z^{\prime})}{\partial Z^{\prime}}\right)^{2}}\textnormal{d}Z^{\prime}$
(5.93)

a normalization factor. This means that on a particular flux surface, the arc
length of each $\Delta\theta$ segment will be equal (see Fig. 5.10). Inserting
this Jacobian definition into Eq. 5.90, the poloidal angle is then given by

$\theta(R,Z)=\frac{1}{s(\Psi(R,Z))}\int_{-Z_{\text{end}}}^{Z}\sqrt{1+\left(\frac{\partial
R(\Psi,Z^{\prime})}{\partial Z^{\prime}}\right)^{2}}\textnormal{d}Z^{\prime}.$
(5.94)

These integrals have no closed form in general and must be evaluated
numerically.

Now we will define the third coordinate, $\alpha$, such that
$\mbox{\boldmath${B}$}=\nabla\alpha\times\nabla\Psi$. To do this, we take
$\alpha$ to be of the form (Kruskal & Kulsrud, 1958)

$\alpha=\phi-q(\Psi)\theta-\nu(\Psi,\theta,\phi),$ (5.95)

where $\nu$ is a to-be-determined function that is periodic in $\theta$ and
$\phi$, and $q(\Psi)$ is the global safety factor, defined as the poloidal
average of the local safety factor $\hat{q}(\Psi,\theta)$,

$q(\Psi)=\frac{1}{2\pi}\int_{0}^{2\pi}\hat{q}(\Psi,\theta)\,\textnormal{d}\theta,$
(5.96)

with

$\hat{q}(\Psi,\theta)=\frac{\mbox{\boldmath${B}$}\boldsymbol{\cdot}\nabla\phi}{\mbox{\boldmath${B}$}\boldsymbol{\cdot}\nabla\theta}=-\mathcal{J}\mbox{\boldmath${B}$}\boldsymbol{\cdot}\nabla\phi=-I(\Psi)\frac{\mathcal{J}}{R^{2}}=\frac{-I(\Psi)s(\Psi)}{R|\nabla\Psi|}.$
(5.97)

It is convenient to define a new toroidal angle $\zeta=\phi-\nu$, so that we
have333Alternatively, we could have defined $\alpha$ in terms of the physical
toroidal angle $\phi$, _i.e._ $\zeta=\phi$, by modifying the poloidal
coordinate to be $\theta^{\prime}=\theta+\nu/q$. In either case, the result is
straight field lines with slope $q$, be it in the $(\theta,\zeta)$ plane or
the $(\theta^{\prime},\phi)$ plane.

$\mbox{\boldmath${B}$}=\nabla\alpha\times\nabla\Psi=\nabla\Psi\times\nabla(q\theta-\zeta).$
(5.98)

Here we can see that the field lines are straight lines with slope $q$ in the
$(\theta,\zeta)$ plane, given by $\alpha=\zeta-q(\Psi)\theta=\text{const}$.

Figure 5.10: Flux surfaces with
$\Psi_{N}=\Psi/\Psi_{\text{sep}}=(1.2,1.4,1.6,1.8,2.0)$ are shown in blue.
Lines of constant $\theta$ are shown in red. Here, $\theta$ is defined to be
an equal-arc-length poloidal angle. This means that on a particular flux
surface, the arc length of each $\Delta\theta=\pi/4$ segment is equal.

To compute $\alpha$, we first note that

$\nabla\alpha=\frac{\partial\alpha}{\partial\theta}\nabla\theta+\frac{\partial\alpha}{\partial\Psi}\nabla\Psi+\frac{\partial\alpha}{\partial\phi}\nabla\phi,$
(5.99)

which gives

$\mbox{\boldmath${B}$}=\nabla\Psi\times\nabla\alpha=\frac{\partial\alpha}{\partial\theta}\nabla\Psi\times\nabla\theta+\frac{\partial\alpha}{\partial\phi}\nabla\Psi\times\nabla\phi.$
(5.100)

Now notice

$\mbox{\boldmath${B}$}\boldsymbol{\cdot}\nabla\phi=\frac{\partial\alpha}{\partial\theta}\nabla\Psi\times\nabla\theta\boldsymbol{\cdot}\nabla\phi=\frac{1}{{J}}\frac{\partial\alpha}{\partial\theta},$
(5.101)

and from Eq. 5.83, we also have

$\mbox{\boldmath${B}$}\boldsymbol{\cdot}\nabla\phi=I(\Psi)|\nabla\phi|^{2}=\frac{I(\Psi)}{R^{2}}.$
(5.102)

Thus we can integrate along $\theta$ (at constant $\Psi$) to find $\alpha$,

$\alpha=C(\Psi,\phi)+\int_{0}^{\theta}\frac{\partial\alpha}{\partial\theta^{\prime}}\textnormal{d}\theta^{\prime}=\phi-I(\Psi)\int_{0}^{\theta}\frac{{J}}{R^{2}}\textnormal{d}\theta^{\prime}=\phi-
RB_{\phi}\int\frac{1}{|\nabla\Psi|R}\textnormal{d}\ell_{p},$ (5.103)

where we have taken the constant of integration to be
$C(\Psi,\phi)=\alpha(\theta=0,\Psi,\phi)=\phi$. As we did for $\theta$, this
integral can be computed by integrating along the contours of $\Psi$
parametrized by $Z$,

$\alpha(R,Z,\phi)=\phi-
RB_{\phi}\int_{0}^{Z}\frac{1}{|\nabla\Psi|R(\Psi,Z^{\prime})}\sqrt{1+\left(\frac{\partial
R(\Psi,Z^{\prime})}{\partial Z^{\prime}}\right)^{2}}\textnormal{d}Z^{\prime}.$
(5.104)

Figure 5.11: Selecting a line of constant $\alpha$ and constant $\Psi$ traces
a field line (blue) from the bottom end-plate (at $\theta=-\pi$ and
$Z=-Z_{\text{end}}=-1.5$) to the top end-plate (at $\theta=\pi$ and
$Z=Z_{\text{end}}=1.5$). A rendering of part of the constant $\Psi_{N}=1.2$
surface is shown in the background.

We now have expressions for $\Psi(\mbox{\boldmath${R}$})$,
$\alpha(\mbox{\boldmath${R}$})$, and $\theta(\mbox{\boldmath${R}$})$. We can
thus define a field-aligned coordinate system with
$(x,y,z)=(\Psi,-\alpha,\theta)$. (The difference of sign between $y$ and
$\alpha$ is a matter of convention, and we follow Beer et al. (1995) here).
The expressions for $\alpha$ and $\theta$ involve integrals that must be
evaluated numerically in most cases. In Fig. 5.10, we show lines of constant
$\theta$ (and $\alpha$) in the poloidal plane for the open-field-line region
of the NSTX-like equilibrium shown in Fig. 5.9, demonstrating the equal-arc-
length poloidal angle. In Fig. 5.11, we show that a line of constant $\alpha$
and constant $\Psi$ traces a field line from the bottom end-plate to the top
end-plate. Further, note that the connection length can be computed from

$L_{c}=\oint\sqrt{g_{zz}}\textnormal{d}\ell_{p}=\int_{-Z_{\text{end}}}^{Z_{\text{end}}}\sqrt{g_{zz}(R,Z^{\prime})}\sqrt{1+\left(\frac{\partial
R(\Psi,Z^{\prime})}{\partial Z^{\prime}}\right)^{2}}\textnormal{d}Z^{\prime}.$
(5.105)

Fig. 5.12 shows the connection length as a function of the poloidal flux
$\Psi$ for the NSTX-like equilibrium.

Figure 5.12: Connection length, $L_{c}$, as a function of normalized poloidal
flux $\Psi_{N}=\Psi/\Psi_{\text{sep}}$ for the NSTX-like equilibrium.

We also need derivatives of the coordinates to compute metric quantities,
which can also be computed numerically via automatic differentiation or finite
differencing; alternatively, integral expressions for the derivatives can also
be derived via the Leibniz integral rule. In Fig. 5.13, we show some of the
resulting geometric quantities. There is significant flux expansion and
magnetic shear for flux surfaces near the separatrix and X-points, which
causes some of the metric quantities to diverge near the top-left and bottom-
left corners of the $(\Psi,\theta)$ domain as $\Psi$ approaches
$\Psi_{\text{sep}}$. Finally, note that the Jacobian of the
$(\Psi,\alpha,\theta)$ coordinate system,

$J=(\nabla\Psi\times\nabla\alpha\boldsymbol{\cdot}\nabla\theta)^{-1}=s(\Psi)\frac{R}{|\nabla\Psi|},$
(5.106)

is equivalent to the Jacobian defined in Eq. 5.92 for the $(\Psi,\theta,\phi)$
magnetic surface coordinate system.

Figure 5.13: Geometric quantities resulting from the field-aligned coordinate
system as derived in this section for the NSTX-like equilibrium shown in Fig.
5.9.

Linear and nonlinear simulations using the shaped SOL geometry presented in
this section are left to future work. An important question is how close the
simulation domain can approach the X-point before the significant magnetic
shear and metric divergence near the X-point become numerically untenable.

###  5.4 Analytical concentric circular equilibria

In the tokamak core, many early results and inter-code benchmarks have used an
_ad hoc_ analytical equilibrium model with circular concentric flux surfaces.
This is commonly given by the popular $s-\alpha$ model with no Shafranov shift
($\alpha=0$). This type of geometry has also been used to model circular SOL
plasmas with a limiter at the inboard midplane, modeling an inner-wall-limited
plasma (Ricci & Rogers, 2013; Halpern et al., 2013; Francisquez et al., 2017;
Zhu et al., 2017). Here we use a slightly modified circular equilibrium model
that is more consistent than the $s-\alpha$ model in the large aspect ratio
approximation (Lapillonne et al., 2009).

As in the previous section, we start from a Solov’ev solution of the form of
Eq. 5.84. We use a toroidal coordinate system $(r,\theta,\phi)$, where
$(r,\theta)$ are the minor radius and poloidal angle coordinates in the
$(R,Z)$ plane such that $R=R_{0}+r\cos\theta$ and $Z=r\sin\theta$, and $\phi$
is the toroidal angle. Taking the large aspect ratio limit $R_{0}/a\gg 1$ and
$\kappa_{0}=1$ for circular flux surfaces, we obtain

$\Psi=\frac{B_{0}}{2\bar{q}}r^{2}.$ (5.107)

The resulting magnetic field is

$\mbox{\boldmath${B}$}=\nabla\phi\times\nabla\Psi+RB_{\phi}\nabla\phi=\frac{R_{0}B_{0}}{R}\left[\mbox{\boldmath${e}$}_{\phi}+\frac{r}{R_{0}\bar{q}}\mbox{\boldmath${e}$}_{\theta}\right].$
(5.108)

While Eq. 5.107 implies $\bar{q}=\text{const}$, we will generalize the
equilibrium to allow $\bar{q}=\bar{q}(r)$, which is related to the true safety
factor via (Lapillonne et al., 2009)

$q(r)=\frac{1}{2\pi}\int_{0}^{2\pi}\frac{{\bf
B}\boldsymbol{\cdot}\nabla\phi}{{\bf B}\boldsymbol{\cdot}\nabla\theta}\
\textnormal{d}\theta=\frac{\bar{q}}{2\pi}\int_{0}^{2\pi}\frac{\textnormal{d}\theta}{1+\epsilon\cos\theta}=\frac{\bar{q}}{\sqrt{1-\epsilon^{2}}},$
(5.109)

where $\epsilon=r/R_{0}$ is the inverse aspect ratio. We instead define the
poloidal flux in terms of its radial derivative,
$\textnormal{d}\Psi/\textnormal{d}r=rB_{0}/\bar{q}$, giving

$\Psi=\int_{0}^{r}\frac{r^{\prime}B_{0}}{\bar{q}(r^{\prime})}\,\textnormal{d}r^{\prime}$
(5.110)

instead of Eq. 5.107.

Here we will choose a straight-field-line poloidal angle $\chi$ defined such
that field lines are straight in the $(\chi,\phi)$ plane with slope $q$. To do
this we take
$({\mbox{\boldmath${B}$}\boldsymbol{\cdot}\nabla\phi})/({\mbox{\boldmath${B}$}\boldsymbol{\cdot}\nabla\chi})=q$,
which leads to
$\textnormal{d}\chi/\textnormal{d}\theta=(\mbox{\boldmath${B}$}\boldsymbol{\cdot}\nabla\phi)/(q\mbox{\boldmath${B}$}\boldsymbol{\cdot}\nabla\theta)$.
Integrating over $\theta$ then gives

$\chi(r,\theta)=\frac{1}{q}\int_{0}^{\theta}\frac{{\bf
B}\boldsymbol{\cdot}\nabla\phi}{{\bf
B}\boldsymbol{\cdot}\nabla\theta^{\prime}}\
\textnormal{d}\theta^{\prime}=\frac{\bar{q}}{q}\int_{0}^{\theta}\frac{\textnormal{d}\theta^{\prime}}{1+\epsilon\cos\theta^{\prime}}=2\arctan\left[\sqrt{\frac{1-\epsilon}{1+\epsilon}}\tan\left(\frac{\theta}{2}\right)\right].$
(5.111)

Now we can define the field-aligned coordinate system $(x,y,z)$ as

$x=r-x_{0},\qquad y=\frac{r_{0}}{q_{0}}\left(q\chi-\phi\right)-y_{0},\qquad
z=\chi,$ (5.112)

where $r_{0}$ is the minor radius of a flux surface of interest, and
$q_{0}=q(r_{0})$.

Figure 5.14: Linear growth rate results for the Cyclone base case, with Gkeyll
results (blue diamonds) compared to several other gyrokinetic codes. Note that
since Gkeyll does not include gyroaveraging we only expect accuracy in the
small $k_{y}\rho_{\mathrm{s}}$ limit. Adapted from Figure 3a of (Görler et
al., 2016).

####  5.4.1 Cyclone base case linear benchmark

Here we perform the now-standard Cyclone base case linear benchmark with
Gkeyll using the circular equilibrium described in the previous section. These
calculations use the electrostatic approximation with adiabatic electrons, as
in the original benchmark (Dimits et al., 2000). Like the KBM calculations in
Section 3.4.2, these simulations use a single narrow cell in the radial
direction, making them effectively “local”, unlike a “global” calculation that
covers some portion of the tokamak minor radius and accounts for the radial
variation of background quantities.444While Gkeyll could in principle be
capable of performing global calculations, one must be careful with the
initial conditions to avoid an equilibrium-scale $n=0$ mode that results from
the neoclassical terms. This can obscure or alter the growth rate of the
$n\neq 0$ mode of primary interest. The so-called “canonical” Maxwellian,
which is formulated in terms of the canonical momentum so that it is an
equilibrium of the full-$f$ gyrokinetic equation including the neoclassical
terms, is used in some full-$f$ gyrokinetic codes to avoid exciting the $n=0$
mode (Angelino et al., 2006). The implementation of a canonical Maxwellian
initial condition is currently in progress, which will enable global
instability calculations. We use an extended domain along the field line with
$-3\pi\leq\chi\leq 3\pi$ and Dirichlet boundary conditions $f(\chi=\pm
3\pi)=F_{0}$ so that no fluctuations are allowed at the domain ends.

Figure 5.15: Linear eigenmode for $n=10$ case, showing ballooning mode
structure along the field line.

We compare our results to the inter-code benchmark of Görler et al. (2016);
relevant parameters for the simulation setup are given in Tables I and II of
that reference. In Fig. 5.14 we have reproduced Figure 3a of Görler et al.
(2016) with the addition of some Gkeyll results. We see good agreement for
$n=5$, where $n=k_{y}r_{0}/q_{0}$ is the toroidal mode number. Since Gkeyll
does not include gyroaveraging, we only expect accuracy in the small
$k_{y}\rho_{\mathrm{s}}$ limit. Consequently, we overpredict the growth rate
as we move to higher mode numbers. We also show the linear eigenmode for the
$n=10$ case in Fig. 5.15, which has the characteristic peak at the outboard
midplane $(\chi=0)$ of a ballooning mode. Here we have removed the $n=0$
component of the potential by subtracting its $y$-average.

Additional work will extend these benchmarks to the two-species electrostatic
cases and the electromagnetic cases that are also examined in Görler et al.
(2016). We will also include global effects in future work.

##  Appendix 5.A Alternative coordinate mappings for helical geometry

One may note that the helical coordinate mapping we used in this chapter, Eq.
5.58, does not match the mapping proposed in Shi et al. (2019) or Chapter 4,
given by

$R=x,\qquad\varphi=\frac{y\sin\vartheta+z\cos\vartheta}{R_{c}},\qquad
Z=z\sin\vartheta,$ (5.113)

where $\sin\vartheta=B_{v}/B=H/L_{c}$, with $\vartheta$ the field line pitch
angle (as shown in Fig. 5.2) so that

$x=R,\qquad z=\frac{Z}{\sin\vartheta}=\frac{L_{c}}{H}Z,\qquad
y=\frac{R_{c}}{\sin\vartheta}\left(\varphi-\frac{Z}{R_{c}}\cot\vartheta\right)=R_{c}\frac{L_{c}}{H}\left(\varphi-\frac{Z}{R_{c}}\frac{B_{\varphi}}{B_{v}}\right).$
(5.114)

After computing the resulting gradient basis vectors, we have

$\nabla x\times\nabla
y=\frac{L_{c}}{H}\left({\frac{B_{\varphi}}{B_{v}}\boldsymbol{\hat{\boldsymbol{\varphi}}}+\frac{R_{c}}{x}\mathbf{\hat{Z}}}\right).$
(5.115)

Taking $\mathcal{C}=HB_{v}/L_{c}$, the resulting background magnetic field in
Clebsch form is

$\mbox{\boldmath${B}$}=\mathcal{C}\nabla x\times\nabla
y=B_{\varphi}\boldsymbol{\hat{\boldsymbol{\varphi}}}+\frac{R_{c}}{x}B_{v}\mathbf{\hat{Z}}.$
(5.116)

We see that this only gives the correct pitch of the magnetic field at
$x=R_{c}$, so _this is not a field-aligned coordinate system_.

We can fix this issue and make the coordinates field-aligned by changing the
$\varphi$ mapping to

$\varphi=\frac{y\sin\vartheta}{R_{c}}+\frac{z\cos\vartheta}{x}=\frac{y}{R_{c}}\frac{H}{L_{c}}+\frac{z}{x}\sqrt{1-\frac{H^{2}}{L_{c}^{2}}},$
(5.117)

so that

$y=\frac{R_{c}}{\sin\vartheta}\left(\varphi-\frac{Z}{R}\cot\vartheta\right)=R_{c}\frac{L_{c}}{H}\left(\varphi-\frac{B_{\varphi}Z}{B_{v}R}\right),$
(5.118)

which is the same as the definition of $y$ from Eq. (5.58) except for a factor
of $1/\sin\vartheta$ (with $x_{0}=R_{c}$ here). This gives

$\nabla x\times\nabla
y=\frac{R_{c}}{x}\left(\cos\vartheta\hat{\boldsymbol{\varphi}}+\sin\vartheta\mathbf{\hat{Z}}\right)=\frac{R_{c}}{x}\frac{L_{c}}{H}\left(\frac{B_{\varphi}}{B_{v}}\boldsymbol{\hat{\boldsymbol{\varphi}}}+\mathbf{\hat{Z}}\right).$
(5.119)

Taking $\mathcal{C}=Bx/R_{c}$, we have

${\bf B}=\mathcal{C}\nabla x\times\nabla
y=B\left(\cos\vartheta\boldsymbol{\hat{\boldsymbol{\varphi}}}+\sin\vartheta\mathbf{\hat{Z}}\right)=B_{\varphi}\boldsymbol{\hat{\boldsymbol{\varphi}}}+B_{z}\mathbf{\hat{Z}},$
(5.120)

which is now the correct form of the magnetic field.

This mapping uses the distance along the field line as the field-aligned $z$
coordinate. For SMT configurations, where the connection length can vary
radially, this choice could mean that the $z$ simulation domain extents must
also vary radially. We could instead normalize the $z$ coordinate to the
connection length, resulting in an equal-arc-length-like parallel coordinate,
but this would result in the parallel coordinate becoming proportional to the
vertical height $Z$. Thus we have used the vertical height $Z$ as the $z$
coordinate in Section 5.2.

##  Appendix 5.B Self-consistently reproducing “simplified” helical geometry

Here we would like to effectively reproduce the “simplified” helical geometry
that we used to produce the results in Section 4.2. As we will see, the actual
geometry that matches the approximations we made in Section 4.1.1 is not
helical, but purely toroidal; it essentially consists of rings of toroidal
field stacked vertically on top of each other, as shown in Fig. 5.16. The
mapping from cylindrical coordinates $(R,\varphi,Z)$ to field-aligned
coordinates $(x,y,z)$ that gives this simplified geometry is

$x=R,\qquad y=Z,\qquad z=R_{c}\varphi.$ (5.121)

The tangent unit vectors are

$\mbox{\boldmath${e}$}_{x}=\mathbf{\hat{R}},\qquad\mbox{\boldmath${e}$}_{y}=\mathbf{\hat{Z}},\qquad\mbox{\boldmath${e}$}_{z}=\frac{x}{R_{c}}\boldsymbol{\hat{\boldsymbol{\varphi}}},$
(5.122)

and the gradient unit vectors are

$\nabla x=\mathbf{\hat{R}},\qquad\nabla y=\mathbf{\hat{Z}},\qquad\nabla
z=\frac{R_{c}}{x}\boldsymbol{\hat{\boldsymbol{\varphi}}}.$ (5.123)

Figure 5.16: Geometry with purely toroidal magnetic field. This is the
geometry that matches the metric-related approximations we made for the
results in Section 4.2. In those simulations, sheath boundary conditions
applied in the $z$ direction are equivalent to a limiter that extends
vertically at some toroidal location (although the “connection length” along
the field lines used in the simulations is not necessarily the circumference
of the toroidal rings). Periodic boundary conditions were applied in the $y$
direction, which is the vertical direction here.

The magnetic field is then given by

$\mbox{\boldmath${B}$}=\mathcal{C}\nabla x\times\nabla
y=\frac{JB}{\sqrt{g_{zz}}}\mathbf{\hat{R}}\times\mathbf{\hat{Z}}=B\boldsymbol{\hat{\boldsymbol{\varphi}}}$
(5.124)

with $B=B_{\text{axis}}(R_{0}/x)$, so that the field is indeed purely
toroidal.

Note that here, the coordinate along the field line, $z$, is defined to be
proportional to the toroidal angle, $\varphi$. Whereas elsewhere in this
chapter we used the poloidal angle as the field-aligned coordinate, this is
not possible here where the field lines have no pitch. In fact, this is
precisely the issue with using a field-aligned poloidal coordinate at the
X-point, where the poloidal magnetic field vanishes and the field is purely
toroidal.

The remaining geometry quantities of interest are

$\displaystyle\texttt{cmag}={\mathcal{C}}=\frac{JB}{\sqrt{g_{zz}}}=B$ (5.125)
$\displaystyle\texttt{b\\_x}=b_{x}=\frac{g_{xz}}{\sqrt{g_{zz}}}=0$ (5.126)
$\displaystyle\texttt{b\\_y}=b_{y}=\frac{g_{yz}}{\sqrt{g_{zz}}}=0$ (5.127)
$\displaystyle\texttt{b\\_z}=b_{z}=\sqrt{g_{zz}}=\frac{x}{R_{c}}$ (5.128)
$\displaystyle\texttt{gxx}=g^{xx}=1$ (5.129)
$\displaystyle\texttt{gxy}=g^{xy}=0$ (5.130)
$\displaystyle\texttt{gyy}=g^{yy}=1$ (5.131)
$\displaystyle\texttt{jacobPhase}=B_{\parallel}^{*}\approx B$ (5.132)
$\displaystyle\texttt{jacobGeo}=J=\sqrt{\det g_{ij}}=\frac{x}{R_{c}}.$ (5.133)

We can see that many of the metric quantities are trivial or vanish, just as
we approximated in Section 4.1.1. The magnetic (curvature and $\nabla B$)
drifts are then

$\displaystyle\mbox{\boldmath${v}$}_{d}\boldsymbol{\cdot}\nabla x=0$ (5.134)
$\displaystyle\mbox{\boldmath${v}$}_{d}\boldsymbol{\cdot}\nabla
y=-\frac{mv_{\parallel}^{2}+\mu B}{qB}\frac{1}{x}$ (5.135)
$\displaystyle\mbox{\boldmath${v}$}_{d}\boldsymbol{\cdot}\nabla z=0,$ (5.136)

which matches Eq. 4.3.

##  Chapter 6 Positivity-preserving discontinuous Galerkin algorithm for
hyperbolic conservation laws without _post-hoc_ diffusion

Physically, the distribution function of particles is a non-negative scalar
function, _i.e._ $f(\mbox{\boldmath${x}$},\mbox{\boldmath${v}$},t)\geq 0$
throughout the phase space. However, there is no guarantee that a numerical
scheme will preserve this property. The discontinuous Galerkin schemes that we
described in Chapter 3 are no exception, as they do not even ensure positivity
of the cell average of the distribution function. In some cases small regions
of negative $f$ do not impact the physics, but in other cases negative regions
can lead to numerical instability. Thus we must develop a method to prevent
negative regions in order for our simulations to be accurate and robust.

There is extensive literature on constructing positivity-preserving (or more
generally, bound-preserving) DG schemes. For example, a widely used and now
standard positivity-limiting scheme presented by Zhang & Shu (2011) works by
limiting the amount of flux leaving a cell surface so that the cell average is
not allowed to become negative. However, this limiter procedure by itself can
produce unphysical steep slopes and higher moments. For this reason, a _post-
hoc_ sub-cell diffusion step is applied whereby the slopes and higher moments
in each cell are adjusted so that the solution remains positive at some
control points. This leads to a robust positivity-preserving algorithm for
many conservation laws, like the Euler equations and the ideal MHD equations.
Generalizations of the Zhang-Shu scheme have also been made (Johnson &
Rossmanith, 2012).

However, the Zhang-Shu (and related) algorithms cannot be used for evolution
of kinetic equations in which the conservation properties are indirect (_i.e._
when there is not a direct equation for the evolution of the energy). The
reason for this is that _post-hoc_ sub-cell diffusion will change the energy
and break energy conservation. One could try to readjust the energy after the
diffusive step to maintain energy conservation, as done in Shi (2017), but
often such adjustments are not possible without disturbing the underlying
physics.

In this chapter we develop a novel positivity-preserving DG scheme without
_post-hoc_ diffusion. After showing a preliminary example of the positivity
issue, we first define what we mean by positivity in the context of the
discontinuous Galerkin representation of the cell. Next we construct the
positivity-preserving scheme and show that conservation properties are
maintained for Hamiltonian systems. Finally we show numerical results in
several dimensionalities and equation systems.

###  6.1 The positivity problem: 1D advection example

Consider an one-dimensional advection equation,

$\frac{\partial f}{\partial t}+v\frac{\partial f}{\partial x}=0.$ (6.1)

Taking $v$ to be constant in this simple example, we know the solution is

$f(x,t)=f(x-vt,0).$ (6.2)

That is, the solution keeps its initial shape as it moves through the domain
with constant velocity $v$. If the domain is periodic with length $L$, the
pulse will return to its initial position at time $t=L/v$.

Figure 6.1: Advection of a square pulse for one period on a 1D periodic
domain, overlaid with the initial condition.

We can easily discretize this system with a discontinuous Galerkin scheme. In
Fig. 6.1, we show the results of a piecewise-linear scheme with 32 cells on a
one-dimensional periodic domain with length $L=1$. The initial condition is a
square pulse centered at $x=0.5$ with width $w=0.5$, shown dash-dotted. After
1 period, the solution (solid lines) returns to its initial position, although
the shape of the pulse has been distorted by numerical artifacts. Notably, we
see unphysical negative overshoot regions at the bottom of the pulse, as well
as positive overshoots at the top of the pulse. Cells with a negative cell
average are marked with points showing the cell center. While in some
applications small negative overshoots may be tolerable, often these
unphysical negative regions can cause severe problems. Our goal in this
chapter will be to devise a positivity-preserving scheme that eliminates these
negative regions.

###  6.2 Defining positivity, in the weak sense

The first challenge is to define what is meant by positivity in the context of
the discontinuous Galerkin representation of the solution. In each cell, the
solution is given by an expansion on some basis set. The goal of this section
will be to develop a method to constrain the expansion coefficients to
maintain positivity of the solution, in some sense. Let’s first consider the
simplest case: a piecewise-linear ($p=1$) representation in one dimension.
Taking an orthonormal basis set $\psi=\\{1/\sqrt{2},\sqrt{3/2}x\\}$ in a cell
$x\in[-1,1]$, we have

$f_{h}=\sum_{k}f_{k}\psi_{k}=\frac{1}{\sqrt{2}}f_{0}+\frac{\sqrt{3}}{\sqrt{2}}xf_{1}.$
(6.3)

How can we constrain the coefficients $f_{0}$ and $f_{1}$ to ensure that the
solution is positive? To start, we should at least ensure that the cell
average is positive, so that $f_{0}\geq 0$. Should the solution be required to
be positive on the whole cell domain $x=[-1,1]$, or can the solution be
negative on some portion of the domain?

One possible way to answer these questions is to define positivity as _weak
equality to a positive-definite function_ (Hakim et al., 2020). For example,
we could consider a non-polynomial positive-definite exponential solution
given by

$g_{h}=g_{0}\exp(g_{1}x).$ (6.4)

Weak equality of $f_{h}$ and $g_{h}$ in the $L_{2}$ sense, which we denote as
$f_{h}\doteq g_{h}$, then requires that the projections of the two
representations onto the basis be equivalent:

$\displaystyle\int_{-1}^{1}\frac{1}{\sqrt{2}}f_{h}\,\textnormal{d}x=\int_{-1}^{1}\frac{1}{\sqrt{2}}g_{h}\,\textnormal{d}x\quad$
$\displaystyle\Rightarrow\quad f_{0}=\frac{\sqrt{2}g_{0}\sinh g_{1}}{g_{1}}$
(6.5)
$\displaystyle\int_{-1}^{1}\frac{\sqrt{3}}{\sqrt{2}}xf_{h}\,\textnormal{d}x=\int_{-1}^{1}\frac{\sqrt{3}}{\sqrt{2}}xg_{h}\,\textnormal{d}x\quad$
$\displaystyle\Rightarrow\quad
f_{1}=\frac{\sqrt{6}g_{0}}{g_{1}^{2}}\left(g_{1}\cosh g_{1}-\sinh
g_{1}\right).$ (6.6)

Note that

$\frac{f_{1}}{\sqrt{3}f_{0}}=\coth g_{1}-\frac{1}{g_{1}}=L(g_{1}),$ (6.7)

where $L(g_{1})$ is the well-known Langevin function in statistical mechanics.
Notably, this function is bounded at $|L(g_{1})|\leq 1$ for all $g_{1}$. This
means that in order for weak-equivalence of the solutions $f_{h}$ and $g_{h}$
to be possible, the coefficients of $f_{h}$ must satisfy

$\frac{|f_{1}|}{\sqrt{3}f_{0}}\leq 1.$ (6.8)

Together with the constraint $f_{0}\geq 0$, we now have positivity constraints
for both coefficients of the piecewise-linear representation of the solution.

We can now express $g_{h}(x)$ in terms of $f_{0}$ and $\bar{x}\equiv
f_{1}/(\sqrt{3}f_{0})$ via

$g_{h}(x)=\frac{f_{0}g_{1}}{\sqrt{2}\sinh g_{1}}e^{g_{1}x},$ (6.9)

where $g_{1}=g_{1}(\bar{x})=L^{-1}(\bar{x})$, and $L^{-1}$ is the inverse
Langevin function.111Although the inverse of the Langevin function does not
have a closed form, a number of Padé approximations for the inverse have been
developed, including (Cohen, 1991)
$L^{-1}(x)\approx\frac{x(3-x^{2})}{1-x^{2}}.$ (6.10) Fig. 6.2 shows an example
of the weak-equivalent linear and exponential solutions with $f_{0}=0.4$ and
$f_{1}=0.8$.

Note that the constraint from Eq. 6.8 does not force the linear solution to be
positive everywhere in the cell. In this case, the linear solution is negative
for $x<-0.5$ but the exponential solution is still realizable. In fact, one
can show that the constraint on $f_{1}$ is equivalent to requiring
$f_{h}(x=\pm 1/3)\geq 0$, so that as long as the linear solution remains
positive at “positivity control nodes” $x=\pm 1/3$, the solution will be
positive in the weak sense. The control nodes are also plotted in Fig. 6.2,
and they are indeed both positive.

Figure 6.2: Weak-equivalent linear and exponential solutions, with $f_{0}=0.4$
and $f_{1}=0.8$. The exponential solution $g_{h}$ is given by Eq. 6.9. The
“positivity control nodes” $f_{h}(x=\pm 1/3)$ are also shown.

####  6.2.1 Generalization to higher dimensionality

It is not immediately clear how to tractably extend the procedure of the
previous section to higher dimensionality. For example, to extend to a two-
dimensional case, one might consider taking the exponential representation to
be of the form $g_{h}=\exp(g_{0}+g_{1}x+g_{2}y+g_{3}xy)$. However, 2D
integrals of this function involve error functions, making it difficult to
evaluate positivity constraints based on weak equality.

We will discuss a more rigorous procedure for generalizing the definition of
the weak-equivalent positive-definite solution to higher dimensionality (and
higher polynomial order) in Appendix 6.B. For now, however, let’s continue
with the idea of “positivity control nodes”. In one dimension, we saw above
that if the piecewise-linear solution is positive at the control nodes at
$x=\pm 1/3$, then the solution can be made weak-equivalent to an exponential
solution. A sensible extension of this idea to higher dimension is to take
tensor products of the control nodes, so that for example in 2D, we have
control nodes at $(-1/3,-1/3),\ (-1/3,1/3),\ (1/3,-1/3)$, and $(1/3,1/3)$.
Thus in the following section, we will consider the solution to be positive
(in the weak sense) if the $N$-dimensional piecewise-linear solution is non-
negative at all of the $2^{N}$ control nodes.

###  6.3 Constructing a positivity-preserving scheme without _post-hoc_
diffusion

Now that we have a definition of positivity, we next focus on how to construct
a discontinuous Galerkin scheme that preserves positivity. In our scheme, we
would like to avoid _post-hoc_ sub-cell diffusion (rescaling slopes or higher
moments of the solution in a cell if they become too extreme and violate
positivity constraints after taking a timestep), which can break conservation
laws involving higher-order moments (such as energy conservation in
Hamiltonian systems like gyrokinetics).

To begin, we will once again consider a generic hyperbolic conservation law of
the form of Eq. 3.11,

$\frac{\partial f}{\partial
t}+\nabla\boldsymbol{\cdot}\mbox{\boldmath${F}$}=0,$ (6.11)

with $\mbox{\boldmath${F}$}(f)$ some arbitrary nonlinear flux. Recall from Eq.
3.12 that the DG discretization of this equation is given by multiplying by a
test function $\psi$ and integrating by parts over a cell $\mathcal{K}_{i}$:

$\int_{\mathcal{K}_{i}}\psi\frac{\partial f_{h}}{\partial
t}\textnormal{d}\mbox{\boldmath${x}$}-\int_{\mathcal{K}_{i}}\mbox{\boldmath${F}$}_{h}\boldsymbol{\cdot}\nabla\psi\,\textnormal{d}\mbox{\boldmath${x}$}+\oint_{\partial\mathcal{K}_{i}}\psi^{-}\mbox{\boldmath${\hat{F}}$}\boldsymbol{\cdot}\textnormal{d}\mbox{\boldmath${s}$}=0.$
(6.12)

Let us first focus on the one-dimensional piecewise-linear ($p=1$) case.
Mapping each cell $\mathcal{K}_{i}$ to the interval $x\in[-1,1]$ using the
transformation $x^{\prime}=2(x-x_{i})/\Delta x$, with $x_{i}$ the cell center
and $\Delta x$ the cell width, this gives (after dropping primes for
simplicity)

$\frac{\Delta x}{2}\int_{-1}^{1}\psi\frac{\partial f_{h}}{\partial
t}\textnormal{d}x-\int_{-1}^{1}F_{h}\frac{\partial\psi}{\partial
x}\textnormal{d}x+\psi(1)\hat{F}(1)-\psi(-1)\hat{F}(-1)=0.$ (6.13)

In our standard discretization scheme, we would substitute the one-dimensional
$p=1$ orthonormal modal basis functions
$\psi_{j}=\\{1/\sqrt{2},\sqrt{3/2}x\\}$ for the test functions, which results
in

$\displaystyle\frac{\partial f_{0}}{\partial t}=-\frac{\sqrt{2}}{\Delta
x}\left[\hat{F}(1)-\hat{F}(-1)\right]$ (6.14) $\displaystyle\frac{\partial
f_{1}}{\partial t}=-\frac{\sqrt{6}}{\Delta
x}\left[\hat{F}(1)+\hat{F}(-1)\right]+\frac{2\sqrt{3}}{\Delta x}F_{0},$ (6.15)

where we have also expanded $f_{h}$ and $F_{h}$ on the orthonormal modal
basis. In terms of these modal coefficients, we learned above that positivity
requires $f_{0}\geq 0$ and $|f_{1}|\leq\sqrt{3}f_{0}$. This is equivalent to
ensuring that control node values at $x=\pm 1/3$ remain non-negative,
$f_{h}(x=\pm 1/3)\geq 0$.

We would like to find a way to limit the surface and volume terms to ensure
that these constraints are not violated as the solution evolves. Existing
positivity-preserving schemes attempt to limit the boundary fluxes so that the
cell-average $f_{0}$ stays positive (the volume term for the cell-average
always vanishes, so the cell-average is only affected by surface terms).
However, it is not immediately clear how to account for an additional
constraint $|f_{1}|\geq\sqrt{3}f_{0}$, since the cell-slope $f_{1}$ is also
affected by the same fluxes; for this reason, existing schemes often rescale
the cell-slope $f_{1}$ _post-hoc_ , which effectively gives sub-cell diffusion
that can break higher-order conservation laws.

It is more convenient to instead consider the evolution of the control nodes,
$f_{\pm}\equiv f_{h}(x=\pm 1/3)=f_{0}/\sqrt{2}\pm f_{1}/\sqrt{6}$. Taking the
appropriate linear combinations of Eqs. 6.14 and 6.15, we have

$\displaystyle\frac{\partial f_{+}}{\partial t}$
$\displaystyle=-\frac{2}{\Delta x}\hat{F}(1)+\frac{\sqrt{2}}{\Delta x}F_{0}$
(6.16) $\displaystyle\frac{\partial f_{-}}{\partial t}$
$\displaystyle=\frac{2}{\Delta x}\hat{F}(-1)-\frac{\sqrt{2}}{\Delta x}F_{0}.$
(6.17)

Unlike in Eqs. 6.14 and 6.15, each of the control nodes is only affected by a
flux from one side of the cell: the left node $f_{-}$ is affected by the flux
on the left boundary $\hat{F}(-1)$, and the right node $f_{+}$ is affected by
the flux on the right boundary $\hat{F}(1)$.

Neglecting the volume terms for now, this means that we can separately limit
the left flux to maintain $f_{-}\geq 0$ and limit the right flux to maintain
$f_{+}\geq 0$. However, note that the neighboring cells are also affected by
these fluxes. To account for this, let us instead examine the evolution of two
cells, $i$ and $i+1$, due to a flux $\hat{F}^{i+1/2}$ at their interface:

$\displaystyle\frac{\partial f_{+}^{i}}{\partial t}$
$\displaystyle=-\frac{2}{\Delta x}\hat{F}^{i+1/2}+\frac{\sqrt{2}}{\Delta
x}F_{0}^{i}$ (6.18) $\displaystyle\frac{\partial f_{-}^{i+1}}{\partial t}$
$\displaystyle=\frac{2}{\Delta x}\hat{F}^{i+1/2}-\frac{\sqrt{2}}{\Delta
x}F_{0}^{i+1}.$ (6.19)

We see that the flux $\hat{F}^{i+1/2}$ is simply exchanging information
between $f_{+}^{i}$ and $f_{-}^{i+1}$, while neither $f_{-}^{i}$ nor
$f_{+}^{i+1}$ is affected by this flux. This also means that only one of
$f_{+}^{i}$ or $f_{-}^{i+1}$ is decreased by the flux, with the other
increasing by the same amount. Upon adopting a forward Euler timestepping
scheme (which can be built into a higher-order Runge-Kutta scheme), and
dropping the volume terms for now, it is easy to see how to limit the flux
$\hat{F}^{i+1/2}$ so that neither $f_{+}^{i}$ nor $f_{-}^{i+1}$ can become
negative after a single timestep. This gives

$\displaystyle f_{+}^{i}-\frac{2\Delta t}{\Delta x}\hat{F}^{i+1/2}$
$\displaystyle\geq 0$ (6.20) $\displaystyle f_{-}^{i+1}+\frac{2\Delta
t}{\Delta x}\hat{F}^{i+1/2}$ $\displaystyle\geq 0.$ (6.21)

The limit on $\hat{F}^{i+1/2}$ to ensure that the flux does not make either
$f_{+}^{i}$ or $f_{-}^{i+1}$ negative in a single step is then

$\displaystyle-f_{-}^{i+1}\frac{\Delta x}{2\Delta t}\leq\hat{F}^{i+1/2}\leq
f_{+}^{i}\frac{\Delta x}{2\Delta t}.$ (6.22)

This is illustrated in the diagram in Fig. 6.3.

Figure 6.3: Illustration of limiter scheme for boundary fluxes based on Eq.
6.22. The flux limit on the flux between cells $i$ and $i+1$,
$\hat{F}^{i+1/2}$, is violated in the left $(a)$ and right $(c)$ panels,
leading to negative control nodes (red points). In the left panel $(a)$,
$\hat{F}^{i+1/2}$ is too large to the right and control node $f_{+}^{i}$
becomes negative (red). In the right panel $(c)$, $\hat{F}^{i+1/2}$ is too
large to the left, leading to control node $f_{-}^{i+1}$ becoming negative
(red). In the middle panel $(b)$, the flux is within the limits and the
control nodes all remain positive (black points). Also note that the control
points $f_{-}^{i}$ and $f_{+}^{i+1}$ remain fixed in all panels because they
are unaffected by the flux (these points are affected by other fluxes, but not
the one depicted here, $\hat{F}^{i+1/2}$).

Once we have limited all fluxes to ensure that the surface terms cannot make
$f_{\pm}$ negative in any cell in the domain, we can limit the volume terms.
Considering Eqs. 6.16 and 6.17, we can see that the volume terms exchange
information between $f_{-}$ and $f_{+}$ within each cell. Thus one way to
limit the volume terms is to simply scale all the volume terms in each cell by
a common factor $0\leq\theta^{i}\leq 1$ to ensure that neither $f_{-}$ or
$f_{+}$ is made negative by the volume terms.

For cell $i$, the final forward-Euler update can be expressed as

$\displaystyle f^{i}_{+}(t_{n}+\Delta t)$
$\displaystyle=f^{i}_{+}(t_{n})-\frac{2\Delta t}{\Delta
x}\hat{F}^{i+1/2}+\frac{\sqrt{2}\Delta t}{\Delta x}\theta^{i}F_{0}^{i}$ (6.23)
$\displaystyle f^{i}_{-}(t_{n}+\Delta t)$
$\displaystyle=f^{i}_{-}(t_{n})+\frac{2\Delta t}{\Delta
x}\hat{F}^{i-1/2}-\frac{\sqrt{2}\Delta t}{\Delta x}\theta^{i}F_{0}^{i},$
(6.24)

with limits on the fluxes given by Eq. 6.22 and the volume scaling factor
given by

$\theta^{i}=\min\left(1,\frac{f_{+}^{i}-\frac{2\Delta t}{\Delta
x}\hat{F}^{i+1/2}}{\frac{-\sqrt{2}\Delta t}{\Delta
x}F_{0}^{i}},\frac{f_{-}^{i}+\frac{2\Delta t}{\Delta
x}\hat{F}^{i-1/2}}{\frac{\sqrt{2}\Delta t}{\Delta x}F_{0}^{i}}\right).$ (6.25)

Here we have prioritized the surface terms over the volume terms in that we
limit the surface terms first. This can, for example, allow a maximal flux out
the left boundary to lower $f_{-}$ to zero. Then if the volume terms wanted to
decrease $f_{-}$ further, the volume terms would be essentially turned off in
this cell. In principle, one could instead prioritize the volume terms,
allowing the maximum flow within the cell and then possibly turning off
boundary fluxes. A comparison of these two approaches is left to future work.

####  6.3.1 Exponential surface extrapolation

While the scheme described above will rigorously preserve positivity of the
solution, we can make an additional improvement involving how the boundary
fluxes are computed. In a standard DG scheme with upwinded fluxes, the flux
between cells $i$ and $i+1$ would be computed as

$\hat{F}^{i+1/2}=\begin{cases}u^{i+1/2}f_{h}^{i}(x^{i+1/2})\qquad&u^{i+1/2}>0\\\
u^{i+1/2}f_{h}^{i+1}(x^{i+1/2})\qquad&u^{i+1/2}<0\end{cases}$ (6.26)

with $f_{h}^{i}(x^{i+1/2})$ and $f_{h}^{i+1}(x^{i+1/2})$ computed using the
piecewise-linear representation of the solution, and $u^{i+1/2}$ the advection
velocity at the cell interface. After mapping to a unit cell on $x\in[-1,1]$,
the boundary values would be given by

$\displaystyle
f_{h}^{i}(1)=\frac{1}{\sqrt{2}}f_{0}+\frac{\sqrt{3}}{\sqrt{2}}f_{1}$ (6.27)
$\displaystyle
f_{h}^{i+1}(-1)=\frac{1}{\sqrt{2}}f_{0}-\frac{\sqrt{3}}{\sqrt{2}}f_{1}.$
(6.28)

Let us consider an extreme case: a cell where the flux from the left boundary
is zero, with advection velocity $u>0$ a constant. In this case, the modal
coefficients of $f_{h}$ in the cell are given by (from Eqs. 6.14 and 6.15)

$\displaystyle\frac{\partial f_{0}}{\partial t}=-\frac{u\sqrt{2}}{\Delta
x}f_{h}^{i}(1)$ (6.29) $\displaystyle\frac{\partial f_{1}}{\partial
t}=-\frac{u\sqrt{6}}{\Delta x}f_{h}^{i}(1)+\frac{2u\sqrt{3}}{\Delta x}f_{0}.$
(6.30)

From above, we know that we need $|f_{1}|/(\sqrt{3}f_{0})<1$ for the solution
to remain positive and realizable. Thus, let us compute the evolution of
$\bar{x}\equiv f_{1}/(\sqrt{3}f_{0})$:

$\frac{\partial\bar{x}}{\partial t}=\frac{1}{\sqrt{3}f_{0}}\frac{\partial
f_{1}}{\partial t}-\frac{\bar{x}}{f_{0}}\frac{\partial f_{0}}{\partial
t}=\frac{u\sqrt{2}}{\Delta
x}\left[2-\sqrt{2}(1-\bar{x})\frac{f_{h}^{i}(1)}{f_{0}}\right].$ (6.31)

If we use the standard linear extrapolation for $f_{h}^{i}(1)$ from Eq. 6.27,
this gives

$\displaystyle\frac{\partial\bar{x}}{\partial t}=\frac{u}{\Delta
x}\left[2-(1-\bar{x})(1+3\bar{x})\right]$ $\displaystyle=\frac{u}{\Delta
x}\left[1-2\bar{x}+3\bar{x}^{2}\right]$ (6.32) $\displaystyle>0\quad\text{for
all $\bar{x}$}.$

This means that without any limiters, $\bar{x}$ always grows without bound in
this extreme case, which would violate the realizability limit $\bar{x}<1$ in
a finite time. Note also that any reduction or limit on the boundary value
$f^{i}_{h}(1)$ only makes the issue worse, so that $\bar{x}$ increases more
quickly and becomes unphysical sooner. In this case, the volume term is
steepening the slope in the cell faster than the boundary flux can flatten it.
In practice, the volume term limiter in our scheme would eventually prevent
$|x|>1$. Nonetheless, perhaps it would help to enhance the extrapolated
boundary flux, in a way that $\partial\bar{x}/\partial t\rightarrow 0$ as
$\bar{x}\rightarrow 1$. This way, perhaps we wouldn’t need to limit the volume
terms as often or as much.

Figure 6.4: $\partial\bar{x}/\partial t$ vs. $\bar{x}$, given by Eq. 6.32, for
different choices of how to evaluate the surface value $f_{h}^{i}(1)$. A
linear extrapolation, Eq. 6.27, gives $\partial\bar{x}/\partial t>0$ for all
$\bar{x}$, which drives the solution to be unrealizable in finite time
(without limiters). Using an exponential extrapolation instead, given by Eq.
6.33, gives $\partial\bar{x}/\partial t\rightarrow 0$ as $\bar{x}\rightarrow
1$, which means that $\bar{x}$ should stop increasing before becoming
unphysical.

One way to _enhance_ the boundary value is to make use of the exponential
reconstruction given by Eq. 6.9. Extrapolating $g_{h}$ to the right edge of
the cell at $x=1$, we have

$g_{h}(1)=\frac{f_{0}g_{1}}{\sqrt{2}\sinh g_{1}}e^{g_{1}}.$ (6.33)

In Fig. 6.4 we plot ${\partial\bar{x}}/{\partial t}$ for three cases: linear
extrapolation, Eq. 6.27; exact exponential extrapolation, Eq. 6.33 with
$g_{1}(\bar{x})=L^{-1}(\bar{x})$; and approximate exponential extrapolation,
Eq. 6.33 with $g_{1}(\bar{x})\approx\bar{x}(3-\bar{x}^{2})/(1-\bar{x}^{2})$,
_i.e._ using the Cohen approximation, Eq. 6.10, for the inverse Langevin
function. We see that the both the exact and the approximate exponential
extrapolation give $\partial\bar{x}/\partial t\rightarrow 0$ as
$\bar{x}\rightarrow 1$, which means that $\bar{x}$ should stop increasing
before becoming unphysical. The approximate exponential extrapolation gives a
region where $\partial\bar{x}/\partial t<0$, which may in fact make the
algorithm more robust because in this case the equilibrium value is
$\bar{x}\approx 0.6$. Meanwhile, the linear extrapolation gives
$\partial\bar{x}/\partial t>0$ for all $\bar{x}$ as we showed above.

####  6.3.2 Extension to higher dimensionality with $p=1$

To extend the scheme to higher dimensionality, we will again track the
evolution of control nodes, which are given by tensor products of the 1D
control nodes. For example, in 2D, we have four control nodes:
$f_{--}=f_{h}(-1/3,-1/3)$, $f_{-+}=f_{h}(-1/3,1/3)$, $f_{+-}=f_{h}(1/3,-1/3)$,
and $f_{++}=f_{h}(1/3,1/3)$. We will illustrate the scheme for the two-
dimensional case, with extension to higher dimensions relatively
straightforward.

In 2D, the DG weak form from Eq. 6.12 mapped to a cell in $x\in[-1,1],\
y\in[-1,1]$ is given by

$\displaystyle\int_{-1}^{1}\textnormal{d}x\int_{-1}^{1}\textnormal{d}y\,\psi\frac{\partial
f_{h}}{\partial
t}-\int_{-1}^{1}\textnormal{d}x\int_{-1}^{1}\textnormal{d}y\left[\frac{2}{\Delta
x}F_{x\,h}\frac{\partial\psi}{\partial x}+\frac{2}{\Delta
y}F_{y\,h}\frac{\partial\psi}{\partial y}\right]$
$\displaystyle\quad+\frac{2}{\Delta
x}\int_{-1}^{1}\textnormal{d}y\left[\psi(1,y)\hat{F}_{x}(1,y)-\psi(-1,y)\hat{F}_{x}(-1,y)\right]$
(6.34) $\displaystyle\quad+\frac{2}{\Delta
y}\int_{-1}^{1}\textnormal{d}x\left[\psi(x,1)\hat{F}_{y}(x,1)-\psi(x,-1)\hat{F}_{y}(x,-1)\right]=0.$
(6.35)

We can then compute the evolution of the four control nodes as

$\displaystyle\frac{\partial f_{--}}{\partial t}$
$\displaystyle=-\frac{1}{\Delta
x}\left(F_{x0}-\frac{1}{\sqrt{3}}F_{x2}\right)-\frac{1}{\Delta
y}\left(F_{y0}-\frac{1}{\sqrt{3}}F_{y1}\right)+\frac{2}{\Delta
x}\hat{F}_{x}(-1,-\frac{1}{3})+\frac{2}{\Delta y}\hat{F}_{y}(-\frac{1}{3},-1)$
(6.36) $\displaystyle\frac{\partial f_{-+}}{\partial t}$
$\displaystyle=-\frac{1}{\Delta
x}\left(F_{x0}+\frac{1}{\sqrt{3}}F_{x2}\right)+\frac{1}{\Delta
y}\left(F_{y0}-\frac{1}{\sqrt{3}}F_{y1}\right)+\frac{2}{\Delta
x}\hat{F}_{x}(-1,\frac{1}{3})-\frac{2}{\Delta y}\hat{F}_{y}(-\frac{1}{3},1)$
(6.37) $\displaystyle\frac{\partial f_{+-}}{\partial t}$
$\displaystyle=\frac{1}{\Delta
x}\left(F_{x0}-\frac{1}{\sqrt{3}}F_{x2}\right)-\frac{1}{\Delta
y}\left(F_{y0}+\frac{1}{\sqrt{3}}F_{y1}\right)-\frac{2}{\Delta
x}\hat{F}_{x}(1,-\frac{1}{3})+\frac{2}{\Delta y}\hat{F}_{y}(\frac{1}{3},-1)$
(6.38) $\displaystyle\frac{\partial f_{++}}{\partial t}$
$\displaystyle=\frac{1}{\Delta
x}\left(F_{x0}+\frac{1}{\sqrt{3}}F_{x2}\right)+\frac{1}{\Delta
y}\left(F_{y0}+\frac{1}{\sqrt{3}}F_{y1}\right)-\frac{2}{\Delta
x}\hat{F}_{x}(1,\frac{1}{3})-\frac{2}{\Delta y}\hat{F}_{y}(\frac{1}{3},1).$
(6.39)

Figure 6.5: Diagram of fluxes affecting each interior control node in 2D. Each
interior control node is affected by only the fluxes at the nearest surface
control nodes (a dotted line between an interior node and a flux indicates
that the flux affects that node).

Notably, each interior control node is affected only by the fluxes at the
nearest surface control nodes, as shown in the diagram in Fig. 6.5. Similar to
above, we can limit each flux Notably, each interior control node is affected
only by the fluxes at the nearest surface control nodes, as shown in the
diagram in Fig. 6.5. Similar to above, we can limit each flux so that the
affected control nodes cannot become negative on a single forward-Euler
timestep. Focusing on the $f_{++}$ control node, if the fluxes
$\hat{F}_{x}(1,1/3)$ and $\hat{F}_{y}(1/3,1)$ are both directed out of the
cell (as depicted in Fig. 6.5), we need to make sure that the combined flux
does not exceed the limit given by

$\frac{2\Delta t}{\Delta x}\hat{F}_{x}(1,\frac{1}{3})+\frac{2\Delta t}{\Delta
y}\hat{F}_{y}(\frac{1}{3},1)\leq f_{++}.$ (6.40)

We can separately limit the $x$ and $y$ fluxes by apportioning a fraction of
$f_{++}$ that is allowed to removed in each direction, which we will denote as
$\eta_{x}$ and $\eta_{y}$, such that $\eta_{x}+\eta_{y}=1$. Now we can limit
the fluxes as

$\displaystyle\hat{F}_{x}(1,\frac{1}{3})$
$\displaystyle\leq\eta_{x}f_{++}\frac{\Delta x}{2\Delta t}$ (6.41)
$\displaystyle\hat{F}_{y}(\frac{1}{3},1)$
$\displaystyle\leq\eta_{y}f_{++}\frac{\Delta y}{2\Delta t}.$ (6.42)

The definition of the $\eta_{d}$ need not be exact; one possible choice is to
use the ratio given by the contribution to the CFL rate from each direction,
$r_{d}$, divided by the total CFL rate $r$, so that $\eta_{d}=r_{d}/r$.

As above, each flux should be limited by the control nodes on each side of the
boundary. Thus the full limit on $\hat{F}_{x}(x^{i+1/2},1/3)$ at the boundary
between cells $i$ and $i+1$ is given by

$\displaystyle-\eta_{x}^{i+1}f_{-+}^{i+1}\frac{\Delta x}{2\Delta
t}\leq\hat{F}_{x}(x^{i+1/2},\frac{1}{3})\leq\eta_{x}^{i}f_{++}^{i}\frac{\Delta
x}{2\Delta t},$ (6.43)

where note the flux fraction $\eta_{x}$ is computed locally in each cell.
Similar limiter expressions can be given for each flux depicted in Fig. 6.5.

Further, the fluxes can be computed with exponential extrapolation as in
Section 6.3.1. We can compute the exponential extrapolation in one direction
at a time, avoiding the need for a multi-dimensional exponential expression.
For example, to compute the exponential extrapolation for
$\hat{F}_{x}(1,1/3)$, we find the exponential $g(x)$ that is weak-equivalent
to $f(x,1/3)$, and then evaluate the exponential expression at the surface.

Once all the surface terms have been limited to ensure that no control point
can become negative, the volume terms can again be limited by scaling all
volume terms by a common factor $0\leq\theta^{i}\leq 1$ in each cell $i$.
Writing the forward-Euler update of each control node $c$ in cell $i$
generically as

$f_{c}^{i}(t_{n}+\Delta t)=f_{c}^{i}(t_{n})+\Delta tS_{c}^{i}+\theta^{i}\Delta
tV_{c}^{i},$ (6.44)

with $S_{c}^{i}$ and $V_{c}^{i}$ the surface and volume terms, respectively,
the volume scaling factor is given by

$\theta^{i}=\min_{c\,|\,V_{c}^{i}<0}\left(1,\frac{f_{c}^{i}+\Delta
tS_{c}^{i}}{-\Delta tV_{c}^{i}}\right).$ (6.45)

###  6.4 Conservation properties for Hamiltonian systems

Although the positivity-preserving scheme presented above could be used to
solve any kind of hyperbolic conservation law of the form of Eq. 6.11, the
primary targets of our scheme are Hamiltonian systems like gyrokinetics. In
Section 3.2 we showed a DG scheme that conserves energy in Hamiltonian
systems. Now let us apply the positivity-preserving limiters and consider how
the conservation properties are modified, if at all.

Starting from Eq. 3.28, the positivity-preserving evolution of the
distribution function $f$ is given by

$\int_{\mathcal{K}_{i}}\psi\frac{\partial(\mathcal{J}f_{h})}{\partial
t}\textnormal{d}\mbox{\boldmath${Z}$}-\int_{\mathcal{K}_{i}}\theta_{i}\mathcal{J}f_{h}\dot{\mbox{\boldmath${Z}$}}_{h}\boldsymbol{\cdot}\frac{\partial\psi}{\partial\mbox{\boldmath${Z}$}}\,\textnormal{d}\mbox{\boldmath${Z}$}+\oint_{\partial\mathcal{K}_{i}}\psi^{-}\Lambda\left[\widehat{\mathcal{J}f_{h}\dot{\mbox{\boldmath${Z}$}}_{h}}\right]\boldsymbol{\cdot}\textnormal{d}\mbox{\boldmath${s}$}=0,$
(6.46)

where $\theta_{i}$ represents the volume term scaling factor in cell $i$, and
the notation $\Lambda[\hat{F}]$ represents limiters applied to surface fluxes.
To check energy conservation, we first insert the discrete Hamiltonian $H_{h}$
for the test function $\psi$ and sum over cells, giving

$\sum_{i}\int_{\mathcal{K}_{i}}H_{h}\frac{\partial(\mathcal{J}f_{h})}{\partial
t}\textnormal{d}\mbox{\boldmath${Z}$}=\sum_{i}\int_{\mathcal{K}_{i}}\theta_{i}\mathcal{J}f_{h}\dot{\mbox{\boldmath${Z}$}}_{h}\boldsymbol{\cdot}\frac{\partial
H_{h}}{\partial\mbox{\boldmath${Z}$}}\,\textnormal{d}\mbox{\boldmath${Z}$}-\sum_{i}\oint_{\partial\mathcal{K}_{i}}H_{h}^{-}\Lambda\left[\widehat{\mathcal{J}f_{h}\dot{\mbox{\boldmath${Z}$}}_{h}}\right]\boldsymbol{\cdot}\textnormal{d}\mbox{\boldmath${s}$}=0.$
(6.47)

Even with the scaling factor $\theta_{i}$, the volume term vanishes as in the
standard case because
$\dot{\mbox{\boldmath${Z}$}}_{h}\boldsymbol{\cdot}\partial
H_{h}/\partial\mbox{\boldmath${Z}$}=\\{H_{h},H_{h}\\}=0$. The surface term
also vanishes just as in the standard case; the fluxes still exactly cancel at
cell boundaries, even with the flux limiters. For Hamiltonian systems written
in canonical form, the remainder of the energy conservation proof from Section
3.2 is unchanged. However, additional complexities arise in our scheme for the
symplectic formulation of the electromagnetic gyrokinetic system, which
requires the inclusion of limiters in some of the field equations. Extension
of the positivity-preserving scheme to EMGK is left to future work, and we
discuss the difficulties briefly in Appendix 6.A.

###  6.5 Results

In this section we implement the positivity-preserving scheme in Gkeyll and
present some numerical results. We first study passive advection and then we
turn to Hamiltonian systems: the incompressible Euler equations and
electrostatic gyrokinetics.

####  6.5.1 1D advection

Figure 6.6: Positivity-preserving advection of a square pulse for one period
on a 1D periodic domain with 32 cells, overlaid with the initial condition.

Let us first return to the one-dimensional advection example from Section 6.1.
Again taking a square pulse on a periodic domain, Fig. 6.6 shows the results
of the new scheme. Compared to Fig. 6.1, we now see no negative overshoots and
no negative cell averages. In fact, not only do the cell averages remain
positive, but the control nodes in each cell also remain positive, which
ensures that slopes do not become unphysically large. Note however that points
on cell boundaries can be negative, so long as the control nodes are positive,
as can be seen at $x=5/32$.

####  6.5.2 2D advection

Figure 6.7: Two-dimensional square pulse after one period of diagonal
advection on a periodic domain with $32\times 32$ cells using standard DG
(left) and positivity-preserving DG (right). In the standard case, cells with
negative cell average are masked in white.

As a first two-dimensional test, we again consider uniform constant advection
of a square pulse, given initially by

$\displaystyle f(x,y,0)=\begin{cases}&1\qquad|x-x_{c}|<1/4\ \mathrm{and}\
|y-y_{c}|<1/4\\\ &0\qquad\mathrm{otherwise}\end{cases}$ (6.48)

with $x_{c}=y_{c}=1/2$. We advect the solution diagonally, with velocity
components $v_{x}=v_{y}=1$, through a periodic domain with $L_{x}=L_{y}=1$.
Fig. 6.7 shows a comparison of the results from the standard DG scheme and the
positivity-preserving scheme. Both cases use piecewise-linear ($p=1$) basis
functions with 32 cells in each direction. In the standard case, cells with
negative cell average are masked in white. The positivity-preserving scheme
successfully eliminates these negative regions.

####  6.5.3 2D vortex waltz

Figure 6.8: Vorticity at $t=100$ for “vortex waltz” problem on a periodic
domain with $128\times 128$ cells using standard DG (left) and positivity-
preserving DG (right). In the standard case, cells with negative cell-average
vorticity are masked in white.

A more stringent test of our positivity-preserving scheme and its conservation
properties is given by the incompressible Euler system. As we saw in Section
3.2.4, this is a Hamiltonian system, with a conserved energy given by

$\mathcal{E}=\int|\nabla_{\perp}\psi|^{2}\textnormal{d}\mbox{\boldmath${Z}$}.$
(6.49)

In the “vortex waltz” problem (Nielsen et al., 1996), we initialize two
Gaussian vortices which merge as they orbit around each other. The domain is
doubly periodic of dimension $10\times 10$ length units. The initial vorticity
given by

$\displaystyle\varpi(x,y,0)=e^{-r_{1}^{2}/0.8}+e^{-r_{2}^{2}/0.8},$ (6.50)

where $r_{i}^{2}=(x-x_{i})^{2}+(y-y_{i})^{2}$ with $(x_{1},y_{1})=(3.5,5.0)$
and $(x_{2},y_{2})=(6.5,5.0)$ the initial locations of the peaks. We
discretize the system with piecewise-linear basis functions $(p=1)$ on a grid
with $128\times 128$ cells. We show a comparison of the vorticity at $t=100$
from the standard DG scheme and the positivity-preserving scheme in Fig. 6.8,
again masking cells in white that have negative cell-average vorticity.

Figure 6.9: Time traces of the total energy in the system $\mathcal{E}(t)$
normalized to the initial energy for three cases: standard DG, the full
positivity-preserving DG scheme, and a non-conservative version of the
positivity-preserving scheme which uses _post-hoc_ diffusion instead of volume
term limiters. In $(a)$, a standard linear extrapolation is used to evaluate
the surface fluxes, while in $(b)$ the exponential extrapolation from Section
6.3.1 is used. Note the difference in scale between $(a)$ and $(b)$; the
‘standard’ trace (solid) is the same in both plots. The negligible difference
between the conservative and non-conservative traces in $(b)$ shows that using
the exponential extrapolation reduces the need for volume limiters.

To verify that the positivity-preserving scheme has not broken energy
conservation, we show in Fig. 6.9 time traces of the total energy, given by
Eq. 6.49, for three cases: standard DG, the full positivity-preserving scheme,
and a non-conservative positivity scheme. In the non-conservative scheme, the
surface term limiters are still applied as in the conservative scheme, but we
do not apply the volume term limiters. This would keep cell averages positive
but could allow unphysical slopes to develop, so we add _post-hoc_ rescaling
of the slopes at the end of the timestep to maintain realizability, which
breaks energy conservation. Indeed, the plot shows that the standard and
conservative positivity-preserving schemes conserve the energy well, while the
non-conservative scheme has energy errors.

####  6.5.4 5D electrostatic gyrokinetics

Our most challenging test of the positivity algorithm is its application to
the 5D electrostatic gyrokinetic system. Without implementing the positivity
algorithm, the standard DG discretization of the electrostatic gyrokinetic
system results in regions of negative distribution function, leading to
regions of negative density and negative temperature. This can lead to
unphysical behavior in the collision operator222To improve robustness, we have
altered the implementation of the collision operator in the standard version
so that collisions are effectively turned off in cells with negative
temperature, which avoids unphysical anti-diffusion in these cells. This
improves robustness but does not completely eliminate positivity-related
issues in the simulations. and the sheath boundary conditions, resulting in
numerical instabilities.

Figure 6.10: Timeseries of electron density at the midplane ($z=0$) for 5D
collisionless electrostatic gyrokinetic seeded blob simulations with the
standard DG scheme (top row) and the positivity-preserving scheme (bottom
row). In the standard case, cells with negative cell-average density are
masked in white.

As a first test of the positivity algorithm in the electrostatic gyrokinetic
system, we perform a collisionless seeded blob test. We initialize a Gaussian
blob in helical NSTX-like geometry with $L_{z}=100$ m and sheath boundary
conditions at $z=\pm L_{z}/2$. As the blob polarizes it begins to advect
radially outwards and also spin due to the Boltzmann spinning effect (Angus et
al., 2012). With the standard DG discretization, this results in cells with
negative cell-average density, as shown in the top row of Fig. 6.10. With the
positivity-preserving algorithm, these negative cells are eliminated. Fig.
6.11 shows that energy conservation is not altered by the positivity
algorithm, with energy still conserved in the system to
$\sim\mathcal{O}(10^{-5})$.

Figure 6.11: Energy conservation in the electrostatic seeded blob simulations
with the standard algorithm and the positivity-preserving algorithm. The
positivity algorithm does not adversely affect energy conservation, as
desired.

An energy-conserving positivity-preserving DG algorithm for the Dougherty
collision operator has also been formulated and will be presented in future
work. This will enable full electrostatic simulations like the ones presented
in Section 4.2 that maintain positivity of the distribution function.
Implementing the positivity-preserving algorithm in the electromagnetic
gyrokinetic system is somewhat more challenging, as we detail in the Appendix
6.A.

###  6.6 Summary

In this chapter we have developed a discontinuous Galerkin scheme for
maintaining positivity of the distribution function. The scheme has been
carefully constructed to avoid _post-hoc_ diffusion so that conservation
properties are preserved for Hamiltonian systems. The results in Section 6.5
show that the scheme is successful in maintaining positivity for passive
advection, incompressible Euler, and collisionless electrostatic gyrokinetic
systems. Extension to include collisions and electromagnetic effects to the
gyrokinetic system is left as future work.

While the simulations in the bulk of this thesis were able to run somewhat
robustly with the standard DG algorithm (without any assurances of positivity
of the distribution function), there were a number of simulations attempted as
part of this thesis that failed due to positivity issues. For example,
simulations failed when we tried to use a collision frequency that varied in
space and time based on local plasma parameters, because negative local values
of density and temperature resulted in an ill-defined collision frequency. We
also expect the positivity problem to only get worse as we move to more
realistic (and complex) simulation setups and geometries. Thus the work of
this chapter is an important and necessary step towards robust, high-fidelity
simulations.

##  Appendix 6.A Difficulties in extending the positivity scheme to
electromagnetic gyrokinetics

Extending the positivity scheme to the electromagnetic gyrokinetics algorithm
presented in Section 3.3 is challenging, in part because one does not have all
the information needed to compute limiters when the limiters themselves are
needed. To illustrate this, first imagine that we (magically) already know
what all the limiters will be at the beginning of the timestep, such that no
terms (surface or volume) can lead to a negative control node. Neglecting
collisions, the DG discretization of the gyrokinetic equation might look
something like

$\displaystyle\int_{\mathcal{K}_{i}}$
$\displaystyle\psi\frac{\partial(\mathcal{J}f_{h})}{\partial
t}\textnormal{d}^{3}\mbox{\boldmath${R}$}\,\textnormal{d}^{3}\mbox{\boldmath${v}$}$
$\displaystyle\quad-\int_{\mathcal{K}_{i}}\theta_{i}^{H}\mathcal{J}f_{h}\dot{\mbox{\boldmath${R}$}}_{h}\boldsymbol{\cdot}\nabla\psi\,\textnormal{d}^{3}\mbox{\boldmath${R}$}\,\textnormal{d}^{3}\mbox{\boldmath${v}$}-\int_{\mathcal{K}_{i}}\mathcal{J}f_{h}\left(\theta_{i}^{H}\dot{v}^{H}_{\parallel
h}-\theta_{i}^{A}\frac{q}{m}\frac{\partial A_{\parallel h}}{\partial
t}\right)\frac{\partial\psi}{\partial
v_{\parallel}}\textnormal{d}^{3}\mbox{\boldmath${R}$}\,\textnormal{d}^{3}\mbox{\boldmath${v}$}$
$\displaystyle\quad+\oint_{\partial\mathcal{K}_{i}}\psi^{-}\Lambda\left[\widehat{\mathcal{J}f_{h}}\dot{\mbox{\boldmath${R}$}}_{h}\right]\boldsymbol{\cdot}\textnormal{d}\mbox{\boldmath${s}$}_{R}\,\textnormal{d}^{3}\mbox{\boldmath${v}$}+\oint_{\partial\mathcal{K}_{i}}\psi^{-}\Lambda\left[\widehat{\mathcal{J}f_{h}}\left(\dot{v}^{H}_{\parallel
h}-\frac{q}{m}\frac{\partial A_{\parallel h}}{\partial
t}\right)\right]\textnormal{d}^{3}\mbox{\boldmath${R}$}\,\textnormal{d}s_{v}$
$\displaystyle\quad=0.$ (6.51)

Here, $\theta_{i}^{H}$ and $\theta_{i}^{A}$ represent volume term scaling
factors in cell $i$ for the Poisson bracket and inductive volume terms,
respectively, and the notation $\Lambda[\hat{F}]$ represents limiters applied
to surface fluxes.

We would like to ensure that this limiter scheme preserves energy
conservation. While the volume term limiters $\theta^{H}$ on the Poisson
bracket terms do not affect energy conservation, we have an additional volume
term outside the bracket involving $\partial A_{\parallel h}/\partial t$. We
also have modifications to the surface terms. In order to maintain energy
conservation, we must account for these limiters in the field equations. To
see this, take $\psi=H_{s\,h}$ to compute

$\displaystyle\sum_{s,i}$
$\displaystyle\int_{\mathcal{K}_{i}}H_{s\,h}\frac{\partial(\mathcal{J}f_{s\,h})}{\partial
t}\,\textnormal{d}^{3}\mbox{\boldmath${R}$}\,\textnormal{d}^{3}\mbox{\boldmath${v}$}=\sum_{s,i}\int_{\mathcal{K}_{i}}\mathcal{J}f_{s\,h}\,\theta_{i}^{H}\left(\dot{\mbox{\boldmath${R}$}}_{h}\mbox{\boldmath${\boldsymbol{\cdot}}$}\nabla
H_{s\,h}+\dot{v}^{H}_{\parallel h}\frac{\partial H_{s\,h}}{\partial
v_{\parallel}}\right)\textnormal{d}^{3}\mbox{\boldmath${R}$}\,\textnormal{d}^{3}\mbox{\boldmath${v}$}$
$\displaystyle-\sum_{s,i}\int_{\mathcal{K}_{i}}\mathcal{J}f_{s\,h}\,\theta^{A}_{i}\frac{q_{s}}{m_{s}}\frac{\partial
A_{\parallel h}}{\partial t}\frac{\partial H_{s\,h}}{\partial
v_{\parallel}}\,\textnormal{d}^{3}\mbox{\boldmath${R}$}\,\textnormal{d}^{3}\mbox{\boldmath${v}$}-\sum_{s,i}\oint_{\partial\mathcal{K}_{i}}H_{s\,h}^{-}\Lambda\left[\widehat{\mathcal{J}f_{h}}\dot{\mbox{\boldmath${R}$}}_{h}\right]\boldsymbol{\cdot}\textnormal{d}\mbox{\boldmath${s}$}_{R}\,\textnormal{d}^{3}\mbox{\boldmath${v}$}$
$\displaystyle\quad-\sum_{s,i}\oint_{\partial\mathcal{K}_{i}}H_{s\,h}^{-}\Lambda\left[\widehat{\mathcal{J}f_{h}}\left(\dot{v}^{H}_{\parallel
h}-\frac{q}{m}\frac{\partial A_{\parallel h}}{\partial
t}\right)\right]\textnormal{d}^{3}\mbox{\boldmath${R}$}\,\textnormal{d}s_{v}.$
(6.52)

As noted above, despite the inclusion of the volume limiter $\theta^{H}$, the
first volume term still vanishes exactly because
$\dot{\mbox{\boldmath${R}$}}_{h}\mbox{\boldmath${\boldsymbol{\cdot}}$}\nabla
H_{h}+\dot{v}^{H}_{\parallel h}\partial H_{h}/\partial
v_{\parallel}=\\{H_{h},H_{h}\\}_{h}=0$. The surface terms also cancel exactly
at cell interfaces, so we are left with

$\displaystyle\sum_{s,i}\int_{\mathcal{K}_{i}}H_{s\,h}\frac{\partial(\mathcal{J}f_{s\,h})}{\partial
t}\,\textnormal{d}^{3}\mbox{\boldmath${R}$}\,\textnormal{d}^{3}\mbox{\boldmath${v}$}$
$\displaystyle=-\sum_{s,i}\int_{\mathcal{K}_{i}}\mathcal{J}f_{s\,h}\,\theta^{A}_{i}\frac{q_{s}}{m_{s}}\frac{\partial
A_{\parallel h}}{\partial t}\frac{\partial H_{s\,h}}{\partial
v_{\parallel}}\,\textnormal{d}^{3}\mbox{\boldmath${R}$}\,\textnormal{d}^{3}\mbox{\boldmath${v}$}$
$\displaystyle=-\int_{\mathcal{T}^{R}}\frac{\partial A_{\parallel h}}{\partial
t}\tilde{J}_{\parallel h}\,\,\textnormal{d}^{3}\mbox{\boldmath${R}$},$ (6.53)

where we will define a _limited_ parallel current, denoted by
$\tilde{J}_{\parallel h}$, as

$\tilde{J}_{\parallel
h}=\sum_{s,j}\frac{q_{s}}{m_{s}}\int_{\mathcal{K}^{v}_{j}}\theta^{A}_{j}\frac{\partial
H_{s\,h}}{\partial
v_{\parallel}}\mathcal{J}f_{s\,h}\,\textnormal{d}^{3}\mbox{\boldmath${v}$}.$
(6.54)

To regain energetic consistency, this limited $\tilde{J}_{\parallel h}$ must
be used in Ampère’s law, which becomes

$\int_{\mathcal{K}^{R}_{i}}\nabla_{\perp}A_{\parallel
h}\mbox{\boldmath${\boldsymbol{\cdot}}$}\nabla_{\perp}\varphi^{(i)}\,\textnormal{d}^{3}\mbox{\boldmath${R}$}-\oint_{\partial\mathcal{K}^{R}_{i}}\varphi^{(i)}\nabla_{\perp}A_{\parallel
h}\boldsymbol{\cdot}\textnormal{d}\mbox{\boldmath${s}$}_{R}=\mu_{0}\int_{\mathcal{K}^{R}_{i}}\varphi^{(i)}\
\tilde{J}_{\parallel h}\,\textnormal{d}^{3}\mbox{\boldmath${R}$}.$ (6.55)

Now we can insert $\varphi^{(i)}=(1/\mu_{0})\partial A_{\parallel h}/\partial
t$ to compute

$\displaystyle\frac{\partial\mathcal{E}_{B\,h}}{\partial t}$
$\displaystyle=\sum_{i}\int_{\mathcal{K}^{R}_{i}}\frac{1}{\mu_{0}}\nabla_{\perp}A_{\parallel
h}\mbox{\boldmath${\boldsymbol{\cdot}}$}\nabla_{\perp}\frac{\partial
A_{\parallel h}}{\partial
t}\,\textnormal{d}^{3}\mbox{\boldmath${R}$}=\int_{\mathcal{T}^{R}}\frac{\partial
A_{\parallel h}}{\partial t}\tilde{J}_{\parallel
h}\,\textnormal{d}^{3}\mbox{\boldmath${R}$},$ (6.56)

which now cancels the term leftover from Eq. 6.53. Now to derive the self-
consistent Ohm’s law, we take the time derivative of Eq. 6.55, giving

$\displaystyle\int_{\mathcal{K}^{R}_{i}}\nabla_{\perp}\frac{\partial
A_{\parallel h}}{\partial
t}\mbox{\boldmath${\boldsymbol{\cdot}}$}\nabla_{\perp}\varphi^{(i)}\,\textnormal{d}^{3}\mbox{\boldmath${R}$}-\oint_{\partial\mathcal{K}^{R}_{i}}\varphi^{(i)}\nabla_{\perp}\frac{\partial
A_{\parallel h}}{\partial
t}\boldsymbol{\cdot}\textnormal{d}\mbox{\boldmath${s}$}_{R}$
$\displaystyle\quad+\mu_{0}\sum_{s}q_{s}\int_{\mathcal{K}_{i}^{R}}\varphi^{(i)}\left[\sum_{j}\oint_{\partial\mathcal{K}^{v}_{j}}\theta^{A}_{j}\bar{v}_{\parallel}^{-}\Lambda\left[\widehat{\mathcal{J}f_{s\,h}}\left({\dot{v}^{H}_{\parallel
h}}-\frac{q_{s}}{m_{s}}\frac{\partial A_{\parallel h}}{\partial
t}\right)\right]\textnormal{d}s_{v}\right]\textnormal{d}^{3}\mbox{\boldmath${R}$}$
$\displaystyle=\mu_{0}\sum_{s}q_{s}\int_{\mathcal{K}^{R}_{i}}\varphi^{(i)}\Bigg{[}\sum_{j}\int_{\mathcal{K}^{v}_{j}}\theta^{A}_{j}\bar{v}_{\parallel}\frac{\partial(\mathcal{J}f_{s\,h})}{\partial
t}^{\star}\textnormal{d}^{3}\mbox{\boldmath${v}$}\Bigg{]}\textnormal{d}^{3}\mbox{\boldmath${R}$},$
(6.57)

where where we have assumed $p=1$, and

$\displaystyle\int_{\mathcal{K}_{i}}\psi\frac{\partial(\mathcal{J}f_{h})}{\partial
t}^{\star}\textnormal{d}^{3}\mbox{\boldmath${R}$}\,\textnormal{d}^{3}\mbox{\boldmath${v}$}=$
$\displaystyle\quad\int_{\mathcal{K}_{i}}\theta^{H}_{i}\mathcal{J}f_{h}\left(\dot{\mbox{\boldmath${R}$}}_{h}\boldsymbol{\cdot}\nabla\psi+\dot{v}^{H}_{\parallel
h}\frac{\partial\psi}{\partial
v_{\parallel}}\right)\textnormal{d}^{3}\mbox{\boldmath${R}$}\,\textnormal{d}^{3}\mbox{\boldmath${v}$}-\oint_{\partial\mathcal{K}_{i}}\psi^{-}\Lambda\left[\widehat{\mathcal{J}f_{h}}\dot{\mbox{\boldmath${R}$}}_{h}\right]\boldsymbol{\cdot}\textnormal{d}\mbox{\boldmath${s}$}_{R}\,\textnormal{d}^{3}\mbox{\boldmath${v}$}.$
(6.58)

This limiter-modified Ohm’s law presents a number of challenges. First, the
surface limiters on the second line in Eq. 6.57 make the problem of solving
for $\partial A_{\parallel h}/\partial t$ a nonlinear one. These limiters act
in phase-space, not real-space, which introduces additional degrees of
freedom. And these are issues even when all the limiters are known at the time
of the solve. In practice, there is the additional complication that all the
limiters (the surface limiters and the volume limiters $\theta^{A}$)
themselves depend on $\partial A_{\parallel h}/\partial t$. Thus we
essentially need to know $\partial A_{\parallel h}/\partial t$ to evaluate
limiters in Ohm’s law in order to solve for $\partial A_{\parallel h}/\partial
t$. Inevitably, one will require an iteration scheme to solve this circular
problem, although it is difficult to say whether such a scheme would converge
quickly, if at all.

##  Appendix 6.B Generalization of positivity constraints to higher polynomial
order and dimensionality

In this section we consider a rigorous procedure for tractably evaluating the
positivity constraints in higher dimensionality and higher polynomial order. A
key result of Section 6.2 was that in 1D, the most-extreme realizable solution
has $f_{1}=\pm\sqrt{3}f_{0}$. The corresponding exponential solution has
$g_{1}\rightarrow\pm\infty$, so that the exponential approaches a delta
function at the cell boundary.

Another way to obtain this result is to project a delta function evaluated at
the cell boundary (or more precisely, just inside the boundary) onto the
piecewise-linear modal basis,

$\displaystyle\delta_{0}=\int_{-1}^{1}\frac{1}{\sqrt{2}}\delta(x-\pm
1)\,\textnormal{d}x=\frac{1}{\sqrt{2}}$ (6.59)
$\displaystyle\delta_{1}=\int_{-1}^{1}\frac{\sqrt{3}}{\sqrt{2}}x\delta(x-\pm
1)\,\textnormal{d}x=\pm\frac{\sqrt{3}}{\sqrt{2}},$ (6.60)

so that the delta function on the boundary has a weak-equivalent piecewise-
linear representation given by

$\delta(x-\pm
1)\doteq\frac{1}{\sqrt{2}}\delta_{0}+\frac{\sqrt{3}}{\sqrt{2}}x\delta_{1}=\frac{1}{2}\pm\frac{3}{2}x\equiv\delta^{\pm}(x).$
(6.61)

Indeed, $\delta_{1}=\pm\sqrt{3}\delta_{0}$, so we have recovered the earlier
result. Also note that the functions $\delta^{\pm}(x)$ have zeros at the
positivity control nodes $x=\pm 1/3$.

Now consider that we could use the functions $\delta^{\pm}(x)$ as basis
functions and expand the solution as

$f_{h}=f_{+}\delta^{+}+f_{-}\delta^{-}=f_{+}\left(\frac{1}{2}+\frac{3}{2}x\right)+f_{-}\left(\frac{1}{2}-\frac{3}{2}x\right).$
(6.62)

This is effectively a nodal basis, with nodal values $f_{\pm}=f_{h}(x=\pm
1/3)$. Note that we can also obtain the coefficients $f_{\pm}$ via projection:

$f_{\pm}=\int_{-1}^{1}\left(\frac{1}{2}\pm\frac{1}{2}x\right)f_{h}\,\textnormal{d}x\equiv\int_{-1}^{1}\delta_{\pm}(x)f_{h}\,\textnormal{d}x.$
(6.63)

Thus we will denote the functions $\delta^{\pm}(x)$ as the positivity
_expansion_ basis functions, and the functions $\delta_{\pm}(x)$ as the
positivity _projection_ basis functions. Finally, as before, the positivity
constraint for the $f_{\pm}$ coefficients is $f_{\pm}\geq 0$.

Note that we have now slightly modified our positivity definition. Instead of
using an exponential basis as the non-polynomial positive-definite basis set
with which we require weak-equality, we will now use delta functions. Thus the
positive-definite representation of the solution is

$g_{h}=g_{+}\delta(x-1)+g_{-}\delta(x+1).$ (6.64)

Requiring this non-polynomial function to be positive on the entire cell
domain gives the constraint $g_{\pm}\geq 0$. Enforcing weak-equality with the
piecewise-linear solution $f_{h}=f_{+}\delta^{+}+f_{-}\delta^{-}$ now simply
gives $f_{\pm}=g_{\pm}$, since by construction the basis functions
$\delta^{\pm}$ are weak-equivalent to the delta functions $\delta(x-\pm 1)$.
This once again gives that the positivity constraints on the coefficients of
$f_{h}$ are $f_{\pm}\geq 0$. Thus despite the change in positive-definite
basis functions from exponentials to delta functions, the result is the same,
which suggests some degree of equivalence between the two choices.

Now consider the one-dimensional piecewise-quadratic case. The orthonormal
modal basis set on the cell $x\in[-1,1]$ is

$\psi=\left\\{\frac{1}{\sqrt{2}},\ \frac{\sqrt{3}}{\sqrt{2}}x,\
\frac{3\sqrt{5}}{2\sqrt{2}}\left(x^{2}-\frac{1}{3}\right)\right\\}.$ (6.65)

Once again, we can project delta functions just inside the cell boundaries,
$\delta(x-\pm 1)$, onto the basis, giving

$\displaystyle\delta^{\pm}=\frac{1}{2}\pm\frac{3}{2}x+\frac{15}{4}\left(x^{2}-\frac{1}{3}\right)\doteq\delta(x-\pm
1).$ (6.66)

We will again use these functions as expansion basis functions. Given the
extra degree of freedom in the piecewise-quadratic case, we need an additional
basis function. We can obtain the final basis function by taking a delta
function at the cell center, $\delta(x-0)$; we choose $x=0$ so that the basis
is symmetric about the cell center. Projecting onto the piecewise-quadratic
basis gives

$\delta^{\circ}=\frac{1}{2}-\frac{15}{8}\left(x^{2}-\frac{1}{3}\right)\doteq\delta(x-0).$
(6.67)

Thus, the piecewise-quadratic positivity expansion basis is given by
$\\{\delta^{+},\delta^{\circ},\delta^{-}\\}$, allowing us to expand the
solution as

$f_{h}=f_{+}\delta^{+}+f_{\circ}\delta^{\circ}+f_{-}\delta^{-}.$ (6.68)

Unlike in the piecewise-linear case, the coefficients
$\\{f_{+},f_{\circ},f_{-}\\}$ do not coincide with nodal values of $f_{h}$.
Instead, they must be obtained by using projection basis functions, which can
be shown to be

$\delta_{+}=\frac{1}{2}x(x+1),\qquad\delta_{\circ}=1-x^{2},\qquad\delta_{-}=\frac{1}{2}x(x-1),$
(6.69)

so that

$\displaystyle f_{+}=\int_{-1}^{1}\delta_{+}f_{h}\,\textnormal{d}x,\qquad
f_{\circ}=\int_{-1}^{1}\delta_{\circ}f_{h}\,\textnormal{d}x,\qquad
f_{-}=\int_{-1}^{1}\delta_{-}f_{h}\,\textnormal{d}x.$ (6.70)

Once again, the requirements
$g_{h}=g_{+}\delta(x-1)+g_{\circ}\delta(x)+g_{-}\delta(x+1)\geq 0$ (for all
$x\in[-1,1]$) and $f_{h}\doteq g_{h}$ give the positivity constraints that the
coefficients of
$f_{h}=f_{+}\delta^{+}+f_{\circ}\delta^{\circ}+f_{-}\delta^{-}$ must be non-
negative: $f_{\pm},f_{\circ}\geq 0$.

We have now successfully generalized the procedure for evaluating positivity
constraints to $p=2$. We can further generalize to arbitrarily high order by
making the following two observations. First, consider that we obtained the
positivity expansion basis functions above by projecting delta functions
centered at $x=-1,1$ for $p=1$ and $x=-1,0,1$ for $p=2$. We can recognize that
these sets of points are Gauss-Lobatto nodes333The Gauss-Lobatto nodes always
include the cell endpoints $x=\pm 1$ in the node set. Variants include the
Legendre-Gauss-Lobatto nodes (commonly referred to as just the Gauss-Lobatto
nodes), where the nodes of order $p$ are given by roots of the polynomial
$P^{\prime}_{p}(x)$, where $P(x)$ is a Legendre polynomial; and the Chebyshev-
Gauss-Lobatto nodes, where the nodes are located at $x_{j}=\cos(\pi j/p)$ for
$j=0,...,p$. For $p\leq 2$, these variants give identical nodes. for $p=1$ and
$p=2$, respectively. Second, note that the positivity projection functions are
the Lagrange basis functions for the same set of Gauss-Lobatto nodes. For
arbitrary $p$, the Lagrange basis functions for nodes $x_{j}$ are given by

$\ell_{j}(x)=\prod_{\begin{smallmatrix}0\leq m\leq p\\\ m\neq
j\end{smallmatrix}}\frac{x-x_{m}}{x_{j}-x_{m}}.$ (6.71)

Thus we now have a general procedure for generating positivity expansion and
projection basis functions in one-dimension for arbitrary polynomial order.
The steps are

  1. 1.

The positivity _expansion_ basis functions can be found by projecting delta
functions centered at Gauss-Lobatto nodes $x_{j}$ onto the orthonormal modal
basis $\mbox{\boldmath${\psi}$}(x)$,

$\displaystyle\ell^{j}(x)=\mbox{\boldmath${\psi}$}(x)\boldsymbol{\cdot}\int_{-1}^{1}\mbox{\boldmath${\psi}$}(x)\
\delta(x-x_{j})\,\textnormal{d}x=\mbox{\boldmath${\psi}$}(x)\boldsymbol{\cdot}\mbox{\boldmath${\psi}$}(x_{j}).$
(6.72)

We can then use this basis to expand $f_{h}$ as

$\displaystyle f_{h}=\sum_{j=0}^{p}g_{j}\ell^{j}(x),$ (6.73)

where we will now use $g_{j}$ to denote the coefficients of the solution
expanded on the positivity basis (to distinguish from modal coefficients
$f_{j}$).

  2. 2.

The positivity _projection_ basis functions are the Lagrange basis functions
for the Gauss-Lobatto nodes:

$\ell_{j}(x)=\prod_{\begin{smallmatrix}0\leq m\leq p\\\ m\neq
j\end{smallmatrix}}\frac{x-x_{m}}{x_{j}-x_{m}}.$ (6.74)

Note that in general, the Lagrange basis functions for a particular set of
nodes $x_{j}$ can be derived by computing the matrix

$M_{jk}=\int_{-1}^{1}\psi_{j}\ell^{k}\,\textnormal{d}x=\psi_{j}(x_{k}),$
(6.75)

where the second equality assumes that the $\psi_{j}$ are orthonormal. Then
the Lagrange basis functions are given by

$\ell_{j}=\sum_{k}(M^{-1})_{jk}\psi_{k}.$ (6.76)

These projection basis functions can be used to find the coefficients

$g_{j}=\int_{-1}^{1}\ell_{j}(x)f_{h}\,\textnormal{d}x,$ (6.77)

since

$\int_{-1}^{1}\ell_{i}\ell^{j}\,\textnormal{d}x=\delta_{ij}.$ (6.78)

The $g_{j}$ are effectively the projection of $f_{h}$ onto a Gauss-Lobatto
nodal (Lagrange) basis. Note however that this is not the same as directly
evaluating $f_{h}$ at the Gauss-Lobatto nodes.

  3. 3.

The piecewise-polynomial solution $f_{h}$ is positive _in the weak sense_ if
all coefficients $g_{j}$ from Eq. 6.77 are non-negative, so that $g_{j}\geq 0$
for all $j$.

The above procedure can then be further generalized to higher dimension by
taking (possibly sparse) tensor products of the one-dimensional positivity
expansion and projection basis functions.

##  Chapter 7 Summary and future work

###  7.1 Summary

The main advance of this thesis was the development of the first capabilities
for simulating electromagnetic gyrokinetic turbulence on open magnetic field
lines. This is an important step towards comprehensive electromagnetic
gyrokinetic simulations of the coupled edge/SOL system. In the past, including
electromagnetic effects in gyrokinetic codes has been challenging, as there
are delicate issues such as the Ampère cancellation problem that must be
handled properly. In our continuum full-$f$ approach, we build on the
successes of continuum $\delta f$ gyrokinetic codes in the core which have
mostly avoided the cancellation problem. The inclusion of electromagnetic
effects in gyrokinetic simulations that can handle the unique challenges of
the boundary plasma (large fluctuations, open and closed field line regions,
_etc._) is critical to the understanding of phenomena such as edge-localized
modes and the pedestal, for which electromagnetic dynamics are expected to
play a key role.

In Chapter 2 we gave a first-principles derivation of the electromagnetic
gyrokinetic system, in the limit of interest for our present work. This
derivation used phase-space-Lagrangian Lie perturbation methods to
systematically derive a self-consistent, energy-conserving, and global
gyrokinetic system, including electromagnetic perturbations. We used the weak-
flow ordering, which simultaneously allows large perturbations $q\Phi/T\sim 1$
at long wavelengths ($k_{\perp}\rho\sim\epsilon_{V}\ll 1$) and small
perturbations $q\Phi/T\sim\epsilon_{V}$ at short wavelengths
($k_{\perp}\rho\sim 1$), along with perturbations at intermediate scales. We
also used the symplectic ($v_{\parallel}$) formulation of electromagnetic
gyrokinetics, which results in the explicit presence of the inductive electric
field in the gyrokinetic equation. After deriving the general formalism
including finite-Larmor-radius (FLR) corrections, we consistently reduced the
system to the long-wavelength limit by neglecting first- and second-order
terms in the single particle Lagrangian to obtain the guiding-center
Lagrangian, which contains no gyroaverages. Variational derivation of the
field equations resulted in a self-consistent, energy-conserving system for
electromagnetic gyrokinetics in the long-wavelength limit. We take this limit
for simplicity of implementation in the Gkeyll code, with extension of the
implementation to include FLR terms left as future work. We summarized the
system implemented in Gkeyll in Section 2.3.

We went to great lengths to ensure that the underlying system is self-
consistent and conservative, so we also needed a robust numerical method with
a discretization scheme that preserves these properties. This was the topic of
Chapter 3. We have employed the discontinuous Galerkin method, a high-order
numerical method that combines attractive features of finite-element and
finite-volume methods. After discussing a discontinuous Galerkin scheme for
general Hamiltonian systems that preserves energy by design (in the
continuous-time limit), we applied the scheme to the electromagnetic
gyrokinetic system. The scheme was then implemented in the gyrokinetic module
of the Gkeyll plasma simulation framework. Linear benchmarks were shown to
verify the implementation. The success of these benchmarks, especially for
cases with high $\beta$ and small $k_{\perp}\rho$, indicated that the Ampère
cancellation problem is avoided. We confirmed this by deriving a semi-discrete
Alfvén wave dispersion relation. As a result, we can handle electromagnetic
fluctuations in a stable, robust, and efficient manner.

The success of the scheme led to the first published simulations of
electromagnetic gyrokinetic turbulence on open field lines, detailed in
Chapter 4. As a rough model of the scrape-off layer in the National Spherical
Torus Experiment (NSTX) experiment at PPPL, we took a simple helical
configuration (like a simple magnetized torus, or SMT) with field lines
wrapping helically around the torus and terminating on conducting plates at
the top and bottom. This model system contains many of the necessary
ingredients for SOL dynamics, including bad curvature and Debye sheath
effects, which are handled via conducting-sheath boundary conditions. Initial
results showed that when electromagnetic effects are included, high $\beta$
blobs can bend and stretch the magnetic field lines as they move radially
outwards in the SOL. Qualitative comparisons to a corresponding electrostatic
simulation showed differences in blob dynamics, with non-adiabatic electron
dynamics playing a key role in the electromagnetic case due to slowing of the
parallel response. We then performed a study of the effects of increasing
$\beta$ on the SOL dynamics. At higher $\beta$, the influence of
electromagnetic effects became stronger, resulting in steepening of pressure
gradients near the source region and flattening of gradients in the remainder
of the domain. We observed a transition from interchange-like modes with
$k_{\parallel}\sim 0$ to ballooning-like modes with finite $k_{\parallel}$ as
pressure gradients $(\alpha^{\mathrm{SMT}})$ increased above the ballooning
stability threshold in the source region. Radially inward magnetic flutter
particle transport off midplane, resulting from parallel motion of electron
along radially-bowed-out field lines, was observed to increase roughly as
$\beta^{2}$. Meanwhile the $E\times B$ component of the radial particle
transport only scaled linearly with $\beta$. This led to slightly reduced
radial transport in the high $\beta$ electromagnetic cases, resulting in
slightly higher peak particle and heat loads on the end plates compared to
corresponding electrostatic cases. These results could have important
implications for the transport of high $\beta$ blobs and ELM filaments.
Further, the electromagnetic mechanism resulting in the steepening of
gradients in the source region could have implications for pedestal formation
and thus deserves more thorough study. Crucially, our electromagnetic
simulations were not significantly more expensive than corresponding
electrostatic simulations, which should allow the routine inclusion of
electromagnetic effects in future results.

We worked on advancing to more realistic SOL geometry in Chapter 5. We adopted
a generalized field-aligned non-orthogonal coordinate system, and expressed
the gyrokinetic system in these coordinates. While these coordinates break
down at the separatrix in diverted geometries due to a singularity at the
X-point, field-aligned coordinates could still be used for efficient
discretization on either side of the separatrix, stitched to a non-aligned
domain in the near vicinity of the separatrix. We then focused on how to
formulate field-aligned coordinate systems for use in flux-tube-like domains
in the SOL. We started with a helical configuration with magnetic shear, which
generalizes the simple geometry used in Chapter 4. Preliminary electrostatic
simulations in this configuration showed that transport is reduced in more
sheared configurations. We then formulated field-aligned coordinate systems
based on an analytical Solov’ev model SOL equilibrium and an analytical
concentric circular equilibrium. The latter is a common geometry used in the
core region, especially for inter-code benchmarking. We presented results from
a preliminary electrostatic ITG benchmark based on the Cyclone base case in
circular core geometry, which compared well with results from other codes in
the long-wavelength limit where our system is valid.

In Chapter 6 we tackled the problem of positivity in our discontinuous
Galerkin scheme. Simulations can suffer from accuracy and robustness issues
because the standard DG scheme does not guarantee that the distribution
function will remain positive (even in the cell-average). We developed a novel
scheme for both defining and preserving positivity in the DG discretization.
Importantly, the scheme was designed without _post-hoc_ diffusion that is used
in many existing positivity-preserving algorithms. This allows the scheme to
preserve energy conservation while maintaining positivity, even in Hamiltonian
systems like gyrokinetics where energy conservation relies on higher-order
moments. We then implemented the scheme in Gkeyll and performed a variety of
numerical tests for advection, the 2D incompressible Euler system, and the 5D
electrostatic gyrokinetic system. The success of the scheme in maintaining
positivity and preserving energy conservation, even in 5D, is a significant
advance.

###  7.2 Future work

While the work of this thesis has advanced the modeling capabilities of the
Gkeyll code, there are a number of areas that remain in order to produce
realistic results for direct comparison with existing experiments or
prediction of future ones. The following list focuses on enhancements
requiring further code and algorithm development, many of which are already in
progress.

  * •

Closed-field-line boundary conditions: The field-aligned geometry formulation
presented in Chapter 5 can also be used for closed-field-line regions. What
remains is the implementation of a boundary condition for closed-field-line
regions. This work is currently underway, led by Mana Francisquez, using the
twist-and-shift approach (Beer et al., 1995). This requires careful
interpolation of mis-aligned sheared grids at the ends of the domain along the
field line. Once ready, this will allow simulations in a limiter configuration
containing both open and closed field line regions. This configuration has
been used by several fluid and gyrofluid codes (Ribeiro & Scott, 2008; Halpern
& Ricci, 2017; Francisquez et al., 2017), and can be used to study SOL flows,
the edge radial electric field, and resulting edge toroidal rotation. These
are all relevant to pedestal formation and the L-H transition. Parra & Catto
(2008, 2010) have stressed the importance of third-order terms in the
Hamiltonian that are required to accurately calculate toroidal rotation. While
these subtleties must be investigated in detail, Parra & Catto’s results are
for the low-flow, up-down symmetric, gyro-Bohm regime. In the edge, the gyro-
Bohm scaling breaks down because eddy sizes are not much smaller than radial
gradient scale lengths, and so including the third-order Hamiltonian terms may
not be required for studying rotation mechanisms in the edge.

  * •

Diverted geometry with X-point: As we have discussed, the X-point is a
significant challenge because field-aligned coordinates are singular on the
separatrix. This has led to a number of new approaches, such as the flux-
coordinate-independent (FCI) approach (Hariri & Ottaviani, 2013; Hariri et
al., 2014; Stegmeir et al., 2016), which move away from the field-aligned
approach. Implementation of FCI or a related approach near the X-point could
allow simulation of diverted geometries. Ideally, one could still use
conventional field-aligned domains in the core and in the SOL, and only use a
non-aligned domain in the immediate vicinity of the separatrix. Stitching
these domains together will require sophisticated interpolation and mapping
schemes, especially if conservation laws are to be preserved.

  * •

Neutral modeling: Neutral interactions play a significant role in plasma-
material interactions that dictate much of the SOL dynamics and evolution. As
such, modeling neutrals is critical to producing experimentally-relevant
results and predictions. Neutral modeling work is underway in Gkeyll, led by
Tess Bernard, leveraging the existing 6D Vlasov kinetic module to produce a
kinetic Boltzmann neutral model. The main interaction mechanisms of electron-
impact ionization, charge exchange, and radiative recombination are modeled.

  * •

Gyroaveraging and higher order terms: While we derived a long-wavelength limit
of the gyrokinetic system, it is important to generalize to shorter
wavelengths $k_{\perp}\rho\sim 1$ within the weak-flow ordering. This involves
gyroaveraging operations in the gyrokinetic equation and the field equations.
Gyroaveraging is relatively simple in the Fourier spectral representation of
many core gyrokinetic codes, where simply multiplying by the Bessel function
$J_{0}(k_{\perp}v_{\perp}/\Omega)$ gives gyroaveraging; however, in real space
implementations, gyroaveraging requires integral operations that sample around
the gyro-orbit. The finite-element implementation of Maurer et al. (2020) is
likely a good starting point for a gyroaveraging implementation in Gkeyll.
Additionally, the second-order $E\times B$ energy term in the Hamiltonian
should be included, so that a time-evolving density can be used in the
polarization term in the Poisson equation instead of the linearized
polarization used in this work. This will be important for cases like pedestal
formation where there is significant evolution of the density profile. With
these additions, we could use the system given in Case 1 from Section 2.2.2.

  * •

More realistic/efficient collision operators: We have used a model Dougherty
collision operator in this work, and we have taken a constant-in-space and
constant-in-time collision frequency. A time- and spatially-varying collision
frequency has been implemented in Gkeyll, but it suffers from robustness
issues likely related to positivity. We have also used an artificially-reduced
collision frequency in this work to avoid severe timestep restrictions; this
issue could be alleviated with an implicit or super-timestepping
implementation of the collision operator. Further, a more realistic collision
operator beyond the simple Dougherty model should be implemented. Preliminary
work on a full nonlinear Fokker-Planck collision model in Rosenbluth potential
form in Gkeyll has been led by Petr Cagas.

  * •

Porting Gkeyll to GPUs: Today, many of the world’s fastest supercomputers all
derive a majority of their computing power from graphics processors (GPUs).
The rise of GPUs in scientific computing over the past decade has been driven
in large part due to their supreme performance for machine learning
applications. To fully leverage the power of these machines, the algorithms in
Gkeyll must be efficiently implemented on GPUs. Work has begun to port the
compute-intensive Gkeyll solver kernels to a CUDA implementation for use on
NVIDIA GPUs, with significant progress made by myself, Ammar Hakim, Jimmy
Juno, Mana Francisquez, and others on the Gkeyll team as part of GPU
hackathons hosted by Princeton.

  * •

Extensions of the positivity algorithm, including collisions, electromagnetic
gyrokinetics, and higher polynomial order: The positivity algorithms detailed
in Chapter 6 are a significant step towards improving robustness of Gkeyll
simulations. Currently, these algorithms are only implemented for the
electrostatic gyrokinetic system. An implementation including collisions has
also been made by Mana Francisquez, but at the time of writing there is some
issue in the implementation with energy conservation. Further extension to the
electromagnetic gyrokinetic system will require additional work, due to the
issues discussed in Appendix 6.A. The algorithm is also formulated in a
general way so that in principle it could be generalized to higher polynomial
order. Finally, the algorithm could also be implemented into the Vlasov-
Maxwell module in Gkeyll.

Along with the further development detailed above, there are a number of
interesting and important physics problems that can leverage the
electromagnetic gyrokinetic capabilities developed in this thesis. An
immediate goal will be to investigate the importance of electromagnetic
effects on SOL dynamics in realistic tokamak geometries at experimental
parameters. Once the additional capability to simultaneously model open- and
closed-field-lines has been developed, we will be able to study the dynamics
of the coupled pedestal/SOL system. This is of critical importance for the
development of a fusion pilot plant, and a major theme of the recent FESAC
Long Range Planning report: “a sustained burning plasma at high power density
is required simultaneously with a solution to the power exhaust challenge:
mitigating the extreme heat fluxes to materials surrounding the plasma”
(Carter et al., 2020).

Electromagnetic effects are expected to play a significant role in the
pedestal region, in part due to large pressure gradients that push the plasma
close to the ideal-MHD stability threshold (Snyder et al., 2011). Further,
while electrostatic turbulence is often suppressed in the pedestal region by
$E\times B$ shear and other effects, the transport can remain above
neoclassical levels due to the presence of electromagnetic instabilities such
as microtearing modes (Hatch et al., 2016). Since turbulence suppression in
the pedestal region plays a key role in pedestal formation and sustenance in
H-mode, understanding the impact of electromagnetic effects on pedestal
transport is of critical importance to the success of current and future
fusion devices such as ITER. Edge-localized modes (ELMs) can also play a key
role in limiting the pedestal pressure gradient, and ELMs are strongly
electromagnetic. Thus self-consistent study of pedestal dynamics requires
modeling the coupled pedestal/SOL system, but previous efforts have relied on
the electrostatic approximation to neglect electromagnetic perturbations
(Idomura et al., 2009; Abiteboul et al., 2013; Churchill et al., 2017; Ku et
al., 2018a). By leveraging and extending the unique capabilities developed in
this thesis to model full-$f$ electromagnetic gyrokinetic turbulence in the
boundary region, we will be able to model the evolution of the coupled
pedestal/SOL system in the presence of electromagnetic microturbulence. This
will enable exciting and impactful research that will be valuable for
understanding current experiments and ensuring the success of future fusion
reactors.

## References

  * Abel et al. (2013) Abel, I. G., Plunk, G. G., Wang, E., Barnes, M., Cowley, S. C., Dorland, W. & Schekochihin, A. A. 2013 Multiscale gyrokinetics for rotating tokamak plasmas: Fluctuations, transport and energy flows. Reports Prog. Phys. 76 (11).
  * Abiteboul et al. (2013) Abiteboul, J., Ghendrih, P., Grandgirard, V., Cartier-Michaud, T., Dif-Pradalier, G., Garbet, X., Latu, G., Passeron, C., Sarazin, Y., Strugarek, A., Thomine, O. & Zarzoso, D. 2013 Turbulent momentum transport in core tokamak plasmas and penetration of scrape-off layer flows. Plasma Phys. Control. Fusion 55 (7), 74001–74012.
  * Angelino et al. (2006) Angelino, P., Bottino, A., Hatzky, R., Jolliet, S., Sauter, O., Tran, T. M. & Villard, L. 2006 On the definition of a kinetic equilibrium in global gyrokinetic simulations. Phys. Plasmas 13 (5), 969.
  * Angus et al. (2012) Angus, J. R., Krasheninnikov, S. I. & Umansky, M. V. 2012 Effects of parallel electron dynamics on plasma blob transport. Phys. Plasmas 19 (8), 82312.
  * Antonsen & Lane (1980) Antonsen, T. M. & Lane, B. 1980 Kinetic equations for low frequency instabilities in inhomogeneous plasmas. Phys. Fluids 23 (6), 1205–1214.
  * Arnold & Awanou (2011) Arnold, D. N. & Awanou, G. 2011 The serendipity family of finite elements. Found. Comput. Math. 11 (3), 337–344.
  * Artun & Tang (1994) Artun, M. & Tang, W. M. 1994 Nonlinear electromagnetic gyrokinetic equations for rotating axisymmetric plasmas. Phys. Plasmas 1 (8), 2682–2692.
  * Bao et al. (2018) Bao, J., Lin, Z. & Lu, Z. X. 2018 A conservative scheme for electromagnetic simulation of magnetized plasmas with kinetic electrons. Phys. Plasmas 25 (2), 22515.
  * Batishcheva et al. (1996) Batishcheva, A. A., Batishchev, O. V., Shoucri, M. M., Krasheninnikov, S. I., Catto, P. J., Shkarofsky, I. P. & Sigmar, D. J. 1996 A kinetic model of transient effects in tokamak edge plasmas. Phys. Plasmas 3 (5), 1634–1639.
  * Beer et al. (1995) Beer, M. A., Cowley, S. C. & Hammett, G. W. 1995 Field–aligned coordinates for nonlinear simulations of tokamak turbulence. Phys. Plasmas 2 (7), 2687–2700.
  * Beer & Hammett (1996) Beer, M. A. & Hammett, G. W. 1996 Toroidal gyrofluid equations for simulations of tokamak turbulence. Phys. Plasmas 3 (11), 4046–4064.
  * Belli & Candy (2010) Belli, E. A. & Candy, J. 2010 Fully electromagnetic gyrokinetic eigenmode analysis of high-beta shaped plasmas. Phys. Plasmas 17 (11), 112314.
  * Bernard et al. (2019) Bernard, T. N., Shi, E. L., Gentle, K. W., Hakim, A., Hammett, G. W., Stoltzfus-Dueck, T. & Taylor, E. I. 2019 Gyrokinetic continuum simulations of plasma turbulence in the Texas Helimak. Phys. Plasmas 26 (4).
  * Bernard et al. (2020) Bernard, T. N., Stoltzfus-Dueck, T., Gentle, K. W., Hakim, A., Hammett, G. W. & Shi, E. L. 2020 Investigating shear flow through continuum gyrokinetic simulations of limiter biasing in the Texas Helimak. Phys. Plasmas 27 (6), 62304\.
  * Boedo et al. (2014) Boedo, J. A., Myra, J. R., Zweben, S., Maingi, R., Maqueda, R. J., Soukhanovskii, V. A., Ahn, J. W., Canik, J., Crocker, N., D’Ippolito, D. A., Bell, R., Kugel, H., Leblanc, B., Roquemore, L. A. & Rudakov, D. L. 2014 Edge transport studies in the edge and scrape-off layer of the National Spherical Torus Experiment with Langmuir probes. Phys. Plasmas 21 (4), 42309.
  * Bohm (1949) Bohm, D. 1949 Minimium ionic kinetic energy for stable sheath. In Charact. Electr. Discharges Magn. Fields (ed. A. Guthrie & K.R. Wakerling), chap. 3. MeGraw-Hill Book Company.
  * Bourdelle et al. (2003) Bourdelle, C., Dorland, W., Garbet, X., Hammett, G. W., Kotschenreuther, M., Rewoldt, G. & Synakowski, E. J. 2003 Stabilizing impact of high gradient of $\beta$ on microturbulence. Phys. Plasmas 10 (7), 2881–2887.
  * Braginskii (1965) Braginskii, S. I. 1965 Transport processes in a plasma; in Reviews of Plasma Physics, M.A. Leontovich (ed.). Rev. Plasma Phys. 1, 205.
  * Brizard (1995) Brizard, A. J. 1995 Nonlinear gyrokinetic Vlasov equation for toroidally rotating axisymmetric tokamaks. Phys. Plasmas 2 (2), 459–471.
  * Brizard & Hahm (2007) Brizard, A. J. & Hahm, T. S. 2007 Foundations of nonlinear gyrokinetic theory. Rev. Mod. Phys. 79 (2), 421–468.