On the importance of parallel magnetic-field fluctuations for electromagnetic instabilities in STEP

This paper discusses the importance of parallel perturbations of the magnetic-field in gyrokinetic simulations of electromagnetic instabilities and turbulence at mid-radius in the burning plasma phase of the conceptual high-β, reactor-scale, tight-aspect-ratio tokamak STEP. Previous studies have revealed the presence of unstable hybrid kinetic ballooning modes (hKBMs) and subdominant microtearing modes at binormal scales approaching the ion Larmor radius. In this STEP plasma it was found that the hKBM requires the inclusion of parallel magnetic-field perturbations to be linearly unstable. Here, the extent to which the inclusion of fluctuations in the parallel magnetic-field can be relaxed is explored through gyrokinetic simulations. In particular, the frequently used MHD approximation (dropping δB∥ and setting the ∇B drift frequency equal to the curvature drift frequency) is discussed and simulations explore whether this approximation is useful for modelling STEP plasmas. It is shown that the MHD approximation can reproduce some of the linear properties of the full STEP gyrokinetic system, but is too stable at low ky and nonlinear simulations using the MHD approximation result in very different transport states. It is demonstrated that the MHD approximation is challenged by the high β′ values in STEP, and that the approximation improves considerably at lower β′ . Furthermore, it is shown that the sensitivity of STEP to δB∥ fluctuations is primarily because the plasma sits close to marginality and it is shown that in slightly more strongly driven conditions the hKBM is unstable without δB∥. Crucially, it is demonstrated that the state of large transport typically predicted by local electromagnetic gyrokinetic simulations of STEP plasmas is not solely due to δB∥ physics.


Introduction
The performance of magnetic-confinement-fusion devices such as spherical tokamaks (STs) is often limited by the presence of turbulent fluctuations that typically dominate the transport losses of heat, particles, and momentum.Therefore, understanding and predicting turbulent transport in next-generation STs is crucial for the optimisation of their performance.The UK STEP programme aims to generate net electric power P el > 100 MW from fusion [1,2], by developing a compact prototype power plant, STEP, based on the ST concept.The first phase of this ambitious programme is to provide a conceptual design of a STEP prototype plant and reference plasma equilibria; in this context, a set of preferred STEP plasma scenarios are being designed [3], using simplified reduced models to model the core transport.
Reduced core-plasma transport models are essential for the integrated modelling of tokamak plasma scenarios, since first principles approaches are generally too computationally expensive to evaluate over confinement or resistive discharge timescales.Conventional-aspect-ratio devices tend to be dominated by electrostatic microturbulence driven by instabilities, such as modes driven by the ion-temperature gradient (ITG), modes driven by the electron-temperature gradient (ETG), and trapped-electron-driven modes (TEMs), and several reduced-complexity models have been developed to model electrostatic microturbulence (see e.g.[4,5]).Such reduced models have been quite successful at predicting turbulent transport in conventional-aspect-ratio present-day devices at low-to-modest β [where β = 2µ 0 p/B 2 is the ratio of the total plasma pressure p to the magnetic-field energy density B 2 /(2µ 0 ), with B the magnetic-field strength and µ 0 the permeability of free space], where the turbulence is predominantly electrostatic: see e.g.[6].However, in STs, the confinement and turbulent transport, reviewed in [7], are both quantitatively and qualitatively different to that in conventional-aspect-ratio devices.ST geometry has favourable stability properties that allow access to higher-β regimes where the turbulence is more electromagnetic in character and the dominant unstable modes include kinetic ballooning modes (KBMs) and microtearing modes (MTMs).These electromagnetic modes, less well understood than their electrostatic counterparts, are not well captured in the most advanced reduced core transport models; transport from electromagnetic turbulence has thus not been reliably captured in scenario modelling for STEP.
Fortunately, most microinstabilities share broad characteristics which are well described by local ‡ linearised δf gyrokinetics (GK) provided that k ⊥ ρ s ∼ 1 and ρ ⋆s ≡ ρ s /a ≪ 1 where ρ s is the gyroradius of species s, a is a typical equilibrium length scale, and k ⊥ is the perpendicular wavenumber of the instability.High-fidelity GK simulations must be exploited: (i) to assess electromagnetic microinstabilities and their associated turbulence in the conceptual power-plant designs developed for STEP; and (ii) to improve the physics basis of reduced models used in scenario design.
In this paper, we concern ourselves with electromagnetic instabilities in STEP, in particular the hybrid-kinetic ballooning mode (hKBM) [8] that GK simulations indicate to dominate turbulent transport [9] in the conceptual designs of STEP.The inclusion of ‡ All simulations in this paper are local simulations; i.e. simulations are performed in a domain of perpendicular size that is infinitesimal in comparison with the length scale over which the equilibrium varies.Plasma equilibrium gradient length scales are taken to be constant across the simulation domain.
compressional magnetic fluctuations δB ∥ was found to be essential in [8] for the hKBM to be unstable in STEP (see also discussion in Section 4).However, δB ∥ physics is often simply missing from many higher-fidelity (e.g.global, full-f ) GK codes and modelling tools.Instead, in GK theory and simulations of microinstabilities, it is common practice to neglect δB ∥ , and to compensate for this by exploiting an approximate cancellation (valid in certain limits) that relies on the form of δB ∥ in gyrokinetics and an exact relationship between the equilibrium ∇B and curvature drifts.In this work, we will refer to this compensation as the MHD approximation § (see e.g.[10][11][12]); so called because it relies on (i) an MHD equilibrium force balance constraint on the relationship between the ∇B and curvature drifts, and (ii) the gyrokinetic form of δB ∥ which excludes high frequency fast compressional Alfvén waves by enforcing perpendicular force balance on fluctuations δj ⊥ × B = ∇ ⊥ • δp ⊥ .The primary focus of this paper is to explore sensitivity of hKBM stability to the inclusion of δB ∥ and the appropriateness of the MHD approximation in STEP plasmas.
The paper is organised as follows: Section 2 briefly reviews the results of [8] and [9] and outlines the motivation for using the MHD approximation in STEP plasmas.Following a brief review of electromagnetic δf gyrokinetics, Section 3 introduces how the MHD approximation fits into the GK framework and its implementation in GK codes.Section 4 introduces the STEP equilibria and associated plasma parameters (with more details available in [8,9]) and presents linear GK simulations (using the GK code GENE [13]) to assess the validity of the MHD approximation in STEP.It is explored whether the performance of the MHD approximation can be improved in conditions where the approximation's underlying assumptions are better satisfied.In Section 5 we demonstrate and explain why the MHD approximation is typically not appropriate for modelling electromagnetic turbulence in the STEP flat-top plasma.We again discuss cases where the MHD approximation can more faithfully model high−β turbulence.Section 6 explains that δB ∥ is only required for instability in STEP because the local equilibrium is close to marginal stability (owing to substantial β ′ stabilisation).We demonstrate that it is in fact possible to drive electromagnetic hKBM turbulence in the absence of δB ∥ in conditions further from marginal stability with or without implementing the MHD approximation on the magnetic drifts.The final summary and conclusions are presented in Section 7. Additional auxiliary material is provided in Appendices, including a brief derivation of the MHD approximation that is applied in GK in Appendix A, a discussion of the why one of two proposed implementations of the MHD approximation is more physical in Appendix B, a useful expression for the magnetic-drift velocity given in Appendix C, and turbulent snapshots in Appendix D from different simulations cited at various points in the main text.§ In Section 3.2, we will describe two separate implementations of the 'MHD approximation' that have been included in some gyrokinetic codes.We will refer to these two different approximations as 'MHD-1' and 'MHD-2'.In all gyrokinetic simulations using the MHD approximation, δB ∥ fluctuations are set to zero and a magnetic drift term is dropped to compensate for this.

