# Magnetic Fluctuations in Gyrokinetic Simulations of Tokamak Scrape-Off Layer
Turbulence

Noah Roth Mandell

###### Abstract

Understanding turbulent transport physics in the tokamak edge and scrape-off
layer (SOL) is critical to developing a successful fusion reactor. The
dynamics in these regions plays a key role in achieving high fusion
performance by determining the edge pedestal that suppresses turbulence in the
high-confinement mode (H-mode). Additionally, the survivability of a reactor
is set by the heat load to the vessel walls, making it important to understand
turbulent spreading of heat as it flows along open magnetic field lines in the
SOL. Large-amplitude fluctuations, magnetic X-point geometry, and plasma
interactions with material walls make simulating turbulence in the edge/SOL
more challenging than in the core region, necessitating specialized
gyrokinetic codes. Further, the inclusion of electromagnetic effects in
gyrokinetic simulations that can handle the unique challenges of the boundary
plasma is critical to the understanding of phenomena such as the pedestal and
edge-localized modes, for which electromagnetic dynamics are expected to be
important.

In this thesis, we develop the first capability to simulate electromagnetic
gyrokinetic turbulence on open magnetic field lines. This is an important step
towards comprehensive electromagnetic gyrokinetic simulations of the coupled
edge/SOL system. By using a continuum full-$f$ approach via an energy-
conserving discontinuous Galerkin (DG) discretization scheme that avoids the
Ampère cancellation problem, we show that electromagnetic fluctuations can be
handled in a robust, stable, and efficient manner in the gyrokinetic module of
the Gkeyll code. We then present results which roughly model the scrape-off
layer of the National Spherical Torus Experiment (NSTX), and show that
electromagnetic effects can affect blob dynamics and transport. We also
formulate the gyrokinetic system in field-aligned coordinates for modeling
realistic edge and scrape-off layer geometries in experiments. A novel DG
algorithm for maintaining positivity of the distribution function while
preserving conservation laws is also presented.

Gregory W. Hammett

###### Acknowledgements.

I have been extremely fortunate to have been supported by an incredible group
of family, friends, mentors, teachers, and collaborators over my academic
career. First I would like to thank my thesis adviser, Greg Hammett. His
physical intuition and insight has been a tremendous resource and inspiration
throughout my time at Princeton. “Physics is about thinking slowly,” Greg told
me on numerous occasions as we puzzled over a problem in his office.111Greg
attributes this quote to Fred Skiff, a very deep thinker. He is masterful at
slowly and carefully thinking through every part of a calculation, model, or
result. This attention to detail, taking the little steps that are needed to
realize big ideas, is among the most important things that I have learned from
Greg. I look forward to many more years of collaboration and friendship. I
would also not be the scientist I am today without the mentorship and
friendship of Bill Dorland. Bill introduced me to the world of fusion,
computational plasma physics, and turbulence when I was a young undergraduate
at the University of Maryland. I was immediately drawn to his passion and his
big, exciting ideas. Thank you for patiently teaching me so much, believing in
me enough to give me such an ambitious project as an undergraduate, and for
the continued mentorship and collaborations while I was at Princeton.
(Remember when he almost stole me for good in year 3, Greg? Haha.) More
importantly, thank you for all the support, guidance, advice, and inspiration
over the years. Thank you to Ammar Hakim for fearlessly leading the
development of the Gkeyll code. Ammar’s computational insights have been
invaluable throughout my thesis work, and there is no question that I have
become a better software developer and physicist because of Ammar’s help.
Working on the Gkeyll code has been incredibly stimulating and rewarding,
mostly thanks to the amazing people in the Gkeyll group. All of this would not
have been possible without all of your hard work. Thank you to Jimmy Juno,
Mana Francisquez, Petr Cagas, Tess Bernard, Rupak Mukherjee, Liang Wang, and
many others for your constant support, often involving late-night Slack chats
and code commits. I also thank Eric Shi, whose thesis paved the path for much
of my work. Thank you to my thesis readers, Matt Kunz and Walter Guttenfelder.
Their careful reading and insightful feedback have greatly improved the text.
Greg Hammett also provided valuable comments. Also, thank you to my thesis
committee members: Matt, Walter, Greg, Ammar, and Sir Steve Cowley. I thank
Stewart Zweben for mentoring me for my first-year project on NSTX gas-puff
imaging, taking a student who had already made up his mind that he wanted to
be a theorist and having the patience to show him the nuances of tokamak
experiments. To all my fellow grad students and post-docs at the lab, thank
you for making my time at Princeton so enjoyable. I especially enjoyed APS
gallivanting with Denis St. Onge and Brian Kraus and others; Denis’ brisket
parties; watching Packers football games with Peter Bolgert, Daniel Ruiz, Jeff
Lestz and others; playing softball with the Tokabats; and of course, playing
ping pong. Some of my fondest memories are from that little room tucked away
behind Science Ed, playing for hours with Peter B. (original commissioner of
the PPPL Ping Pong League, or PPPLPPL), Daniel R., Jonathan Ng, Jeff L.,
Vasily Geyko, Brian K., Lee Ellison, Jack Matteucci, Hongxuan Zhu, Joaquim
Loizu, Vinicius Duarte, David Pfefferle, Jacob Schwartz, Charles Swanson, Ian
Ochs, Alex Glasser, Elijah Kolmes, Andy Alt, Nick McGreivy, and many
others/left-handed-alter-egos. The Bolgert Open tournament was an annual
summer tradition (could never get past Hongxuan in the finals…), and Brian
even developed a computerized ratings system. We probably spent way too much
time playing ping pong when we should’ve been working, but we all got pretty
good and it was a lot of fun. I am also honored to have been a member of the
party office, which has a rich history of housing outstanding Princeton
graduates. Thank you to Dara Lewis and Beth Leman for all their help with
navigating the administrative details of being a graduate student at PPPL.
There are many more teachers and professors to thank for helping me along the
way, including the outstanding professors and lecturers in the Princeton
Program in Plasma Physics, but I would like to say a special thank you to my
high school physics teacher, Pieter Kreunen. My love of physics started back
in Mr. Kreunen’s classroom, getting shocked by Van der Graff generators and
building robots. Thank you for inspiring me and encouraging me to pursue
physics those many years ago. I have been blessed with an amazingly loving and
supportive family. Mom and Dad, I am forever grateful for every opportunity
you gave me growing up, and the sacrifices you made to give them to me.
Madison, you mean more to me than you know, and I am so proud to have you as
my sister. Thank you for supporting me unconditionally every step of the way.
I am also so grateful for the invaluable months I have been able to spend over
the past year222I do not want to say this without solemnly acknowledging the
COVID-19 pandemic that has devastated the country and the world for over a
year. We have been fortunate to keep our health, but millions have not been so
lucky. In these strange times, the Stoppard quote at the beginning, originally
included as a quip about mathematical consistency, has truly taken on a new
meaning… back home in Gurnee with Mom, Dad, and Madison, and in Maine with
Calla, Arthur, Mariah, Stephen and Jack. Lastly, most importantly: Haley, we
made it! You’ve been with me every step of this journey, and I wouldn’t have
gotten here without you. This is yours as much as it is mine. I am inspired by
you every day. Thank you for always lovingly supporting me and believing in
me. I wouldn’t trade a single day with you for the world.
**************************************************************************

I was very fortunate to be funded for four years by the Department of Energy
Computational Science Graduate Fellowship, provided under DOE grant DE-
FG02-97ER25308. Thank you for the generous support, and for all the friends I
made as part of the fellowship program. Additional funding came from DOE
contract DE-AC02-09CH11466. Simulations were performed on the Perseus and Eddy
clusters at Princeton University/PPPL, and the Cori system at the National
Energy Research Scientific Computing Center. to Haley, with love

##  Chapter 1 Introduction

###  1.1 Motivation: the promise of fusion energy

After Einstein first discovered the relationship between energy and mass,
governed by the iconic equation $E=mc^{2}$, it was soon realized that this
relationship was the key to the process that produces the energy of the Sun
and stars: nuclear fusion. Fifteen years after Einstein’s discovery, British
astrophysicist Arthur Eddington was the first to describe how the Sun and
similarly-sized stars create their energy by fusing hydrogen atoms into
helium. Eddington realized that the tiny difference in mass between a helium
atom and its constituent hydrogen parts, as had been recently shown by Aston,
meant that the ‘missing’ mass is converted into energy via Einstein’s
equation. “The store is well nigh inexhaustible, if only it could be tapped,”
Eddington said in a lecture at the annual meeting of the British Association
for the Advancement of Science in Cardiff (Eddington, 1920). Thus began the
promise of man-made fusion as a terrestrial energy source.

One of the main allures of fusion power is the abundance of the fuel. Unlike
fossils fuels, which at current energy-consumption rates would be burned
through in less than 1,000 years (causing catastrophic global warming in the
process), the fuel for fusion is virtually limitless because it can be
extracted from seawater. The most promising fusion reaction for use on Earth
is not the proton-proton reaction that powers the Sun, but an easier-to-
initiate reaction between deuterium (2H) and tritium (3H):

${}^{2}\mathrm{H}+{}^{3}\mathrm{H}\rightarrow{}^{4}\mathrm{He}\ (3.5\
\mathrm{MeV})+\mathrm{n}\ (14.1\ \mathrm{MeV}).$ (1.1)

Deuterium is a naturally abundant isotope of hydrogen that can be readily
extracted from seawater at minimum cost, with each liter of seawater
containing $\sim 0.02$ g of deuterium. Tritium is not naturally abundant due
to its relatively short half-life of 12.3 years. However, fusion reactors can
use the energetic neutron from the D-T reaction to breed their own tritium via
lithium (6Li) blankets via the reaction

$\mathrm{n}+{}^{6}\mathrm{Li}\rightarrow{}^{4}\mathrm{He}\ (2.1\
\mathrm{MeV})+{}^{3}\mathrm{H}\ (2.7\ \mathrm{MeV}).$ (1.2)

Current world lithium supplies are approximately 13.5 million tons, but
lithium is also contained in seawater at a concentration of 0.2 mg per liter.
Thus there is enough fusion fuel readily available in the oceans to power the
Earth for millions of years, several orders of magnitude longer than other
terrestrial fuel sources other than solar energy (Cowley, 2016).

Other major benefits of fusion are its minimal environmental impact and
operational safety. Fusion would be clean and virtually carbon-neutral,
emitting no greenhouse gases and not contributing to climate change. While
this benefit is also shared by nuclear fission (where energy is produced by
splitting the nuclei of heavy elements like uranium), fusion has the
additional advantage that it has no long-lived radioactive byproducts. Helium
is an inert gas, and while the energetic neutron from the D-T reaction can
transmute the materials in the walls of a reactor and make them radioactive
over time, the use of low-activation wall materials would make the waste
substantially safer than fission waste. Further, a fusion power plant would be
safer to operate than a fission power plant, as there is no runaway meltdown
scenario. Unlike fission reactions, fusion reactions immediately shut down
when the fuel is removed or cooled.

Unfortunately, a fusion reaction is very difficult to get started. To produce
fusion, the positively-charged fuel nuclei must have enough energy to overcome
the repulsive Coulomb force between them; only then can they get close enough
to fuse together via the nuclear strong force and release energy. Unlike in
the Sun, where immense gravitational pressure creates conditions necessary for
fusion, terrestrial fusion must achieve fusion conditions via other methods.
The most promising approach is to heat a gas of deuterium and tritium to very
high temperature, so that particles have enough energy that random collisions
can overcome the Coulomb repulsion. The energies required, of order $10$ keV,
are well above the electron binding energy, resulting in the fuel gases
ionizing fully and becoming a plasma.

At these extreme temperatures, the fuel cannot simply be contained by material
walls. Instead, we can take advantage of the fact that a plasma is composed of
charged particles. In the presence of magnetic fields, charged particles
spiral helically around the field lines, providing a way to control the
particle motion. Particles are still free to move parallel to the magnetic
field lines, so one way to confine them is to wrap the field lines into a
torus shape, creating a ‘magnetic bottle’. However, a simple ring
configuration leads to vertical particle drifts, making the configuration
intrinsically unstable. This problem can be overcome by twisting the field
lines into a helical shape wrapping around the torus. The twisting magnetic
field guides the particles up or down, counteracting the drift motion and
enhancing confinement. This configuration is the basis of both the tokamak and
stellarator concepts. In a tokamak, the twist in the field is produced by a
current driven toroidally through the plasma, whereas in a stellarator the
twist is produced by shaped helical field coils. These two configurations are
shown schematically in Fig. 1.1. We will focus on tokamaks in this work.

Figure 1.1: Schematic diagrams for a tokamak (left) and a stellarator (right).
In a tokamak, the poloidal component of the magnetic field that gives the
helical twist is produced by current driven toroidally through the plasma. In
a stellarator, specially-shaped helical field coils produce the poloidal
component of the magnetic field. (Source: CEA)

###  1.2 Turbulent transport in fusion plasmas

With the plasma confined (to lowest order), the plasma can be heated without
direct contact with the vessel walls. The goal is then to keep the plasma hot
enough and dense enough for long enough for fusion to occur. This is the idea
behind the fusion triple product, $nT\tau_{E}$, where $n$ is the plasma
density, $T$ is the mean temperature, and $\tau_{E}$ is the energy confinement
time. Lawson’s criterion gives the condition for ‘break-even’, at which point
the plasma’s self-heating from fusion exceeds its losses (Wesson, 2005),

$nT\tau_{E}\geq 10^{21}\ \mathrm{keV\boldsymbol{\cdot}s/m^{3}}.$ (1.3)

In practice, the energy confinement time has proved to be the most challenging
component to maximize. It is defined as

$\tau_{E}=\frac{W}{P_{\mathrm{loss}}},$ (1.4)

where $W$ is the energy content of the plasma and $P_{\mathrm{loss}}$ is the
energy loss rate. While the particles are well-confined along the magnetic
field lines so that the parallel (with respect to the field lines) confinement
time is very large, particles can also diffuse radially outward, perpendicular
to the magnetic field. As a result, the perpendicular diffusion rate is the
limiting factor on the energy confinement time. The transport was originally
thought to be dominated by collisional processes (yielding “classical” and
geometry-modified “neoclassical” transport), but these processes were found to
greatly under-predict the transport seen in tokamaks, with most of the
measured transport denoted “anomalous”. It is now recognized that plasma
turbulence is responsible for this anomalous component, so that tokamak plasma
confinement is dominated by turbulent transport.

In the tokamak core, turbulence is driven by small-scale, low-frequency
“micro-instabilities”. These instabilities feed off the density and
temperature gradients that inherently result from the requirement that the
temperature must be low ($\sim 10^{3}$ K) near the walls of the device but
very hot in the core ($\sim 10^{8}$ K). Despite fluctuation levels of only
order $1\%$, core turbulence leads to significant transport of particles,
momentum, and heat. The fluctuations typically have length scales
perpendicular to the background magnetic field on the order of the ion
gyroradius $\rho_{i}=v_{ti}/\Omega_{i}$ and frequencies (and growth rates) on
the order of the diamagnetic drift frequency,
$\omega_{\ast}=k_{\theta}\rho_{i}v_{ti}/L_{n}$, where $k_{\theta}$ is a
typical poloidal wavenumber, $v_{ti}=\sqrt{T_{i}/m_{i}}$ is the ion thermal
speed, $\Omega_{i}=ZeB/m_{i}$ is the ion cyclotron frequency, and
$L_{n}=-(\textnormal{d}\ln n/\textnormal{d}r)^{-1}$ is the density scale
length. Given these length and time scales, we can make a simple mixing-length
estimate of the diffusivity,

$D\sim\frac{(\Delta x)^{2}}{\Delta
t}\sim\rho_{i}^{2}\omega_{\ast}\sim\rho_{i}^{2}k_{\theta}\rho_{i}\frac{v_{ti}}{L_{n}}.$
(1.5)

Taking $k_{\theta}\rho_{i}\sim 1$ yields the so-called gyro-Bohm diffusivity,
$D_{\mathrm{gB}}\sim\rho_{i}^{2}v_{ti}/L_{n}$. While this gives a rough
scaling of the transport, additional theory and numerical simulation are
required for meaningful understanding and quantitative prediction of turbulent
transport in tokamaks.

###  1.3 The boundary plasma

Figure 1.2: Diagram of plasma profiles in L- and H-mode regimes. In H-mode
profiles are elevated by the pedestal, which is the result of a transport
barrier that forms near the last closed flux surface. Figure from (Lang et
al., 2013).

While the core attracted much of the focus in the early days of fusion
research, it was soon realized that the edge and scrape-off layer (SOL), which
we together refer to as the boundary plasma, greatly affect the device
performance and dynamics. Performance is strongly determined by the edge
profiles because core profiles of density and temperature are relatively stiff
(Doyle et al., 2007; Kinsey et al., 2011). A primary example of this is the
high-confinement mode (H-mode), first discovered by Wagner et al. (1982),
where a steep-gradient transport barrier region called the pedestal forms in
the edge and raises the core profiles (as if they were standing on a
pedestal), as shown in Fig. 1.2. Strong sheared poloidal flows are observed in
this region, correlated with a reduction in turbulent fluctuation levels and
fluxes. Understanding pedestal formation and predicting the pedestal height
are of great current interest (Snyder et al., 2011), and a major motivator for
first-principles modeling of the boundary plasma.

The scrape-off layer (SOL) is the region outside the last closed flux surface
(LCFS) where the field lines are open and terminate on material walls. Charged
particles move freely along the field lines and are lost when they strike the
walls (until they recombine and reenter the plasma as cold neutrals, a process
called recycling). The dynamics in the SOL is primarily set by the interplay
between particles and heat crossing the LCFS from the edge, parallel losses to
the walls, cross-field turbulent transport, and plasma surface interactions
(PSIs), including recycling and impurity fluxes. As a result of these
processes, the SOL plasma is rather cold, with $T_{e}\sim 10-100$ eV.

Figure 1.3: Diagram of limiter $(a)$ and divertor $(b)$ configurations of
tokamak operation. Figure from Wesson (2005).

The termination points of the open field lines in the SOL are determined by
whether the tokamak is operated in a limiter or divertor configuration, as
shown in the diagram in Fig. 1.3. In the former, material limiters are placed
at various locations on the first wall. The field lines that intersect the
limiters then define the SOL. While the limiter configuration is operational,
the divertor configuration is generally preferred in high-performance devices.
In the divertor configuration, an external current in the direction of the
plasma current is applied at the top and/or bottom of the device, resulting in
the formation of X-point nulls. This moves the plasma-wall interactions onto
the divertor targets, which are much further away from the main core plasma
than limiter plates. This is beneficial since neutrals and impurities released
from the divertor plates cannot directly enter the core plasma. Divertor
configurations are also preferable for handling the heat exhaust requirements
of the SOL and removing impurities and fusion ash via pumping (Wesson, 2005).

