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

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.

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

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/m² in some present devices and expected to be $\sim 1$ GW/m² in ITER (Loarte et al., 2007). The maximum heat load for present materials with active cooling is typically $10$ MW/m² normal to the surface in steady state and $20$ MW/m² 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).

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

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

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,

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

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

As a result, the parallel force balance equation becomes

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)

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

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

The electric field is given by

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

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

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

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

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

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

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

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.

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

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

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,

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

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)

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,

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

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)

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,

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

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

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.

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

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.²²2A 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

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

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

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

the Lagrangian can be transformed as

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

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

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

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

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,

where we have defined

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

Here,

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

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

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 by³³3The 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.

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

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

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

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.