Nonlinear second order electromagnetic gyrokinetic theory for a tokamak plasma

The steep plasma pressure gradient that forms at the edge of the high confinement, H-mode regime of tokamak operation provides free energy to drive electromagnetic micro-instabilities that are widely believed to influence the transport processes in this so-called pedestal region. This high pressure gradient also provides a high current density (bootstrap current), known to influence ballooning mode stability and to be important for driving kink modes in the ideal magneto-hydrodynamic plasma model (so-called peeling-ballooning modes). Furthermore, efficient, steady state future tokamak power plants must operate with a large bootstrap current in the core and especially concerning spherical tokamaks, confinement will be influenced by electromagnetic turbulence. To accommodate these important situations, conventional electromagnetic gyrokinetic theory is extended to incorporate neoclassical effects in the equilibrium drives, allowing Bϑ∼B0 (B 0 is the confining magnetic field, and B ϑ is its poloidal component). This provides a global gyrokinetic model that self-consistently captures the consequences of large bootstrap current fractions on the equilibrium distribution functions.


Introduction
The tokamak is the most advanced fusion energy system, confining the fuel in plasma that is formed by a toroidal magnetic field geometry. The fusion performance is strongly influenced by the effectiveness of the confinement system, which is usually degraded by turbulence. This turbulence arises as a consequence of a zoo of drift instabilities with short wavelengths in directions perpendicular to the magnetic field, typically at * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. or below the ion Larmor radius; the wavelength along the magnetic field lines is comparable to the system size. Under certain conditions, including increasing the heating power through a well-defined threshold, a tokamak plasma bifurcates to a state of high confinement, called the H-mode [1]. This is a consequence of the suppression of turbulence in the outer few centimetres of the plasma providing an insulating transport barrier and consequent high pressure gradient. This transport barrier also called the pedestal region, results in high core pressure and is attractive for high fusion energy output. The H-mode therefore forms the basis for the baseline operational regime for ITER [2].
The pressure at the core of an H-mode plasma is highly sensitive to the properties of the pedestal: a wide pedestal with steep pressure gradient results in high core pressure and therefore high fusion power in a tokamak like ITER. For this reason, optimising the fusion power in ITER is likely to require the optimisation of the pedestal. There are two key pieces of physics that influence these pedestal properties [3].
First, intermediate toroidal mode number (n ∼ 10) ideal magneto-hydrodynamic (MHD) instabilities that tap the combined effects of the plasma pressure gradient and current density, called peeling-ballooning modes, trigger plasma eruptions called edge localised modes (or ELMs) [4][5][6]. Each ELM results in a crash in the pedestal pressure gradient, and with many ELMs per second they therefore have a strong impact on the time-averaged pedestal properties. Peeling-ballooning modes are well-understood, with quantitative predictions available for their trigger conditions. There are three main pedestal parameters that influence this trigger in a tokamak plasma of given shape-the pressure gradient, the current density and the pedestal width. The onset of peeling-ballooning modes provides our first condition, which constrains the pedestal structure.
The second piece of physics that has an impact on the pedestal structure is the local drift instabilities that influence turbulence and transport there [7][8][9][10]. These typically have a much shorter wavelength than the peeling-ballooning modes, and are often characterised by wavenumbers perpendicular to the magnetic field k ⊥ ≳ 1/ρ i , where ρ i is the ion Larmor radius. As mentioned above, the transport barrier forms as one class of turbulence is suppressed. This enables the pressure gradient to rise until another class of micro-instabilities is triggered, driving residual turbulence and related transport that controls the pressure gradient in the pedestal. A number of features differentiate the pedestal from the core plasma in a conventional tokamak, resulting in a different nature to the turbulence: (a) the high pressure gradient compared to the magnetic energy density drives stronger fluctuations in the magnetic field; (b) strong electric field shear influences the micro-stability and turbulence, and (c) the high pressure gradient drives strong neoclassical flows and bootstrap current. Understanding this turbulence and related transport provides another constraint on the pedestal parameters which, when combined with the peeling-ballooning mode constraint, provides a minimal model for the pedestal structure [3].
Spherical tokamak plasmas achieve much higher ratios of thermal to magnetic energy (β) than conventional aspect ratio tokamaks. This then makes the core plasma vulnerable to electromagnetic instabilities, such as the kinetic ballooning mode and microtearing modes. This is especially true in reactor scenarios, such as those envisaged for STEP [11,12], where a high β is beneficial for commercial viability. Benefits include achieving a high fraction of bootstrap current and creating a spherical tokamak core plasma turbulence regime with several similarities to the pedestal plasma of conventional tokamaks, described above.
Consider the bootstrap current, which is driven in a tokamak plasma by the density and temperature gradients that the magnetic field structure supports; it flows parallel to the magnetic field [13]. In the pedestal, the high pressure gradient drives a high bootstrap current in a relatively narrow region, resulting in steep current density gradients. Such high current density gradients are well known to drive long wavelength kink-modes and, as mentioned above, play an important role in the intermediate n peeling-ballooning modes. However, because this drive falls as 1/n, it is not usually thought to be important for short wavelength drift modes at ion Larmor radius length scales. It is thus ordered out of the conventional gyrokinetic theory, which provides a reasonable approach in the core of a conventional tokamak. However, the steep current density gradient in the pedestal challenges this assumption. Recent work [14] has extended the ELITE ideal MHD code [6,15] to quantify the relative sizes of the kink and curvature drive terms, and the stabilising effect of field line bending. It is found that for a model equilibrium, characteristic of mid-sized tokamaks, the kink drive dominates over the curvature term out to n ∼ 40, which corresponds to k ⊥ ρ i ∼ 1 (see figure 4.19 of [14]). While this kink drive is predominantly external (associated with a vacuum rational surface), this provides significant evidence for the importance of the equilibrium bootstrap current density for pedestal plasma stability, which then raises the question of how to effectively incorporate it into microinstability calculations.
Gyrokinetic theory provides the basis for most microstability analysis in tokamak plasmas. This includes full-f kinetic simulations (e.g. [16][17][18][19]) where the equilibrium and fluctuation scales are merged, and a δf approach (global, e.g. [20] or [21,22] based on [23]; local theory, e.g. [24][25][26][27] or [28] for drift cyclotron kinetics, and simulations [29][30][31][32], and hybrid, e.g. [33][34][35], where the radial profile variation is introduced perturbatively). One of the complications associated with the global full-f techniques arises from the gyrokinetic quasineutrality equation that is not derived to sufficiently high order to determine the electrostatic potential fully consistent with the fluctuating distribution function (this problem is due to a coupling between the radial electric field [36] and the toroidal plasma rotation and was addressed in [37] in the limit of B ϑ /B 0 ≪ 1). In the δf formalism, one separates the distribution function into two parts: (a) an 'equilibrium' distribution function, F j (j is the particle species index), which evolves on slow, transport timescales and has length scales of order the system size, and (b) a fluctuating part, δf j , which has short wavelength across magnetic field lines (comparable to the Larmor radius) and long wavelength along them (comparable to the system size).
One impact of the bootstrap current can be understood by incorporating it into the MHD equilibrium solution, which can then be analysed for micro-instabilities using a gyrokinetic code. However, such analyses will typically adopt a Maxwellian equilibrium distribution function, which is then not consistent with the bootstrap current used to evaluate the MHD equilibrium. The neoclassical effects on the equilibrium distribution function can be introduced by expanding the drift kinetic equation that describes it in powers of ρ ϑj /L, where ρ ϑj is the Larmor radius of species j evaluated using only the poloidal component of the magnetic field (i.e. the poloidal Larmor radius), and L characterises the length scale of equilibrium variations. The leading order provides a Maxwellian solution without flows, while the next order provides the physics of the neoclassical flows, including the bootstrap and Pfirsch-Schluter currents. To ensure consistent ordering, one would also need to retain features of the second order correction that depend on the gyrophase angle. This equilibrium solution is then used to derive the gyrokinetic equation for the fluctuating quantities in response to that equilibrium. Solving the drift kinetic equation for the equilibrium plasma eliminates the finite orbit width effects and is equivalent to working in terms of the magnetic moment extended to higher order in the Larmor radius relative to the system size. The corresponding gyrokinetic equation is described by the generalised gyrokinetics of [38] in the absence of nonlinear effects; by [39] when the nonlinear effects are retained; or by [37,[40][41][42] that focuses on purely electrostatic fluctuations to address a theory for toroidal angular momentum transport. This generalisation implicitly [39,43] or explicitly [37,40] implies B ϑ /B 0 ≪ 1. This assumption significantly simplifies the gyrokinetic derivation and allows one to solve for the gyrophase-independent fluctuating distribution function only. B ϑ /B 0 ∼ ε/q, where ε is the inverse aspect ratio of the flux surface and q is the corresponding safety factor. B ϑ /B 0 ≪ 1 is typically a good approximation for conventional tokamaks, including for tokamak pedestals (provided ε ≪ 1, and q is large at the edge), but becomes questionable for spherical tokamak equilibria. To avoid this limitation and to allow general magnetic field configurations, in this paper we relax the B ϑ /B 0 ∼ ε/q ≪ 1 assumption, which (as we shall see in the sections below) ultimately results in solving a system of gyrokinetic equations, including for the gyrophase angle dependent piece of the fluctuating distribution function, and also requires one to retain the finite Larmor radius effects in the equilibrium distribution function.
Previous work also considered the impact of higher order equilibrium effects on the gyrokinetic physics. The emphasis of [44], for example, was to make a connection between gyrokinetic theory and ideal MHD, and therefore employed a number of approximations relevant for that comparison, which we will not make here. An electromagnetic theory was developed in [45], retaining only the parallel component of the vector potential; they adopted a model for an equilibrium with an electron flow to provide the current density, but dropped the finite orbit width effects, which we keep in order to capture the bootstrap current. Other works [46][47][48] provide a phase space Lagrangian to address the second order corrections in the Larmor radius over the system size; they restricted consideration to the electrostatic limit, which therefore provides a benchmark in the Lagrangian formulation of the electromagnetic theory presented here.
The reminder of the paper is organised as follows. In section 2 we describe the magnetic topology and coordinate system. In section 3 we introduce the general kinetic equation in guiding centre coordinates. This is followed by the calculation in section 4 of the equilibrium distribution function, expanding it in powers of ρ ϑ /L (ρ ∼ ρ ϑ , where ρ is the particle Larmor radius and ρ ϑ is the poloidal Larmor radius) 4 up to second order to ensure consistent ordering. The latter also ensures that the equilibrium solution 4 The particle species indices are omitted for simplicity and to be introduced only when it is necessary to distinguish ions and electrons. is consistent with the MHD equilibrium in situations when the bootstrap current is large. In section 5 we proceed to the calculation of the fluctuating distribution function and derive a set of gyrokinetic equations to ensure, first, consistent ordering, second, incorporating current density gradient effects, third, B ϑ /B 0 ∼ 1. The set of equations presented in section 5 is integro-differential. In section 5 we also consider the limit of B ϑ /B 0 ≪ 1 to recover the Frieman-Chen nonlinear gyrokinetics [39]. In section 6 we discuss limitations associated with adopting a local (flux tube) approximation. This is followed by conclusions and a summary in section 7.

Magnetic topology and coordinates
We consider an axisymmetric tokamak plasma with an equilibrium magnetic field given by where I = RB φ depends on the poloidal current, R is the tokamak major radius, and B φ is the toroidal component of the magnetic field. Here, ψ is the poloidal flux function introduced to label nested magnetic flux surfaces, and φ denotes the toroidal angle. We choose the orthogonal coordinate system with the basis formed by the unit vectors (b, e 2 , e 3 with B ϑ being the poloidal component of the magnetic field and B 0 = |B 0 |. Note, b = e 2 × e 3 , e 2 · de 3 /dt = −e 3 · de 2 /dt and e 2 · de 2 /dt = 0 and so on, where d/dt denotes the time derivative following the actual particle orbit.
To distinguish between equilibrium and perturbed quantities, we write where Y denotes an arbitrary field. Here Y 0 describes the equilibrium piece, and Y is the fluctuating piece with A measure of the size of the perturbation is then For the distribution function we write where Considering rapid spatial variations of the field/distribution function perturbations, in gyrokinetics all equilibrium quantities evolve on transport time scales, which we neglect, i.e.
We note if ∇ acts on unperturbed quantities and if ∇ acts on perturbed quantities with ∆ = ρ/L ≪ 1, in accordance with the conventional gyrokinetic ordering. Here ∥/ ⊥ denotes the vector component along/across the equilibrium magnetic field lines. The second small parameter we define is where Θ = B ϑ /B 0 . We allow Θ ∼ 1 and hence The E × B and magnetic drifts are estimated as respectively, where V t is the particle thermal velocity. We write the velocity as where is parallel and s = s (e 2 cos α + e 3 sin α) ≡ sŝ (9) is perpendicular to the equilibrium magnetic field. Here α denotes the gyrophase angle. The Larmor radius vector is then where ω c = eZB 0 /m is the cyclotron frequency, eZ and m are the particle charge and mass, respectively; ρ = s/ω c . The position of the guiding centre is then defined as where x is the actual position of the particle. We define the magnetic moment (per unit mass) as µ = s 2 /2B 0 . The total energy (per unit mass) is defined as U = µB 0 + u 2 /2 + eZΦ 0 /m, where Φ 0 is the equilibrium electrostatic potential. The equation of motion for these particles then reads where eZV × B 0 = −mω c sρ, B is the perturbed magnetic field and E = E 0 + E is the electric field.

General kinetic equation in {t, X, U, µ, α} space
The general kinetic equation for the particle distribution function, f, in {t, X, U, µ, α} space reads The collisional effects are described by the collision operator, C (unspecified, its form depends on the particular problem of interest). The left hand side contains the Vlasov operator written in terms of the guiding centre coordinates, X. From the equation of motion, equation (12), it can be shown that and The derivation is similar to that presented in [49]. Here we have introduced B ∥ = B · b. Equations (14), (15) and (17) can be further simplified to and Here j ∥ is the component of the equilibrium current density along the equilibrium magnetic field lines, and the drift velocity is defined as C X , C µ and C α are defined as and Here Equation (13) wherē and

Equilibrium
In accordance with equation (4), we write the particle distribution function in a form: where F describes the equilibrium piece (in gyrokinetic ordering) and δf is the fluctuating piece. F satisfies the equation: Note, to include turbulence effects on the neoclassical transport, one should add ⟨Lδf ⟩ X ⊥ to the left hand side of equation (27), . ⟨Ỹ⟩ X ⊥ = 0 and X ⊥ corresponds to rapid spatial variations [39]. However, it can be demonstrated that (1) ⟨Lδf ⟩ X ⊥ is only required to calculate the α independent part of the O(δ 2 F 0 ) equilibrium solution, which does not contribute to the fluctuation equation and (2) To solve equation (27) for F, we expand it in powers of δ (e.g. see [50][51][52]), allowing equation (6). This perturbative treatment of equation (27) provides the following straightforward result: whereC C µ and C X are given by equations (21) and (22), respectively. Note that the collision frequency is assumed to be small compared to the cyclotron frequency to ensure that the leading order equilibrium distribution function is independent of α at fixed X. We also anticipate that the difference between the collision operator and its gyro-averaged form can be neglected 5 . The first term in equation (28) is , and therefore only the first term is to be retained in F 2 in the limit of B ϑ ≪ B 0 . An equivalent approach to obtain equation (28) at B ϑ ≪ B 0 would be to consider the magnetic moment extended to higher order in ∆, where ⟨μ⟩ is independent of the gyrophase angle (e.g. see [37,38,40,49,53] etc), and then to expand F around µ: which is equivalent to equation (28) up to O(δ∆) with B ϑ ≪ B 0 and the following replacement: F (X, U, µ) = F 0 (X, U) + F 1 (X, U, µ). Note, ⟨μ⟩ is typically found to ensure thatμ is an adiabatic invariant: [54]. As we shall see below, at B ϑ ≪ B 0 this extension to the equilibrium distribution function results in an extended source term, e.g. equation (35) of [38] for the linear or equation (48) of [39] for the nonlinear case 6 . Based on the O(∆ 0 δ 1 Vt L F 0 ) equilibrium equation [50], averaged over α at fixed X, we can estimate the size of the deviation from Maxwellian, F M , employed in the conventional gyrokinetic theory for F 0 , by considering the edge pedestal region of an H-mode plasma, which we know has a pressure gradient close to the ideal MHD ballooning boundary, α MHD ∼ (Rq 2 /B 2 0 )dp/dr ∼ 1, where q is the safety factor, p is the plasma pressure, and r is the varying tokamak minor radius. Interestingly, the plasma scenarios being considered for STEP are in a similar regime [11]. This then yields where ε is the tokamak inverse aspect ratio and β ϑ is the ratio of the thermal energy density of the plasma to the magnetic energy density in the poloidal component of the magnetic field. This is a factor order B 0 /B ϑ larger than the finite Larmor radius correction to the Maxwellian. There are some important consequences from incorporating the correction F 1 into the equilibrium distribution function: (a) it introduces an equilibrium flow and current density, which are strong in regions of high density/temperature gradient; (b) it introduces some additional anisotropy into the equilibrium distribution function, and (c) the equilibrium distribution function is not constant on a flux surface. While F 1 captures the neoclassical physics, including the bootstrap current associated with the finite banana orbit widths of trapped particles, the gyro-angle dependent part of the second order correction to the equilibrium distribution, F 2 , must be retained as well to ensure the ordering is fully consistent.

Perturbations: gyrokinetic theory
The plasma response to the fluctuating electromagnetic field is described by δf: 6 A similar approach can be applied to reconstruct the second term on the right hand side of equation (28) to allow B ϑ ∼ B 0 . It would require an iterative procedure to obtain the X dependence of ρ in equation (11), i.e.
Blue/green rectangles correspond to conventional (e.g. [27])/ generalised [38,39] gyrokinetics. Blue/yellow dashed rectangles indicate partial contributions to the O(∆ 1 Vt provided F satisfies equation (27). To capture turbulence effects in the equilibrium solution, one should add ⟨Lδf ⟩ X ⊥ to the left hand side of equation (27), and hence to ensure consistency, −⟨Lδf ⟩ X ⊥ should contribute to the left hand side of equation (30). However, it can be demonstrated that ⟨⟨Lδf ⟩ X α ⟩ X ⊥ = 0, and hence ⟨Lδf ⟩ X ⊥ will be omitted. The where we learn that δf 0 is independent of α at fixed X. δf 0 is to be determined from the solubility condition on the O(∆ 1 Vt The Equation (32), together with the solubility condition from equation (33), provides δf 1 . Here L 1 and L 2 are given by equation (25) with ( E, B) being replaced with (E ′ , B ′ ) and (E ′ ′ , B ′ ′ ), respectively, and δf 2 = O(∆ 3 F 0 ). Note that equation (33) contains certain nonlinear O(∆ 1 Vt L F 0 ) terms that contribute to equation (32), and L 2 δf 0 + L 1 δf 1 contains nonlinear O(∆ 2 Vt L F 0 ) corrections that contribute to equation (33) (see figure 1). Splitting the fluctuating distribution function into an adiabatic response and a resonant piece: where and Here we have defined g Note that the distinction between g and h is that g is driven by the fields and higher order nonlinear terms associated with Y ′ and g.
WritingLδf 0 + L 1 F explicitly in equation (32) (see appendix A) and expanding g in powers of δ, we obtain (33) provide an equation for g 0 and the α dependent part of h 0 + g 1 : Here −iω is to be understood as ∂/∂t when acting on fluctuating quantities. Note that, in the absence of nonlinear terms, the gyro-averaging of equation (38) provides the conventional linear gyrokinetic equation (e.g. [27] for electrostatic fluctuations). Gyro-averaging equation (38) over α at fixed X eliminates the first two terms on the left hand side of equation (38) and provides a solubility condition to be solved for g 0 = g 0 (t, X, U, µ). Substituting the obtained g 0 into equation (38) and integrating equation (38) over α at fixed X provides the α dependent part of h 0 + g 1 ≡ H = H(α) +H, denoted byH. The α independent part of h 0 + g 1 , i.e.H =H (t, X, U, µ), is to be determined from a combination of the O ∆ 1 δ 1 Vt L F 0 and O ∆ 2 Vt L F 0 equations, in accordance with the perturbation theory and figure 1.
Note, in equation (39) we have omitted the O(∆ 2 δ 1 F 0 ) term in δf 1 as it is not relevant for the subsequent derivations. Averaging equation (39) over α at fixed X annihilates the first term on the left-hand side of equation (39) and provides the first part of the O ∆ 2 Vt L F 0 equation to be solved for the α independent part of h 0 + g 1 .

The
In this subsection, we focus on the third line in figure 1. WritingLδf 1 + L 1 δf 0 + L 2 F explicitly in equation (33) and neglecting the O ∆ 2 δ Vt L F 0 corrections (or higher), we find where h 0 is related to δf 1 via equation (36). In equation (40) we have taken into account that L 3 F + L 2 δf 0 + L 1 δf 1 contains the O ∆ 2 Vt L F 0 corrections that have to be retained in the O ∆ 2 Vt L F 0 equation as well (see appendix B). Equation (40) contains the O ∆ 1 δ 1 Vt L F 0 and O ∆ 2 Vt L F 0 corrections and provides the second part of the O ∆ 2 Vt L F 0 equation to be solved for the gyro-angle independent part of h 0 + g 1 . Here we have defined∇ as ∇ acting on equilibrium quantities, e.g.

Source term
Averaging equations (38) and (39) and equation (40) over α at fixed X, we obtain: with the right hand side (source term) given by (see appendix C) and Here S 10 L contains the O ∆ 1 δ 0 Vt L F 0 terms which correspond to the source term of the conventional linear gyrokinetic equation. S 11 L contains the O ∆ 1 δ 1 Vt L F 0 corrections associated with the first term on the right hand side of equation (28) and provides the generalised gyrokinetics of [38]. S 20 L is driven by the O ∆ 2 F 0 contribution to the equilibrium solution, i.e. the second term on the right hand side of equation (28). T 20 L is the O ∆ 2 Vt L F 0 contribution provided by equation (40). S 1 NL , S 2 NL and S 3 NL are the nonlinear corrections from equation (40), in accordance with figure 1. The parallel streaming operator in guiding centre coordinates, (b · ∇) X , is related to b · ∇ via equation (D.3) (see appendix D for a detailed derivation).
In the limit of low collisions, a combination of equations (42)- (44) can be further simplified. Indeed, assuming that the particle collision frequency is much smaller compared to the free streaming along the equilibrium magnetic field lines/characteristic drift frequency in the equilibrium solution and differentiating the gyro-averaged, O(∆ 0 δ 1 Vt L F 0 ) equation for F 1 [50] with respect to µ and U, one can find and hence Equation (37) is the O ∆ 0 Vt L F 0 equation, where we learn that g 0 is gyrophase independent. Proceeding to O ∆ 1 Vt L F 0 in equation (41), which is equivalent to equation (38) averaged over α at fixed X, we obtain Equation (51) reproduces the leading order, conventional gyrokinetics and can be solved for g 0 = g 0 (t, X, U, µ). Proceeding to O ∆ 2 Vt L F 0 in equation (41), we write Note, S 1 NL on the right hand side of equation (52) should exclude the nonlinear O(∆ 1 Vt L F 0 ) correction already employed in equation (51). Substituting g 0 obtained as a solution of equation (51) into equation (38) and integrating equation (38) over α at fixed X providesH with the constant of integration,H =H (t, X, U, µ), being determined from equation (52) as follows: The B ϑ ∼ B 0 ordering is equivalent to ∆ ∼ δ and thus ultimately requires solving for the gyro-angle dependent piece of the fluctuating distribution function when retaining neoclassical physics in the equilibrium plasma. The latter is relevant, for example, to the edge pedestal region of a tokamak where pressure gradients are strong or to a spherical tokamak core plasma. We note that in the absence of collisions, the integration over α is straightforward for the right hand side of equation (38) and allows one to express ⟨LH⟩ X α on the right hand side of equation (53) in terms of g 0 and Y ′ analytically. One can then combine equations (51) and (53) to obtain a single equation for g 0 +H which is gyrophase independent. The latter is addressed in section 5.4 for electrostatic and in section 5.5 for electromagnetic fluctuations.

The nonlinear electrostatic case
In the absence of collisions and the fluctuating vector potential, integrating equation (38) over α at fixed X provides which can then be substituted into equation (53). A sum of equations (51) and (53) then provides where we have defined G = g 0 +H. As discussed in section 4, an alternative approach to obtain equation (55) is to retain the O(∆µ) and O(∆ρ) corrections in the magnetic moment and the guiding centre variable, respectively. For example, Imposing a straight, constant magnetic field approximation and taking into account equations (35) and (36), it can be demonstrated that equations (55) and (56) contain those equations provided by the recursive approach of [55] for nonlinear electrostatic fluctuations. Note that the equation presented in appendix A of [55] for the gyrophase angle independent piece of the fluctuating distribution function can be further simplified to the form of equation (55) when the second order, gyrophase angle dependent correction, i.e.g i2 of [55], is substituted explicitly. As discussed in [55], the main difference between equation (19) of [55] or equation (55) of this paper and the electrostatic limit of the nonlinear Frieman-Chen gyrokinetics [39] is that [39] contains the leading order, E × B nonlinearity only. As discussed in section 4, the latter is due to the assumptions associated with the form of the equilibrium distribution function employed in [39].

The nonlinear electromagnetic case
Retaining the fluctuating vector potential, we obtaiñ and therefore for a sum of equations (51) and (53) in the absence of collisions. Here, we have defined In the limit of low collisions in the equilibrium solution, the F 0 contribution to the last line of equation (57) can be rearranged using equation (49) to provide Note, the g 0 contribution to the last line of equation (57) can be rewritten in a similar way if one split equation (57) into a system of gyrokinetic equations for g 0 = O(∆F 0 ) (equation (51)) andH = O(∆ 2 F 0 ). Therefore, equation (57) has the structure of equation (55) with Φ ′ and Φ ′ ′ being replaced with the generalised potential, χ ′ and χ ′ ′ , respectively.
In appendix E, we construct the particle and energy balance equations associated with equation (57) in the absence of collisions to illustrate how equation (57), retaining the equilibrium current density gradient physics, extends the conventional 'gyrokinetic-fluid' model associated with the conventional gyrokinetic equation of [27].

The B ϑ ≪ B 0 limit
For benchmarking purposes, in this section we discuss how certain existing extended gyrokinetic models, including the Frieman-Chen nonlinear gyrokinetics of [39], are contained in equation (57). In the limit of B ϑ ≪ B 0 , ∆ ≪ δ and therefore only the first term on the right hand side of equation (28) is to be retained. While the magnetic moment variation is captured in the equilibrium solution, no second order finite Larmor radius effects are required when defining the guiding centre coordinate, and hence X = x − ρ(X, V). Equations (51) and (53) (or equation (57) in the absence of collisions) then reduce to (58) (2) [37,40,41] in the absence of the fluctuating vector potential. Equation (58) can be further simplified to where Note, the first line of equation (44) reduces from O(∆δ Vt L F 0 ) to O(∆ 2 Vt L F 0 ) in the limit of low collisions (see equation (49)) and thus is not captured in [38,39].
The plasma turbulent transport is driven by microinstabilities that are typically described by the gyrokinetic simulations, and can be sensitive to the collisional effects. The latter requires the gyro-averaged collision operator with the retained finite Larmor orbit width effects. As discussed above, for the second order gyrokinetics, an additional complication arises from the fact that X = x − ρ (X + ρ, V). In this paper, we do not discuss the gyrokinetic collision operator. There are a number of works and reviews [56] (and references therein) that address this, e.g. [57][58][59] for multiple ion species linearised model collisions, for the gyrokinetic linearised Landau/Fokker-Planck collision operator [60,61] etc.

Maxwell's field equations in gyrokinetics
To ensure a self-consistent description of the nonlinear gyrokinetic system, equations (38), (51) and (53) [39] can be significantly simplified if (1) b · ∇ is rewritten in terms of the guiding centre variable, as shown in equation (58), and (2) the O(∆µ) term in the magnetic moment is written explicitly via VD. Note, Cµ/B 0 is equivalent to µ 1 − ⟨µ 1 ⟩ in equation (23) of [38] and equation (22) of [39]. typically reduces to the quasi-neutrality condition) consistent with equation (57) reads where ε 0 is the permittivity of free space, and Z denotes the phase space integral [55]. j is the particle species index. For a system of equations (38), (51) and (53), we must separate the first and second order contributions and hence equation (61) reduces to coupled to equations (51) and (53), respectively. In a similar way, Ampère's law consistent with equation (57) can be written as where µ 0 is the magnetic permeability of free space. The parallel and perpendicular components of equation (63) then provide A ′ ∥ , A ′ ′ ∥ and B ′ ∥ , B ′ ′ ∥ , respectively. In the absence of the B ∥ effects that are generally weak for gyrokinetic simulations [20], equation (63) reduces to For a system of equations (38), (51) and (53), equation (64) then provides similar to equation (62). Equations (61) and (63) are general and match the Maxwell system of [62] for truncated gyrokinetics and the Poisson equation in the recursive approach of [55] etc. When the flux tube approximation is employed, equations (61) and (63) reduce to a system of local gyrokinetic Maxwell equations, e.g. equations (10)-(12) of [63] (in the absence of neoclassical physics). In the following section, we address limitations of the local, flux tube approximation illustrated by MHD simulation results.

Limitations of the ballooning transform for tokamak pedestal
Since the perturbed fields are allowed to vary significantly over the Larmor radius length scale, equations (38), (51) and (53) (or equations (55) and (57) for collisionless plasmas) form a system of integro-differential equations. Conventionally, such an integro-differential form can be reduced to a differential system by employing an eikonal approximation for the fast spatial variation of fluctuations (e.g. see [27] or [64] for the extended Fourier-ballooning transform that allows one to reconstruct the radial mode structure). A local form of equation (58), for example, which is valid in a low B ϑ limit, is provided by equation (50) of [39] or by equation (41) of [38] for a linearised model. Employing the eikonal approximation for the higher order gyrokinetic theory is associated with certain complications. For example, the eikonal function, S, defined such that b · ∇S(x) = 0 is typically expanded around the guiding centre coordinate up to O(∆S), i.e. S(x) = S(X) + ρ · ∇S, which does not provide a sufficient accuracy in the limit of ∆ ∼ δ. Furthermore, while the local form allows internal modes, it is insufficient to capture modes that couple to a resonant surface outside the plasma ('external'). Indeed, let us illustrate the structure of the kink mode in the pedestal by considering ideal MHD. One then finds that instabilities exist that involve a coupling between the curvature drive (pressure gradient) and the kink drive (current gradient). These are called peeling-ballooning modes. The most robust kink drive occurs when the instability can tap into the rational surface just outside the plasma. This situation is shown in figure 2 (left) for a n = 20 peeling-ballooning mode derived using ELITE [65], which shows that when the pressure gradient peaks near the plasma edge, the eigenfunction has a large amplitude there. The plasma radial displacement is decomposed into poloidal Fourier harmonics, labelled by poloidal mode number, m, and the radial variation of each is plotted in figure 2 (right). Note that there are a significant number of harmonics, each peaking at their rational surface, where m = nq. They are all in phase, resulting in a poloidal mode structure that peaks on the outboard (low field) side of the tokamak plasma. Furthermore, each of the harmonics resonant inside the plasma has qualitatively the same shape, and only differs from its neighbours by the relative amplitude. This (ballooning) symmetry is exploited in 'local theories' (e.g. [66] in ideal MHD or [25] in gyrokinetics), which to leading order describes the radial shape of the Fourier amplitudes, and then the higher order provides their relative size. However, the Fourier harmonics, which are resonant outside the plasma (i.e. those increasing towards the plasma edge) do not have the same shape as the harmonics resonant inside the plasma. Therefore in this situation, the ballooning symmetry breaks  down and local theory must be replaced with a general set of equations (38), (51) and (53) (or equation (55)/(57) for collisionless plasmas) coupled to Poisson's equation or plasma quasi-neutrality and Ampère's law for perturbed fields to close the system. Such a gyrokinetic-Maxwell system can be further simplified, e.g. employing a full spectral decomposition of the fluctuating quantities [65,67].
We consider a special case in figure 3, where we show the effect of moving the peak in the normalised pressure gradient, α MHD , further into the plasma. The eigenfunction now peaks away from the plasma edge, close to where α MHD is a maximum, and we see that the largest amplitude Fourier harmonics now do have qualitatively the same shape. Thus, in this situation we do expect to be able to describe the dominant instability characteristics with a local theory. Although the coupling to the Fourier harmonics resonant in the vacuum is much reduced in this case, the kink drive is still significant -looking at the contributions to the perturbed energy, about 60% arises from the curvature drive compared to 40% from the kink.

Conclusions
The principal motivation for the calculations presented in this paper is to ensure that the effects associated with sharp pressure gradients and the consequential high bootstrap and Pfirsch-Schluter currents are fully captured in the nonlinear electromagnetic gyrokinetic theory, while allowing arbitrary magnetic field configurations and finite orbit width effects and ensuring consistent ordering. Bootstrap current can readily be incorporated into the ideal MHD Grad-Shafranov solution for the equilibrium, which can then be analysed for microinstabilities using standard gyrokinetic tools. However, such gyrokinetic analysis usually explore fluctuations by expanding about an equilibrium distribution function that is Maxwellian, which is then not consistent with a large bootstrap current (and hence the physics of kink modes).
In ideal MHD, the kink term represents the instability drive from current density gradients, and therefore originates from the contribution to force balance j ∥ b × B ′ , where j ∥ b is the equilibrium current density which flows parallel to the magnetic field. In the gyrokinetic formalism, this requires the treatment of the term V × B ′ · ∇ V F which contributes to the force balance (i.e. the first velocity moment of the kinetic equation) from the part of the equilibrium distribution function, F, which carries flow. While the flow is carried out by the first order correction in a δ expansion of the equation, describing the equilibrium distribution function, this equation must be expanded to second order in δ to be self-consistent, and capture the kink drive. To see this, consider where we have employed the local approximation for illustration. It is immediately apparent that an isotropic distribution function F = F(U) does not contribute to this term, and therefore one must include physics in the equilibrium beyond the Maxwellian, as captured by the first order equilibrium solution, F 1 . The first term gyro-averages to zero when acting on F 1 , which is independent of α and so will not contribute to the current. Thus, only the α-dependent parts of F, embodied in the second order correction, F 2 , can contribute to this drive term. This reinforces the importance of capturing the physics of F 2 . Some previous works (e.g. [37][38][39]43]) extend the gyrokinetic equation to incorporate these effects; however, implicitly or explicitly requiring an extra expansion parameter, B ϑ /B 0 ≪ 1, which is typically a good approximation for conventional tokamaks but must be relaxed for a spherical tokamak plasma. The O(∆ 2 Vt L F 0 ) terms whose estimates contain B ϑ /B 0 are also anticipated to play a role for a full treatment of the current density gradient physics in situations when the equilibrium profile gradients are steep (specifically in transport barrier regions). In this paper, we have extended a nonlinear electromagnetic gyrokinetic formalism to a set of equations valid for any tokamak geometry and plasma beta, allowing B ϑ ∼ B 0 . This gyrokinetic system is provided by a set of equations (38), (51) and (53) with no restrictions on the plasma collisionality. For collisionless plasmas, this reduces to equation (55) for nonlinear electrostatic and equation (57) for electromagnetic fluctuations. We have demonstrated that to allow B ϑ ∼ B 0 , a gyrokinetic formalism correct to third order in ρ/L must be derived which ultimately requires solving for the gyroangle dependent part of the fluctuating distribution function (or alternatively retaining the O(∆µ) and O(∆ρ) corrections in the magnetic moment and the guiding centre variable). Note that current gyrokinetic codes are typically valid up to the first order in ρ/L.
To explore internal modes when the fraction of bootstrap current is large, one can exploit a local version of equations (38), (51) and (53). However, in making this simplification, we restrict the applicability of our formalism. For example, such a local approach requires equilibrium quantities to vary slowly from one rational surface to the next, and that may not be compatible with the strong gradients required to drive the kink dynamics. The classical local approximation (e.g. [26]) is commonly applied to reduce the global integro-differential problem to a differential system. The extended local/hybrid models (e.g. [34,64]) allow one to retain slow radial profile variations of the density/temperature gradients. Both incorporate the Fourier-ballooning transform, whose symmetry breaks down at the boundary that separates the plasma from the vacuum. Thus, the local version of gyrokinetics is only suitable to capture modes that are contained inside the plasma (i.e. internal modes). To capture external modes, a global version must be employed. Global effects are expected to be important for capturing all of the key dynamics associated with kink mode physics, so it is necessary to retain the equilibrium profile effects. Note that, while the system of global equations derived in section 5 is accurate, it remains computationally expensive, especially when collisional effects are retained. As part of future work, we will address possible ways for its reduction, while keeping the formulation global and capturing the plasma-vacuum interface.
The approach presented here is iterative. Therefore, despite consistent ordering, it is not generally obvious that the gyrokinetic equation obtained by the iterative method can be rewritten in a phase space conserving form or a form that exactly conserves an energy-like quantity. This problem is well illustrated in [55]. In contrast, it might be necessary to retain certain higher order corrections to ensure that the equations have an exact energy-like invariant [46,55,68] which generally violates consistent ordering in the gyrokinetic-Maxwell system. In section 5.4 we briefly discussed how equation (55) (coupled to equation (61)) for electrostatic fluctuations is related to the nonlinear electrostatic gyrokinetic equation of [55] for straight field lines, which is energy and phase space conserving. We leave an investigation of the conservative properties (and the higher order gyrokinetic field theory in particular 8 ) of equation (55) for electrostatic fluctuations in general geometry and equation (57) for electromagnetic fluctuations for future work.
We started our motivation by exploring kink physics in the pedestal, showing that the kink drive can be significant up to large/intermediate toroidal mode numbers in ideal MHD calculations (e.g. see figure 4.19 of [14]). It is therefore expected to be important for electromagnetic drift instabilities, and the associated turbulence, as found in the pedestal or in reactorgrade spherical tokamak plasmas (e.g. STEP [11]). In tokamak plasmas, gyrokinetic theory is employed to describe these instabilities. Equations (38), (51) and (53) extend conventional gyrokinetics, capturing the kink/peeling (current density gradient) physics for any equilibrium magnetic field configuration. While equations (38), (51) and (53) (or equation (57) for collisionless plasmas) can readily be incorporated into the existing gyrokinetic community codes (e.g. ORB5 [21,22]), the problem with the boundary conditions for external modes remains. Indeed, coupling to the fluctuating magnetic field in the vacuum is important for those MHD instabilities, and understanding how gyrokinetic plasma instabilities couple to the vacuum will need further work and is important for a full understanding of pedestal physics.

Data availability statement
No new data were created or analysed in this study.

Appendix A. Derivation ofLδf
where a is the acceleration provided by the equilibrium component of the Lorentz force and ∇ V is the velocity space gradient at fixed x,Lδf 0 + L 1 F can be reduced tō and the O ∆δ 2 Vt L F 0 corrections or higher have been neglected. Note, here u = ub(x). Equation (A.1) is to be substituted into equation (32). In the limit of low collisions, equation (A.1) can be further simplified. One can write provided the perturbative treatment of collisions is allowed in the equilibrium solution 9 . Equilibrium quantities on the right Appendix B. The O ∆ 2 Vt L F 0 corrections from higher order gyrokinetics that must be added to the O ∆ 2 Vt L F 0 equation.

Appendix C. Source term
To obtain equations (42)-(44), we find it convenient to do the following splitting based on the O(δ∆F 0 ) and O(∆ 2 F 0 ) contributions in equation (28):

Appendix E. Particle and energy balance
To simplify the analysis, in this section we neglect the B ∥ effects that are typically anticipated to be weak at low beta values, characteristic of conventional tokamaks, and in certain finite beta cases [20]. The particle and energy balance equations are obtained by integrating equation (57) over velocity space with appropriate scalar weights [68]. Based on equations (35) and (36), we anticipate a balance of (a) turbulent-driven contributions associated with fluctuating quantities, g 0 and g 1 + h 0 , (b) neoclassical contributions provided by F 1 and the gyro-phase angle dependent part of F 2 and (c) adiabatic contributions associated with the extended adiabatic pieces of equations (35) and (36), as well as the nonlinear contributions provided by the nonlinear drive on the right hand side of equation (57). For the particle balance, we obtain ∂ ∂t n (0) ad,j + n (1) ad,j + n (2) ad,j + n ′ j + n ′ ′ Here the adiabatic contributions are defined as ad,j = −n eqm ad,j∥ = Φ ′ I B 0 dn eqm dψ with n eqm and T j being the equilibrium plasma density and temperature of species j. The neoclassical contribution is provided by where p j is the equilibrium plasma pressure of species j. n ′ j + n ′ ′ j are the turbulent-driven contributions to density, associated with the O(∆F 0 ) and O(∆ 2 F 0 ) parts of G, respectively. Γ ′ j + Γ ′ ′ j is defined similar to equations (76), (77) and (80) of [68] and includes for each particle species. Here the parallel drift velocity is defined as u D = µ 0 s 2 /(2B 0 ω c )j ∥ b [49]. Note, u D · ∇H = O(β∆ 3 Vt L F 0 ) and thus the corresponding contribution has been neglected in Γ ′ ′ L . B ′ ⊥ is the perturbed magnetic field component across the equilibrium magnetic field lines. The energy balance equation is obtained by multiplying equation (57) by mV 2 /2 and integrating over velocity space, Here Γ ′ j and Γ ′ ′ j are defined as in equations (E.2) and (E.3). Note, in equation (E.4) we have taken into account that U is based on the total equilibrium electrostatic potential, Φ 0 , in contrast to [49], where U is based on a total Φ = Φ 0 + Φ. ad,j T j .
Note, 5/2 in equation (E.4) includes the heat flux associated with the neoclassical contribution, (3/2)n eqm T j u neo j∥ , and the work associated with the pressure gradient, n eqm T j u neo j∥ . The energy densities, Q ′ j and Q ′ ′ j , are provided by g 0 and H respectively integrated over velocity space with m j V 2 /2. Equations (E.1) and (E.4) obtained by direct integration of equation (57) over velocity space provide the anticipated balance of adiabatic, neoclassical and turbulent-driven contributions. The B ′ ⊥ · ∇ ⊥ (n eqm u neo j∥ ), omitted in conventional (Maxwellian based) gyrokinetics, is typically referred to as the kink drive term in the fluid-kinetic plasma description.