The importance of parallel magnetic-field perturbations
Recent GK analyses for mid-radius of the preferred flat-top operating point, STEP-EC-HD [3], revealed the presence of unstable hybrid KBMs∥ (hKBM) and subdominant MTMs at binormal scales approaching the ion Larmor radius [8,9] (note that similar MTMs were reported for an earlier proposed variant of STEP [14]).Nonlinear local GK simulations for the same local equilibrium in [9] find that the hKBM turbulence can drive heat fluxes that exceed the available heating power by orders of magnitude in the absence of equilibrium flow shear, and that this large transport state is characterised by highly radially extended turbulent eddies.Further simulations in [9] indicate that a transport steady state can be reached if equilibrium flow shear and/or β ′ stabilisation are sufficient at the flat-top to regulate the turbulence.(A transport steady state for a STEP flat-top was recently reported [15] using a physics-based model of transport from hKBM turbulence.)At lower β ′ and/or in the absence of equilibrium flow shear, however, hKBMs drive very large turbulent transport in all channels that may challenge access to the burning flat top through lower β ′ states.
A natural extension to nonlinear simulations in [9] is to attempt higher-fidelity global gyrokinetic simulations for conditions where local GK predicts large turbulent fluxes.δB ∥ , indicated by local GK to be essential for the hKBM to be unstable in STEP-EC-HD [8], is unfortunately neglected in most global GK codes.To progress towards a global description of hKBM-driven turbulence, we must both (i) implement δB ∥ in a global-capable code; ¶ and (ii) strive to better understand the sensitive dependence of the hybrid mode on δB ∥ .This paper focuses on the latter.