####  1.3.1 Intermittent SOL transport and blob dynamics

The cross-field transport in the SOL is highly intermittent. Unlike in the
core, where the transport is dominated by small fluctuations, fluctuations in
the SOL can be comparable to the equilibrium quantities. This is primarily due
to the convective transport of coherent structures of enhanced density and
temperature called blobs or filaments. These structures propagate quasi-
ballistically, moving radially outwards and resulting in significant particle
and heat transport. Blobs are highly extended along the field line with
parallel lengths $\sim 1-10$ m and much smaller scales $\sim 1-10$ cm
perpendicular to the field (Zweben et al., 2017). The intermittent nature of
blob transport suggests that a simple picture of diffusive transport is
inadequate (Naulin, 2007). Instead, the transport is avalanche-like,
suggesting that the system gets pushed up against some critical gradient
threshold and then intermittently releases bursts of transport when the
threshold is exceeded (LaBombard et al., 2005; Labombard et al., 2008). An
expansive review of experimental evidence and theoretical understanding of
intermittent edge turbulence and blobs is given by D’Ippolito et al. (2011).

The basic mechanism of blob transport is plasma polarization due to magnetic
drifts. On the outboard side of the tokamak, the curvature and $\nabla B$
drifts are vertical, with ions drifting in one direction and electrons
drifting in the other. The resulting charge polarization produces a vertical
electric field across the blob, giving a radially outward $E\times B$ drift
(Krasheninnikov, 2001). This is shown schematically in Fig. 1.4.

Figure 1.4: Diagram of basic blob transport mechanism from charge
polarization. Figure from (Krasheninnikov et al., 2008). Figure 1.5: Diagram
of equivalent blob circuit. Charge polarization in the blob provides a
constant current source. The current flows along the field lines in a dipolar
structure. The current can be closed through the sheath or through the plasma,
depending on the resistivities of the respective current paths. Figure from
(Krasheninnikov et al., 2008). Note that the circuit element through which the
polarization current flows may be more appropriately characterized as a
capacitor (Xu et al., 2010).

The magnitude of the blob electric field, and thereby the blob speed, is
affected by the balance between the polarization current and parallel
currents. To explain this, it is useful to visualize the currents in the blob
via a blob equivalent circuit (Myra & D’Ippolito, 2005; Krasheninnikov et al.,
2008; Xu et al., 2010), as shown in the circuit diagram in Fig. 1.5 from
Krasheninnikov et al. (2008). (Note that the circuit element through which the
polarization current flows may be more appropriately characterized as a
capacitor due to plasma inertia (Xu et al., 2010).) The magnetic drifts act as
a local current source. At constant current, the potential drop across the
blob is determined by the resistance in the circuit. If the plasma has low
resistivity ($\eta_{\parallel}$), the current flows freely along the field
lines to the sheath, and the effective sheath resistance
($\eta_{\mathrm{sheath}}$) will determine the blob potential and thereby the
blob velocity. This is known as the sheath-limited regime, and it can lead to
reduced blob speed and transport as the blob polarization current can be
effectively shorted out by the current closure through the sheath. Conversely,
if the plasma resistivity $\eta_{\parallel}$ is larger due to increased
collisionality, the effective resistivity in the circuit will increase and
lead to larger blob velocity. At large enough resistivity, parallel currents
are hindered enough that cross-field current closure happens away from the
sheath via ion polarization currents or collisional currents, with complete
disconnection giving the inertial or resistive-ballooning regime. Magnetic
shear (especially near the X-point) can have the opposite effect on the blob
velocity, as it can lead to a thin, elongated region of the blob where
magnetic shear is strong. This makes it easier for cross-field currents to
close the circuit through the thin sheared part of the blob, reducing the
resistivity of the current loop and thereby slowing the blob. Current closure
through regions of high magnetic shear can thus also effectively disconnect
the blob from the sheath. However, notice that sheath disconnection due to
increased collisionality gives the opposite effect on blob velocity than
sheath disconnection via magnetic shear; the former results in increased
effective blob circuit resistivity and larger velocities, while the latter
decreases resistivities and slows the blobs (D’Ippolito et al., 2011;
Krasheninnikov et al., 2008).

####  1.3.2 SOL heat exhaust problem

Particles and heat from the core are transported across the LCFS and exhausted
in the SOL. The heat flows quickly along the open field lines to the walls,
with the parallel heat flux in the SOL reaching above 500 MW/m2 in some
present devices and expected to be $\sim 1$ GW/m2 in ITER (Loarte et al.,
2007). The maximum heat load for present materials with active cooling is
typically $10$ MW/m2 normal to the surface in steady state and $20$ MW/m2 for
transients (Loarte et al., 2007). Thus the heat load must be reduced below
these material limits in order to avoid damage to the wall plates and the
introduction of impurities that degrade fusion performance. The heat load can
be reduced in part by making the incidence angle of the magnetic field lines
on the walls very shallow $(\sim 2-5^{\circ})$ to reduce the component of the
flux normal to the walls, but this still leaves a significant portion of the
heat to be dissipated via other means. The width of the heat flux channel
becomes an important parameter, since spreading the heat over a larger area
reduces the peak heat load. Here, cross-field turbulent transport is
beneficial as it can widen the heat flux width. An empirical scaling of the
heat flux width, $\lambda_{q}$, computed from a multi-machine database has
shown that the heat-flux width, mapped to the outboard midplane, varies
strongest with the inverse of the plasma current (or equivalently, the inverse
of the poloidal magnetic field strength) (Eich et al., 2013). Simply
extrapolating the scaling from present-day experiments to the upcoming ITER
experiment suggests that the heat flux width for the ITER $Q=10$ baseline
could be $\approx 1$ mm (Eich et al., 2013), much smaller than the $3-3.5$ mm
result from the ITER physics basis based mostly on JET ELM-averaged data
(Loarte et al., 2007). SOLPS transport modeling has suggested
$\lambda_{q}=3.6$ mm (Kukushkin et al., 2013). The validity of these empirical
scalings for the ITER heat flux width is an important issue that must be
addressed by first-principles modeling. A recent XGC1 electrostatic
gyrokinetic simulation predicted $\lambda_{q}\approx 5.9$ mm (Chang et al.,
2017), with the width widened due to electron turbulence. Additional analysis
of XGC1 data has suggested that trapped electron mode (TEM) turbulence in
particular is responsible for increased SOL heat transport. While $E\times B$
shear suppresses TEM in the SOL of present devices, $E\times B$ shear is
predicted to be weaker in ITER, allowing TEM to drive transport. These results
have suggested a new scaling of $\lambda_{q}\sim
1/B_{\mathrm{pol}}(a/\rho_{i,\mathrm{pol}})$, with the new parameter
$a/\rho_{i,\mathrm{pol}}$ related to the neoclassical $E\times B$ shearing
rate (Chang et al., 2020).

###  1.4 Electromagnetic effects in the boundary plasma

In this thesis we will focus in particular on electromagnetic effects in the
plasma boundary. The edge/SOL region features steep pressure gradients,
especially in the H-mode transport barrier and SOL regions, which contribute
to the importance of electromagnetic effects. Experimental evidence has
indicated that the edge plasma state is controlled by electromagnetic drift
wave dynamics (LaBombard et al., 2005; Labombard et al., 2008). In this
regime, the parallel electron dynamics is no longer fast relative to the drift
turbulence, so electrons can no longer be treated adiabatically (Scott, 1997).
This leads to coupling of the perpendicular vortex motions and kinetic shear
Alfvén waves, which results in field-line bending (Xu et al., 2010). The
slowing of parallel electron dynamics can also add impedance along the field
line, leading to blobs becoming electrically disconnected from the sheath and
resulting in enhanced blob velocities. While in the electrostatic case the
sheath potential is communicated to the upstream plasma rapidly on the order
of the electron transit time ($\tau_{e}=L_{\parallel}/v_{te}$, with
$v_{te}=\sqrt{T_{e}/m_{e}}$), in the electromagnetic case Alfvén waves
communicate the potential on the order of the Alfvén time,
$\tau_{A}=L_{\parallel}/v_{A}=L_{\parallel}/(\sqrt{2/\beta}c_{\mathrm{s}})=L_{\parallel}\sqrt{\mu_{0}nm_{i}}/B$,
with $c_{\mathrm{s}}=\sqrt{T_{e}/m_{i}}$ the sound speed. Thus a basic
condition for electromagnetic effects to alter sheath connection is
$v_{A}<v_{te}$, or $\hat{\beta}\equiv(\beta/2)m_{i}/m_{e}>1$. If in the time
$\tau_{A}$ the blob is able to move more than its width across the field, the
information about the sheath will never reach it. Thus the blob will move as
if the sheath did not exist if $\tau_{A}\gtrsim L_{\perp}/v_{\perp}$, or
${\beta}\gtrsim(L_{\perp}/L_{\parallel})^{2}(c_{\mathrm{s}}/v_{\perp})^{2}$,
where $L_{\perp}$ is the typical length scale of the potential of the blob,
and $v_{\perp}$ is the blob radial velocity at the midplane (Lee et al.,
2015a; Hoare et al., 2019). Given these conditions, electromagnetic effects
could especially be important for the high beta filaments found in edge
localized modes (ELMs), which are large-scale magnetohydrodynamic (MHD) modes
that result in large, high pressure filaments originating from the pedestal.
ELM filaments also carry a large unidirectional current, distinguishing them
from standard blobs and further enhancing the electromagnetic effects of ELMs
by inducing magnetic field perturbations (Myra, 2007; Kirk et al., 2005, 2006;
Migliucci & Naulin, 2010; Vianello et al., 2011). Additionally, experiments
have found correlations between (non-ELM) large blobs and MHD modes (Zweben et
al., 2020).

The following subsections briefly illustrate the role of magnetic induction in
determining the parallel electron dynamics and producing field-line bending,
mostly following Xu et al. (2010) and Scott (1997).

####  1.4.1 Parallel electron dynamics and the role of magnetic induction

The strong mobility of electrons along the field line makes it important to
understand how the parallel current responds to forces due to parallel
gradients in the density $n$, electron temperature $T_{e}$, and the
electrostatic potential $\Phi$. The linear response determines the propagation
of wave-like disturbances along the field line. The dynamics is governed by
the electron parallel force balance equation, also known as the parallel
component of the generalized Ohm’s law. In the electrostatic limit, we have

$\frac{m_{e}}{ne^{2}}\frac{\textnormal{d}J_{\parallel}}{\textnormal{d}t}+\eta_{\parallel}J_{\parallel}=\frac{1}{ne}{\nabla}_{\parallel}p_{e}-{\nabla}_{\parallel}\Phi.$
(1.6)

On the right-hand side, we have the balance between the parallel pressure and
electric forces, where $p_{e}=nT_{e}$ is the electron pressure,
$n=n_{i}=n_{e}$ is the plasma density (assuming a quasi-neutral plasma with
singly charged ions), and $E_{\parallel}=-\nabla_{\parallel}\Phi$ is the
parallel electric field in the electrostatic limit. Here,
$\nabla_{\parallel}=\mathbf{\hat{b}}\boldsymbol{\cdot}\nabla$ denotes a
derivative in the direction of the background magnetic field,
$\mathbf{\hat{b}}=\mbox{\boldmath${B}$}_{0}/B$. On the left-hand side, the
first term is electron inertia, which gives finite-electron-mass ($m_{e}$),
collisionless effects. Here, $\textnormal{d}/\textnormal{d}t=\partial/\partial
t+\mbox{\boldmath${v}$}_{E}\boldsymbol{\cdot}\nabla$ is the total time
derivative, with
$\mbox{\boldmath${v}$}_{E}=(1/B)\mathbf{\hat{b}}\times\nabla\Phi$ the $E\times
B$ velocity. The second term on the left-hand side is resistive friction, with
$J_{\parallel}\approx-enu_{\parallel e}$ the parallel current (dominated by
electron parallel flow $u_{\parallel e}$) and
$\eta_{\parallel}=0.51m_{e}\nu_{ei}/(ne^{2})$ the parallel resistivity, which
is proportional to the electron collision frequency. The electrons are said to
be “adiabatic” when the forces on the right-hand side balance. After
linearizing and assuming the electrons are sufficiently fast to isothermalize
along the field line so that $\nabla_{\parallel}T_{e}=0$, we have

$T_{e0}\nabla_{\parallel}n=n_{0}e\nabla_{\parallel}\Phi,$ (1.7)

which results in the adiabatic electron density response, given by the
Boltzmann distribution $n=n_{0}e\Phi/T_{e0}$, with subscript $0$ denoting
background quantities.

Now we will introduce electromagnetic (finite $\beta$) effects. We will
consider only perpendicular magnetic fluctuations of the form
$\mbox{\boldmath${B}$}_{1}=\nabla\times(A_{\parallel}\mathbf{\hat{b}})\approx-\mathbf{\hat{b}}\times\nabla
A_{\parallel}$, where $A_{\parallel}$ is the parallel component of the
magnetic vector potential. This is related to the parallel current via the
parallel component of Ampère’s law,

$-\nabla_{\perp}^{2}A_{\parallel}=\mu_{0}J_{\parallel}.$ (1.8)

Electromagnetic effects enter into parallel force balance in two ways. First,
the parallel gradient must be taken along perturbed field lines, resulting in
an additional “magnetic flutter” component due to $\mbox{\boldmath${B}$}_{1}$:

$\tilde{\nabla}_{\parallel}\equiv\frac{1}{B}(\mbox{\boldmath${B}$}_{0}+\mbox{\boldmath${B}$}_{1})\boldsymbol{\cdot}\nabla=\nabla_{\parallel}-\mathbf{\hat{b}}\times\nabla
A_{\parallel}\boldsymbol{\cdot}\nabla.$ (1.9)

Second, magnetic induction adds to the parallel electric field, which is now
given by

$E_{\parallel}=-\tilde{\nabla}_{\parallel}\Phi-\frac{\partial
A_{\parallel}}{\partial t}.$ (1.10)

As a result, the parallel force balance equation becomes

$\frac{\partial A_{\parallel}}{\partial
t}+\frac{m_{e}}{ne^{2}}\frac{\textnormal{d}J_{\parallel}}{\textnormal{d}t}+\eta_{\parallel}J_{\parallel}=\frac{1}{ne}\tilde{\nabla}_{\parallel}p_{\parallel}-\tilde{\nabla}_{\parallel}\Phi.$
(1.11)

Balancing the first two terms on the left-hand side, we can see that induction
is dominant over inertia at perpendicular scales larger than the collisionless
skin depth $d_{e}=\sqrt{m_{e}/(n_{0}e^{2}\mu_{0})}$, so that
$k_{\perp}d_{e}<1$. Balancing the first and third terms on the left-hand side
gives that induction is dominant over resistivity at perpendicular scales
larger than the collisional skin depth, so that
$k_{\perp}d_{e}\sqrt{\nu_{ei}/\omega}<1$, with $\omega$ some characteristic
frequency so that $\partial/\partial{t}\sim\omega$. Any imbalance of the
forces on the right-hand side will result in non-adiabatic electrons,
providing a channel to exchange the internal particle energy with the magnetic
energy of field-line bending (via induction), or producing irreversible
dissipation of magnetic energy (via resistivity). Defining the parallel
electromotive force (emf) as (Hinton et al., 2003; Xu et al., 2010)

$\psi\equiv\int(\eta_{\parallel}J_{\parallel}-E_{\parallel})\textnormal{d}\ell,$
(1.12)

with $\ell$ the length along the perturbed field line, we can rewrite parallel
force balance compactly as

$\frac{\partial A_{\parallel}}{\partial
t}+\eta_{\parallel}J_{\parallel}=\tilde{\nabla}_{\parallel}(\psi-\Phi).$
(1.13)

####  1.4.2 Field-line bending

To compute the evolution (bending) of magnetic field lines, we use Faraday’s
law,

$\frac{\partial\mbox{\boldmath${B}$}_{1}}{\partial
t}=-\nabla\times\mbox{\boldmath${E}$}.$ (1.14)

The electric field is given by

$\mbox{\boldmath${E}$}=-\nabla\Phi-\mathbf{\hat{b}}\frac{\partial
A_{\parallel}}{\partial
t}=-\nabla\phi-\tilde{\nabla}_{\parallel}(\psi-\Phi)\mathbf{\hat{b}}+\eta\mbox{\boldmath${J}$}=-\nabla_{\perp}\Phi-\tilde{\nabla}_{\parallel}\psi\mathbf{\hat{b}}+\eta\mbox{\boldmath${J}$},$
(1.15)

where we have used parallel force balance and dropped the parallel subscripts
on the resistivity term
($\eta_{\parallel}J_{\parallel}\gg\eta_{\perp}J_{\perp}$). Note that
$\nabla_{\perp}\Phi=\mbox{\boldmath${v}$}_{E}\times(\mbox{\boldmath${B}$}_{0}+\mbox{\boldmath${B}$}_{1})$.
Substituting into Faraday’s law, we have

$\frac{\partial\mbox{\boldmath${B}$}_{1}}{\partial
t}=-\nabla\times\left[\mbox{\boldmath${v}$}_{E}\times(\mbox{\boldmath${B}$}_{0}+\mbox{\boldmath${B}$}_{1})\right]+\nabla\times(\tilde{\nabla}_{\parallel}\psi\mathbf{\hat{b}})+D_{m}\nabla^{2}\mbox{\boldmath${B}$}_{1}-\nabla\eta\times\mbox{\boldmath${J}$}.$
(1.16)

From left to right, on the right-hand side we have a frozen-in term, a drift
term, a magnetic diffusion term (with $D_{m}=\eta/\mu_{0}$) and a resistivity
gradient term. In the limit of small resistivity, we can write

$\frac{\partial\mbox{\boldmath${B}$}_{1}}{\partial
t}\approx\nabla\times\left[\tilde{\nabla}_{\parallel}(\psi-\Phi)\mathbf{\hat{b}}\right].$
(1.17)

This shows that the net parallel gradient force (the right-hand side of Eq.
1.13) drives line bending via non-adiabatic electrons exchanging energy with
the magnetic field (Xu et al., 2010).

Note that we can also write the electric field as

$\mbox{\boldmath${E}$}=(\mbox{\boldmath${B}$}_{0}+\mbox{\boldmath${B}$}_{1})\times\mbox{\boldmath${v}$}_{F}-\nabla\psi+\eta\mbox{\boldmath${J}$},$
(1.18)

where
$\mbox{\boldmath${v}$}_{F}\equiv(1/B)\mathbf{\hat{b}}\times\nabla(\Phi-\psi)$
is the velocity of field lines (neglecting resistive magnetic diffusion). From
this, Faraday’s law becomes

$\frac{\partial\mbox{\boldmath${B}$}_{1}}{\partial
t}=-\nabla\times\left[\mbox{\boldmath${v}$}_{F}\times(\mbox{\boldmath${B}$}_{0}+\mbox{\boldmath${B}$}_{1})\right]+D_{m}\nabla^{2}\mbox{\boldmath${B}$}_{1}-\nabla\eta\times\mbox{\boldmath${J}$}.$
(1.19)

