Quantifying the role of higher order neoclassical corrections to gyrokinetics in tokamak plasmas

We implement the higher order gyrokinetic theory developed in Dudkovskaia et al (2023 Plasma Phys. Control. Fusion 65 045010), reduced to the limit of Bϑ/B0≪1 , where B 0 is the tokamak equilibrium magnetic field, and B ϑ is its poloidal component, in the local gyrokinetic turbulence code, GS2. The principal motivation for this extension is to quantify the importance of neoclassical flows in electromagnetic gyrokinetics, with a particular interest in sharp pressure gradient regions where the bootstrap current becomes dominant. To incorporate neoclassical equilibrium physics, GS2 is coupled to NEO, a multi-species drift kinetic solver. It is found that the regions where microinstabilities are most likely to be influenced by neoclassical equilibrium effects are in a pedestal plasma and a spherical tokamak core plasma.


Introduction
Gyrokinetic theory is typically employed to describe turbulent fluctuations in tokamak plasmas. Yet, simulation codes are based on a variety of alternative formulations. Most generally, there are full-f formulations [1,2] where the equilibrium and fluctuation scales are merged. More common are δf-type formulations [3][4][5][6] where one separates the total distribution function into an equilibrium part, F, that evolves on slow transport timescales and has length scales of order the system * 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. size, and a fluctuating part, δf, which has a short perpendicular wavelength and δf ≪ F. More recently, hybrid δf schemes have been proposed [7,8] that treat the radial profile variation perturbatively.
Conventional local gyrokinetics [3] typically adopt a Maxwellian equilibrium, which does not include neoclassical flow and current density. There are a number of extended models that capture deviations from Maxwellian [4,5,[9][10][11], deriving the gyrokinetic equation in response to this extended non-Maxwellian equilibrium. However, the majority of these models are valid in certain simplified cases: these typically are (1) B ϑ /B 0 ≪ 1, where B 0 is the equilibrium magnetic field and B ϑ is its poloidal component [4,5,10], (2) rare collisions compared to particle parallel streaming along the equilibrium magnetic field lines and/or the characteristic drift frequency, and (3) electrostatic fluctuations. When applied to conventional tokamaks, B ϑ /B 0 ∼ ε/q ≪ 1 is a relatively reasonable approximation: (1) the inverse aspect ratio ε ≪ 1 in the core plasma and (2) the safety factor, q, is large at the plasma edge. However, B ϑ /B 0 ≪ 1 can be violated in spherical tokamaks (STs). Note that, as discussed in [12], the approach of [4,5,10] is based on a difference between the two expansion parameters: ∆ = ρ/L and δ = ρ ϑ /L in the limit of B ϑ /B 0 ≪ 1, where ρ is the particle Larmor radius, ρ ϑ ∼ ρB 0 /B ϑ is the particle poloidal Larmor radius, and L is the system size. References [4,5,10] expand the equilibrium distribution function in powers of δ (which is also equivalent to retaining the O(∆µ) corrections in the magnetic moment, µ) up to O(∆δF 0 ), neglecting the O(∆ 2 F 0 ) corrections, where F 0 is Maxwellian. The B ϑ /B 0 ≪ 1 assumption is relaxed in [12].
While there are a number of higher order gyrokinetic theories, only a few have been implemented in codes. For example, [13][14][15] extend GS2 (a nonlinear local gyrokinetic code designed to study low frequency drift wave turbulence in magnetised plasmas) [16][17][18] to incorporate higher order gyrokinetics based on the diamagnetic/neoclassical corrections to Maxwellian equilibria. Hornsby et al [19] provides a similar implementation of neoclassical effects but in GKW (a nonlinear gyrokinetic flux tube code). All of them employ the B ϑ /B 0 ≪ 1 approximation and are electrostatic. Pusztai et al [20] extends electrostatic low flow gyrokinetics of [10] to a simplified, low beta electromagnetic case, and implements this in GS2 to study high mode number kink modes.
In the present paper, we also exploit the limit of B ϑ /B 0 ≪ 1. This assumption is convenient for numerical implementation and testing and is sufficient to quantify the impact of higher order gyrokinetic corrections associated with neoclassical equilibrium current density effects. We implement the theory developed in [12] (reduced to the limit of small B ϑ for simplicity) in GS2. In contrast to the works listed above, (1) our implementation is electromagnetic, (2) it allows for a finite/high beta plasma and (3) it also allows for finite collisions in the equilibrium solution. To ensure consistency, we also extend the GS2 Maxwell's field equations accordingly.
The paper is structured as follows: in the following section we introduce the tokamak geometry and coordinate system we employ. In section 3 we describe the calculation of the equilibrium distribution function, retaining the higher order, neoclassical corrections. The corresponding generalised gyrokinetic-Maxwell system is presented in section 4. To quantify the importance of these neoclassical corrections, we adopt a local (flux tube) approximation in section 4.1. The impact of the neoclassical corrections is quantified in section 5, and the obtained results are discussed in section 6.

Magnetic geometry and coordinate system
We employ an axisymmetric tokamak geometry with an equilibrium magnetic field written as: where φ is the toroidal angle and ψ is the poloidal flux coordinate. Constant ψ surfaces are assumed to form a set of nested toroidal magnetic flux surfaces, within which the equilibrium magnetic field lines lie. The function I = I(ψ), which is constant on these flux surfaces, is related to the toroidal component of the magnetic field, I = RB φ , where R is the major radius coordinate measuring the distance from the axis of symmetry.
To describe the distribution of ions and electrons in a magnetised plasma, we introduce a velocity coordinate system resolved parallel and perpendicular to the equilibrium magnetic field: with and s = s (e 2 cos α + e 3 sin α) ≡ sŝ.
Here α is the gyrophase angle, and b, e 2 and e 3 are the unit vectors defined as 3 and B 0 = |B 0 |. The gyro-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 . A convenient representation of the velocity coordinates employs α together with the magnetic moment per unit mass, µ, and the total particle energy per unit mass, U: µ = s 2 /2B 0 and U = µB 0 + u 2 /2 + eZΦ 0 /m ≡ V 2 /2 + eZΦ 0 /m, where Φ 0 is the total equilibrium electrostatic potential and V is the particle speed.
We work with spatial coordinates which describe the position of a particle's guiding centre. This guiding centre coordinate, X, is related to the actual position of the particle, x, as it gyrates around the magnetic field lines as

Neoclassical equilibrium solution
To distinguish equilibrium and fluctuating quantities, we separate the particle distribution function, f, into an equilibrium piece, F, and fluctuating piece, δf: 3 Here δf /F ∼ ∆ ≪ 1, where ∆ = ρ/L with the gyro-radius, ρ, and L characterising the system size. F here is the solution of the conventional drift kinetic equation [21,22] that is generally obtained by expanding the kinetic equation in powers of ∆ and subsequent gyro-averaging in the absence of the finite Larmor orbit width effects, i.e. ⟨χ (x, t)⟩ The leading order equilibrium distribution function, F 0 , is provided by a Maxwellian solution that does not carry flows, while the first order solution, F 1 , contains the physics of the neoclassical flows. In particular, F 1 is the solution of a constraint equation for the second order equilibrium distribution function, F 2 : where V D is the drift velocity and C on the right hand side of equation (8) is the gyro-averaged collision operator. Note that the spatial differential operators are defined at fixed U, µ and α, unless stated otherwise. If written in terms of V, µ and α, equation (8) becomes A detailed recursive derivation of the drift kinetic equation can be found in [21], and its numerical solution with the Fokker-Planck collisions consistent with the Poisson equation is described in detail in [22][23][24]. To ensure ordering consistency, we must also retain certain terms in the gyro-angle dependent piece of the second order equilibrium distribution function, F 2 , that can be expressed in terms of F 0 and F 1 (see equation (28) of [12]). As discussed in [12], retaining the gyrophase dependent piece of F 2 in the limit of B ϑ /B 0 ≪ 1 is equivalent to expanding F around µ: where the correctionμ = O(∆µ) (e.g. see [25] or equations (23) and (27) of [4]); this is the approach adopted in [4,26,27] etc. In the present paper, we employ two models for F 1 . The first model is employed for comparative purposes: it adopts a large aspect ratio tokamak approximation and provides F 1 as the solution of the constraint equation based on the momentum-conserving, pitch angle scattering collision operator (see appendix A). In the second model, F 1 is provided by NEO (a multi-species drift kinetic solver). Specifically, it is found as the solution of the first order drift kinetic -Poisson system (i.e. equation (9) for each particle species, coupled by the Poisson equation) with the full linearised Fokker-Planck collisions [22][23][24]. The latter allows consideration of a generalised, finite aspect ratio tokamak geometry.

Generalised gyrokinetic equation for fluctuations
δf provides the fluctuating response to the fluctuating electromagnetic fields. We separate δf into an extended adiabatic piece and a resonant piece to find in the absence of nonlinear effects 4 . (11) and (12) is provided by F 0 + F 1 , as described in the previous section. The gyro-phase angle dependent piece of F 2 is expressed in terms of F 0 and F 1 and incorporated in equation (12). The generalised perturbed potential is defined as (12) can be found in [12]. Equation (12) captures the neoclassical equilibrium physics while ensuring consistent ordering, provided F 1 = O(δF 0 ), where δ = ρ ϑ /L. Indeed, in addition to conventional gyrokinetics [3], the right hand side of equation (12) retains the O(∆δ Vt L F 0 ) corrections associated with the neoclassical effects, i.e. neoclassical flows and bootstrap currents. Note that the focus of this paper is to quantify the impact of these neoclassical corrections, and therefore the O(∆ 2 Vt L F 0 ) terms have been neglected in equation (12), which is reasonable, provided B ϑ /B 0 ≪ 1. Generally, all the second order terms must be retained for consistency -these are provided by equation (57) of [12]. Collisions are described by the operator C on the left hand side of equation (12).
In the limit of rare collisions in the equilibrium solution, equation (12) can be further simplified. Indeed, differentiating equation (8) with respect to µ and U, to leading order in collisions, the last two lines of equation (12) can be rewritten as: providing the O(∆ 2 Vt L F 0 ) contribution, and therefore equation (12) reduces to: 4 At O(∆δ Vt L F 0 ), the nonlinear effects are provided by 1 B0 ⟨∇χ ′ × b⟩ X α · ∂g ∂X added to the right hand side of equation (12).
(b · ∇) X is the parallel streaming operator in guiding centre coordinates. In equations (12) and (14), it is related to b · ∇ via [12] ( is the component of the perturbed vector potential perpendicular to B 0 . Equation (14) matches the linearised Frieman-Chen gyrokinetics of [5] 5 , as well as the 'generalised' gyrokinetics of [4]. In the absence of electromagnetic effects, equation (14) also matches the electrostatic 'low flow' gyrokinetics employed in [10,27,28]. It is therefore apparent that [4,5] etc are only consistent, provided (1) ρ ≪ ρ ϑ (which is equivalent to assuming B ϑ /B 0 ≪ 1) and (2) collisions are rare (the particle collision frequency is much smaller than the particle parallel streaming along the equilibrium magnetic field lines). Note, equation (14) can be further reduced to the standard gyrokinetic equation: where g is related to H via Conventional gyrokinetics based on a Maxwellian equilibrium does not allow for the neoclassical, parallel current density effects to be incorporated in the solution, and is therefore inconsistent with the ideal magneto-hydrodynamic (MHD) model. Indeed, by constructing the force balance equation from Maxwellian based gyrokinetics, one would recover (1) the electric component of the Lorentz force associated with the time variation of the perturbed magnetic field and (2) the pressure gradient related piece in the force balance equation (see appendix B). However, the magnetic component of the Lorentz force and hence the parallel current density gradient effects require the physics of F 1 (F 1 carries neoclassical equilibrium flows) and F 2 (the gyro-angle dependent piece of F 2 is required to ensure consistent ordering and, being expressed in terms of F 1 , also contributes to flows in equations (12) and (14)). A simplified case is considered in appendix B to illustrate how the equilibrium parallel current density gradient physics is retained in ideal MHD associated with equation (14).

Local formulation
In this section, we consider a local theory (suitable for implementation in GS2), representing the flute-like nature of fluctuating quantities. A local version of the gyrokinetic equation is provided by the eikonal approximation [9,21,[29][30][31] employed for fluctuations. This simplifies the integrodifferential form of equations (12) and (14), reducing it to a local differential equation. We employ the Fourier-ballooning transform for fluctuating fields separating the fast spatial variation and the slowly varying amplitude,Φ, and ∇S. Note, Φ ′ here can be replaced by any fluctuating field. Here n is the toroidal mode number and S denotes the eikonal function, S(x, p) = φ − q(ϑ − p) with ϑ being the poloidal angle, defined such that b · ∇S(x) = 0. Expanding S around the guiding centre coordinate S(x) = S(X) + ρ · ∇S and noting that the perpendicular wave number Here ξ is defined as (17) provides the gyro-angle dependence of fluctuating fields at fixed position of the guiding centre. g in equations (12) and (14) is independent of α when expressed in guiding centre coordinates. Therefore, adopting an eikonal form for g, i.e.
g =ˆRg (t, x, p) e inS(x) dp, we must ensure that the slow variation withing is of the form and therefore g =ˆRĝ (t, X, p) e inS(X) dp.
Substituting equations (17) and (20) in equation (12), we obtain Here we have defined the drift frequency as ω D = V D · k ⊥ and . J k is the Bessel function of the first kind, of order k. The collisional term on the left hand side of equation (21) depends on a specific form of the collision operator and therefore has been left unspecified. V D · ∇ ⊥ in equation (21) is to be replaced with V D · i k ⊥ when acting on fluctuating quantities.

Maxwell's field equations in local gyrokinetics
The system of equation (21) for each particle species is closed via the plasma quasi-neutrality condition: and Ampère's law: and Here λ = 2µ/V 2 is the pitch angle and u = σV(1 − λB 0 ) 1/2 with σ = u/|u|; n eqm,j and T j are the equilibrium density and temperature of species j. ρ j is the Larmor radius of species j. µ 0 is the magnetic permeability of free space. According to equation (24), B ′ ∥ ∼ β∆B 0 , where β is the ratio of plasma pressure to magnetic field pressure, and therefore B ′ ∥ is typically anticipated to be more important at higher beta. Based on equation (11), equations (22)-(24) extend a conventional set of the time-independent local gyrokinetic Maxwell equations (e.g. equations (10)-(12) of [32]). For the purposes of this paper, equations (22)- (24) are sufficient, provided B ϑ /B 0 ≪ 1. Generally, when all the second order corrections are retained (see equation (57) of [12]), equations (22)-(24) must retain higher order terms in ∆, including corrections associated with ∇ acting on equilibrium quantities. For example, limitations associated with the conventional gyrokinetic quasi-neutrality equation are discussed in [33] for equilibria that deviate from the Maxwellian.
Equations (21)-(24) provide a relatively compact system, valid in the limit of B ϑ /B 0 ≪ 1, that can be readily incorporated into existing local gyrokinetic codes. Note that typically current gyrokinetic codes are accurate up to O(∆ Vt L F 0 ). To quantify the impact of these neoclassical, higher order corrections, we have extended GS2, incorporating the effects of neoclassical tokamak plasma equilibria in accordance with equations (21)- (24). The extension is compatible with the NEO code for the equilibrium distribution function (an Eulerian code that solves the drift kinetic -Poisson system) [22][23][24]. This therefore allows for a general, D-shaped tokamak geometry, complete with experimental plasma profiles. Furthermore, equation (21) removes the low collisionality limitation in the equilibrium solution. Note that in the limit of rare collisions in the equilibrium solution and in the absence of the fluctuating magnetic field, equation (21) reduces to the low flow gyrokinetics of [10,27]. In the following section, we quantify the importance of the neoclassical, higher order corrections of equations (21)-(24).

Numerical assessment of neoclassical corrections
We extend GS2 to capture the neoclassical physics associated with F 1 in equations (21)- (24). Two models are employed to calculate F 1 . In the first (analytic) model, F 1 is based on a momentum-conserving, pitch angle scattering collision operator in a large aspect ratio tokamak (see appendix A). In this analytic model, collisions are assumed to be much smaller than the particle free streaming along the equilibrium magnetic field lines/drift frequencies and hence can be treated perturbatively. The second model is more general: F 1 is provided by NEO (to be referred to as NEO GS2) with no restrictions on the particle collision frequency, apart from the fact that the collision frequency should be much smaller than the cyclotron frequency to ensure F 0 is gyro-phase independent. In NEO, the drift kinetic equation for the first order non-adiabatic distribution function is solved by expanding it in a Fourier-Legendre series in pitch angle space and a Fourier-Chebyshev [22,23]/Laguerre series [24] in energy space. The non-adiabatic piece of the neoclassical distribution function, F 1 , is then where P a is a Legendre polynomial of degree a ⩾ 0 and L k(a)+1/2 b is the generalised Laguerre polynomial of degree b ⩾ 0. The index k (a) can, in principle, be arbitrary and allows one to test different kernel functions in the energy expansion. In this paper, k (0) = 0 and k (a) = 1 at a > 0. ξ (λ) = u/V with ξ 2 = 1 − λB 0 , ϑ is the poloidal angle, and V t is the particle thermal speed. The first order electrostatic potential associated with F 1 is Φ 1 0 = O(∆Φ 0 ) and eZΦ 0 /T ∼ 1 with T being the leading order particle equilibrium temperature 6 . NEO then finds the generalised Fourier coefficients c ab and the corresponding electrostatic potential. This solution is then provided as an input in GS2 and transformed back into real space onto the GS2 velocity grid in accordance with equations (25) and (26). This 'real space' F 1 then contributes to a system of equations (21)-(24), incorporating neoclassical equilibrium physics. Note that the velocity space transforms in equation (25) are performed at fixed ψ, which therefore allows one to calculate the radial derivatives of F 1 on the right hand side of equation (21) in inverse, Fourier-Legendre/Laguerre space and then perform the transformation into 'real space', i.e. (ψ, ϑ, λ, V). 6 Note To ensure that the converted F 1 satisfies the neoclassical equilibrium of equation (9), we (1) benchmark F 1 against its analytic form obtained for a large aspect ratio, low collisional plasma (as shown in appendix A) and (2) compare the orbitaveraged moments of this 'real space' distribution function against the NEO results. As an example, in figures 1(a) and (b) we show F 1 plotted as a function of λ for different collision operators. 'Connor', 'HS0' and 'FP' correspond to the 'real space' F 1 calculated with the Connor (i.e. the pitch angle scattering operator of appendix A that conserves momentum and number of particles) [34], zeroth-order Hirshman-Sigmar [35] and full linearised Fokker-Planck [24] collision operators. 'Analytic' represents the analytic F 1 calculated in appendix A. As anticipated, the analytic F 1 is in agreement with the 'Connor' F 1 away from the trapped-passing boundary. In the vicinity of the trapped-passing boundary, gradients in pitch angle become steep, and hence collisional dissipation becomes comparable to parallel streaming. This, in turn, creates a thin dissipation layer that surrounds the trapped-passing boundary (e.g. see section 5.2 of [36]), invalidating the perturbative treatment of collisions in the 'analytic' F 1 and explaining the difference between the analytic solution and 'Connor' F 1 in figure 1(b) near λ = 0.8. Since the analytic model F 1 also requires the large aspect ratio tokamak limit, below it will be used only in certain examples for comparative purposes. Note that the zeroth-order Hirshman-Sigmar collision operator contains slowing down effects in addition to pitch angle scattering, providing results in close agreement with the full linearised Fokker-Planck collisional model (figures 1(a) and (b)). This deceleration effect is important to accurately calculate the bootstrap current.
To illustrate the second benchmark test, in figures 1(c) and (d) we plot the NEO and NEO GS2 orbit-averaged parallel current density and its radial derivative as functions of plasma collisionality. The parallel current density is constructed as where only the σ dependent part of the equilibrium distribution function carries parallel flows. Orbit-averaging is defined as an average over ϑ at fixed ψ for passing particles and an integral between bounce points with a sum over σ for trapped particles. The agreement between NEO and NEO GS2 is remarkably good in both cases.
We note that there is no equilibrium electric field in conventional GS2 (i.e. when the conventional, Maxwellian based gyrokinetic-Maxwell system is being solved), i.e. dΦ 0 0 /dψ = 0, with the equilibrium plasma quasi-neutrality being satisfied automatically by the Maxwellian. Inclusion of the higher order corrections in the equilibrium solution makes the equilibrium distribution function not constant on the flux surface. The neoclassical, first order electrostatic potential, Φ 1 0 , is then  Example of the neoclassical, first order electrostatic potential, Φ 1 0 , (green solid line) and its second radial derivative, ∂ 2 Φ 1 0 /∂ε 2 , divided by the radial electric field, − ∂Φ 1 0 /∂ε , (purple dash-dot line) plotted as a function of poloidal angle, ϑ, at the minor radius r = 0.5a (a is the value of r at the last closed flux surface). The inverse aspect ratio ε = 0.2, elongation κ = 1 and triangularity △ = 0; ν i = 10 −3 (Fokker-Planck collisions) ρ i /a = 10 −2 and the equilibrium density and temperature gradient length scales are RL −1 n = 2.2 and RL −1 T = 6.9, respectively. Here Φ 1 0 is normalised to Te/e, where Te is the leading order equilibrium electron temperature.
calculated to ensure that the equilibrium plasma is quasineutral to O(∆ Vt L F 0 ): j eZ j´F1,j dV = 0 (e.g. see figure 2 for Φ 1 0 and its derivative calculated by NEO and figure C1 of appendix C for how it is transformed onto the GS2 grid).
Note that in addition to the limitations listed above, the analytic equilibrium solution of appendix A does not allow for the poloidal angle variation in the equilibrium electrostatic potential, imposing Φ 0 = Φ 0 (ψ). Below, we quantify the impact of these higher order, neoclassical corrections on the linear gyrokinetic stability, based on three different equilibrium models.

Cyclone base case
The Cyclone Base Case (CBC) provides a set of local equilibrium parameters from an ITER relevant DIII-D H-mode discharge [37,38] at minor radius r = 0.5a, where a is the r value at the last closed flux surface. It is known to be unstable to kinetic-ballooning modes (KBMs) [39] at finite beta, as well as to ion temperature gradient (ITG) and trapped electron modes (TEM) in the electrostatic limit [38]. The equilibrium parameters are as follows: s-alpha large aspect ratio tokamak geometry, the safety factor q = 1, magnetic shear s = 0.79, inverse aspect ratio at the surface of interest ε = 0.2, k y ρ i = 0.3 (GS2 notation) 7 with the ratio of the ion and electron temperatures, T i /T e = 1, and their mass ratio, m e /m i = 2.74 × 10 −4 . The second term on the right hand side of equation (21) also requires the radial derivative of F 1 . This is calculated based on equilibrium parameters at r/a = (0.49, 0.5, 0.51).
The neoclassical terms of equations (21)-(24) are anticipated to become more significant at higher ρ and/or at larger equilibrium density or temperature gradients. These conditions are relevant to a plasma, which is characterised by high fractions of bootstrap currents, such as the tokamak pedestal. To probe the impact of the neoclassical corrections on the local linear gyrokinetic stability, below we perform a series of artificial parameter scans around the experimental equilibria, and the s-alpha approximation exploited in this subsection allows one to easily switch between the equilibrium parameters. Note that the s-alpha approximation is to be replaced by a Miller model in the following subsection. In figure 3 we provide an illustrative example of how the complex mode frequency, 7 GS2 operates in terms of kx and ky with k 2 ⊥ = k 2 x + k 2 y that can be directly related to the form employed in section 4.1.
ω + i γ, changes with the ion Larmor radius, ρ i . The scan in figure 3 is obtained at fixed k y ρ i , noting that then the conventional GS2 solution (i.e. with the Maxwellian F) is independent of ρ i . The neoclassical corrections start playing a role at ρ i /R ≳ 10 −3 . Figure 3 is obtained at a low beta, representative of conventional tokamaks. In figure 4 we show how the complex mode frequency changes with increasing β. In addition to ρ, the second factor that influences the impact of the neoclassical terms is the equilibrium density/temperature gradient. In figures 5 and 6 we show the complex mode frequency as a function of the density and temperature gradient length scales. This clearly demonstrates a higher impact of neoclassical terms in steep gradient regions. Note that in the pedestal region, we expect typical values of RL −1 n,T ∼ (15-50). In the initial value solver mode, GS2 provides the solution for the most unstable mode. For the conventional GS2 CBC case shown in figures 5 and 6, these are KBMs that are characterised by a rapid increase in the mode growth rate with increasing β ( figure 4). In addition, their growth rate increases with the density and temperature gradients, also as expected for KBMs. The mode frequency shown in figures 5 and 6 for the NEO GS2 case smoothly matches the conventional ω at lower gradient length scales. This therefore can still be identified as KBMs, modified by the neoclassical corrections. However, it is interesting to note that the growth rate saturates with increasing RL −1 n for the neoclassical case of figure 6 at ρ * = 10 −2 .
To investigate this saturation further, in figure 7 we present the effect of the radial electric field shear, defined here via − ∂ 2 Φ 1 0 /∂ε 2 / ∂Φ 1 0 /∂ε , on the mode growth rate. In accordance with figure 7(a), it slightly decreases with the inverse temperature gradient length scale, and so the growth rate increases at larger RL −1 T compared to the conventional GS2 prediction ( figure 7(b)). In contrast, − ∂ 2 Φ 1 0 /∂ε 2 / ∂Φ 1 0 /∂ε grows rapidly with RL −1 n in figure 7(a), correlating with the rapid mode growth rate saturation at a larger RL −1 n in figures 6 and 7(c). Note that the growth rate saturation is not observed for the analytic F 1 (see figure 7(c)), where Φ 0 is a function of ψ only and hence does not impact F 1 nor the fluctuating solution.     This therefore highlights the importance of retaining the poloidal variation in the equilibrium electrostatic potential, Φ 1 0 = Φ 1 0 (ψ, ϑ), calculated consistently with plasma quasineutrality (see appendix D for the effects of different collision operators in the equilibrium drift kinetic equation on the KBM growth rate). Stabilising effects of the radial electric field (and the toroidal flow shear) have previously been reported in [40] to explain the turbulence reduction at higher gradients, in [41] for ITGs or [42] in flux tube gyrokinetic turbulence simulations with the Maxwellian background, for example.
In figure 8 we perform a density gradient length scale scan in the absence of the electron temperature gradient and a significantly reduced ion temperature gradient to show that the growth rate saturation still persists in this case and is therefore not related to the ITG-KBM mode coupling. The mode transition can be seen as a jump in the mode frequency, ω. For the conventional GS2 CBC case, this can be identified as the transition from TEM to KBM at the normalised density gradient length scale, RL −1 n ≈ 27 (β = 0.5%) and RL −1 n ≈ 8 (β = 2%). Note that while the KBM growth rate is known to increase rapidly with beta, as well as with the density and temperature gradients, the growth rate characteristic of TEM typically has a much weaker dependence on beta and density gradient 8 . As one would expect, inclusion of neoclassical corrections provides a higher impact at higher equilibrium gradients. In particular, the equilibrium density gradient value when the neoclassical corrections start playing a role is lower at higher beta values: RL −1 n ≈ 10 (β = 0.5%), RL −1 n ≈ 8 (β = 2%), RL −1 n ≈ 5 (β = 20%) for the equilibrium parameters of figure 8, for example 9 . For the NEO GS2 case shown in figure 8, the mode can be identified as TEM at β = 0.5% with slowly decreasing growth rate with increasing density gradients at larger RL −1 n . Similar to the conventional GS2 case, there is the mode transition at RL −1 n ≈ 8 at β = 2% from TEM to KBM, modified by the neoclassical contributions (its growth rate increases with β and with the equilibrium gradients). Similar to figure 6, the mode growth rate is again significantly suppressed at β = 2% for larger RL −1 n . It is also interesting to note that the impact of the neoclassical corrections is much less, when the mode is a TEM than when it transitions to a KBM.
In figure 9 we consider the situation when both RL −1 n and RL −1 T have pedestal-relevant values. In figures 9(a) and (b) we plot the frequency and growth rate of the most unstable mode as a function of β at k y ρ i = 0.3. The jump in the GS2 mode frequency at β ≈ 0.3% represents a transition from the electrostatic limit with γ being relatively insensitive to β to 8 The TEM core plasma cases typically also include the ion temperature gradient mode (ITG) for the same range of the kyρi variation. In the case considered in figure 8, as can be anticipated, the observed TEM disappears in the absence of the electron gyrokinetic component. 9 The density gradient length scale scan at β = 20% is not presented here. a KBM at higher beta. In this case, the electrostatic limit is provided by a TEM driven by the density gradient and with the mode frequency being in the electron diamagnetic direction (ω < 0). Inclusion of the neoclassical effect smooths this mode switching, significantly reducing the KBM growth rate. The latter is mainly the consequence of the radial electric field shear effect, as well as F 1 being calculated consistently with plasma quasi-neutrality. Indeed, the value of RL −1 n in figure 9 is large enough to deviate the radial electric field shear from its conventional gyrokinetic value characteristic of the low density gradient (as shown in figure 7(a)). The corresponding growth rate spectrum at large density and temperature gradients is presented in figure 9(c) for two different beta values.
All the provided examples demonstrate that the neoclassical corrections influence the gyrokinetic-Maxwell solution even for the CBC (conventional tokamak core plasma) case. Their impact grows with the equilibrium gradients, ρ and with β. Therefore, we would expect their contribution to become non-negligible for the tokamak pedestal region and STs. An example of the latter is considered in the following section.

ST core plasma
The equilibrium used in this section is based on the finite beta ST equilibrium of [43] (for a detailed description of the Miller parametrisation process, see [44]). The equilibrium parameters are as follows: the Miller equilibrium tokamak geometry [45], the normalised radius of the surface of interest, r/a = 0.69, normalised major radius of the centre of the surface of interest, R 0 /a = 1.85, the safety factor q = 3.5, magnetic shear s = 0.54, plasma elongation κ = 2.82 and triangularity △ = 0.33, elongation shear s κ = −0.21 and triangularity shear s △ = −0.1, β = 0.16 and β ′ = −1.15 (where here a prime denotes the derivative with respect to r) defined via the gradient of the plasma pressure, ignoring the gradient of the magnetic field, and ρ i /a = 4 × 10 −3 . For simplification, we assume a model based on two kinetically modelled species (bulk Z i = 1 ions and electrons) with the ratio of the  In figures 10 and 11 we show the artificial density and temperature gradient length scale scans based on the solution of the conventional gyrokinetic-Maxwell system and the solution of equations (21)- (24) in the presence of the neoclassical equilibrium effects. At this relatively low value of the normalised Larmor radius, the impact of the neoclassical corrections on the stability is relatively small, except at the highest gradients. Similar to the CBC case, the dominant instability in figures 10 and 11 can be identified as a KBM (k y ρ i = 0.21): its growth rate increases with both density and temperature gradient length scales, as well as with beta ( figure 12(a)). From figure 12, at the reference equilibrium ion/electron density and temperature gradient length scales aL −1 n = 0.43 and aL −1 T = 3.35, and β consistent with the equilibrium parameters at the flux surface of interest, r/a = 0.69, the mode growth rate almost matches the conventional GS2 γ in the absence of neoclassical corrections. However, switching to a more central flux surface, r/a = 0.53, results in an increase of β, and the differences between the NEO GS2 solution and the conventional one are more visible. In particular, higher order gyrokinetics is found to decrease the KBM threshold beta, increasing the KBM growth rate, compared to the conventional gyrokinetic predictions (see figure 12(b)).
In this example, at k y ρ i > 1, there is a mode transition from KBMs to micro-tearing modes (MTMs) for the most unstable modes (see figure 13) with the mode frequencies being in the electron diamagnetic direction ω < 0. A discontinuity in the mode frequency in figure 13(a) at k y ρ i ≈ 3 represents a transition between two different MTMs. Like KBMs, MTMs strongly depend on beta. Indeed, both these modes require a threshold beta to be exceeded to drive them unstable [46,47]. However, in contrast to KBMs, the impact of the higher order  theory on the MTM growth rate is found to be minimal in this case.

Discussion
The simulations presented here have been motivated by the importance of the bootstrap current (and associated pressure gradients) for ideal MHD stability (e.g. high n kink-ballooning modes), suggesting that it is likely to be important for electromagnetic drift instabilities and the associated turbulence. We have therefore extended gyrokinetic theory to incorporate neoclassical flows and the bootstrap current in particular. To simplify the gyrokinetic Maxwell field equations, we have only considered the B ϑ /B 0 ≪ 1 limit of [12], i.e. equation (12) (global formulation) or the system of equations (21)-(24) (local formulation). The latter has been implemented in the local gyrokinetic code, GS2, to quantify the impact of the higher order, neoclassical corrections in gyrokinetics and can be considered as an extension of the previously developed, electrostatic, low flow version of GS2 [10,[13][14][15]] to a finite/high beta, electromagnetic case.
As can be anticipated, there are three main factors that determine the importance of the neoclassical corrections to gyrokinetics: (1) the particle Larmor radius, (2) the equilibrium density and temperature gradient length scales and (3) beta. The conditions of the tokamak pedestal are most likely to be where the neoclassical corrections are important, especially when the KBM plays a role. While the value of ρ i /L ≳ 10 −3 is already sufficient for the neoclassical corrections to play a role, the equilibrium gradient values must be carefully selected to ensure consistency with the perturbative treatment of the equilibrium solution and F 1 in particular.
Based on the test cases considered in section 5, we summarise • The KBM growth rate is found to be significantly suppressed by the inclusion of neoclassical effects at large density gradients, representative of the tokamak pedestal values, at low (figure 8(d))/moderate (figure 6)/large (figure 9(b)) temperature gradients. This growth rate suppression primarily results from the effect of the radial electric field shear at large density gradients. • The radial electric field shear response to increasing temperature gradients is opposite to that associated with increasing density gradients. However, in contrast to RL −1 n , increasing RL −1 T only slightly impacts the radial electric field shear (especially at lower RL −1 n ), which then explains why no growth rate saturation is observed at a larger RL −1 T in figure 5 or in section 5.2.
• While then the main factor responsible for suppressing/enhancing the mode growth rate is the neoclassical electrostatic potential, Φ 1 0 = Φ 1 0 (ψ, ϑ), calculated consistently with plasma quasi-neutrality, we also note that it is important to accurately resolve the collisional boundary layer that surrounds the trapped-passing boundary in pitch angle space. This layer introduces the physics of trapped particles into the dynamics, influencing the form of Φ 1 0 (ψ, ϑ) and its radial derivatives.
• In certain cases, inclusion of neoclassical effects is found to prevent mode switching (e.g. TEM-KBM in figure 8(b)). Note that in the initial value solver mode, GS2 provides the result for the most unstable mode, i.e. other, more stable modes might still be present. • The main impact of the neoclassical corrections is found to be in the range of 0 < k y ρ i ≲ 1. • While the dominant impact of the neoclassical equilibrium physics is found to be on KBMs (both the KBM growth rate and threshold beta), electrostatic modes (k y ρ i ≲ 1) are influenced as well at large density and temperature gradients (figure 9).

Future work
Here we have limited the numerical simulations to linear theory to provide an initial assessment of the importance of neoclassical corrections, ensuring consistent ordering by considering the limit of B ϑ /B 0 ≪ 1. At B ϑ /B 0 ≪ 1, the only nonlinear effect will be provided by the conventional nonlinear term, i.e. 1 B0 ⟨∇χ ′ × b⟩ X α · ∂g ∂X to be added to the right hand side of equation (12). Indeed, it can be demonstrated that turbulent effects on the neoclassical transport (that must be included in the equilibrium equation when the problem is nonlinear) contribute only to the gyro-angle independent piece of the second order equilibrium distribution function, F 2 , and the gyro-phase average of the corresponding term is zero [5]. The nonlinear effects will be addressed in our future publication. We also note that the B ϑ /B 0 ≪ 1 assumption might be inaccurate in certain cases (especially those relevant to STs), and therefore the O(∆ 2 Vt L F 0 ) terms of [12] (neglected here) are likely to be important for a full treatment of the current density gradient effects (and kink/peeling physics in particular) in regions where the equilibrium profile gradients are strong (specifically in STs). The O(∆ 2 Vt L F 0 ) terms of [12] will be incorporated in GS2 as part of future work. Note that in the pedestal, the pedestal width, ∆ ped ∼ (10-40)ρ i , [48][49][50][51] is typically chosen for the characteristic length scale, which therefore might impose additional limitations on ∆ = ρ i /L ∼ ρ i /∆ ped (and therefore δ) as the expansion parameter in the δf-type kinetic formulations. In particular, the latter might limit the applicability of the presented higher order gyrokinetic theory (and generally any δf-type approach based on ∆ ≪ 1) in certain tokamak pedestal cases, i.e. depending on how ∆ ped scales with ρ i , as well as B ϑ /B 0 in the pedestal. In situations when δ ≲ 1 but ∆ remains small, one can employ ∆ as an expansion parameter in the equilibrium solution, which will then require one to retain the finite Larmor radius corrections in F 2 , as discussed in [12].
Finally, while local gyrokinetic simulations are relatively well understood, they are limited by the assumptions of the ballooning transform. The eikonal approximation of section 4.1 (and hence the set of equations (21)-(24)) eliminates formulational problems of nonlocal approaches (e.g. absence of unique sources or radial boundary conditions), which makes it convenient for the initial assessment of the impact of higher order, neoclassical equilibrium effects. Furthermore, this is fully justified, provided there is no strong functional dependence of the fluctuating quantities on B ϑ /B 0 . Generally, to allow for the radial profile variation around the flux surface of interest, future work is required to implement the gyrokinetic theory of [12] in a hybrid gyrokinetic code. In addition, while the equations derived in [12] are generally global, they require accurate boundary conditions/generalisation of the ballooning transform to be employed to capture coupling to the vacuum solution and therefore to allow for the external modes. This physics is likely important for a full understanding of the impact of the current density (i.e. kink drive) in pedestal gyrokinetics, and will be addressed in future work.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https://doi.org/10.5281/ zenodo.4643844.

Appendix A. Example neoclassical equilibrium solution for large aspect ratio
For benchmarking purposes, in this section we employ a number of simplifying assumptions and present F 1 , based on a specific two species model suitable for analytic treatment. F 1 is provided by the gyro-averaged O(∆ 1 Vt Here the drift velocity is defined as with E 0 being the equilibrium electric field. Note, we assume that the particle collision frequency ∼V t /L, i.e. is ∆ times smaller than the cyclotron frequency, ω c ; and V D = |V D | ∼ ∆V t . C on the right hand side of equation (A.1) is the gyroaveraged collision operator. In the drift kinetic approach, the difference between C and its gyro-averaged version is neglected. It can be demonstrated that for the like-species Fokker-Planck collision operator, this difference is zero [21].
where Φ 0 is the equilibrium electrostatic potential, and prime denotes its derivative with respect to ψ. Treating collisions perturbatively in equation (A.3) (e.g. relevant to the banana collisionality regime), we write for the leading order solution in collisions.F 1 is the constant of integration provided by the solubility condition at next order: where q is the safety factor. ⟨. . .⟩ ψ ϑ represents the orbit averaging operator at fixed ψ, i.e. an average over poloidal angle, ϑ, for passing particles and an integral between bounce points with a sum over σ = u/|u| for trapped particles. In both cases the integrals in poloidal angle are on the flux surface of interest, ψ. Equation (A.5) requires a specific model for the collision operator. Here we employ a momentum-conserving, pitch angle scattering collision operator (equations (62)-(64) of [52]): (A.6) for ions and (A.8) for electrons. Here u ∥j = (3π 1/2 V 3 tj )/(2n 0 )´F j u/V 3 dV and u ∥j = (1/n 0 )´F j udV (j = i for ions and j = e for electrons); V tj is the thermal velocity of species j and n 0 is the equilibrium density. F M j is the Maxwellian distribution of species j, and ν jj and ν ei denote the corresponding collision frequency. λ = 2µ/V 2 is the pitch angle. We also adopt a large aspect ratio, circular cross section tokamak equilibrium for simplicity. The system of equations (A.5) with equations (A.6)- (A.8) is an integro-differential system coupled via the electron-ion collisions. It is solved with the following boundary conditions in λ space: (1) we require a mixed boundary condition at the deeply passing/trapped region to ensure the solution is finite and (2) the solution and its first derivatives must be matched at the trapped-passing boundary. After one iteration in the momentum-conserving term, we find  Figure C1. Example of the neoclassical, first order electrostatic potential, Φ 1 0 , and its radial, ∂Φ 1 0 /∂ε, and poloidal, ∂Φ 1 0 /∂ϑ, derivatives plotted as functions of ϑ at r = 0.5a (a is the minor radius, r, at the last closed flux surface). The plasma equilibrium parameters and normalisation are chosen as in figure 2. 'NEO' denotes results calculated by NEO on the NEO theta (ϑ) grid. 'GS2' corresponds to Φ 1 0 and its derivatives reconstructed from the NEO Φ 1 0 on the GS2 theta (ϑ) grid.

Appendix D. Impact of collisions on the growth rate saturation at larger density gradients
In this section, we consider how significant the impact of different collision operators is on the growth rate saturation observed in figure 6. In figure D1 we show the KBM growth rate, γ, as a function of the equilibrium density gradient length scale, RL −1 n , at RL −1 T = 6.9 for different collisional models. 'GS2' corresponds to the solution of the conventional gyrokinetic-Maxwell system with the Maxwellian background in the absence of the neoclassical effects. 'NEO GS2, FP', 'NEO GS2, Connor', 'NEO GS2, HS0' and 'NEO GS2, HS' correspond to the solution of equations (21)- (24) with F 1 obtained with the full linearised Fokker-Planck [24], pitch angle scattering (Connor) [34], zeroth-order Hirshman-Sigmar [35] and full Hirshman-Sigmar collisional model [55], respectively, for s-alpha geometry. 'Analytic' denotes the solution of equations (21)-(24) based on the analytic F 1 obtained in appendix A. The main features of each model are summarised in table D1.
In figure D1 the growth rate saturation at larger density gradients is observed in NEO GS2 even for a simplified pitch angle scattering collision operator in the absence of friction and diffusion-type operators. This therefore suggests that the main factor responsible for suppressing the KBM growth rate is the self-consistent neoclassical electrostatic potential, Φ 1 0 = Φ 1 0 (ψ, ϑ). We also highlight the importance of accurately capturing the collisional trapped-passing boundary layer physics as this provides collisional dissipation even in the limit of rare collisions via matching the passing and trapped particle distribution functions in pitch angle space, as well as introduces the trapped particle physics into the dynamics. The latter then influences the form of Φ 1 0 (ψ, ϑ) via quasi-neutrality.  Figure D1. The KBM growth rate, γ, plotted as a function of the equilibrium density gradient length scale, RL −1 n , at RL −1 T = 6.9 for the equilibrium parameters of section 5.1. The normalised ion collision frequency is ν i = 0.19Vt/a.