The MHD approximation for the magnetic-drift velocity
As remarked in Section 1, in GK theory and simulations of microinstabilities it is common to neglect the (typically destabilising [20,21]) + parallel magnetic perturbation δB ∥ at low β, and to compensate for this by artificially adding a correction term proportional to ∇p to the magnetic drift velocity (see e.g.[10,12,22] and references therein).This idea of modifying the magnetic drifts, hereinafter referred to as the MHD approximation, is recommended in many GK codes on the neglect of δB ∥ .The merits, or lack thereof, of the MHD approximation have been tested in different GK simulations by various authors (see e.g.[20], [12]).Whilst proper treatment of δB ∥ physics is often crucial [20], the cancellation of δB ∥ and the ∇p term in the ∇B drift has been demonstrated to be a good approximation when calculating the growth rates of ∥ The term 'hybrid KBM' was introduced in [8] to describe the dominant instabilities seen in STEP plasmas; typically these modes share attributes associated with KBM, ITG, and TEM.¶ Exciting progress is being made in this direction (see e.g.[16][17][18]), and there are several parallel efforts working towards implementing δB ∥ globally in both Eulerian codes such as GENE [13] and Particle-In-Cell codes such as ORB5 [19].+ It was demonstrated in [21] that neglecting δB ∥ typically is stabilising for MHD-like modes because this omits the "magnetic compression" term in δW (the potential energy of the plasma).If δB ∥ is ignored then δW is minimised incorrectly which has a stabilising effect on pressure-driven modes.
KBMs at low β and long wavelengths [12].Here we explore whether this approximation may be appropriate for use in STEP plasmas * .
To reiterate, the sensitivity of hKBM stability to δB ∥ poses difficulties for exploiting existing codes for STEP.If the MHD approximation were valid for STEP, it would be useful for two key reasons: • Global gyrokinetic codes often neglect δB ∥ .
Testing high turbulent fluxes predicted, but potentially not well resolved, by local GK in some local equilibria is of high priority for STEP design.Progress requires global gyrokinetic simulations including profile variation that should better resolve the box-scale streamers observed in local simulations [9].However, most globalcapable codes simply do not include δB ∥ at present because of the associated algorithmic complexity.The GK equation is closed by self-consistently solving the gyrokinetic form of Maxwell's equations to evaluate the perturbed electromagnetic fields.Retaining the perpendicular part of Ampère's law introduces a coupling between the electrostatic potential and the perpendicular magnetic potential δB ∥ .Including δB ∥ consistently requires the solution of a potentially poorly-conditioned system of equations involving δϕ and δB ∥ , as well as derivatives and gyro-averages of these quantities.In a local code, these derivatives and gyro-averages are trivial to handle via spectral methods.In a global code, such an approach is not possible and one has to instead introduce complex "gyrodisk averages" to obtain a system which is numerically well behaved.(see [17,18]).The MHD approximation has been used in a global code to study electromagnetic turbulence in NSTX plasmas [23].
• Integrated modelling tools often neglect δB ∥ Quasilinear turbulence models such as TGLF [4] provide fast transport predictions that are essential for the integrated modelling codes used in scenario design.These models perform a fast computation of the linear-mode properties, and feed these into saturation rules developed to parameterise the turbulent fluxes.As such, these reduced models rely on both the fidelity of the saturation rules, and the accuracy of the linear physics calculation, which requires capturing δB ∥ effects in STEP.
Sections 4-6 of this paper are devoted to: (i) exploring the degree to which the MHD approximation is appropriate for use in STEP plasmas; and (ii) seeking to understand conditions where the requirement to include δB ∥ fluctuations can be relaxed.