From this it follows that magnetic flux is conserved in the limit of no
resistivity or diffusion, with field lines advected with velocity
$\mbox{\boldmath${v}$}_{F}$. The difference between the $E\times B$ velocity
and the velocity of the field lines is a function of the parallel emf:
$\Delta\mbox{\boldmath${v}$}=\mbox{\boldmath${v}$}_{E}-\mbox{\boldmath${v}$}_{F}=(1/B)\mathbf{\hat{b}}\times\nabla\psi$
(Xu et al., 2010).

###  1.5 Modeling the boundary plasma

The boundary of a tokamak is a complicated nonlinear system. As such,
numerical modeling is a critical tool for helping to understand the physics of
the boundary plasma. As detailed below, several approaches have produced
valuable results and insights at varying levels of complexity and
computational expense. Overviews of some of the numerical modeling approaches
and associated simulation codes for the boundary plasma are given by Ricci
(2015); Loarte et al. (2007); Shi (2017), and briefly detailed below.

####  1.5.1 Empirical modeling

Empirical extrapolation of data obtained from present devices serves as the
basis for much of the modeling and design of tokamak divertors and wall
systems. Codes solve a simplified set of transport equations based on the
Braginskii fluid equations in two dimensions, assuming axisymmetry. Since
plasma turbulence is not captured directly in these models, _ad hoc_ cross-
field anomalous diffusion coefficients are used, with the parameters adjusted
to fit existing experimental data. This system is often coupled to a Monte-
Carlo neutral particle model so that pumping, fueling, and plasma-wall
interactions can be modeled. Codes using this approach include SOLPS (formerly
B2-EIRENE) (Reiter et al., 1991; Schneider et al., 1992), UEDGE (Rognlien et
al., 1994), EDGE2D (Simonini et al., 1994), and SOLDOR (Shimizu et al., 2003).
SOLPS has been used extensively as the SOL simulation code for the ITER
divertor design (Pitts et al., 2009; Kukushkin et al., 2011).

####  1.5.2 Fluid modeling

Given the relatively low temperatures and high collisionalities of the scrape-
off layer, a fluid approach is reasonable to reduce the computational cost of
global turbulence simulations. As such, models based on the drift-reduced
Braginskii equations (Braginskii, 1965; Zeiler et al., 1997) have provided
valuable results and insights on boundary plasma phenomenon. Since these
models only evolve the first three moments of the distribution function, they
rely on high collisionality to provide fluid closure. This implicitly assumes
that the distribution function is close to thermal equilibrium. Codes
employing the drift-reduced Braginskii approach include BOUT++ (Xu et al.,
2008), GBS (Ricci et al., 2012; Halpern et al., 2016), TOKAM3X (Tamain et al.,
2010), GDB (Zhu et al., 2018), and GRILLIX (Stegmeir et al., 2018).

There have also been efforts to extend the validity of the moment-based
approach to more kinetic regimes by using gyrofluid models (Ribeiro & Scott,
2008; Madsen, 2013; Held et al., 2016), based on earlier work on gyrofluid
models for core turbulence (Hammett & Perkins, 1990; Dorland & Hammett, 1993;
Beer & Hammett, 1996; Snyder & Hammett, 2001). Another recent approach uses an
Hermite-Laguerre formulation to allow the use of an arbitrary number of
moments, although this approach has not yet produced numerical results (Jorge
et al., 2017; Frei et al., 2020).

####  1.5.3 Gyrokinetic modeling

Despite the high collisionality of the plasma boundary, kinetic treatments
will inevitably be necessary for reliable quantitative predictions in some
cases (Jenko & Dorland, 2001; Cohen & Xu, 2008). Significant deviations from
thermal equilibrium can occur due to transient events such as ELMs
(Batishcheva et al., 1996). Kinetic treatments are also required if one wishes
to cross the LCFS and model the coupled dynamics of the pedestal and the SOL
within a single framework, since the fluid approximations break down in the
hot pedestal.

While the most general approach would involve solving the full six-dimensional
Vlasov-Maxwell or Fokker-Planck-Maxwell system to model the plasma, this is
impractical due to the high dimensionality and wide range of timescales
involved, including the fast cyclotron motion. Instead, we can take advantage
of the fact that the characteristic turbulent modes have frequencies much
lower than the cyclotron frequency, allowing us to average over the cyclotron
motion and eliminate one of the velocity dimensions (the gyrophase angle). The
result is the gyrokinetic model, which describes the evolution of particle
guiding centers in a reduced five-dimensional phase space. Gyrokinetic theory
and direct numerical simulation have become important tools for studying
turbulence and transport in fusion plasmas, especially in the core region
(Dimits et al., 2000). This includes the simulation codes GEM (Parker et al.,
1993a), GS2 (Kotschenreuther et al., 1995; Dorland et al., 2000), GTC (Lin et
al., 2000), GENE (Jenko, 2000), EUTERPE (Jost et al., 2001), GYRO (Candy &
Waltz, 2003), GT3D (Idomura et al., 2003), GKV (Watanabe & Sugama, 2005), GTS
(Wang et al., 2006), ORB5 (Jolliet et al., 2007; Lanti et al., 2019), GT5D
(Idomura et al., 2008), GKW (Peeters et al., 2009), CGYRO (Candy et al.,
2016), and GX (Mandell et al., 2018). In the edge and SOL, gyrokinetic
simulations are particularly challenging because the large, intermittent
fluctuations in the SOL make assumptions of scale separation between
equilibrium and fluctuations not strongly valid. This necessitates a full-$f$
approach that self-consistently evolves the full distribution function, $f$
(as opposed to the $\delta f$ approach commonly used in the core, where one
assumes $f=F_{0}+\delta f$ with a fixed background $F_{0}$ so that only
$\delta f$ perturbations must be evolved, and the parallel electric field
nonlinearity is frequently neglected). Additional complications of the
edge/SOL region include: open field line regions requiring sheath boundary
conditions and models of plasma-wall interactions; X-point geometry in
diverted configurations, which makes the use of efficient field-aligned
coordinate systems challenging; a wide range of collisionality regimes, from
the hot pedestal top to the cold SOL; and atomic physics and neutral
interactions. Major extensions to existing core gyrokinetic codes or
altogether new efforts are required to meet these challenges. To this end,
steady progress in gyrokinetic boundary plasma modeling has been made with
both particle-in-cell (PIC) and continuum methods. Codes employing the PIC
method in the plasma boundary include XGC1 (Ku et al., 2009, 2016) and ELMFIRE
(Korpilo et al., 2016). Continuum methods are used by the codes COGENT (Dorf
et al., 2016), Gkeyll (Shi et al., 2017, 2019; Mandell et al., 2020) and a
modified version of GENE (Pan et al., 2018). Both PIC and continuum methods
have their own advantages and disadvantages, as we detail briefly below. XGC1
is currently the most sophisticated code for the plasma boundary, capable of
simulating electrostatic gyrokinetic turbulence in realistic diverted
geometries, including neutral and atomic physics. It is critical to have at
least a few successful codes that can cross-check against each other on the
difficult problems in the edge/SOL, so as to give more confidence to the
predictions.

##### Particle-in-cell (PIC) approach

The first gyrokinetic simulation algorithms used particle-in-cell (PIC)
methods (Lee, 1983; Dimits & Lee, 1993; Parker & Lee, 1993; Denton &
Kotschenreuther, 1995; Dimits et al., 1996). In the PIC approach, the 5D phase
space is sampled with an ensemble of $N_{p}$ markers or ‘superparticles’,
representing some clump of physical particles with given position and
velocity. The markers are advanced through the domain according to the
characteristics of the gyrokinetic equation, while the electromagnetic fields
are evaluated and solved on a fixed three-dimensional grid. Communication
between the markers and fields requires interpolation: in order to solve the
field equations, the markers must be interpolated onto the grid positions so
that charges and currents can be computed; likewise, the effects of the fields
must be interpolated onto the marker positions to advance the particles. Since
the PIC method is essentially a Monte Carlo sampling technique, sampling noise
(which scales as $1/\sqrt{N_{p}}$) arises in moment calculation and can be
problematic in some cases (Nevins et al., 2005; Krommes, 2007; Wilkie &
Dorland, 2016). The sampling noise can be reduced but not completely
eliminated by $\delta f$ methods (Denton & Kotschenreuther, 1995). Various
other techniques have been used to reduce noise (Chen & Parker, 2007; Garbet
et al., 2010). Noise-related issues have contributed to the challenges of
handling electromagnetic fluctuations in PIC codes due to the Ampère
cancellation problem, as we discuss below. In general, PIC methods benefit
from being rather intuitive, fairly efficient, and easily parallelizable, with
straight-forward generalization to higher dimensionality. The lack of a need
for a velocity-space grid is attractive. Further, PIC methods automatically
guarantee positivity of the distribution function. PIC methods also have a
longer history to draw on than continuum methods.

##### Continuum (grid-based) approach

The first continuum gyrokinetic codes were developed some years later
(Kotschenreuther et al., 1995; Dorland et al., 2000; Jenko, 2000; Jenko &
Dorland, 2001; Candy & Waltz, 2003). In the continuum method, the full five-
dimensional gyrokinetic distribution function is discretized on a 5D phase-
space grid. Conventional numerical methods for solving partial differential
equations are then used to advance the distribution function according to the
gyrokinetic equation, including finite-difference, finite-volume,
(pseudo)spectral, finite-element, and discontinuous Galerkin (DG) methods.
Since the electromagnetic fields are discretized on the configuration-space
subset of the grid, no interpolation is required to solve the field equations,
only moment calculations. Continuum methods do not suffer from statistical
noise issues, which has contributed to the success of continuum codes in
including electromagnetic effects where some PIC codes have failed.
Discretization on a five-dimensional phase-space grid presents some additional
challenges for parallelization and memory handling, but these issues can still
be handled efficiently with well-designed schemes. In particular, continuum
schemes can make use of high-order methods that perform more calculations per
grid point and potentially enable faster convergence. One key disadvantage of
continuum methods is the strict Courant-Friedrichs-Lewy (CFL) stability limit
placed on the time step for explicit time-advance schemes, which can be
especially restrictive for electrostatic simulations (Lee, 1987) and highly-
collisional regimes. Another disadvantage of the continuum approach is that
the typical numerical methods used do not guarantee positivity of the
distribution function, which can cause numerical stability issues.

##### Including electromagnetic effects

Including electromagnetic effects in gyrokinetic simulations has proved
numerically and computationally challenging, both in the core and in the edge.
The so-called Ampère cancellation problem is one of the main numerical issues
that has troubled primarily PIC codes (Reynders, 1993; Cummings, 1994).
Various $\delta f$ PIC schemes to address the cancellation problem have been
developed and there are interesting recent advances in this area (Chen &
Parker, 2003; Mishchenko et al., 2004; Hatzky et al., 2007; Mishchenko et al.,
2014; Startsev & Lee, 2014; Bao et al., 2018). Meanwhile, some continuum
$\delta f$ core codes avoided the cancellation problem completely (Rewoldt et
al., 1987; Kotschenreuther et al., 1995), while others had to address somewhat
minor issues resulting from it (Jenko, 2000; Candy & Waltz, 2003). With
respect to the cancellation problem, one possible reason for the differences
might be that in continuum codes the fields and particles are discretized on
the same grid, whereas in PIC codes the particle positions do not coincide
with the field grid. Because particle positions are randomly located relative
to the field grid, one might need to be more careful in some way when treating
the interaction of the particles and electromagnetic fields.

Prior to the work described in this thesis, all published nonlinear
electromagnetic gyrokinetic results had focused on the core region, mostly
within the $\delta f$ formulation neglecting the $E_{\parallel}$ nonlinearity
(although the ORB5 PIC code includes the $E_{\parallel}$ nonlinearity and is
effectively full-$f$ (Lanti et al., 2019)). The XGC1 code is also full-$f$ and
is focused on both the core and the edge/SOL; it has an option for a
gyrokinetic ion/drift-fluid massless electron hybrid model (Hager et al.,
2017), with a fully kinetic implicit electromagnetic scheme based on Chen et
al. (2015) recently implemented and under further development (Ku et al.,
2018b). The GENE-X code is a recent extension of the core gyrokinetic code
GENE (Jenko, 2000) to a full-$f$ electromagnetic formulation similar to the
one presented in this work. GENE-X has now produced preliminary (but not yet
published) global electromagnetic gyrokinetic simulations including the SOL
and X-point. Other gyrokinetic codes working on the SOL are not yet
electromagnetic. To our knowledge, the results presented here were the first
nonlinear electromagnetic full-$f$ gyrokinetic turbulence simulations on open
field lines. The demonstration of full-$f$ electromagnetic capabilities,
handled in stable and efficient manner that does not significantly increase
the computational cost, is a major contribution of this thesis.

##### Handling diverted geometries with X-point

Another challenge is the magnetic geometry of the edge/SOL region, which
requires treatment of open and closed magnetic field-line regions and the
resulting plasma interactions with material walls on open field lines. The
X-point in a diverted geometry is an additional complication which makes the
use of field-aligned coordinates challenging.

Core gyrokinetic codes typically use such a field-aligned coordinate system,
which allows one to take advantage of the elongated nature of the turbulence
along the field line. This reduces the computational demands by allowing a
coarse grid along the field line. Unfortunately, field-aligned coordinate
systems are singular at the separatrix in diverted geometries due to the
presence of the X-point (Stegmeir et al., 2016). This has lead some codes to
abandon field-aligned coordinates altogether, opting instead for simpler
cylindrical coordinates. XGC1 uses a cylindrical coordinate system for the
particle motion and an unstructured field-following triangular mesh for the
field solver (Ku et al., 2016). BOUT++ uses multiple blocks, each with
separate field-aligned coordinates systems, that conform to the X-point but
still avoid it (Leddy et al., 2017). Recent interest has focused on ideas like
the flux-coordinate independent (FCI) approach, which abandons field- and
flux-aligned coordinates in the poloidal plane but retains a field-line-
following discretization of the parallel gradient operator to regain some of
the advantages of field-aligned domains (Hariri & Ottaviani, 2013; Hariri et
al., 2014; Stegmeir et al., 2016). This approach has been pioneered by the
GRILLIX fluid code (Stegmeir et al., 2016, 2018), and recently adopted by
several codes, including GDB, GBS, and GENE-X. Another recent approach by the
COGENT code uses a flux-aligned poloidal grid with controlled dealignment near
the X-point (McCorquodale et al., 2015; Dorf et al., 2016; Dorf & Dorr, 2020).
After breaking the toroidal direction into several blocks (wedges), a local
field-aligned coordinate system is used in each block. Interpolation (similar
to what is done in the FCI approach) is required to compute the parallel
derivatives between blocks.

Among gyrokinetic codes, currently only XGC1 (Ku et al., 2016) has published
results simulating turbulence in a three-dimensional diverted geometry with an
X-point. As mentioned above, the recently-developed GENE-X code is also
capable of including the X-point in global gyrokinetic turbulence simulations.

###  1.6 Thesis overview

First-principles modeling is crucial for understanding the dynamics in the
boundary plasma. In particular, there is a need for comprehensive gyrokinetic
simulations including electromagnetic effects. To this end, our efforts in
this thesis are focused on demonstrating and advancing the capabilities of the
gyrokinetic modules of the Gkeyll plasma simulation framework (which also
includes solver modules for the Vlasov–Maxwell system (Cagas et al., 2017;
Juno et al., 2018) and multi-moment fluid equations (Wang et al., 2015)).
Gkeyll was the first successful continuum gyrokinetic code on open field lines
due to the pioneering work of Shi (2017); Shi et al. (2017, 2019). In this
work we make the critical step of including electromagnetic fluctuations in
Gkeyll and demonstrating that this additional physics can be handled in a
stable and efficient manner. A primary goal is then to investigate how
electromagnetic effects can influence SOL turbulence and transport dynamics.

In Chapter 2 we derive the 5D full-$f$ electromagnetic gyrokinetic system in
Hamiltonian form using the symplectic ($v_{\parallel}$) formulation. In
Chapter 3 we describe an energy-conserving high-order discontinuous Galerkin
discretization scheme for the EMGK system that has been implemented in Gkeyll,
building on the electrostatic scheme of Shi (2017). In Chapter 4 we leverage
Gkeyll’s new electromagnetic capabilities to produce the first published
electromagnetic gyrokinetic results on open field lines. These simulations use
a simple helical scrape-off layer as a model of the SOL of the National
Spherical Torus Experiment (NSTX) experiment at PPPL, extending the
electrostatic results of Shi et al. (2019). Chapter 5 moves towards more
realistic geometry by describing and formulating field-aligned coordinate
systems for use in SOL geometries with magnetic shear and shaping. Chapter 6
develops a novel positivity-preserving DG scheme which can improve robustness
and accuracy of the simulations while maintaining critical conservation laws.
Finally, we conclude in Chapter 7 by reviewing the main results and describing
important areas for future work.

##  Chapter 2 Theoretical background: the full-$f$ electromagnetic gyrokinetic
system

Turbulence in strongly magnetized plasma is characterized by frequencies much
smaller than the ion cyclotron frequency $(\omega\ll\Omega_{i})$ and strong
anisotropy, with correlation lengths along the background field much longer
than perpendicular to it $(k_{\parallel}\ll k_{\perp})$. These two properties
are the basis for gyrokinetic theory, which reduces the full six-dimensional
(three position dimensions and three velocity dimensions) kinetic phase space
to five dimensions by averaging over the cyclotron motion. This eliminates a
velocity coordinate (the gyrophase angle) and results in a kinetic description
of the dynamics of charged gyro-rings. The first derivations of gyrokinetics
used a recursive procedure to generate an order-by-order asymptotic expansion,
yielding the local $\delta f$ gyrokinetic equation with the distribution
function separated into equilibrium ($F_{0}$) and perturbed ($\delta f$) parts
(Catto, 1978; Antonsen & Lane, 1980; Frieman & Chen, 1982; Abel et al., 2013).

Alternative approaches were later presented which derived (global, full-$f$)
gyrokinetics via Lagrangian and Hamiltonian Lie-transform perturbation methods
(Dubin et al., 1983; Hahm et al., 1988; Brizard & Hahm, 2007; Sugama, 2000).
We will take this latter approach in this chapter, using phase-space-
Lagrangian Lie perturbation methods (Littlejohn, 1983) to systematically
derive self-consistent, energy-conserving, global gyrokinetic equations. We
primarily follow Brizard & Hahm (2007) and references therein, but we have
also found a series of Ph.D. dissertations from the GENE group (Dannert, 2005;
Pueschel, 2009; Görler, 2009; Lapillonne, 2010; Told, 2012)111(Dannert, 2005)
is written in German but Google Translate does an admirable job of parsing it.
to be helpful for understanding parts of the derivation at a more introductory
level.

###  2.1 Gyrokinetic single-particle dynamics

The goal of this first section is to obtain gyrokinetic equations of motion
for single particles. We start from the description of a non-relativistic
charged particle with charge $q$, mass $m$, and velocity ${v}$ at position
${x}$ in the presence of an electrostatic potential
$\Phi(\mbox{\boldmath${x}$})$ and magnetic potential
$\mbox{\boldmath${A}$}(\mbox{\boldmath${x}$})$. The single-particle phase-
space Lagrangian is

$L(\mbox{\boldmath${x}$},\mbox{\boldmath${v}$},\dot{\mbox{\boldmath${x}$}},\dot{\mbox{\boldmath${v}$}},t)=\left(m\mbox{\boldmath${v}$}+q\mbox{\boldmath${A}$}\right)\boldsymbol{\cdot}\dot{\mbox{\boldmath${x}$}}-\left(\frac{1}{2}m|\mbox{\boldmath${v}$}|^{2}+q\Phi\right)\equiv\mbox{\boldmath${p}$}\boldsymbol{\cdot}\dot{\mbox{\boldmath${x}$}}-H,$
(2.1)

where $\mbox{\boldmath${p}$}=m\mbox{\boldmath${v}$}+q\mbox{\boldmath${A}$}$ is
the canonical momentum of the particle (in SI units),
$\mbox{\boldmath${v}$}=\dot{\mbox{\boldmath${x}$}}\equiv\textnormal{d}\mbox{\boldmath${x}$}/\textnormal{d}t$,
and $H$ is the Hamiltonian (for an introduction to the phase-space Lagrangian
formulation of mechanics, see section II of Cary & Brizard, 2009). This
Lagrangian will now be subjected to a series of coordinate transformations
that will separate the fast gyromotion from the guiding-center and gyrocenter
dynamics.

####  2.1.1 Ordering assumptions

The fundamental ordering requirement that allows us to effectively ignore the
fast gyromotion of a charged particle in a magnetic field is

$\frac{\omega}{\Omega}\sim\epsilon\ll 1,$ (2.2)

where $\omega$ is a typical frequency of interest, and $\Omega=qB/m$ is the
gyrofrequency for some species of interest. We will focus on drift wave
turbulence, which has characteristic frequencies $\omega\sim\omega_{*}\sim
v_{t}/L_{p}$, where $v_{t}$ is the thermal speed and $L_{p}$ is a
characteristic macroscopic scale length over which profile quantities vary.
This implies that

$\frac{\rho}{L_{p}}\sim\epsilon\ll 1$ (2.3)

is another small parameter, where $\rho=v_{t}/\Omega$ is the thermal
gyroradius. This ordering is valid in many tokamaks over a wide range of
experimental conditions, including the edge and scrape-off layer. The
frequency and length-scale orderings from Eqs. 2.2 and 2.3 comprise the
primary ordering in $\epsilon$. We must then decide how to deal with
fluctuations, flows, and magnetic geometry within the model, resulting in
additional parameters ordered with $\epsilon$ and $\epsilon^{2}$.

The standard nonlinear gyrokinetic ordering (Frieman & Chen, 1982) also
assumes small fluctuations,

$\frac{\delta f}{F_{0}}\sim\frac{q\delta\Phi}{T}\sim\epsilon,$

where $\delta f$ and $\delta\Phi$ are perturbations of the distribution
function and potential, respectively, and $F_{0}$ is the equilibrium
distribution function. Typical wavenumbers (of the perturbations) are then
ordered as

$k_{\parallel}\rho\sim\epsilon,\qquad k_{\perp}\rho\sim 1.$

This is the “$\delta f$” ordering, which is usually well-satisfied in the core
of tokamak plasmas, and has been used successfully to study core
microturbulence for many years (Parker et al., 1993a; Kotschenreuther et al.,
1995; Lin et al., 2000; Dimits et al., 2000; Dorland et al., 2000; Jenko,
2000; Jost et al., 2001; Candy & Waltz, 2003; Idomura et al., 2003; Watanabe &
Sugama, 2005; Jolliet et al., 2007; Idomura et al., 2008; Peeters et al.,
2009; Lanti et al., 2019). In the edge region, however, the $\delta f$
ordering is not strongly valid due to the presence of large fluctuations, even
though the fundamental frequency ordering is still satisfied.

The ordering can be generalized to allow larger perturbations by instead
taking a drift ordering (Dimits et al., 1992; Parra & Catto, 2008; Dimits,
2012)

$\epsilon_{V}\equiv\frac{v_{E}}{v_{t}}\simeq
k_{\perp}\rho\frac{q\Phi}{T}\sim\epsilon\ll 1,$ (2.4)

where $v_{E}$ is the $E\times B$ drift velocity from $\Phi$. Here, we have
defined a new ordering parameter $\epsilon_{V}$, which we take to be
$\mathcal{O}(\epsilon)$. This is commonly referred to as the “weak-flow”
ordering since it takes $E\times B$ flows to be small compared to the thermal
speed; this is generally satisfied in the edge and scrape-off layer (Brower et
al., 1987; Gohil et al., 1994; Zweben et al., 2015). By constraining gradients
of $\Phi$ instead of $\Phi$ itself, the weak-flow ordering 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. Another way to think about this is
that one can use $\Phi(\mbox{\boldmath${R}$})$ at the center of a gyro-orbit
(denoted ${R}$) as a reference point, and then one can require that the
variation of the potential energy around a gyro-orbit be small compared to the
kinetic energy,

$q\Phi(\mbox{\boldmath${R}$}+\mbox{\boldmath${\rho}$})-q\Phi(\mbox{\boldmath${R}$})\approx
q\mbox{\boldmath${\rho}$}\boldsymbol{\cdot}\nabla\Phi\ll T,$ (2.5)

which leads to the same criterion as above (Hammett, 2016). Here ${\rho}$ is
the gyroradius vector which points from the center of the gyro-orbit ${R}$ to
the particle location ${x}$ (it will be defined more precisely below). The
ordering has also been extended further to allow strong $E\times B$ flows of
order the thermal speed ($\epsilon_{V}\sim 1$) (Artun & Tang, 1994; Brizard,
1995; Hahm, 1996; Qin et al., 2007; Hahm et al., 2009; Dimits, 2010; Sharma &
McMillan, 2015; McMillan & Sharma, 2016; Sharma & McMillan, 2020); we will not
consider this here.

We also need an ordering parameter pertaining to the magnetic geometry. For
this, we introduce the equilibrium magnetic field scale length
$L_{B}\sim|\nabla_{\perp}\ln B|^{-1}$ $\sim R$, where $B$ is the background
magnetic field and $R$ is the major radius. In the core we expect $L_{B}\sim
L_{p}$, but the edge features much stronger gradients so that $L_{p}/L_{B}\sim
L_{p}/R\lesssim\rho/L_{p}\sim\epsilon$ is small (Gohil et al., 1994; Burrell
et al., 1994; Zweben et al., 2007). This leads us to a strong-gradient
ordering, in which we define an additional ordering parameter

$\epsilon_{B}\equiv\frac{\rho}{L_{B}}\sim\epsilon^{2}.$ (2.6)

This ordering has been employed in several edge gyrokinetic models (Hahm et
al., 2009; Dimits, 2012; Frei et al., 2020), as it allows the gyrokinetic
derivation to proceed fully consistently (up to second order in $\epsilon$) in
a two-step process, first by using the $\epsilon_{B}$ ordering to derive the
guiding-center motion, and then subsequently using the $\epsilon_{V}$ ordering
to introduce electromagnetic perturbations. Without the strong-gradient
ordering, $\epsilon_{B}\sim\epsilon_{V}\sim\epsilon$ are of the same order,
which is usually the case in the core. Even though many derivations still use
the two-step procedure in this case (see _e.g._ Brizard & Hahm, 2007), Parra &
Calvo (2011) have shown that the two-step procedure does not yield fully
consistent results at second order (and higher), as it misses terms of order
$\mathcal{O}(\epsilon_{V}\epsilon_{B})$ that involve both geometric effects
and field perturbations. The strong-gradient ordering
$\epsilon_{B}\sim\epsilon^{2}$ eliminates these concerns, since the geometry
only enters at even order in $\epsilon$.

With these ordering assumptions, we take the background magnetic field to be
$\mbox{\boldmath${B}$}_{0}=\nabla\times\mbox{\boldmath${A}$}_{0}$, with
$\mbox{\boldmath${A}$}_{0}\sim\mathcal{O}(1/\epsilon_{B})$ the background
vector potential. We do not include an $\mathcal{O}(1/\epsilon_{B})$
background electrostatic potential because this would violate the weak-flow
ordering. We then consider electromagnetic perturbations
$\mbox{\boldmath${A}$}_{1}$ and $\Phi_{1}$ of the form (Dimits et al., 1992)

$\displaystyle\mbox{\boldmath${A}$}_{1}(\mbox{\boldmath${x}$},t)=\mbox{\boldmath${A}$}_{1}(\mbox{\boldmath${R}$},t)+\epsilon_{V}\,\delta\mbox{\boldmath${A}$}_{1}(\mbox{\boldmath${R}$},\mbox{\boldmath${\rho}$},t)$
(2.7)
$\displaystyle\Phi_{1}(\mbox{\boldmath${x}$},t)=\Phi_{1}(\mbox{\boldmath${R}$},t)+\epsilon_{V}\,\delta\Phi_{1}(\mbox{\boldmath${R}$},\mbox{\boldmath${\rho}$},t),$
(2.8)

where the guiding-center component of the perturbation is $\mathcal{O}(1)$,
and the $\mathcal{O}(\epsilon_{V})$ part of the perturbation is the deviation
of the potential around the gyro-orbit, effectively giving the finite-Larmor-
radius (FLR) correction to the potential,

$\displaystyle\delta\mbox{\boldmath${A}$}_{1}$
$\displaystyle\equiv\mbox{\boldmath${A}$}_{1}(\mbox{\boldmath${x}$},t)-\mbox{\boldmath${A}$}_{1}(\mbox{\boldmath${R}$},t)\approx\mbox{\boldmath${\rho}$}\boldsymbol{\cdot}\nabla\mbox{\boldmath${A}$}_{1}\sim\frac{T}{q}\frac{B_{1\perp}}{B}\sim\frac{T}{q}\frac{v_{f}}{v_{t}}\sim\mathcal{O}(\epsilon_{V})$
(2.9) $\displaystyle\delta\Phi_{1}$
$\displaystyle\equiv\Phi_{1}(\mbox{\boldmath${x}$},t)-\Phi_{1}(\mbox{\boldmath${R}$},t)\approx\mbox{\boldmath${\rho}$}\boldsymbol{\cdot}\nabla\Phi_{1}\sim\frac{T}{q}\frac{v_{E}}{v_{t}}\sim\mathcal{O}(\epsilon_{V}),$
(2.10)

with $v_{f}=v_{\parallel}B_{1\perp}/B$ the magnetic flutter velocity. Thus the
total electromagnetic potentials can be written as

$\displaystyle\mbox{\boldmath${A}$}(\mbox{\boldmath${x}$})$
$\displaystyle=\frac{1}{\epsilon_{B}}\mbox{\boldmath${A}$}_{0}(\mbox{\boldmath${x}$})+\mbox{\boldmath${A}$}_{1}(\mbox{\boldmath${x}$})$
$\displaystyle=\left[\frac{1}{\epsilon_{B}}\mbox{\boldmath${A}$}_{0}(\mbox{\boldmath${R}$})+\mbox{\boldmath${\rho}$}\boldsymbol{\cdot}\nabla\mbox{\boldmath${A}$}_{0}(\mbox{\boldmath${R}$})+\mathcal{O}(\epsilon_{B})\right]+\mbox{\boldmath${A}$}_{1}(\mbox{\boldmath${R}$})+\epsilon_{V}\delta\mbox{\boldmath${A}$}_{1}(\mbox{\boldmath${R}$},\mbox{\boldmath${\rho}$})$
(2.11) $\displaystyle\Phi(\mbox{\boldmath${x}$})$
$\displaystyle=\Phi_{1}(\mbox{\boldmath${x}$})=\Phi_{1}(\mbox{\boldmath${R}$})+\epsilon_{V}\delta\Phi_{1}(\mbox{\boldmath${R}$},\mbox{\boldmath${\rho}$}),$
(2.12)

where we have Taylor expanded the background magnetic potential
$\mbox{\boldmath${A}$}_{0}(\mbox{\boldmath${x}$})=\mbox{\boldmath${A}$}_{0}(\mbox{\boldmath${R}$}+\mbox{\boldmath${\rho}$})$
to first order in $\epsilon_{B}$ around the guiding-center position ${R}$.
With these definitions, the Lagrangian from Eq. 2.1 can be written as

$\displaystyle{L}={L}_{0}+\epsilon_{V}\delta{L}$
$\displaystyle=\left[\frac{q}{\epsilon_{B}}\mbox{\boldmath${A}$}_{0}(\mbox{\boldmath${R}$})+q\mbox{\boldmath${\rho}$}\boldsymbol{\cdot}\nabla\mbox{\boldmath${A}$}_{0}(\mbox{\boldmath${R}$})+q\mbox{\boldmath${A}$}_{1}(\mbox{\boldmath${R}$})+m\mbox{\boldmath${v}$}+\mathcal{O}(\epsilon_{B})\right]\boldsymbol{\cdot}\dot{\mbox{\boldmath${x}$}}$
$\displaystyle\quad-\left[\frac{1}{2}mv^{2}+q\Phi_{1}(\mbox{\boldmath${R}$})\right]+\epsilon_{V}\left[q\,\delta\mbox{\boldmath${A}$}_{1}(\mbox{\boldmath${R}$},\mbox{\boldmath${\rho}$})\boldsymbol{\cdot}\dot{\mbox{\boldmath${x}$}}-q\,\delta\Phi_{1}(\mbox{\boldmath${R}$},\mbox{\boldmath${\rho}$})\right].$
(2.13)

Note that we will hereafter drop the subscript $0$ in the background magnetic
field $\mbox{\boldmath${B}$}_{0}$, unless it is needed for clarity, and simply
write ${B}$. For magnetic perturbations we will retain the subscript in
$\mbox{\boldmath${B}$}_{1}$, but these will frequently be expressed instead in
terms of the perturbed vector potential. The total magnetic field, including
background and perturbation, will be written as
$\mbox{\boldmath${B}$}_{0}+\mbox{\boldmath${B}$}_{1}$.

Writing the potentials in the form of Eqs. 2.11 and 2.12 has the advantage
that we can clearly see that FLR corrections are higher order. Thus, we could
self-consistently neglect FLR corrections by taking only the zeroth-order
terms. While we will proceed with the general derivation up to order
$\epsilon\sim\epsilon_{V}$ in this chapter, for simplicity we have elected to
neglect FLR corrections in the current implementation of the system in the
Gkeyll code. Extension to the more general system including FLR corrections is
left as important future work. The system that we solve in the current version
of Gkeyll is summarized in Section 2.3.

####  2.1.2 Transformation to guiding-center coordinates

Following Littlejohn (1983) and Cary & Brizard (2009), we first transform the
zeroth order Lagrangian $L_{0}$ to guiding-center coordinates,
$\mbox{\boldmath${Z}$}=(t,\mbox{\boldmath${R}$},v_{\parallel},\mu,\vartheta)$,
where ${R}$ is the guiding-center position, $v_{\parallel}$ is the guiding-
center velocity along the background magnetic field, $\mu=mv_{\perp}^{2}/(2B)$
is the lowest-order magnetic moment, and $\vartheta$ is the gyrophase angle.
In terms of the guiding-center coordinates, the particle coordinates
$\mbox{\boldmath${z}$}=(t,\mbox{\boldmath${x}$},\mbox{\boldmath${v}$})$
coordinates can be expressed (with the time coordinate $t$ staying the same)
as

$\displaystyle\mbox{\boldmath${x}$}=\mbox{\boldmath${R}$}+\epsilon_{B}\rho(\mbox{\boldmath${R}$})\mathbf{\hat{a}}(\mbox{\boldmath${R}$},\vartheta)=\mbox{\boldmath${R}$}+\epsilon_{B}\mbox{\boldmath${\rho}$}$
(2.14)
$\displaystyle\mbox{\boldmath${v}$}=v_{\parallel}\mathbf{\hat{b}}(\mbox{\boldmath${R}$})+v_{\perp}(\mbox{\boldmath${R}$})\mathbf{\hat{c}}(\mbox{\boldmath${R}$},\vartheta)+\epsilon_{V}\mbox{\boldmath${u}$}_{\perp}(\mbox{\boldmath${R}$},v_{\parallel})=v_{\parallel}\mathbf{\hat{b}}+\sqrt{\frac{2\mu
B}{m}}\mathbf{\hat{c}}+\epsilon_{V}\mbox{\boldmath${u}$}_{\perp},$ (2.15)

where
$\rho(\mbox{\boldmath${R}$})=v_{\perp}/\Omega(\mbox{\boldmath${R}$})=\sqrt{2m\mu/(q^{2}B(\mbox{\boldmath${R}$}))}$
is the gyroradius, $\mathbf{\hat{b}}=\mbox{\boldmath${B}$}/B$ is the unit
vector along the background field, and
$\mbox{\boldmath${u}$}_{\perp}(\mbox{\boldmath${R}$},v_{\parallel})$ is a to-
be-defined $\mathcal{O}(\epsilon_{V})$ velocity perpendicular to
$\mathbf{\hat{b}}$, taken to be the velocity of the reference frame; note that
$\mbox{\boldmath${u}$}_{\perp}$ is evaluated at the guiding-center position
${R}$ and is assumed to be gyrophase-independent. From standard guiding-center
motion, we might expect this reference frame velocity to be something like the
$E\times B$ drift velocity.222A moving reference frame is typically used in
strong-flow derivations of gyrokinetics, where $\epsilon_{V}\sim 1$ so that
all terms in Eq. 2.15 are the same order. In most weak-flow derivations (Frei
et al., 2020, is an exception), the reference frame is assumed to be
stationary ($\mbox{\boldmath${u}$}_{\perp}=0$), but here we allow for a slowly
moving frame. We also define

$\displaystyle\mathbf{\hat{a}}(\mbox{\boldmath${R}$},\vartheta)=\cos\vartheta\,\mbox{\boldmath${e}$}_{1}(\mbox{\boldmath${R}$})-\sin\vartheta\,\mbox{\boldmath${e}$}_{2}(\mbox{\boldmath${R}$})$
(2.16)
$\displaystyle\mathbf{\hat{c}}(\mbox{\boldmath${R}$},\vartheta)=\frac{\partial\mathbf{\hat{a}}(\mbox{\boldmath${R}$},\vartheta)}{\partial\vartheta}=-\sin\vartheta\,\mbox{\boldmath${e}$}_{1}(\mbox{\boldmath${R}$})-\cos\vartheta\,\mbox{\boldmath${e}$}_{2}(\mbox{\boldmath${R}$})$
(2.17)