Linear electromagnetic gyrokinetics and the MHD approximation
We are interested in plasmas that are well described by the GK framework (see e.g., [24,25]): i.e. we are concerned with fluctuations ∼ e i(k•r−ωt) , having characteristic frequency ω and wavenumbers k ∥ and k ⊥ parallel and perpendicular to the equilibrium where Ω s = q s B/m s is the cyclotron frequency of species s with charge q s , equilibrium density and temperature n 0s and T 0s , respectively, mass m s and thermal speed v ths = 2T 0s /m s , ν ss ′ is the typical collision frequency.The perturbed electrostatic potential is represented by δϕ.Electromagnetic perturbations enter GK through δB ∥ and δB ⊥ , the fluctuations of the magnetic-field parallel and perpendicular to the equilibrium magnetic-field direction.Electromagnetic effects are most conveniently described in local GK by writing the fluctuating magnetic field δB = ∇ × (δA It is convenient to write the GK guiding-centre distribution function in the form Here, the GK distribution function consists of a Maxwellian piece F 0s and an order-ρ ⋆s -small perturbation δf s , with g s being the non-adiabatic part of δf s .Under the above ordering, the collisionless, linear, GK equation is given in Fourier space by where: b s = k ⊥ v ⊥ /Ω s ; J 0 and J 1 are Bessel functions of the first kind that arise from gyroaveraging; the drive frequency ω T ⋆s = ω ⋆s [1 + η s (v 2 /v 2 ths − 3/2)] involves η s = d ln T 0s /dn 0s and the diamagnetic frequency ω ⋆s = (k ⊥ T 0s /q s B)d ln n 0s /dr is defined in terms of the radial coordinate r.The magnetic drift frequency is given by where the ∇B drift frequency ∝ ω ∇B = k • b × ∇B/B, and the curvature drift frequency The equilibrium magnetic field can be written as B 0 = ∇ψ×∇α, where ψ is poloidal flux, and α = ξ − q(ψ)θ − ν(ψ, θ) is the binormal coordinate, with ξ the toroidal angle, θ the poloidal angle, and ν(ψ, θ) a periodic function of θ that depends on flux-surface shaping.Perpendicular gradients can be expressed in terms of binormal and radial wavenumbers as The fluctuating field quantities appearing in the GK equation are determined through the field equations.The perturbed electrostatic potential δϕ is determined through quasineutrality The parallel magnetic vector potential δA ∥ is determined by parallel Ampère's law ♯ Note that ω κ and ω ∇B [appearing in (2)] are not themselves frequencies.
The magnetic fluctuation δB ∥ is determined by perpendicular Ampère's law, leading to: The MHD approximation is obtained in the long wavelength limit, where (5) reduces to: involving the perturbed gyrokinetic perpendicular pressure that is defined in this limit as: Equation 6 implies that perpendicular force balance is satisfied to leading order in gyrokinetics, which is essential to exclude fast compressional Alfvén modes with frequencies ∼ Ω/ √ β that lie outside of the gyrokinetic approximation (as discussed in Appendix A of [26]).

The MHD approximation in GK
The historical precedent for the MHD approximation in GK comes from the fact that, in the one fluid MHD limit (k y ρ s → 0), the compressional magnetic fluctuations cancel the contribution of the ∇B drift in a low-β expansion in the absence of magnetic curvature [10].Our first goal in this paper is to understand how this cancellation manifests itself in the GK framework.
Linear GK force balance in 6 requires finite δB ∥ to sustain a perpendicular pressure perturbation.It is therefore only appropriate to drop δB ∥ if δP ⊥ can also be neglected.In Appendix A we calculate δP ⊥ from the perturbed gyrokinetic distribution function, exploiting an exact relationship (given in Appendix C) between curvature and ∇B drifts that follows straightforwardly from equilibrium force balance: whereby ω κ and ω ∇B differ by a term proportional to the pressure gradient.It is demonstrated in Appendix A that δP ⊥ vanishes in the long-wavelength, low-β limit with k ∥ v ths ≪ ω ds ≪ ω, if: Enforcing equation 9 can therefore be used to implement the MHD approximation and drop δB ∥ from GK in the appropriate limits.Two approaches have been adopted to enforce equation 9 in gyrokinetics.Denoting the physical values (i.e.those set by the actual equilibrium magnetic geometry) of ω κ and ω ∇B by ω act κ and ω act ∇B , we can satisfy 9 by either: • (MHD-1) setting the ∇B drift equal to the curvature drift Both have been implemented in GK codes.Appendix B explains why MHD-1 is better motivated physically for the study of a particular equilibrium, which is consistent with arguments given in Appendix A of [27] and Appendix F.5 of [28].

Implementation of the MHD approximation in gyrokinetic codes
It is helpful to connect the MHD approximation to the magnetic-drift velocity v d rather than the magnetic-drift frequency ω ds = k ⊥ • v ds .The exact magnetic drift velocity is given by where it is shown in Appendix C that: We can now understand both implementations of the MHD approximation in terms of adding or subtracting corrections to this definition.

The MHD-1 approximation
The historical precedent for setting the ∇B drift equal to the curvature drift comes from a result of [11] showing that, in the one fluid MHD limit (k y ρ D → 0), there is a cancellation of terms such that δB ∥ can be dropped as long as the ∇B drift is corrected by a term proportional to ∇p so that it is equal in magnitude to the curvature drift and points in the same direction.This will be referred to as "MHD-1": In MHD-1 a term proportional to the pressure gradient has been added to the right hand side of Equation 10.

The MHD-2 approximation
In MHD-2 the drift velocity is written as: As discussed in Appendix C, this effectively corresponds to changing the equilibrium magnetic field curvature by setting b • ∇b = ∇B/B, and again removing the pressure gradient contribution to the ∇B drift in this different local equilibrium.MHD-2 allows B ∥ to be dropped from the GK equation ( 1), but modifies the local equilibrium.
MHD-1 is the more physical implementation because it simply exploits a cancellation between two terms without changing the local equilibrium curvature; henceforth where we use the term "MHD approximation" and this will imply MHD-1.† †

Linear gyrokinetic simulations using the MHD approximation
The gyrokinetic analysis in this paper focuses primarily on a single equilibrium fluxsurface taken from close to mid-radius (q = 3.5, Ψ n = 0.49) in the STEP reference scenario STEP-EC-HD † (where EC stands for Electron Cyclotron heating and current drive, and HD stands for High Density), which is designed to deliver a fusion power P fus = 1.8 GW.We note that this is precisely the same flux surface as examined in [9], a choice that has been made intentionally to allow benchmarking between our results against those reported in [9].It should also be noted that whereas [9] focused on accurate investigation of the turbulent transport in STEP; the focus of this paper is instead on a more conceptual study of the hKBM.
A Miller parameterisation [29] was used to model the local magnetic equilibrium geometry, and the shaping parameters were fitted to the chosen surface using pyrokinetics [30], a Python library developed to facilitate pre-and post-processing of gyrokinetic analysis performed using a range of different GK codes.pyrokinetics also contains an ideal-ballooning solver which has been used throughout this work to ensure that all equilibria remain MHD-stable when the local equilibrium parameters are varied.Table 1 provides the values of various local equilibrium quantities at the flux surface examined in this paper, including magnetic shear ŝ; safety factor q; normalised minor radius ρ/a; elongation κ and its radial derivative κ ′ ; triangularity δ and its radial derivative δ ′ ; the radial derivative of the Shafranov shift ∆ ′ ; and the normalised inverse density and temperature gradient scale lengths of species s, a/L ns and a/L T s , respectively.Our simulations evolve three species: electrons, deuterium, and tritium and neglect entirely impurities and fast particles.The interested reader is referred to [8,9] for more details on the equilibrium and on the setup of the computational grids, which are identical to those used in the aforementioned works.
Previous linear analysis shows that the hKBM is the dominant ion-scale instability on this surface, with a subdominant MTM also found to be unstable on a subset of these binormal scales (see Fig 19 and Fig 20 of [8]).No unstable microinstabilities are observed at the electron Larmor radius scale.The dominant hKBM and the subdominant MTM can both be recovered physically; that is, one can recover the subdominant mode by either by exploiting the up-down symmetry in the local equilibrium and forcing the parity of the perturbed distribution function in an initial-value calculation, or by using an eigenvalue solver to return the unstable linear spectrum.However, importantly for our work, it was also shown that for this equilibrium it is possible to recover the   number codes and modelling tools) then we recover the previously subdominant MTM (note the change in frequency) as the fastest growing unstable mode in the system.Succinctly, the hKBM is linearly stable on this surface (along with many other surfaces in the STEP equilibrium) without δB ∥ .