to be unit vectors in the radial and tangential directions to the gyro-orbit
that rotate with $\vartheta$, where $\mbox{\boldmath${e}$}_{1}$ and
$\mbox{\boldmath${e}$}_{2}$ are some arbitrary pair of perpendicular unit
vectors in the plane perpendicular to the background field such that
$\mbox{\boldmath${e}$}_{1}\times\mbox{\boldmath${e}$}_{2}=\mathbf{\hat{b}}$.
Here and in the following, we will keep track of the order of various terms in
$\epsilon_{B}$ and $\epsilon_{V}$, but formally these parameters are equal to
unity so that the expressions retain the same dimensional form after taking
$\epsilon_{B}=\epsilon_{V}=1$.

Inserting the guiding-center coordinate transformations into Eq. 2.13, the
zeroth order Lagrangian in guiding-center coordinates is

${L}_{0}=\left[\frac{q}{\epsilon_{B}}\mbox{\boldmath${A}$}_{0}+q\mbox{\boldmath${A}$}_{1}+q(\mbox{\boldmath${\rho}$}\boldsymbol{\cdot}\nabla)\mbox{\boldmath${A}$}_{0}+mv_{\parallel}\mathbf{\hat{b}}+mv_{\perp}\mathbf{\hat{c}}\right]\boldsymbol{\cdot}\left(\dot{\mbox{\boldmath${R}$}}+\dot{\mbox{\boldmath${\rho}$}}\right)-H_{0},$
(2.18)

where $H_{0}=mv_{\parallel}^{2}/2+\mu B+q\Phi_{1}$ is the zeroth order
Hamiltonian, and spatially-varying quantities are evaluated at the guiding-
center position ${R}$ unless otherwise noted. Note that although the
gyroradius vector ${\rho}$ is $\mathcal{O}(\epsilon_{B})$, its time derivative
$\dot{\mbox{\boldmath${\rho}$}}$ is $\mathcal{O}(1)$, as it is given by

$\displaystyle\dot{\mbox{\boldmath${\rho}$}}$
$\displaystyle=\epsilon_{B}\left[(\dot{\mbox{\boldmath${R}$}}\boldsymbol{\cdot}\nabla)\mbox{\boldmath${\rho}$}+\dot{\mu}\frac{\partial\mbox{\boldmath${\rho}$}}{\partial\mu}\right]+\dot{\vartheta}\frac{\partial\mbox{\boldmath${\rho}$}}{\partial\vartheta}=\epsilon_{B}\left[\frac{\mbox{\boldmath${\rho}$}}{2B}\dot{\mbox{\boldmath${R}$}}\boldsymbol{\cdot}\nabla
B+\frac{1}{qv_{\perp}}\mathbf{\hat{a}}\dot{\mu}\right]+\frac{v_{\perp}}{\Omega}\mathbf{\hat{c}}\dot{\vartheta}$
$\displaystyle=\frac{v_{\perp}}{\Omega}\mathbf{\hat{c}}\dot{\vartheta}+\mathcal{O}(\epsilon_{B}),$
(2.19)

since $\dot{\vartheta}=\Omega\sim\epsilon_{B}^{-1}$ (Cary & Brizard, 2009).

We then make a series of gauge transformations to eliminate the dependence on
the gyrophase to lowest order in $\epsilon_{B}$, following Littlejohn (1983).
These gauge transformations take the form of adding a total time derivative to
the Lagrangian, ${L}\rightarrow{L}+\dot{S}$, which does not affect the
equations of motion. Taking
$S=-q\mbox{\boldmath${\rho}$}\boldsymbol{\cdot}(\mbox{\boldmath${A}$}_{0}+\epsilon_{B}\mbox{\boldmath${A}$}_{1})$,
so that

$\displaystyle\dot{S}$
$\displaystyle=-\left(\frac{q}{\epsilon_{B}}\mbox{\boldmath${A}$}_{0}+q\mbox{\boldmath${A}$}_{1}\right)\boldsymbol{\cdot}\dot{\mbox{\boldmath${\rho}$}}-q\dot{\mbox{\boldmath${R}$}}\boldsymbol{\cdot}\nabla\mbox{\boldmath${A}$}_{0}\boldsymbol{\cdot}\mbox{\boldmath${\rho}$}+\mathcal{O}(\epsilon_{B}),$
(2.20)

the Lagrangian can be transformed as

$\displaystyle{L}_{0}\rightarrow{L}_{0}+\dot{S}$
$\displaystyle=\left(\frac{q}{\epsilon_{B}}\mbox{\boldmath${A}$}_{0}^{*}+q\mbox{\boldmath${A}$}_{1}\right)\boldsymbol{\cdot}\dot{\mbox{\boldmath${R}$}}+\left[q\mbox{\boldmath${\rho}$}\boldsymbol{\cdot}\nabla\mbox{\boldmath${A}$}_{0}+mv_{\parallel}\mathbf{\hat{b}}+mv_{\perp}\mathbf{\hat{c}}\right]\boldsymbol{\cdot}\dot{\mbox{\boldmath${\rho}$}}-H_{0},$
(2.21)

where cancellations resulted from noting that
$q\left[\mbox{\boldmath${\rho}$}\boldsymbol{\cdot}\nabla\mbox{\boldmath${A}$}_{0}-(\nabla\mbox{\boldmath${A}$}_{0})\boldsymbol{\cdot}\mbox{\boldmath${\rho}$}\right]=-q\mbox{\boldmath${\rho}$}\times(\nabla\times\mbox{\boldmath${A}$}_{0})=-q\mbox{\boldmath${\rho}$}\times\mbox{\boldmath${B}$}=-mv_{\perp}\mathbf{\hat{c}}$,
and we have defined the modified vector potential

$\mbox{\boldmath${A}$}_{0}^{*}\equiv\mbox{\boldmath${A}$}_{0}+\epsilon_{B}\frac{mv_{\parallel}}{q}\mathbf{\hat{b}}.$
(2.22)

Recognizing $\mbox{\boldmath${A}$}_{0}^{*}$ as the gyroaveraged canonical
momentum from the background field, we can see that this gauge transformation
has effectively gyroaveraged the first term in Eq. 2.13.

We then make an additional gauge transformation with
$S=-\epsilon_{B}(q/2)(\mbox{\boldmath${\rho}$}\boldsymbol{\cdot}\nabla)\mbox{\boldmath${A}$}_{0}\boldsymbol{\cdot}\mbox{\boldmath${\rho}$}$,
which gives

$\displaystyle\dot{S}=-\frac{q}{2}\left[(\nabla\mbox{\boldmath${A}$}_{0})\boldsymbol{\cdot}\mbox{\boldmath${\rho}$}+\left(\mbox{\boldmath${\rho}$}\boldsymbol{\cdot}\nabla\right)\mbox{\boldmath${A}$}_{0}\right]\boldsymbol{\cdot}\dot{\mbox{\boldmath${\rho}$}}-\epsilon_{B}\frac{q}{2}\mbox{\boldmath${\rho}$}\boldsymbol{\cdot}\frac{\textnormal{d}(\nabla\mbox{\boldmath${A}$}_{0})}{\textnormal{d}t}\boldsymbol{\cdot}\mbox{\boldmath${\rho}$},$
(2.23)

where we will drop the last term because it is higher order. The Lagrangian is
then transformed as

$\displaystyle{L}_{0}\rightarrow{L}_{0}+\dot{S}$
$\displaystyle=\left(\frac{q}{\epsilon_{B}}\mbox{\boldmath${A}$}_{0}^{*}+q\mbox{\boldmath${A}$}_{1}\right)\boldsymbol{\cdot}\dot{\mbox{\boldmath${R}$}}+\frac{m\mu}{q}\dot{\vartheta}-\left[\frac{1}{2}mv_{\parallel}^{2}+\mu
B+q\Phi_{1}\right]+\mathcal{O}(\epsilon_{B}).$ (2.24)

This Lagrangian describes the motion of charged particles in a strong
background magnetic field and slowly varying electromagnetic potentials (with
no variation at the gyroradius scale), and it could be used to derive the
drift-kinetic Vlasov equation (up to zeroth order in $\epsilon_{B}$). Since
Eq. 2.24 is independent of gyrophase, Noether’s theorem gives that the
quantity $\partial{L}_{0}/\partial\dot{\vartheta}=m\mu/q$ is a constant, which
is confirmation that $\mu$ is an adiabatic invariant in the absence of
electromagnetic perturbations on the scale of the gyroradius (to lowest
order). For an alternative derivation of the guiding-center Lagrangian, which
eliminates the gyrophase dependence via gyroaveraging instead of gauge
transformations, see Helander & Sigmar (2002).

####  2.1.3 Transformation to gyrocenter coordinates

Now we must account for the variations of the electromagnetic fields on the
scale of the gyroradius, which are contained in $\delta L$ from Eq. 2.13,

$\displaystyle\delta{L}$
$\displaystyle=\left[q\,\delta\mbox{\boldmath${A}$}_{1}(\mbox{\boldmath${R}$},\mbox{\boldmath${\rho}$})+m\mbox{\boldmath${u}$}_{\perp}\right]\boldsymbol{\cdot}(\dot{\mbox{\boldmath${R}$}}+\dot{\mbox{\boldmath${\rho}$}})-\left[q\,\delta\Phi_{1}(\mbox{\boldmath${R}$},\mbox{\boldmath${\rho}$})+m\mbox{\boldmath${v}$}_{\perp}\boldsymbol{\cdot}\mbox{\boldmath${u}$}_{\perp}+\frac{\epsilon_{V}}{2}mu_{\perp}^{2}\right]$
$\displaystyle=q\,\delta\mbox{\boldmath${A}$}_{1}^{*}\boldsymbol{\cdot}\dot{\mbox{\boldmath${R}$}}+\frac{m}{B}\mbox{\boldmath${v}$}_{\perp}\boldsymbol{\cdot}\delta\mbox{\boldmath${A}$}_{1}^{*}\,\dot{\vartheta}-\left[q\,\delta\Phi_{1}+m\mbox{\boldmath${v}$}_{\perp}\boldsymbol{\cdot}\mbox{\boldmath${u}$}_{\perp}+\frac{\epsilon_{V}}{2}mu_{\perp}^{2}\right]+\mathcal{O}(\epsilon_{B}),$
(2.25)

where we have defined

$\delta\mbox{\boldmath${A}$}_{1}^{*}\equiv\delta\mbox{\boldmath${A}$}_{1}(\mbox{\boldmath${R}$},\mbox{\boldmath${\rho}$})+(m/q)\mbox{\boldmath${u}$}_{\perp}(\mbox{\boldmath${R}$},v_{\parallel}).$
(2.26)

Since the perturbations $\delta\mbox{\boldmath${A}$}_{1}$ and $\delta\Phi_{1}$
depend on ${\rho}$, we have reintroduced gyrophase dependence in the
Lagrangian and broken $\mu$ conservation. Unlike in the lowest-order case
above, we cannot simply use gauge transformations to eliminate the gyrophase
dependence at this order because the perturbations depend non-trivially on
${\rho}$. Instead, we use another kind of coordinate transformation known as a
Lie transform (for details, see Cary, 1981; Littlejohn, 1982; Cary &
Littlejohn, 1983). The Lie transform offers a systematic method for making
perturbative coordinate transformations and computing the resulting changes in
functions of those coordinates.

For these transformations, it will be convenient to adopt the Poincaré-Cartan
one-form formalism (see _e.g._ Cary & Littlejohn, 1983), where the one-form
$\gamma(\mbox{\boldmath${Z}$})$ is defined via the action integral

$\mathcal{I}=\int
L(\mbox{\boldmath${Z}$})\,\textnormal{d}t=\int\gamma(\mbox{\boldmath${Z}$}).$
(2.27)

Here,

$\displaystyle\gamma$
$\displaystyle=\gamma_{Z^{\alpha}}\textnormal{d}Z^{\alpha}=L\,\textnormal{d}t$
(2.28)

for $Z^{\alpha}$ (with $\alpha=0,1,\dots,6$) the components of the extended
phase-space coordinates that include time as the zeroth element so that
$\gamma_{t}=-H$, the Hamiltonian. The remaining components, $\gamma_{Z^{i}}$
with $i=1,\dots,6$, are together called the _symplectic_ component of the one-
form.

We will define $\gamma=\gamma_{0}+\epsilon_{V}\gamma_{1}$, with

$\displaystyle\gamma_{0}$ $\displaystyle\equiv
L_{0}\,\textnormal{d}t=q\left(\mbox{\boldmath${A}$}_{0}^{*}+\mbox{\boldmath${A}$}_{1}\right)\boldsymbol{\cdot}\textnormal{d}\mbox{\boldmath${R}$}+\frac{m\mu}{q}\textnormal{d}\vartheta-\left[\frac{1}{2}mv_{\parallel}^{2}+\mu
B+q\Phi_{1}\right]\textnormal{d}t$ (2.29) $\displaystyle\gamma_{1}$
$\displaystyle\equiv\delta
L\,\textnormal{d}t=q\,\delta\mbox{\boldmath${A}$}_{1}^{*}\boldsymbol{\cdot}\textnormal{d}\mbox{\boldmath${R}$}+\frac{m}{B}\mbox{\boldmath${v}$}_{\perp}\boldsymbol{\cdot}\delta\mbox{\boldmath${A}$}_{1}^{*}\,\textnormal{d}\vartheta-\left[q\,\delta\Phi_{1}+m\mbox{\boldmath${v}$}_{\perp}\boldsymbol{\cdot}\mbox{\boldmath${u}$}_{\perp}+\frac{\epsilon_{V}}{2}mu_{\perp}^{2}\right]\textnormal{d}t,$
(2.30)

where all quantities are evaluated at ${R}$ (although
$\delta\mbox{\boldmath${A}$}_{1}$ and $\delta\Phi_{1}$ still also depend on
${\rho}$), and we have dropped the $1/\epsilon_{B}$ ordering parameter on
$\mbox{\boldmath${A}$}_{0}^{*}$. The Lie transformation that yields the
gyrocenter one-form
$\Gamma=\Gamma_{0}+\epsilon_{V}\Gamma_{1}+\epsilon_{V}^{2}\Gamma_{2}+\dots$ is
given (up to second order) by

$\displaystyle\Gamma_{0}$ $\displaystyle=\gamma_{0}+\textnormal{d}S_{0}$
(2.31) $\displaystyle\Gamma_{1}$
$\displaystyle=\gamma_{1}-\mathscr{L}_{1}\gamma_{0}+\textnormal{d}S_{1}$
(2.32) $\displaystyle\Gamma_{2}$
$\displaystyle=\gamma_{2}-\mathscr{L}_{2}\gamma_{0}-\frac{1}{2}\mathscr{L}_{1}(\gamma_{1}+\Gamma_{1})+\textnormal{d}S_{2},$
(2.33)

where $\mathscr{L}_{n}$ denotes the $n$th order Lie derivative and $S_{n}$ is
an arbitrary $n$th order scalar gauge function. The Lie derivative acting on a
one-form $\gamma$ is given by333The expression in Eq. 2.34 is only part of the
formal Lie derivative. There is another part of the form
$\textnormal{d}(G_{n}^{\beta}\gamma_{Z^{\beta}})$, but this part can be
absorbed into $\textnormal{d}S$ in Eqs. 2.32 and 2.33 because $S$ can be
chosen arbitrarily.

$\mathscr{L}_{n}\gamma=G_{n}^{\beta}\omega_{\alpha\beta}\textnormal{d}Z^{\alpha}=G_{n}^{\beta}\left(\frac{\partial\gamma_{Z^{\alpha}}}{\partial
Z^{\beta}}-\frac{\partial\gamma_{Z^{\beta}}}{\partial
Z^{\alpha}}\right)\textnormal{d}Z^{\alpha},$ (2.34)

where the functions $G_{n}^{\beta}$ are the components of the $n$th order
generating vector field of the Lie transform, and $\omega_{\alpha\beta}$ are
the elements of the Lagrange tensor. With this, the first order gyrocenter
one-form can be rewritten as

${\Gamma}_{1,Z^{\alpha}}=\gamma_{1,Z^{\alpha}}-G_{1}^{\beta}\left(\frac{\partial\gamma_{0,Z^{\alpha}}}{\partial
Z^{\beta}}-\frac{\partial\gamma_{0,Z^{\beta}}}{\partial
Z^{\alpha}}\right)+\frac{\partial S_{1}}{\partial Z^{\alpha}}.$ (2.35)

The goal now is to find generating functions $\mbox{\boldmath${G}$}_{n}$ and
gauge functions $S_{n}$ such that the gyrocenter one-form no longer depends on
the gyrophase at each order. At zeroth order, we will simply take $S_{0}=0$,
so that

${\Gamma}_{0}=\gamma_{0},$ (2.36)

since we have already removed gyrophase dependence from the zeroth order one-
form, Eq. 2.29.

At first order, we can use Eq. 2.35 to compute

$\displaystyle\Gamma$
$\displaystyle=q\,\delta\mbox{\boldmath${A}$}_{1}^{*}+qG_{1}^{\mathbold{R}}\times\mbox{\boldmath${B}$}^{*}-mG_{1}^{v_{\parallel}}\mathbf{\hat{b}}+\nabla
S_{1}$ (2.37) $\displaystyle\Gamma_{1,v_{\parallel}}$
$\displaystyle=m\mathbf{\hat{b}}\boldsymbol{\cdot}G_{1}^{\mathbold{R}}+\frac{\partial
S_{1}}{\partial v_{\parallel}}$ (2.38) $\displaystyle\Gamma_{1,\mu}$
$\displaystyle=\frac{m}{q}G_{1}^{\vartheta}+\frac{\partial
S_{1}}{\partial\mu}$ (2.39) $\displaystyle\Gamma_{1,\vartheta}$
$\displaystyle=\frac{m}{B}\mbox{\boldmath${v}$}_{\perp}\boldsymbol{\cdot}\delta\mbox{\boldmath${A}$}_{1}^{*}-\frac{m}{q}G_{1}^{\mu}+\frac{\partial
S_{1}}{\partial\vartheta}$ (2.40) $\displaystyle\Gamma_{1,t}$
$\displaystyle=-q\,\delta\Phi_{1}-\epsilon_{V}qG_{1}^{\mathbold{R}}\boldsymbol{\cdot}\mbox{\boldmath${E}$}_{1}^{*}+G_{1}^{v_{\parallel}}mv_{\parallel}+G_{1}^{\mu}B+\frac{\partial
S_{1}}{\partial t},$ (2.41)

where
$\mbox{\boldmath${E}$}_{1}^{*}=-\epsilon_{V}(\nabla\Phi_{1}+\partial\mbox{\boldmath${A}$}_{1}/\partial
t)-\epsilon_{B}(\mu/q)\nabla
B=\epsilon_{V}\mbox{\boldmath${E}$}_{1}-\epsilon_{B}(\mu/q)\nabla
B\sim\mathcal{O}(\epsilon_{V})$ and
$\mbox{\boldmath${B}$}^{*}=\nabla\times\mbox{\boldmath${A}$}_{0}^{*}+\epsilon_{V}\nabla\times\mbox{\boldmath${A}$}_{1}=\mbox{\boldmath${B}$}_{0}^{*}+\epsilon_{V}\mbox{\boldmath${B}$}_{1}$.
We have also taken $G_{1}^{t}=0$ since we do not need to make a coordinate
transformation in time.

We now have some freedom to choose the $\mbox{\boldmath${G}$}_{n}$ and $S_{n}$
to simplify the form of the gyrocenter Lagrangian. To this end, we choose to
enforce $\Gamma_{1,v_{\parallel}}=\Gamma_{1,\mu}=\Gamma_{1,\vartheta}=0$,
which gives

$\displaystyle\mathbf{\hat{b}}\boldsymbol{\cdot}G_{1}^{\mathbold{R}}=-\frac{1}{m}\frac{\partial
S_{1}}{\partial v_{\parallel}},\qquad
G_{1}^{\vartheta}=-\frac{q}{m}\frac{\partial S_{1}}{\partial\mu},\qquad
G_{1}^{\mu}=\frac{q}{B}\mbox{\boldmath${v}$}_{\perp}\boldsymbol{\cdot}\delta\mbox{\boldmath${A}$}_{1}^{*}+\frac{q}{m}\frac{\partial
S_{1}}{\partial\vartheta}.$ (2.42)

Dotting Eq. 2.37 with $\mbox{\boldmath${B}$}^{*}$ gives

$\displaystyle\mbox{\boldmath${B}$}^{*}\boldsymbol{\cdot}\Gamma_{1,\mathbold{R}}=q\mbox{\boldmath${B}$}^{*}\boldsymbol{\cdot}\delta\mbox{\boldmath${A}$}_{1}^{*}-mG_{1}^{v_{\parallel}}B_{\parallel}^{*}+\mbox{\boldmath${B}$}^{*}\boldsymbol{\cdot}\nabla
S_{1},$ (2.43)

while crossing Eq. 2.37 with $\mathbf{\hat{b}}$ gives

$\displaystyle\mathbf{\hat{b}}\times\Gamma_{1,\mathbold{R}}$
$\displaystyle=q\mathbf{\hat{b}}\times\delta\mbox{\boldmath${A}$}_{1}^{*}+qB_{\parallel}^{*}G_{1}^{\mathbold{R}}+\frac{q}{m}\mbox{\boldmath${B}$}^{*}\frac{\partial
S_{1}}{\partial v_{\parallel}}+\mathbf{\hat{b}}\times\nabla S_{1},$ (2.44)

so that we have

$\displaystyle G_{1}^{v_{\parallel}}$
$\displaystyle=\frac{\mbox{\boldmath${B}$}^{*}}{mB_{\parallel}^{*}}\boldsymbol{\cdot}\left(-\Gamma_{1,\mathbold{R}}+\nabla
S_{1}+q\delta\mbox{\boldmath${A}$}_{1}^{*}\right)$ (2.45) $\displaystyle
G_{1}^{\mathbold{R}}$
$\displaystyle=-\frac{\mbox{\boldmath${B}$}^{*}}{mB_{\parallel}^{*}}\frac{\partial
S_{1}}{\partial
v_{\parallel}}-\frac{\mathbf{\hat{b}}}{qB_{\parallel}^{*}}\times\left(-\Gamma_{1,\mathbold{R}}+\nabla
S_{1}+q\delta\mbox{\boldmath${A}$}_{1}^{*}\right),$ (2.46)

where
$B_{\parallel}^{*}\equiv\mathbf{\hat{b}}\boldsymbol{\cdot}\mbox{\boldmath${B}$}^{*}\approx
B+\mathcal{O}(\epsilon_{B})$. The first order gyrocenter Hamiltonian is then

$\displaystyle H_{1}$
$\displaystyle=-\Gamma_{1,t}=q\,\delta\Phi_{1}+m\mbox{\boldmath${v}$}_{\perp}\boldsymbol{\cdot}\mbox{\boldmath${u}$}_{\perp}+\epsilon_{V}\frac{1}{2}mu_{\perp}^{2}-\epsilon_{V}\frac{q}{mB_{\parallel}^{*}}\frac{\partial
S_{1}}{\partial
v_{\parallel}}\mbox{\boldmath${B}$}^{*}\boldsymbol{\cdot}\mbox{\boldmath${E}$}_{1}^{*}$
$\displaystyle\qquad-\frac{\epsilon_{V}\mbox{\boldmath${E}$}_{1}^{*}\times\mathbf{\hat{b}}+v_{\parallel}\mbox{\boldmath${B}$}^{*}}{B_{\parallel}^{*}}\boldsymbol{\cdot}\left(-\Gamma_{1,\mathbold{R}}+\nabla
S_{1}+q\,\delta\mbox{\boldmath${A}$}_{1}^{*}\right)-q\mbox{\boldmath${v}$}_{\perp}\boldsymbol{\cdot}\delta\mbox{\boldmath${A}$}_{1}^{*}-\Omega\frac{\partial
S_{1}}{\partial\vartheta}-\frac{\partial S_{1}}{\partial t}$
$\displaystyle=q\left(\delta\Phi_{1}-\mbox{\boldmath${v}$}_{\perp}\boldsymbol{\cdot}\delta\mbox{\boldmath${A}$}_{1}\right)+\epsilon_{V}\frac{1}{2}mu_{\perp}^{2}-\frac{\epsilon_{V}\mbox{\boldmath${E}$}_{1}^{*}\times\mathbf{\hat{b}}+v_{\parallel}\mbox{\boldmath${B}$}^{*}}{B_{\parallel}^{*}}\boldsymbol{\cdot}\left(-\Gamma_{1,\mathbold{R}}+q\,\delta\mbox{\boldmath${A}$}_{1}^{*}\right)-\frac{\textnormal{d}S_{1}}{\textnormal{d}t}.$
(2.47)

We now choose the gauge function $S_{1}$ to cancel the gyrophase dependence in
Eq. 2.47. We will leave $\Gamma_{1,\mathbold{R}}$ unspecified for now; by
construction, the gyrocenter one-form will be gyrophase-independent, so the
choice of $\Gamma_{1,\mathbold{R}}$ will not affect the choice of $S_{1}$. We
can define $S_{1}$ as the solution to

$\displaystyle\frac{\textnormal{d}S_{1}}{\textnormal{d}t}\approx\Omega\frac{\partial
S_{1}}{\partial\vartheta}+\mathcal{O}(\epsilon_{V})=q\left(\widetilde{\delta\Phi_{1}}-\widetilde{\mbox{\boldmath${v}$}_{\perp}\boldsymbol{\cdot}\delta\mbox{\boldmath${A}$}_{1}}\right)-\frac{\epsilon_{V}\mbox{\boldmath${E}$}_{1}^{*}\times\mathbf{\hat{b}}+v_{\parallel}\mbox{\boldmath${B}$}^{*}}{B_{\parallel}^{*}}\boldsymbol{\cdot}q\,\widetilde{\delta\mbox{\boldmath${A}$}_{1}},$
(2.48)

where

$\widetilde{A}=A(\mbox{\boldmath${R}$}+\mbox{\boldmath${\rho}$})-\langle
A\rangle$ (2.49)

is the gyrophase-dependent part of a quantity
$A(\mbox{\boldmath${R}$}+\mbox{\boldmath${\rho}$})$, and
$\langle\boldsymbol{\cdot}\rangle$ denotes a gyroaverage, defined by

$\langle
A(\mbox{\boldmath${R}$}+\mbox{\boldmath${\rho}$})\rangle\equiv\frac{1}{2\pi}\int_{0}^{2\pi}A(\mbox{\boldmath${R}$}+\mbox{\boldmath${\rho}$})\,\textnormal{d}\vartheta.$
(2.50)

Noting that $\widetilde{\delta\Phi_{1}}=\widetilde{\Phi_{1}}$,
$\widetilde{\delta\mbox{\boldmath${A}$}_{1}}=\widetilde{\mbox{\boldmath${A}$}_{1}}$,
and
$\widetilde{\mbox{\boldmath${v}$}_{\perp}\boldsymbol{\cdot}\delta\mbox{\boldmath${A}$}_{1}}=\widetilde{\mbox{\boldmath${v}$}_{\perp}\boldsymbol{\cdot}\mbox{\boldmath${A}$}_{1}}-\mbox{\boldmath${v}$}_{\perp}\boldsymbol{\cdot}\mbox{\boldmath${A}$}_{1}$,
the gauge function $S_{1}$ becomes

$\displaystyle
S_{1}\approx\frac{q}{\Omega}\int^{\vartheta}\left[\widetilde{\Phi_{1}}-\widetilde{\mbox{\boldmath${v}$}_{\perp}\boldsymbol{\cdot}\mbox{\boldmath${A}$}_{1}}+\mbox{\boldmath${v}$}_{\perp}\boldsymbol{\cdot}\mbox{\boldmath${A}$}_{1}-\frac{\epsilon_{V}\mbox{\boldmath${E}$}_{1}^{*}\times\mathbf{\hat{b}}+v_{\parallel}\mbox{\boldmath${B}$}^{*}}{B_{\parallel}^{*}}\boldsymbol{\cdot}\widetilde{\mbox{\boldmath${A}$}_{1}}\right]\textnormal{d}\vartheta^{\prime},$
(2.51)

where we will take the solution with $\langle S_{1}\rangle=0$ (this is
required to prevent $S_{1}$ from becoming unbounded; see Cary & Littlejohn
(1983)). The Hamiltonian then becomes

$\displaystyle H_{1}$
$\displaystyle=q\left[\langle\Phi_{1}-\mbox{\boldmath${v}$}_{\perp}\boldsymbol{\cdot}\mbox{\boldmath${A}$}_{1}\rangle-\Phi_{1}\right]+\epsilon_{V}\frac{1}{2}mu_{\perp}^{2}$
$\displaystyle\quad-\frac{\epsilon_{V}\mbox{\boldmath${E}$}_{1}\times\mathbf{\hat{b}}+\epsilon_{B}(\mu/q)\mathbf{\hat{b}}\times\nabla
B+v_{\parallel}\mbox{\boldmath${B}$}^{*}}{B_{\parallel}^{*}}\boldsymbol{\cdot}\left(-\Gamma_{1,\mathbold{R}}+q\langle\mbox{\boldmath${A}$}_{1}\rangle-q\mbox{\boldmath${A}$}_{1}+m\mbox{\boldmath${u}$}_{\perp}\right),$
(2.52)

so that the gyrocenter one-form is

$\Gamma=\Gamma_{0}+\epsilon_{V}\Gamma_{1}=q\left(\mbox{\boldmath${A}$}_{0}^{*}+\mbox{\boldmath${A}$}_{1}+\epsilon_{V}\Gamma_{1,\mathbold{R}}\right)\boldsymbol{\cdot}\textnormal{d}\mbox{\boldmath${R}$}+\frac{m\mu}{q}\textnormal{d}\vartheta-\left[\frac{1}{2}mv_{\parallel}^{2}+\mu
B+q\Phi_{1}+\epsilon_{V}H_{1}\right]\textnormal{d}t.$ (2.53)

Now we must choose $\Gamma_{1,\mathbold{R}}$. One option is to use
$\Gamma_{1,\mathbold{R}}=-q\mbox{\boldmath${A}$}_{1}$. This would eliminate
$\mbox{\boldmath${A}$}_{1}$ from the symplectic part of the one-form, opting
instead to move all dependence on field perturbations to the Hamiltonian. This
is known as the “Hamiltonian” formulation of electromagnetic gyrokinetics
(Brizard & Hahm, 2007), so-named because all field perturbations reside in the
Hamiltonian. This approach has the advantage that the equations of motion do
not contain explicit time derivatives of the magnetic potential, which can be
advantageous in some discretization schemes. As a result, however, the
parallel momentum coordinate becomes the canonical momentum, which depends on
the perturbed magnetic potential. In other discretization schemes (namely, in
the one that we pursue in Chapter 3) having the perturbed magnetic potential
in the Hamiltonian can be disadvantageous.

Thus we will take another approach, known as the “symplectic” formulation, so-
named because the symplectic part of the gyrocenter one-form is allowed to
retain gyrophase-independent parts of the perturbed fields. In this approach,
we eliminate the dependence on $\mbox{\boldmath${A}$}_{1}$ from the
Hamiltonian at first order. This results in the parallel momentum coordinate
remaining the kinetic momentum. Thus we take

$\displaystyle\Gamma_{1,\mathbold{R}}=q\langle\delta\mbox{\boldmath${A}$}_{1}\rangle=q\langle\mbox{\boldmath${A}$}_{1}(\mbox{\boldmath${R}$}+\mbox{\boldmath${\rho}$})-\mbox{\boldmath${A}$}_{1}(\mbox{\boldmath${R}$})\rangle=q\langle\mbox{\boldmath${A}$}_{1}\rangle-q\mbox{\boldmath${A}$}_{1},$
(2.54)

where note that
$\Gamma_{1,\mathbold{R}}\neq-\widetilde{\mbox{\boldmath${A}$}}_{1}$, since the
non-averaged term in Eq. 2.54 is evaluated at ${R}$, not
$\mbox{\boldmath${R}$}+\mbox{\boldmath${\rho}$}$. With this choice, the
gyrocenter one-form is given by

$\Gamma=q\left(\mbox{\boldmath${A}$}_{0}^{*}+\langle\mbox{\boldmath${A}$}_{1}\rangle\right)\boldsymbol{\cdot}\textnormal{d}\mbox{\boldmath${R}$}+\frac{m\mu}{q}\textnormal{d}\vartheta-\mathcal{H}\textnormal{d}t,$
(2.55)

with the total gyrocenter Hamiltonian

$\displaystyle\mathcal{H}$ $\displaystyle=\frac{1}{2}mv_{\parallel}^{2}+\mu
B+q\langle\Phi_{1}-\mbox{\boldmath${v}$}_{\perp}\boldsymbol{\cdot}\mbox{\boldmath${A}$}_{1}\rangle+\epsilon_{V}^{2}\frac{1}{2}mu_{\perp}^{2}$
$\displaystyle\qquad-\left[\frac{v_{\parallel}\mbox{\boldmath${B}$}_{0}^{*}+\epsilon_{V}\left(\mbox{\boldmath${E}$}_{1}\times\mathbf{\hat{b}}+v_{\parallel}\mbox{\boldmath${B}$}_{1\perp}\right)+\epsilon_{B}(\mu/q)\mathbf{\hat{b}}\times\nabla
B}{B_{\parallel}^{*}}\right]\boldsymbol{\cdot}\epsilon_{V}m\mbox{\boldmath${u}$}_{\perp},$
(2.56)

and
$\mbox{\boldmath${B}$}_{1\perp}=(\nabla\times\mbox{\boldmath${A}$}_{1})_{\perp}$.
By taking the velocity $\mbox{\boldmath${u}$}_{\perp}$ to be the
$\mathcal{O}(\epsilon_{V})$ part of the term in square brackets above, which
is equivalent to the sum of the guiding-center $E\times B$ velocity,
$\mbox{\boldmath${v}$}_{E}$, and the “magnetic flutter” component of the
parallel velocity perpendicular to the background field,
$\mbox{\boldmath${v}$}_{f}$, so that

$\mbox{\boldmath${u}$}_{\perp}(\mbox{\boldmath${R}$},v_{\parallel})\equiv\frac{\mbox{\boldmath${E}$}_{1}\times\mathbf{\hat{b}}}{B_{\parallel}^{*}}+v_{\parallel}\frac{\mbox{\boldmath${B}$}_{1\perp}}{B_{\parallel}^{*}}=\mbox{\boldmath${v}$}_{E}+\mbox{\boldmath${v}$}_{f},$
(2.57)

we can reduce the Hamiltonian to

$\mathcal{H}=\frac{1}{2}mv_{\parallel}^{2}+\mu
B+q\langle\Phi_{1}-\mbox{\boldmath${v}$}_{\perp}\boldsymbol{\cdot}\mbox{\boldmath${A}$}_{1}\rangle-\epsilon_{V}^{2}\frac{m}{2}\left|\mbox{\boldmath${v}$}_{E}+\mbox{\boldmath${v}$}_{f}\right|^{2}-\epsilon_{V}\epsilon_{B}\,m\mbox{\boldmath${v}$}_{d}\boldsymbol{\cdot}(\mbox{\boldmath${v}$}_{E}+\mbox{\boldmath${v}$}_{f}).$
(2.58)

Here,

$\mbox{\boldmath${v}$}_{d}\equiv\frac{mv_{\parallel}^{2}}{qB_{\parallel}^{*}}\mathbf{\hat{b}}\times(\mathbf{\hat{b}}\boldsymbol{\cdot}\nabla\mathbf{\hat{b}})+\frac{\mu}{qB_{\parallel}^{*}}\mathbf{\hat{b}}\times\nabla
B\sim\mathcal{O}(\epsilon_{B})$ (2.59)

is the combined curvature and $\nabla B$ drifts, with
$(\nabla\times\mathbf{\hat{b}})_{\perp}=\mathbf{\hat{b}}\times(\mathbf{\hat{b}}\boldsymbol{\cdot}\nabla\mathbf{\hat{b}})$.
Thus the final form of the gyrocenter one-form is

$\displaystyle\Gamma$
$\displaystyle=q\left(\mbox{\boldmath${A}$}_{0}^{*}+\langle\mbox{\boldmath${A}$}_{1}\rangle\right)\boldsymbol{\cdot}\textnormal{d}\mbox{\boldmath${R}$}+\frac{m\mu}{q}\textnormal{d}\vartheta-\left[\frac{1}{2}mv_{\parallel}^{2}+\mu
B+q\langle\Phi_{1}-\mbox{\boldmath${v}$}_{\perp}\boldsymbol{\cdot}\mbox{\boldmath${A}$}_{1}\rangle\right.$
$\displaystyle\qquad-\left.\epsilon_{V}^{2}\frac{1}{2}m\left|\mbox{\boldmath${v}$}_{E}+\mbox{\boldmath${v}$}_{f}\right|^{2}-\epsilon_{V}\epsilon_{B}\,m\mbox{\boldmath${v}$}_{d}\boldsymbol{\cdot}(\mbox{\boldmath${v}$}_{E}+\mbox{\boldmath${v}$}_{f})\right]\textnormal{d}t.$
(2.60)