Linear gyrokinetic simulations using the MHD approximation in STEP
The MHD approximation introduced in Section 3.1 may allow us to study the hKBM in the absence of δB ∥ fluctuations.In Section 3.1, we introduced two implementations: (i) MHD-1 (v d given by equation 12); and (ii) MHD-2 (v d given by equation 13), that were designed to compensate for the effect of excluding δB ∥ .We now seek to explore whether either of these approximations are appropriate for modelling STEP plasmas.(note the change in frequency) as the fastest growing unstable mode in the system.‡ However, using the MHD-1 approximation [setting the ∇B drift parallel to the curvature drift (Fig 2 , green)] does find the hKBM but with a strongly reduced growth rate.These results indicate that the MHD-1 approximation is the more appropriate treatment of the magnetic drift when δB ∥ is neglected in STEP plasmas, as we might expect from the fact that MHD-1 does not change the local equilibrium.However, we note that even though this approximation recovers the hKBM, accurately capturing this mode clearly requires a proper treatment of δB ∥ and exact magnetic drifts (particularly at low k y ρ D where MHD-1 suggests hKBMs are stable at wavenumbers where full physics finds the mode to be unstable).

Improving the fidelity of the MHD approximation
This section explores whether it is possible to study the hKBM in the absence of δB ∥ by modifying the local equilibrium parameters such that the MHD approximation is better satisfied.There are three key assumptions in the derivation of the MHD approximation presented in Appendix A: (ii) An MHD-like ordering of the frequencies k ∥ v ths ≪ ω ds ≪ ω.
These assumptions are explored in the following subsections through parameter scans.
4.2.1.Assumption (i) k y ρ D ≪ 1 The long-wavelength assumption in STEP is well satisfied since we are primarily concerned with modes with k y ρ D ≪ 1.

Assumption (ii)
k ∥ v ths ≪ ω ds ≪ ω and Assumption (iii) β ≪ 1 Here we show that the difficulty in improving the fidelity of the MHD approximation lies in satisfying Assumption (ii) and Assumption (iii) simultaneously.
The argument leading to the MHD approximation, given in Appendix A and Appendix B, requires δP ⊥ /(n s T s ) ≪ q s δϕ/(n s T s ).From Equation B.6: it follows that we require:

Now for an Alfvénic instability
‡ In fact, the simulation using the MHD-2 approximation (red) returns the exact same linear spectrum as the simulation that excludes δB ∥ but uses the exact form of the magnetic-drift velocity (orange).This arises because (i) the hKBM is more stable at reduced curvature in the MHD-2 local equilibrium, and (ii) MTMs are typically only very weakly sensitive to ∇p (see Fig 15b of [14]).
Equation 14, together with β ≪ 1, describes the parameter regime where we expect the MHD approximation to hold, and shows that reducing β on its own (i.e. at fixed β ′ ) can make satisfying ( 14) more difficult!Fig 3 shows the results of simulations using different values of β (at fixed β ′ ) and both including, f B = 1 (blue), and excluding, f B = 0 (green, red), δB ∥ fluctuations.Simulations without δB ∥ fluctuations employ either the MHD-1 (green) or MHD-2 (red) approximation.As expected, reducing β at fixed β ′ gives no noticeable improvement in the performance of the MHD approximation.
On the other hand, equation 14 shows that reducing β ′ at fixed β should always improve the fidelity of the MHD approximation, as will be discussed further in Section 6.