We can now see that the first-order correction to the gyrocenter one-form,
$\Gamma_{1}$, has effectively replaced the guiding-center potentials
$\mbox{\boldmath${A}$}_{1}(\mbox{\boldmath${R}$})$ and
$\Phi_{1}(\mbox{\boldmath${R}$})$ that appeared in $\Gamma_{0}$ with
gyroaveraged versions. This also resulted in additional higher-order terms in
the Hamiltonian. While we could (and should) develop the Lie transform to
second order using Eq. 2.33 to obtain FLR corrections to these second-order
terms, and possibly other second-order terms, the guiding-center (long-
wavelength) versions of these terms as they appear in Eq. 2.60 will be
sufficient for our current purposes, since we will only use the first-order
Hamiltonian to compute the equations of motion.

Proper treatment of the second-order $E\times B$ energy term
$-m/2|\mbox{\boldmath${v}$}_{E}|^{2}$ is necessary for deriving an
energetically-consistent gyrokinetic Poisson equation, as we will see in
Section 2.2. We have obtained this term without needing to compute the next-
order Lie transform by making a convenient choice for the reference frame
velocity $\mbox{\boldmath${u}$}_{\perp}$ in Eq. 2.15, effectively guessing an
$\mathcal{O}(\epsilon_{V})$ correction to the velocity ${v}$ that one could
also find from continuing the Lie transform. We can compare the second-order
terms in Eq. 2.60 to Eq. (54) from Brizard & Hahm (2007), which gives the
Hamiltonian that results from computing the Lie transform to second order,
with the second-order terms given in the long-wavelength limit. We see that
indeed we have recovered some of the second-order terms, but missed a term of
the form $\mu|\mbox{\boldmath${B}$}_{\perp 1}|^{2}/(2B)$.

####  2.1.4 Gyrocenter equations of motion

Now that we have the gyrocenter one-form given by Eq. 2.60, we can derive the
gyrocenter Poisson bracket and the gyrocenter equations of motion. At this
point, we will simplify the system by assuming
$\mbox{\boldmath${A}$}_{1}=A_{\parallel}\mathbf{\hat{b}}$, so that

$\mbox{\boldmath${B}$}_{1}=\nabla\times(A_{\parallel}\mathbf{\hat{b}})=-\mathbf{\hat{b}}\times\nabla
A_{\parallel}+A_{\parallel}\nabla\times\mathbf{\hat{b}}.$ (2.61)

This results in the neglect of most compressional fluctuations of the magnetic
field, although even in this form there remains a small compressional
component,
$\mathbf{\hat{b}}\boldsymbol{\cdot}\mbox{\boldmath${B}$}_{1}=A_{\parallel}\mathbf{\hat{b}}\boldsymbol{\cdot}\nabla\times\mathbf{\hat{b}}$,
which may vanish or be finite depending on the particular magnetic geometry.
Note that the second term in Eq. 2.61 is frequently dropped since it is
smaller than the first term by $\mathcal{O}(\epsilon_{B})$, but we will choose
to keep it, in part so that
$\nabla\boldsymbol{\cdot}\mbox{\boldmath${B}$}_{1}=0$ exactly. Future work
will include the full compressional fluctuations $\delta B_{\parallel}$, which
can influence microinstabilities not only at large $\beta\sim 1$ but also when
gradients of $\beta$ are large, particularly in spherical torus machines like
NSTX (Bourdelle et al., 2003; Joiner et al., 2010; Belli & Candy, 2010; Zocco
et al., 2015).

We will also drop some second-order terms in the Hamiltonian, but for now we
will leave the exact form of the Hamiltonian unspecified, since the
Hamiltonian does not affect the form of the Poisson bracket. Thus we will
write the gyrocenter Lagrangian as

$\mathcal{L}=q\left(\mbox{\boldmath${A}$}_{0}^{*}+\langle{A}_{\parallel}\rangle\mathbf{\hat{b}}\right)\boldsymbol{\cdot}\dot{\mbox{\boldmath${R}$}}+\frac{m\mu}{q}\dot{\vartheta}-H=\mbox{\boldmath${\Lambda}$}\boldsymbol{\cdot}\dot{\mbox{\boldmath${Z}$}}-H,$
(2.62)

where ${\Lambda}$ denotes the symplectic part of the Lagrangian.

The phase-space Euler-Lagrange equations are then given by

$\frac{\textnormal{d}}{\textnormal{d}t}\left(\frac{\partial\mathcal{L}}{\partial\dot{Z}^{i}}\right)=\frac{\partial\mathcal{L}}{\partial
Z^{i}},$ (2.63)

which yields

$\dot{\Lambda}_{i}=\frac{\partial\Lambda_{j}}{\partial
Z^{i}}\dot{Z}^{j}-\frac{\partial H}{\partial Z^{i}},$ (2.64)

where $Z^{i}\ (i=1,\dots,6)$ are the phase-space coordinates ${Z}$ not
including time. We can expand the total time derivative $\dot{\Lambda}_{i}$
and rearrange terms to obtain

$\omega_{ij}\dot{Z}^{j}=\frac{\partial H}{\partial
Z^{i}}+\frac{\partial\Lambda_{i}}{\partial t},$ (2.65)

where the Lagrange tensor $\boldsymbol{\omega}$ is defined by

$\omega_{ij}\equiv\frac{\partial\Lambda_{j}}{\partial
Z^{i}}-\frac{\partial\Lambda_{i}}{\partial Z^{j}}.$ (2.66)

Assuming $\det\boldsymbol{\omega}\neq 0$, we can define the Poisson tensor
$\mathbf{\Pi}$ to be the inverse of the Lagrange tensor, (i.e.,
$\Pi^{ik}\omega_{kj}=\delta^{i}_{j}$), so that the Euler-Lagrange equations
from Eq. 2.65 can be inverted to give the equations of motion as

$\dot{Z}^{i}=\Pi^{ij}\left(\frac{\partial H}{\partial
Z^{j}}+\frac{\partial\Lambda_{j}}{\partial t}\right).$ (2.67)

Defining the (non-canonical) Poisson bracket as

$\\{f,g\\}\equiv\frac{\partial f}{\partial Z^{i}}\Pi^{ij}\frac{\partial
g}{\partial Z^{j}},$ (2.68)

and recognizing that $\Pi^{ij}=\\{Z^{i},Z^{j}\\}$, we can also write the
equations of motion as

$\dot{Z}^{i}=\\{Z^{i},H\\}+\\{Z^{i},Z^{j}\\}\frac{\partial\Lambda_{j}}{\partial
t}.$ (2.69)

Inserting the gyrocenter phase-space Lagrangian from Eq. 2.60 into Eq. 2.66,
the non-zero tensor elements are (Cary & Brizard, 2009)

$\displaystyle\omega_{R_{i}v_{\parallel}}=-\omega_{v_{\parallel}R_{i}}=-m\,b_{i}$
(2.70)
$\displaystyle\omega_{R_{i}R_{j}}=-\omega_{R_{j}R_{i}}=q\,\epsilon_{ijk}\bar{B}^{*k}$
(2.71)
$\displaystyle\omega_{\mu\vartheta}=-\omega_{\vartheta\mu}=-\frac{m}{q},$
(2.72)

where $\epsilon_{ijk}$ is the Levi-Civita tensor and $\bar{B}^{*k}$ are the
components of
$\mbox{\boldmath${\bar{B}}$}^{*}\equiv\langle\mbox{\boldmath${B}$}^{*}\rangle=\mbox{\boldmath${B}$}_{0}^{*}+\nabla\times(\langle
A_{\parallel}\rangle\mathbf{\hat{b}})$, so that the full tensor takes the form

$\mbox{\boldmath${\omega}$}=\begin{pmatrix}0&q\bar{B}^{*3}&-q\bar{B}^{*2}&-mb_{1}&0&0\\\
-q\bar{B}^{*3}&0&q\bar{B}^{*1}&-mb_{2}&0&0\\\
q\bar{B}^{*2}&-q\bar{B}^{*1}&0&-mb_{3}&0&0\\\ mb_{1}&mb_{2}&mb_{3}&0&0&0\\\
0&0&0&0&0&-\frac{m}{q}\\\ 0&0&0&0&\frac{m}{q}&0\end{pmatrix}.$ (2.73)

We can then invert this to obtain the Poisson tensor,

$\mathbf{\Pi}=\mbox{\boldmath${\omega}$}^{-1}=\frac{1}{\bar{B}_{\parallel}^{*}}\begin{pmatrix}0&-\frac{1}{q}b_{3}&\frac{1}{q}b_{2}&\frac{1}{m}\bar{B}^{*1}&0&0\\\
\frac{1}{q}b_{3}&0&-\frac{1}{q}b_{1}&\frac{1}{m}\bar{B}^{*2}&0&0\\\
-\frac{1}{q}b_{2}&\frac{1}{q}b_{1}&0&\frac{1}{m}\bar{B}^{*3}&0&0\\\
-\frac{1}{m}\bar{B}^{*1}&-\frac{1}{m}\bar{B}^{*2}&-\frac{1}{m}\bar{B}^{*3}&0&0&0\\\
0&0&0&0&0&\frac{q}{m}\\\ 0&0&0&0&-\frac{q}{m}&0\end{pmatrix},$ (2.74)

with
$\bar{B}_{\parallel}^{*}=\mathbf{\hat{b}}\boldsymbol{\cdot}\mbox{\boldmath${\bar{B}}$}^{*}$.
The Poisson bracket is then given by Eq. 2.68,

$\displaystyle\\{f,g\\}$
$\displaystyle=\frac{\mbox{\boldmath${\bar{B}}$}^{*}}{m\bar{B}_{\parallel}^{*}}\boldsymbol{\cdot}\left(\nabla
f\frac{\partial g}{\partial v_{\parallel}}-\frac{\partial f}{\partial
v_{\parallel}}\nabla
g\right)-\frac{\mathbf{\hat{b}}}{q\bar{B}_{\parallel}^{*}}\boldsymbol{\cdot}\nabla
f\times\nabla g+\frac{q}{m}\left(\frac{\partial
f}{\partial\vartheta}\frac{\partial g}{\partial\mu}-\frac{\partial
f}{\partial\mu}\frac{\partial g}{\partial\vartheta}\right).$ (2.75)

We can also compute the Jacobian of the transformation from particle
coordinates to gyrocenter coordinates
$\mbox{\boldmath${z}$}=(\mbox{\boldmath${x}$},\mbox{\boldmath${v}$})\rightarrow\mbox{\boldmath${Z}$}=(\mbox{\boldmath${R}$},v_{\parallel},\mu,\vartheta)$,

$\mathcal{J}=\frac{1}{m^{3}}\sqrt{\det\boldsymbol{\omega}}=\frac{\bar{B}_{\parallel}^{*}}{m}=\frac{1}{m}\left[B+\left(\frac{m}{q}v_{\parallel}+\langle
A_{\parallel}\rangle\right)\mathbf{\hat{b}}\boldsymbol{\cdot}\nabla\times\mathbf{\hat{b}}\right]=\frac{B}{m}+\mathcal{O}(\epsilon_{B}),$
(2.76)

where we will neglect
$\mathbf{\hat{b}}\boldsymbol{\cdot}\nabla\times\mathbf{\hat{b}}\lesssim\mathcal{O}(\epsilon_{B})$
in the Jacobian. This approximation breaks the exact equivalence of
$\mbox{\boldmath${\omega}$}^{-1}$ and $\mathbf{\Pi}$, but otherwise does not
affect conservation properties. Note that the factor of $1/m^{3}$ comes from
the Jacobian of the transformation from canonical to non-canonical particle
coordinates
$(\mbox{\boldmath${x}$},\mbox{\boldmath${p}$})\rightarrow(\mbox{\boldmath${x}$},\mbox{\boldmath${v}$})$.444In
some texts this $1/m^{3}$ factor does not appear and the Jacobian is given as
$m^{2}B_{\parallel}^{*}$, which is the Jacobian of the transformation
$(\mbox{\boldmath${x}$},\mbox{\boldmath${p}$})\rightarrow\mbox{\boldmath${Z}$}$;
an additional factor of $1/m$ can also appear when the parallel momentum
$p_{\parallel}$ is used as a gyrocenter coordinate instead of $v_{\parallel}$.

Now we can use Eq. 2.69 to obtain the gyrocenter equations of motion,

$\displaystyle\dot{\mbox{\boldmath${R}$}}$
$\displaystyle=\\{\mbox{\boldmath${R}$},H\\}+\frac{\mathbf{\hat{b}}}{q\bar{B}_{\parallel}^{*}}\times\frac{\partial\mbox{\boldmath${\Lambda}$}}{\partial
t}=\frac{\mbox{\boldmath${\bar{B}}$}^{*}}{m\bar{B}_{\parallel}^{*}}\frac{\partial
H}{\partial
v_{\parallel}}+\frac{\mathbf{\hat{b}}}{q\bar{B}_{\parallel}^{*}}\times\nabla
H$ (2.77) $\displaystyle\dot{v}_{\parallel}$
$\displaystyle=\\{v_{\parallel},H\\}-\frac{\mbox{\boldmath${\bar{B}}$}^{*}}{m\bar{B}_{\parallel}^{*}}\boldsymbol{\cdot}\frac{\partial\mbox{\boldmath${\Lambda}$}}{\partial
t}=-\frac{\mbox{\boldmath${\bar{B}}$}^{*}}{m\bar{B}_{\parallel}^{*}}\boldsymbol{\cdot}\nabla
H-\frac{q}{m}\frac{\partial\langle A_{\parallel}\rangle}{\partial t}.$ (2.78)

Finally, note that the zeroth-order (guiding-center) equations of motion are
the same, except with all gyroaverages replaced by evaluation of the quantity
at the guiding-center position.

###  2.2 Gyrokinetic field theory

In the previous section, we derived the phase-space Lagrangian and equations
of motion for a single charged particle in the presence of electromagnetic
fields. Now we describe the collective behavior of a system of many such
particles and the interactions between the particles and the fields.

The system Lagrangian is given by integrating the single-particle Lagrangian
over phase space, weighted by the distribution function $f_{s}$ and summed
over all species $s$, plus an additional field term:

$\mathcal{L}=\sum_{s}\int\mathcal{J}f_{s}\mathcal{L}_{s}\,\textnormal{d}^{6}\mbox{\boldmath${Z}$}+\int\mathcal{L}_{f}\,\textnormal{d}^{3}\mbox{\boldmath${x}$}.$
(2.79)

Here,
$\textnormal{d}^{6}\mbox{\boldmath${Z}$}\equiv\textnormal{d}^{3}\mbox{\boldmath${R}$}\,\textnormal{d}v_{\parallel}\,\textnormal{d}\mu\,\textnormal{d}\vartheta$
is the gyrocenter phase-space volume element with $\mathcal{J}=B/m$ the
Jacobian (dropping the $\mathcal{O}(\epsilon_{B})$ terms in Eq. 2.76),
$\mathcal{L}_{s}$ is the single-particle Lagrangian for species $s$, and
$\mathcal{L}_{f}$ is the field Lagrangian. Note that formally the Jacobian
might be included in the definition of
$\textnormal{d}^{6}\mbox{\boldmath${Z}$}$, but we instead opt to have the
Jacobian appear explicitly in the expressions.

####  2.2.1 The gyrokinetic Vlasov equation

In the absence of sources and collisions (which we address later), the
evolution of the distribution function $f$ is governed by the gyrokinetic
Vlasov equation. This takes the form of Liouville’s equation, which states
that the distribution function is conserved along the nonlinear phase-space
characteristics. This is expressed by

$\frac{\textnormal{d}f(\mbox{\boldmath${Z}$},t)}{\textnormal{d}t}=\frac{\partial
f}{\partial t}+\dot{\mbox{\boldmath${Z}$}}\boldsymbol{\cdot}\frac{\partial
f}{\partial\mbox{\boldmath${Z}$}}=\frac{\partial f}{\partial
t}+\dot{\mbox{\boldmath${R}$}}\boldsymbol{\cdot}\nabla
f+\dot{v}_{\parallel}\frac{\partial f}{\partial v_{\parallel}}=0,$ (2.80)

with the phase-space characteristics given by Eqs. 2.77 and 2.78. From Eq.
2.69, this can also be written in terms of the Poisson bracket as

$\frac{\partial f}{\partial
t}+\\{f,H\\}+\\{f,\mbox{\boldmath${Z}$}\\}\boldsymbol{\cdot}\frac{\partial\mbox{\boldmath${\Lambda}$}}{\partial
t}=\frac{\partial f}{\partial t}+\\{f,H\\}-\frac{q}{m}\frac{\partial\langle
A_{\parallel}\rangle}{\partial t}\frac{\partial f}{\partial v_{\parallel}}=0.$
(2.81)

Together with Liouville’s theorem, which states that phase-space volume is
conserved, as expressed by

$\frac{\partial\mathcal{J}}{\partial
t}+\frac{\partial}{\partial\mbox{\boldmath${Z}$}}\boldsymbol{\cdot}\left(\mathcal{J}\dot{\mbox{\boldmath${Z}$}}\right)=\frac{\partial\mathcal{J}}{\partial
t}+\nabla\boldsymbol{\cdot}\left(\mathcal{J}\dot{\mbox{\boldmath${R}$}}\right)+\frac{\partial}{\partial
v_{\parallel}}\left(\mathcal{J}\dot{v}_{\parallel}\right)=0,$ (2.82)

the gyrokinetic Vlasov equation can also be written in conservative form as

$\frac{\partial(\mathcal{J}f)}{\partial
t}+\frac{\partial}{\partial\mbox{\boldmath${Z}$}}\boldsymbol{\cdot}\left(\mathcal{J}f\dot{\mbox{\boldmath${Z}$}}\right)=\frac{\partial(\mathcal{J}f)}{\partial
t}+\nabla\boldsymbol{\cdot}\left(\mathcal{J}f\dot{\mbox{\boldmath${R}$}}\right)+\frac{\partial}{\partial
v_{\parallel}}\left(\mathcal{J}f\dot{v}_{\parallel}\right)=0.$ (2.83)

####  2.2.2 Variational derivation of the gyrokinetic field equations

We follow Sugama (2000); Scott & Smirnov (2010), to derive the gyrokinetic
field equations. The field equations are derived directly from the Lagrangian
by requiring variations of the action,
$\mathcal{I}=\int\mathcal{L}\,\textnormal{d}t$, to vanish with respect to the
fields $\Phi_{1}$ and $\mbox{\boldmath${A}$}_{1}$. In this way, approximations
and simplifications can be made at the level of the Lagrangian, and then the
resulting field equations will be consistent with those approximations, so
that momentum and energy conservation are preserved.

Thus we must first specify the form of the Lagrangian. We consider three
different cases: (1) keeping second-order terms in the single-particle
Hamiltonian; (2) dropping second-order terms in the Hamiltonian; (3) dropping
first- and second-order terms in the single-particle Lagrangian, resulting in
the guiding-center Lagrangian.

##### Case 1: Single-particle Hamiltonian with second-order terms

In the first case, we take the single-particle Lagrangian to be the gyrocenter
Lagrangian from Eq. 2.62. For the Hamiltonian, we will keep second-order
terms, but we will neglect all second-order terms involving magnetic
fluctuations in Eq. 2.58, so that we are left with

$H=\frac{1}{2}mv_{\parallel}^{2}+\mu
B+q\langle\Phi\rangle-\frac{1}{2}mv_{E}^{2}=\frac{1}{2}mv_{\parallel}^{2}+\mu
B+q\langle\Phi\rangle-\frac{m}{2B^{2}}|\nabla_{\perp}\Phi|^{2}.$ (2.84)

Note that here and after we will drop the subscript on $\Phi_{1}$, since there
is no $\Phi_{0}$ to confuse it with.

The field Lagrangian $\mathcal{L}_{f}$ comes from the standard electrodynamic
field term $(E^{2}-B^{2})/(2\mu_{0})$, but we assume quasineutrality, which
eliminates the electric field term. After neglecting parallel fluctuations of
the magnetic field, we have

$\mathcal{L}_{f}=-\frac{B_{1\perp}^{2}}{2\mu_{0}}\approx-\frac{|\nabla_{\perp}A_{\parallel}|^{2}}{2\mu_{0}}.$
(2.85)

The system Lagrangian for this case is now

$\displaystyle\mathcal{L}$
$\displaystyle=\sum_{s}\int\mathcal{J}f_{s}\left[q_{s}\left(\mbox{\boldmath${A}$}_{0}^{*}+\langle
A_{\parallel}\rangle\mathbf{\hat{b}}\right)\boldsymbol{\cdot}\dot{\mbox{\boldmath${R}$}}+\frac{m_{s}\mu}{q_{s}}\dot{\vartheta}\right.$
$\displaystyle\qquad\left.-\left(\frac{1}{2}m_{s}v_{\parallel}^{2}+\mu
B+q\langle\Phi\rangle-\frac{m_{s}}{2B^{2}}|\nabla_{\perp}\Phi|^{2}\right)\right]\textnormal{d}^{6}\mbox{\boldmath${Z}$}-\int\frac{|\nabla_{\perp}A_{\parallel}|^{2}}{2\mu_{0}}\textnormal{d}^{3}\mbox{\boldmath${x}$}.$
(2.86)

The field equation for the electrostatic potential $\Phi$ is found from the
requirement that variations of the action with respect to $\Phi$ vanish. This
gives the condition $\delta\mathcal{I}/\delta\Phi(\mbox{\boldmath${x}$})=0$,
where the functional derivative is given by

$\displaystyle\frac{\delta\mathcal{I}}{\delta\Phi(\mbox{\boldmath${x}$})}$
$\displaystyle=\sum_{s}\int\left[-\mathcal{J}f_{s}\frac{\partial
H_{s}}{\partial\Phi(\mbox{\boldmath${x}$})}+\nabla\boldsymbol{\cdot}\left(\mathcal{J}f_{s}\frac{\partial
H_{s}}{\partial\nabla\Phi(\mbox{\boldmath${x}$})}\right)\right]\textnormal{d}^{6}\mbox{\boldmath${Z}$}$
$\displaystyle=\sum_{s}\int\left[-q_{s}\mathcal{J}f_{s}\delta(\mbox{\boldmath${x}$}-\mbox{\boldmath${R}$}-\mbox{\boldmath${\rho}$})-\nabla\boldsymbol{\cdot}\left(\mathcal{J}f_{s}\delta(\mbox{\boldmath${x}$}-\mbox{\boldmath${R}$})\frac{m_{s}}{B^{2}}\nabla_{\perp}\Phi\right)\right]\textnormal{d}^{6}\mbox{\boldmath${Z}$}.$
(2.87)

Requirement that this quantity vanish yields an equation for
$\Phi(\mbox{\boldmath${x}$})$ that takes the form of the quasineutrality
condition,

$\sigma_{gy}+\sigma_{pol}=\sigma_{gy}-\nabla\boldsymbol{\cdot}\mbox{\boldmath${P}$}=0,$
(2.88)

where the gyrocenter charge density is

$\sigma_{gy}=\sum_{s}q_{s}\int\langle\mathcal{J}f_{s}\rangle^{\dagger}\textnormal{d}^{3}\mbox{\boldmath${v}$}\equiv\sum_{s}q_{s}\bar{n}_{s},$
(2.89)

with $\textnormal{d}^{3}\mbox{\boldmath${v}$}\equiv
2\pi\textnormal{d}v_{\parallel}\textnormal{d}\mu$. We will continue writing
$\textnormal{d}^{3}\mbox{\boldmath${v}$}$ in integrals defined this way
throughout, even though there are only two evolved velocity dimensions; the
factor of $2\pi$ comes from a trivial integration over the gyroangle, the
third velocity dimension. This factor cancels the factors of $(2\pi)^{-1}$
included in the defintions of the gyroaveraging operations, where integration
over the gyroangle does appear explicitly. Again, we also choose to leave the
phase-space Jacobian out of our definition of
$\textnormal{d}^{3}\mbox{\boldmath${v}$}$ so that it appears explicitly in the
expressions. The notation

$\langle
f\rangle^{\dagger}\equiv\frac{1}{2\pi}\int_{0}^{2\pi}f(\mbox{\boldmath${x}$}-\mbox{\boldmath${\rho}$})\,\textnormal{d}\vartheta$
(2.90)

denotes a gyroaverage taken at constant ${x}$ (as opposed to constant ${R}$).
Note that this operator is the adjoint of the gyroaverage taken at constant
${R}$ defined in Eq. 2.50, _i.e._ it satisfies the property

$\int\langle f\rangle\,g\,\textnormal{d}^{3}\mbox{\boldmath${R}$}=\int
f\langle g\rangle^{\dagger}\,\textnormal{d}^{3}\mbox{\boldmath${x}$}.$ (2.91)

The polarization charge density $\sigma_{pol}$ is given as the divergence of
the polarization vector

$\mbox{\boldmath${P}$}=-\sum_{s}\int\mathcal{J}f_{s}\frac{m_{s}}{B^{2}}\nabla_{\perp}\Phi\,\textnormal{d}^{3}\mbox{\boldmath${v}$}=-\sum_{s}\frac{m_{s}n_{s}}{B^{2}}\nabla_{\perp}\Phi.$
(2.92)

The quasineutrality condition can then be written as the gyrokinetic Poisson
equation,

$-\nabla\boldsymbol{\cdot}\sum_{s}\frac{m_{s}n_{s}}{B^{2}}\nabla_{\perp}\Phi=\sum_{s}q_{s}\bar{n}_{s}.$
(2.93)

Finally, note that the second order $E\times B$ energy term in the Hamiltonian
was required to obtain the polarization charge density. Without it, the
quasineutrality condition would not give us an equation for the potential.

The field equation for the parallel vector potential $A_{\parallel}$ is found
from the requirement that variations of the action with respect to
$A_{\parallel}$ vanish. This gives the condition $\delta\mathcal{I}/\delta
A_{\parallel}(\mbox{\boldmath${x}$})=0$, where the functional derivative is
given by

$\displaystyle\frac{\delta\mathcal{I}}{\delta
A_{\parallel}(\mbox{\boldmath${x}$})}$
$\displaystyle=-\nabla\boldsymbol{\cdot}\frac{\partial\mathcal{L}_{f}}{\partial\nabla
A_{\parallel}(\mbox{\boldmath${x}$})}+\sum_{s}\int\mathcal{J}f_{s}\frac{\partial\mathcal{L}_{s}}{\partial
A_{\parallel}(\mbox{\boldmath${x}$})}\,\textnormal{d}^{6}\mbox{\boldmath${Z}$}$
$\displaystyle=\frac{1}{\mu_{0}}\nabla_{\perp}^{2}A_{\parallel}+\sum_{s}q_{s}\int\mathbf{\hat{b}}\boldsymbol{\cdot}\dot{\mbox{\boldmath${R}$}}\,\mathcal{J}f_{s}\,\delta(\mbox{\boldmath${x}$}-\mbox{\boldmath${R}$}-\mbox{\boldmath${\rho}$})\,\textnormal{d}^{6}\mbox{\boldmath${Z}$}$
$\displaystyle=\frac{1}{\mu_{0}}\nabla_{\perp}^{2}A_{\parallel}+\sum_{s}q_{s}\int\frac{1}{m_{s}}\frac{\partial
H_{s}}{\partial
v_{\parallel}}\mathcal{J}f_{s}\,\delta(\mbox{\boldmath${x}$}-\mbox{\boldmath${R}$}-\mbox{\boldmath${\rho}$})\,\textnormal{d}^{6}\mbox{\boldmath${Z}$}$
$\displaystyle=\frac{1}{\mu_{0}}\nabla_{\perp}^{2}A_{\parallel}+\sum_{s}q_{s}\int
v_{\parallel}\mathcal{J}f_{s}\,\delta(\mbox{\boldmath${x}$}-\mbox{\boldmath${R}$}-\mbox{\boldmath${\rho}$})\,\textnormal{d}^{6}\mbox{\boldmath${Z}$},$
(2.94)

where we have used Eq. 2.77 to substitute for
$\mathbf{\hat{b}}\boldsymbol{\cdot}\dot{\mbox{\boldmath${R}$}}=v_{\parallel}$.
The requirement that this quantity vanish results in the parallel component of
Ampère’s law,

$-\nabla_{\perp}^{2}A_{\parallel}=\mu_{0}\bar{J}_{\parallel},$ (2.95)

with

$\bar{J}_{\parallel}\equiv\sum_{s}q_{s}\int
v_{\parallel}\langle\mathcal{J}f_{s}\rangle^{\dagger}\,\textnormal{d}^{3}\mbox{\boldmath${v}$}=\sum_{s}q_{s}\overline{n_{s}u}_{\parallel
s}$ (2.96)

the gyrocenter parallel current density.

##### Case 2: Single-particle Hamiltonian without second-order terms

In this case, we will drop all second-order terms in the Hamiltonian, so that
we are left with

$H=\frac{1}{2}mv_{\parallel}^{2}+\mu B+q\langle\Phi\rangle.$ (2.97)

As we noted above, the second order $E\times B$ energy term (which we have now
dropped) was needed to obtain the polarization charge density in the
quasineutrality equation. Without the second-order term, we need another way
to obtain the polarization term.

Following Sugama (2000); Scott & Smirnov (2010), we can instead cast the
higher-order $E\times B$ energy term into a field term, _i.e._ as part of the
field Lagrangian $\mathcal{L}_{f}$. To do this, we replace the distribution
function multiplying this term in the system Lagrangian by a time-independent
background distribution function $f_{0}$, giving

$\displaystyle\mathcal{L}$
$\displaystyle=\sum_{s}\int\mathcal{J}f_{s}\left(\mathcal{L}_{s0}+\mathcal{L}_{s1}\right)\textnormal{d}^{6}\mbox{\boldmath${Z}$}+\sum_{s}\mathcal{J}f_{0s}\mathcal{L}_{s2}\,\textnormal{d}^{6}\mbox{\boldmath${Z}$}-\int\frac{|\nabla_{\perp}A_{\parallel}|^{2}}{2\mu_{0}}\textnormal{d}^{3}\mbox{\boldmath${x}$}$
(2.98)
$\displaystyle=\sum_{s}\int\mathcal{J}f_{s}\left[q_{s}\left(\mbox{\boldmath${A}$}_{0}^{*}+\langle
A_{\parallel}\rangle\mathbf{\hat{b}}\right)\boldsymbol{\cdot}\dot{\mbox{\boldmath${R}$}}+\frac{m_{s}\mu}{q_{s}}\dot{\vartheta}-\left(\frac{1}{2}m_{s}v_{\parallel}^{2}+\mu
B+q_{s}\langle\Phi\rangle\right)\right]\textnormal{d}^{6}\mbox{\boldmath${Z}$}$
$\displaystyle\qquad+\sum_{s}\int\mathcal{J}f_{0s}\frac{m_{s}}{2B^{2}}|\nabla_{\perp}\Phi|^{2}\,\textnormal{d}^{6}\mbox{\boldmath${Z}$}-\int\frac{|\nabla_{\perp}A_{\parallel}|^{2}}{2\mu_{0}}\textnormal{d}^{3}\mbox{\boldmath${x}$},$
(2.99)

with the entire second line of Eq. 2.99 now comprising $\mathcal{L}_{f}$.

After following the same steps as in Case 1 to derive the quasineutrality
condition from $\delta\mathcal{I}/\delta\Phi(\mbox{\boldmath${x}$})=0$, we
obtain the same expression for the gyrocenter charge density in Eq. 2.89, but
the polarization vector is modified as

$\mbox{\boldmath${P}$}=-\sum_{s}\int\mathcal{J}f_{0s}\frac{m_{s}}{B^{2}}\nabla_{\perp}\Phi\,\textnormal{d}^{3}\mbox{\boldmath${v}$}=-\sum_{s}\frac{m_{s}n_{0s}}{B^{2}}\nabla_{\perp}\Phi,$
(2.100)

where the density in Eq. 2.92 has been replaced by some time-independent
background density $n_{0}$. This result is sometimes referred to as
_linearized polarization_. From a computational perspective, it is helpful
that in this case the kernel that must be inverted in the gyrokinetic Poisson
equation does not change in time, allowing for parts of the inversion (_e.g._
matrix factorization) to be done only at the beginning of the calculation.
Thus the linearized polarization is commonly used for computational efficiency
(Idomura et al., 2008; Ku et al., 2018a; Shi et al., 2019), even in the
edge/SOL where large density fluctuations could lead to questions of the
validity of replacing the full density with a background density.

It is important to note that using the linearized polarization approximation
requires neglecting the second order $E\times B$ energy term in the
Hamiltonian, and vice versa, in order for the resulting system to be
energetically consistent. This will be shown explicitly in Section 2.2.3,
where we show conservation properties of the system.

Finally, the parallel Ampère equation for $A_{\parallel}$ remains unchanged
from Eq. 2.95; since the Hamiltonian does not depend on $A_{\parallel}$ (in
the symplectic formulation), dropping second-order terms in the Hamiltonian
has no effect on variations of the action with respect to $A_{\parallel}$.

##### Case 3: Guiding-center single-particle Lagrangian

We finally consider the case where we drop first- and second-order terms in
the single-particle Lagrangian, resulting in the guiding-center Lagrangian
(with no gyroaverages). Similar to Case 2, we can cast the first- and second-
order terms into field terms multiplying a background distribution function:

$\displaystyle\mathcal{L}$
$\displaystyle=\sum_{s}\int\mathcal{J}f_{s}\mathcal{L}_{s0}\,\textnormal{d}^{6}\mbox{\boldmath${Z}$}+\sum_{s}\int\mathcal{J}f_{0s}\left(\mathcal{L}_{s1}+\mathcal{L}_{s2}\right)\textnormal{d}^{6}\mbox{\boldmath${Z}$}-\int\frac{|\nabla_{\perp}A_{\parallel}|^{2}}{2\mu_{0}}\textnormal{d}^{3}\mbox{\boldmath${x}$}$
(2.101)
$\displaystyle=\sum_{s}\int\mathcal{J}f_{s}\left[q_{s}\left(\mbox{\boldmath${A}$}_{0}^{*}+A_{\parallel}\mathbf{\hat{b}}\right)\boldsymbol{\cdot}\dot{\mbox{\boldmath${R}$}}+\frac{m_{s}\mu}{q_{s}}\dot{\vartheta}-\left(\frac{1}{2}m_{s}v_{\parallel}^{2}+\mu
B+q_{s}\Phi\right)\right]\textnormal{d}^{6}\mbox{\boldmath${Z}$}$
$\displaystyle\qquad+\sum_{s}\int\mathcal{J}f_{0s}\left(q_{s}[\langle
A_{\parallel}\rangle-
A_{\parallel}]\mathbf{\hat{b}}\boldsymbol{\cdot}\dot{\mbox{\boldmath${R}$}}-q_{s}[\langle\Phi\rangle-\Phi]+\frac{m_{s}}{2B^{2}}|\nabla_{\perp}\Phi|^{2}\right)\textnormal{d}^{6}\mbox{\boldmath${Z}$}$
$\displaystyle\qquad-\int\frac{|\nabla_{\perp}A_{\parallel}|^{2}}{2\mu_{0}}\textnormal{d}^{3}\mbox{\boldmath${x}$}.$
(2.102)

The functional derivative of the action with respect to variations of
$\Phi(\mbox{\boldmath${x}$})$ gives

$\displaystyle\frac{\delta\mathcal{I}}{\delta\Phi(\mbox{\boldmath${x}$})}$
$\displaystyle=\sum_{s}\int\bigg{[}-q_{s}\mathcal{J}f_{s}\delta(\mbox{\boldmath${x}$}-\mbox{\boldmath${R}$})-q_{s}\mathcal{J}f_{0s}[\delta(\mbox{\boldmath${x}$}-\mbox{\boldmath${R}$}-\mbox{\boldmath${\rho}$})-\delta(\mbox{\boldmath${x}$}-\mbox{\boldmath${R}$})]\Bigg{.}$
$\displaystyle\qquad\left.+\nabla\boldsymbol{\cdot}\left(\mathcal{J}f_{0s}\delta(\mbox{\boldmath${x}$}-\mbox{\boldmath${R}$})\frac{m_{s}}{B^{2}}\nabla_{\perp}\Phi\right)\right]\textnormal{d}^{6}\mbox{\boldmath${Z}$}$
$\displaystyle=-\sum_{s}q_{s}\int\mathcal{J}\left(f_{s}+[\langle
f_{0s}\rangle^{\dagger}-f_{0s}]\right)\textnormal{d}^{3}\mbox{\boldmath${v}$}-\nabla\boldsymbol{\cdot}\left(\sum_{s}\frac{m_{s}}{B^{2}}\nabla_{\perp}\Phi\int\mathcal{J}f_{0s}\,\textnormal{d}^{3}\mbox{\boldmath${v}$}\right)$
$\displaystyle=-\sum_{s}q_{s}(n_{s}+\bar{n}_{0s}-n_{0s})-\nabla\boldsymbol{\cdot}\sum_{s}\frac{m_{s}n_{0s}}{B^{2}}\nabla_{\perp}\Phi.$
(2.103)

If we assume that the background density $n_{0}$ varies slowly on the
gyroradius scale (consistent with the ordering $\rho/L_{p}\ll 1$), we can
approximate $\bar{n}_{0s}-n_{0s}\approx 0$ (and to be consistent, we should
also drop the $f_{0}\mathcal{L}_{1}$ field term in Eq. 2.101), so that the
gyrokinetic Poisson equation becomes

$-\nabla\boldsymbol{\cdot}\sum_{s}\frac{m_{s}n_{0s}}{B^{2}}\nabla_{\perp}\Phi=\sum_{s}q_{s}{n}_{s}.$
(2.104)