Nonlinear GK simulations using the MHD approximation (MHD-1)
From here we proceed using the more physical implementation of the MHD approximation, MHD-1, which still captures hKBMs as the fastest growing linear instability in STEP, albeit with reduced growth rates especially at low k y (see Section 4.1).This suggests that it might be possible to simulate hKBM driven turbulence in STEP using the MHD approximation in a global code, where the full inclusion of δB ∥ is not yet routinely available.As a first step towards this goal, we perform a local nonlinear simulation using the MHD approximation; this simulation is otherwise identical to the hKBM simulation performed in [9] in the absence of equilibrium flow shear.(10) for the drift velocity (this simulation is equivalent to that shown in Fig 3a of [9] and will henceforth be referenced as the 'STEP baseline' simulation).The second simulation (green) shown in Fig 4 does not include δB ∥ and instead uses the MHD approximation (MHD-1).Although the linear physics is broadly similar (both linear spectra have hKBMs as the dominant unstable mode), the nonlinear physics produced by these two simulations is strikingly different.The simulation using the MHD approximation appears to achieve a robustly steady † saturated state at values of the heat flux around two orders of magnitude smaller than those reached by the simulation that includes full δB ∥ physics.It is interesting to note † In this paper, we use the same definition of robustly steady as was introduced in [9].The augmented Dickey-Fuller statistical test [31] is applied to the time trace of the total heat flux in the saturated or pseudo-saturated phase of each simulation (i.e. after the linear growth phase) to determine whether the saturation corresponds to a robust stationary state.The null hypothesis of the augmented Dickey-Fuller test is the presence of a unit root, while the alternative hypothesis is the stationarity of the time series.The null hypothesis is rejected when the p-value returned by the statistical test is below a threshold value that is taken to be p= 0.1 in this paper.The only simulation in this paper that does not saturate by this criterion is the reduced β ′ calculation shown in Figure 5.  that the lower flux MHD approximation simulation does not exhibit the high-frequency oscillations, seen in the full δB ∥ physics case, that appear to be typical of such simulations (see [9]); this is under investigation.Fig 4 also shows the time trace of the k y spectrum of δϕ and δA ∥ for the two simulations.In both cases, the zonal mode dominates over the non-zonal modes and it is likely that zonal flows and fields play a role in saturating the hKBM instability ‡.Despite the difference between the two simulations, it is important to remark that the heat flux predicted using the MHD approximation is still orders of magnitude larger than the heat flux driven by the subdominant MTM (see Fig 14 of From this figure, we can see that long wavelength modes, which are unstable with δB ∥ , are stable with the MHD approximation, and it is precisely these wavenumbers at k y ρ D ≪ 1, that are associated with the very large turbulent fluxes.To be explicit, it appears to be possible to find a lower-flux saturated state even when hKBMs are unstable across a wide range of the linear spectrum provided that they are not unstable up to very long wavelengths.Further evidence for this claim is demonstrated in Fig 5.This figure compares the total heat and particle fluxes from the baseline simulation (blue) to a second simulation (purple) that also includes the full-physics model but on a reduced grid of poloidal wavenumbers such that the only unstable modes evolved are those which are also unstable when the MHD-1 approximation is used.Removing the low-k y modes results in saturation at much smaller levels even when all of the retained modes have identical linear properties.

The MHD approximation does not accurately reproduce the large-transport state predicted by full-physics simulations
The results in Fig 4 appear to support the conclusions of [8,9] that δB ∥ fluctuations play an essential role in hKBM driven turbulence in these STEP local equilibria, and in particular in the transition to a large-transport state.However, it has been shown [8,9] that it is possible to vary the local equilibrium in such as way as to change the hKBM threshold.In Section 6 we study whether it is possible to extend this result by posing the question; can we find local equilibria where the hKBM is linearly unstable in the absence of δB ∥ with exact drifts?

The importance of parallel magnetic perturbations to the hKBM
In [9] it was shown that the total heat and particle fluxes from nonlinear simulations vary with β when β ′ is varied consistently (see Fig. 11 of [9]).However, the heat and particle fluxes remained very large even at much smaller values of β than the STEP baseline.‡ The exact saturation mechanism of the hKBM in the absence of equilibrium flow shear demands further study which will be the focus of future work.It was argued that reduced β ′ stabilisation of the hKBM is largely responsible for such large turbulent fluxes at lower β.In the STEP baseline case, β ′ stabilisation is very important since it pushes the hKBM towards marginality.The sensitive dependence of the hKBM on β ′ is made most obvious when β ′ is varied holding all other local equilibrium parameters fixed (albeit inconsistent).It is interesting to note, as remarked in Section 4.2.2, that we should expect the MHD approximation to improve as |β ′ | is reduced.
Fig 6 shows the linear growth rate (a, c) and frequency (b, d) as functions of the binormal wavenumber for simulations with reduced β ′ with, f B = 1 (blue), and without, f B = 0 (orange, green, red), δB ∥ fluctuations.In all cases the hKBM is much more unstable than at the nominal value of β ′ (see Fig. 2) due to the much decreased role of β ′ stabilisation.In each case, the hKBM is likely to drive much more transport.
There are several key observations from Fig 6 that we wish to stress.
• We note that the MHD approximation does a much better job of capturing the linear growth rates of the hKBM at smaller values of |β ′ |, in agreement with the discussion in Section 4.2.2.
• The hKBM can be unstable up to n = 1 in the absence of δB ∥ with exact drifts, and this is certainly a no-go zone for any actual device.
In more strongly driven regimes at lower β ′ in STEP, it is clearly possible to study hKBM-driven turbulence using global GK codes that do not include B ∥ .STEP (with all other parameters are identical).As expected from the linear spectrum (Fig 6), the nonlinear physics of this new case qualitatively mirrors the STEP baseline case even without δB ∥ fluctuations.The heat fluxes rise to very large levels and no robustly steady saturated state is reached over the period of this simulation: (note that this is the only simulation in this paper that does not reach a robustly steady No approximations are made concerning the magnetic-drift velocity.The lower β ′ simulation without δB ∥ (magenta) exhibits a qualitatively similar electromagnetic runaway (transition to very large heat fluxes) to that seen with full physics since the hKBM has been driven unstable at long wavelength, but it is the only simulation in this paper that does not meet the statistical definition of saturation introduced in [9].saturated state).It is currently unknown whether this simulation would obtain a steady state if it were continued and such a question is very difficult to answer numerically due to the computational cost.Simulations such as those discussed here motivate the pressing need for a rigorous theory of electromagnetic turbulence saturation in the fluxtube limit.Snapshots of the turbulent fields Fig D1 (g,h,i) reveal the presence of boxscale radially elongated streamers.
The fact that at lower β ′ we obtain electromagnetic runaway fluxes without δB ∥ is an important result; it reduces the number of ingredients needed to observe electromagnetic runaway fluxes in local gyrokinetic simulations.Put more simply, we observe qualitatively similar runaway fluxes in local simulations when there is electromagnetic instability at the longest-wavelengths in the system.We take great care to emphasise that the results in this paper suggests that there is an important distinction to be made between (A) runaway fluxes that eventually saturate at very large values such as those discussed in this paper, and (B) fluxes that never saturate such as the runaways reported by [33] for cyclone-base-case [34] geometry without δB ∥ .All nonlinear results for the STEP equilibrium at the nominal β ′ in this paper definitively fall into category (A).It remains unclear which category the lower β ′ simulation falls into.Identification of these two categories motivates the need for a rigorous theory of electromagnetic turbulence saturation that can distinguish between the two types of behaviour.

Conclusions
In this paper, we have examined the importance of including parallel magnetic perturbations (δB ∥ ) in gyrokinetic simulations of electromagnetic turbulence at midradius in the burning plasma phase of the conceptual high-β reactor-scale, tight-aspectratio tokamak STEP.It has previously been found [8] that δB ∥ is essential for the hKBM to be unstable in local equilibria taken from the STEP burning flat top.
7.1.The MHD approximation is typically not appropriate for use in high-β e reactor-scale plasmas In Section 4.2.2 it was argued that the MHD approximation § should perform well in the limit β ′ ≪ √ β ≪ 1.This limit is not well satisfied in STEP-EC-HD (where β e = 0.09 and |β ′ | = 0.48 at ρ/a = 0.64), and it may also be questionable in existing experiments at internal/external transport barriers with high pressure gradients, and in other reactor-scale plasmas.
In Section 4, it was shown that inclusion of δB ∥ and a rigorous treatment of the magnetic drift velocity is essential to capture the correct linear spectrum of the hKBM in the STEP flat-top.The MHD approximation roughly captures the linear spectrum of the hKBM, but predicts lower growth rates and is significantly more stable at long wavelengths for this STEP local equilibrium.
Even though the MHD approximation reproduces qualitatively similar linear physics, Section 5 reveals that the corresponding nonlinear simulations do not accurately track nonlinear simulations for STEP with full physics that were reported in [9].Typically, simulations which attempt to resolve the hKBM with full physics (including δB ∥ ) yield fluxes rising to very large values and require substantial computational resource to reach saturation.When the MHD approximation is used, the turbulence saturates at a level around two orders of magnitude lower than that reached in a comparable time in the full physics simulations.The reason for this difference is that the longest wavelength modes (i.e., corresponding to n < 10) are stable under the MHD approximation, but unstable with full physics.It was argued that these linearly unstable modes at very long-wavelength are responsible for the transition to very large fluxes (see Fig 5 and discussion).7.2.δB ∥ is not always essential for studying hKBM driven turbulence.However, δB ∥ is essential for studying turbulence in the high β ′ conditions of the STEP flat-top.
In Section 6, strong β ′ stabilisation [32] was identified as being responsible for the sensitivity of hKBM stability to δB ∥ fluctuations, and it was noted that this plasma sits close to marginality.It was also shown, however, that the hybrid mode becomes unstable in the absence of δB ∥ if β ′ is reduced to push the local equilibrium further from marginality.In such local equilibria, which are more unstable at long wavelength, electromagnetic runaway fluxes emerge in nonlinear simulations that neglect δB ∥ (without, and likely also with, the MHD approximation).This suggest that turbulence can be simulated in some high-β ST plasmas using presently available tools and codes that neglect δB ∥ .
In the present design of the STEP flat-top, however, it is essential to include δB ∥ because the modes expected to dominate turbulence are both: (i) unstable only with δB ∥ and; (ii) unstable up to very long wavelengths.Global gyrokinetic codes currently being developed to include δB ∥ rigorously will be important and timely for STEP.
Finally this work motivates the need for a rigorous theory of electromagnetic turbulence saturation, and an elucidation of the role of δA ∥ in particular.
which usefully relates the gradient in magnetic-field strength to the curvature, and can be expressed as: Using (C.3) to substitute for ∇ ln B in (C.1), the magnetic drift frequency can be reexpressed as:

Figure 1 :
Figure 1: Growth rate (a) and mode frequency (b) as functions of the binormal wavenumber from linear simulations of the dominant instability in STEP-EC-HD on a mid-radius flux surface.Simulations are shown both with, f B = 1 (blue), and without, f B = 0 (orange), δB ∥ .The two simulations are otherwise identical and the exact magnetic-drift velocity is used in each simulation.

Fig 1
Fig 1 shows the linear growth rate (a) and frequency (b) (normalised to the deuterium sound speed, c D = T e /m D , divided by the minor radius of the last closed flux surface) as functions of the binormal wavenumber k y ρ D = nρ * D dρ/dΨ n .Table 1 also includes the binormal wavenumber k n=1 y ρ D corresponding to the toroidal mode number n = 1.Simulations are shown both including, f B = 1 (blue), and neglecting, f B = 0 (orange), δB ∥ fluctuations (the two simulations are otherwise identical) and no approximations are made concerning the magnetic drift velocity.Importantly, we see from Fig 1 that if δB ∥ is artificially excluded from calculations (as routinely assumed in a

Figure 2 :
Figure 2: Growth rate (a) and mode frequency (b) as functions of the binormal wavenumber from linear simulations of the dominant instability in STEP-EC-HD on a mid-radius flux surface.Simulations are shown both with, f B = 1 (blue), and without, f B = 0 (orange, green, red), δB ∥ .For the simulations without δB ∥ , different treatments of the drift velocity are also shown: (i) the full drift velocity [v ds given by equation 10 (orange)]; (ii) MHD-1 [v ds given by equation 12 (green)]; and MHD-2 [v ds given by equation 13 (red)].Note that the MHD-2 curve (red) lies on top of the curve that excludes δB ∥ but uses the exact form of the magnetic-drift velocity (orange).

Fig 2
Fig 2 shows the linear growth rate (a) and frequency (b) as functions of the binormal wavenumber for simulations both including, f B = 1 (blue), and excluding, f B = 0 (orange, green, red) δB ∥ fluctuations.The simulations without δB ∥ fluctuations employ different prescriptions of the magnetic drift velocity: (i) the full drift velocity [v d given by equation 10 (orange)]; (ii) MHD-1 [v d given by equation 12 (green)]; and MHD-2 [v d given by equation 13 (red)].Using the MHD-2 approximation (Fig 2, red) we see that the hKBM is linearly stable and we once again recover the previously subdominant MTM

Fig 4
Fig 4  shows time traces of the total heat flux from two nonlinear GENE simulations.The first simulation (blue) includes δB ∥ and uses(10) for the drift velocity (this simulation is equivalent to that shown in Fig3aof[9] and will henceforth be referenced as the 'STEP baseline' simulation).The second simulation (green) shown in Fig4does not include δB ∥ and instead uses the MHD approximation (MHD-1).Although the linear physics is broadly similar (both linear spectra have hKBMs as the dominant unstable mode), the nonlinear physics produced by these two simulations is strikingly different.The simulation using the MHD approximation appears to achieve a robustly steady † saturated state at values of the heat flux around two orders of magnitude smaller than those reached by the simulation that includes full δB ∥ physics.It is interesting to note

Figure 3 :
Figure 3: Growth rates (a, c) and mode frequencies (b, d) as functions of the binormal wavenumber from linear simulations of the dominant instability in STEP-EC-HD on a mid-radius flux surface.The value of β e has been varied keeping β ′ fixed (all other local equilibrium quantities are the same as in Table 1).Results are shown for β e = 0.15 (a, b) and β e = 0.06 (c, d).Simulations are shown both with, f B = 1 (blue), and without, f B = 0 (green, red), δB ∥ .For the simulations without δB ∥ , different treatments of the drift velocity are also shown: (i) MHD-1 (green); and (ii) MHD-2 (red).The absolute error made by the MHD-1 approximation for the growth rate and the mode frequency are shown in panels (e) and (f ) respectively.The error is computed only for modes which are unstable and propagating in the ion-direction.

Figure 4 :
Figure 4: Time traces of the total heat flux (a) and total particle flux (b) from GENE simulations.Simulations are shown both with, f B = 1 (blue), and without, f B = 0 (green), δB ∥ .The case without δB ∥ fluctuations uses the MHD approximation (MHD-1) with v d given by equation 12. Also shown is the time trace of the k y spectrum of δϕ (c, e) and δA ∥ (d, f ) from the same simulations.Only the first decade of poloidal wavenumbers is shown.
[9]), showing the strong effect of the hKBM instability when the MHD approximation is used to drop δB ∥ .Snapshots of the turbulence are shown in Appendix D, Fig D1.The considerably reduced transport fluxes with the MHD approximation, compared to with full physics, are associated with lower amplitude turbulent fluctuations shown in Fig.4, and less radially extended turbulent structures in the perturbed fields (see Fig.D1(a-f)).The reason for saturation at much lower fluxes under the MHD approximation in Fig 4 is revealed by closer inspection of the MHD-1 linear results in Fig 2a.

Figure 5 :
Figure 5: Time traces of the total heat flux (a) and total particle flux (b) from GENE simulations.Both simulations include δB ∥ and use the exact prescription of the magnetic-drift velocity.The value n min is the toroidal mode number corresponding to the value of k y,min .The simulation with the larger value of n min only evolves those unstable modes that are also unstable when the MHD-1 approximation is used.Removing the lowest k y modes (purple) in this way results in saturation at values two orders of magnitude lower than the baseline (blue).

Figure 6 :
Figure 6: Growth rates (a, c) and mode frequencies (b, d) as functions of the binormal wavenumber from linear simulations of the dominant instability in STEP-EC-HD on a mid-radius flux surface.Results are shown with reduced β ′ = 0.4β ′ STEP (a, b) and β ′ = 0.1β ′ STEP (c, d).All other parameters are kept fixed.Simulations are shown both with, f B = 1 (blue), and without, f B = 0 (orange, green, red), δB ∥ fluctuations.For the simulations without δB ∥ , different treatments of the drift velocity are also shown: (i) MHD-1 (green); and (ii) MHD-2 (red).In all cases, the hKBM is much more unstable than the reference case (see Fig 2) due to reduced β ′ stabilisation [32].

Fig 7
Fig 7  shows time traces of the total heat flux and particle flux from two nonlinear GENE simulations: (i) with the nominal value of β ′ that includes δB ∥ fluctuations (blue); and (ii) excluding δB ∥ fluctuations and using exact drifts (magenta) at a lower value of β ′ = 0.4β ′ STEP (with all other parameters are identical).As expected from the linear spectrum (Fig6), the nonlinear physics of this new case qualitatively mirrors the STEP baseline case even without δB ∥ fluctuations.The heat fluxes rise to very large levels and no robustly steady saturated state is reached over the period of this simulation: (note that this is the only simulation in this paper that does not reach a robustly steady

Figure 7 :
Figure 7: Time traces of the total heat flux (a) and total particle flux (b) from GENE simulations.Simulations are shown both with δB ∥ at the nominal β ′ (f B = 1, blue), and without at β ′ = 0.4β ′ STEP (f B = 0, magenta).No approximations are made concerning the magnetic-drift velocity.The lower β ′ simulation without δB ∥ (magenta) exhibits a qualitatively similar electromagnetic runaway (transition to very large heat fluxes) to that seen with full physics since the hKBM has been driven unstable at long wavelength, but it is the only simulation in this paper that does not meet the statistical definition of saturation introduced in[9].

Table 1 :
Local equilibrium quantities at the STEP mid-radius flux surface considered in this paper.
subdominant mode simply by artificially suppressing δB ∥ (thus stabilising the hKBM).This latter point is exemplified in Fig 1.