On electromagnetic turbulence and transport in STEP

In this work, we present first-of-their-kind nonlinear local gyrokinetic (GK) simulations of electromagnetic turbulence at mid-radius in the burning plasma phase of the conceptual high-β , reactor-scale, tight-aspect-ratio tokamak Spherical Tokamak for Energy Production (STEP). A prior linear analysis in Kennedy et al (2023 Nucl. Fusion 63 126061) reveals the presence of unstable hybrid kinetic ballooning modes (KBMs), where inclusion of the compressional magnetic field fluctuation, δ B ∥ , is crucial, and subdominant microtearing modes (MTMs) are found at binormal scales approaching the ion-Larmor radius. Local nonlinear GK simulations on the selected surface in the central core region suggest that hybrid KBMs can drive large turbulent transport, and that there is negligible turbulent transport from subdominant MTMs when hybrid KBMs are artificially suppressed (through the omission of δ B ∥ ). Nonlinear simulations that include perpendicular equilibrium flow shear can saturate at lower fluxes that are more consistent with the available sources in STEP. This analysis suggests that hybrid KBMs could play an important role in setting the turbulent transport in STEP, and possible mechanisms to mitigate turbulent transport are discussed. Increasing the safety factor or the pressure gradient strongly reduces turbulent transport from hybrid KBMs in the cases considered here. Challenges of simulating electromagnetic turbulence in this high-β regime are highlighted. In particular the observation of radially extended turbulent structures in the absence of equilibrium flow shear motivates future advanced global GK simulations that include δ B ∥ .


Introduction
Understanding and predicting turbulence in the core of next-generation spherical tokamaks (STs) is critical for the optimisation of their performance.In order to be economically competitive, ST power plants require a high ratio, β, of thermal pressure to magnetic pressure in the tokamak core.At high β, turbulence is expected to become more electromagnetic in nature, and to be influenced by kinetic ballooning modes (KBMs) and microtearing modes (MTMs), as frequently reported for STs (see [1] and references therein).Whilst gyrokinetic (GK) simulations have thus far proven to be a reliable tool for modelling turbulent transport in predominantly electrostatic regimes at low β in conventional aspect ratio tokamaks, obtaining saturated nonlinear simulations of plasmas at higher β with unstable KBMs and MTMs has proved more challenging [2,3].
The UK Spherical Tokamak for Energy Production (STEP) programme is developing a compact prototype power plant design based on a high-β ST, with the primary goal of generating net electric power P el > 100 MW from fusion [4][5][6][7].The first phase of this ambitious programme is to provide a conceptual design of the STEP prototype plant and a reference equilibrium for the preferred flat-top operating point.Various reference plasma scenarios for STEP [5,8] have been developed using the integrated modelling suite JINTRAC [9], with plasma transport calculated using the Bohm-gyro-Bohm model [10] with transport coefficients rescaled (i) to give an energy confinement time enhancement factor H 98 = 1.3, over the ITER98,y2 [11] confinement scaling law, and (ii) to give dominant electron heat transport (as frequently observed in STs [1]).
We will focus here on 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.The global parameters for this plasma flat-top operating point are reported in [12], while details on the JINTRAC modelling of the chosen STEP flat-top operating point are reported in [5,8].At this STEP flat-top operating point, the heat source is mainly from fusion α-particles and from electron-cyclotron heating, while the particle source is from pellet injection that provides 7.4 ×10 21 particles/s §.
This paper seeks to progress the important objective of testing the transport and confinement assumptions used to develop STEP scenarios against first-principlesbased local GK simulations; these simulations are expected to be considerably more challenging than for conventional aspect ratio tokamaks at low β, because of the expected electromagnetic nature of the turbulence in high-β STs.The main goals of this work are (a) to model plasma turbulence in a STEP flat-top operating point with the best currently available tools; (b) to assess the compatibility of the predicted transport with the anticipated sources in STEP; (c) to explore strategies for optimising transport in the conceptual STEP design; and (d) to identify modelling challenges and suggest future directions for tackling them.‡ More specifically we use STEP-EC-HD=v5 with SimDB UUID: 2bb77572-d832-11ec-b2e3-679f5f37cafe, Alias: smars/jetto/step/88888/apr2922/seq-1 § In modelling this flat-top operating point the particle confinement time τ p ∼ 4.5τ E is assumed consistent with typical findings in JET.Assuming the same confinement time for helium ash gives a saturated helium abundance of about 9%.Higher helium confinement would degrade fusion performance due to the core accumulation of helium ash.
The nonlinear simulations presented here are carried out using the local gyrokinetic code CGYRO [13] (commit 399deb4c).All simulations include three species (electrons, deuterium and tritium) and neglect the impacts of impurities and fast α particles, which will be assessed in future work.
The paper is organised as follows.Firstly the reference local equilibrium used throughout is introduced in §2 and its local microstability properties are summarised: this is the q = 3.5 surface close to mid-radius in STEP-EC-HD.Then §3 presents fully electromagnetic (i.e.including δB ∥ ) local nonlinear gyrokinetic simulations of the hybrid-KBM dominated turbulence, and simulations are performed for various assumed values of the perpendicular equilibrium flow shearing rate.In §4, we explore possible mechanisms to mitigate large turbulent fluxes that are driven by the hybrid-KBMs and we assess sensitivity of the turbulent transport to local pressure gradient, β and safety factor, in order to inform the design of future STEP flat-top operating points.Simulations of isolated turbulence from the subdominant MTM instability are performed in §5 in a situation where hybrid-KBMs are artificially suppressed by removing δB ∥ .Finally, a summary of findings is presented in §6.

Local equilibrium and microstability at q = 3.5 in STEP-EC-HD
The local gyrokinetic analysis in this paper is carried out at the q = 3.5 surface of the STEP-EC-HD flat-top equilibrium.This surface is highlighted in red in Figure 1, which also shows radial profiles of the three plasma species included in the gyrokinetic simulations; electrons, deuterium and tritium∥.A Miller parameterisation [14] was used to model the local magnetic equilibrium geometry, and the shaping parameters were fitted to the chosen surface using pyrokinetics [15], which is a python library developed to facilitate pre-and post-processing of gyrokinetic analysis.Table 1 provides the values of various local equilibrium quantities on this surface, including: magnetic shear, ŝ; normalised minor radius, ρ/a; elongation and its radial derivative, κ and κ ′ ; triangularity and its radial derivative, δ and δ ′ ; radial derivative of the Shafranov shift, ∆ ′ ; temperature T and density n for electrons, deuterium and tritium, and their normalised inverse gradient scale lengths.
Local linear GK analysis on this q = 3.5 surface of STEP-EC-HD is reported in [12], and reveals the presence, at ion Larmor radius scales, of a dominant hybrid-KBM instability that has significant linear drive contributions from trapped electrons and electromagnetic terms.A subdominant MTM is also unstable at ion binormal scales, and no unstable mode is observed at electron Larmor radius scales.Figure 2 shows the linear growth rate, γ, and mode frequency, ω, of the dominant and sub-dominant modes as functions of the normalised binormal wave-number k y ρ s , where ρ s denotes the deuterium Larmor radius (evaluated using the electron temperature) ρ s = c s /Ω D , ∥ The source reference equilibrium includes impurities and helium, but our GK calculations neglect these species and scale the deuterium and tritium density by a factor close to unity to restore quasineutrality.c s = T e /m D is the deuterium sound speed, Ω D = eB 0 /m D is the deuterium cyclotron frequency, m D is the deuterium mass, and B 0 is the toroidal magnetic field at the centre of the flux surface.Full details on the linear physics are reported in [12], where the linear properties of the hybrid-KBM are characterised in detail.The work in [12] points out the critical role played by compressional magnetic fluctuations, δB ∥ .If δB ∥ is artificially excluded from calculations the hybrid-KBM is entirely stabilised, while the linear properties of the MTM are unaffected.Here we exploit this in order to model the turbulent transport that might be expected to be driven by these two classes of mode if they could be modelled independently.
Whilst in this paper the analysis is restricted to q = 3.5, we expect general results and trends to be relevant across a range of surfaces in STEP-EC-HD, which have been shown in [12] to share qualitatively similar linear microstability properties.Further details on the local equilibrium are reported in [12].

Nonlinear simulations including δB ∥ : hybrid-KBM turbulence
The focus here is to analyse turbulence driven by hybrid-KBMs.This requires fully electromagnetic nonlinear gyrokinetic simulations that include perturbations in electrostatic potential δϕ, parallel magnetic vector potential δA ∥ , and parallel magnetic field δB ∥ .
Table 1: Local parameters at q = 3.5 in the STEP-EC-HD flat-top.Parameters not defined in the text include: β e = 2µ 0 n e T e /B 2 0 ; minor radius, a[m] = 2, is the halfdiameter at the equatorial midplane of the last closed flux-surface; local minor radius, r; ν ee is the electron collision frequency; X ′ = dX dρ denotes the derivative of quantity X with respect to the flux label ρ = r/a; major radius R; surface area A surf ; total heating power crossing the surface P surf ; density n s and temperature T s for species s; and the normalised inverse gradient scale length for quantity X, a/L X = a X dX dρ .Figure 2: Growth rate (a) and mode frequency (b) as a function of k y ρ s from linear simulations of the dominant (hybrid-KBM) and subdominant (MTM) instabilities of STEP-EC-HD at the radial surface corresponding to q = 3.5 (see [12] for details on the linear analysis).Stable modes are shown as open markers.

Numerical resolution
The CGYRO simulations have generally been performed using grid sizes ¶ in the parallel, binormal, radial, energy and pitch-angle dimensions, (n θ , n ky , n kx , n v , n ξ ) = (32,96,128,8,32), and using the Sugama collision model [16].The binormal k y grid spans 0 ≤ k y ρ s ≤ 0.95, so includes the full range of linear instability and linearly stable grid-points at both high and low k y (see Figure 2).The radial box size L x = j/(ŝk y,min ) is set to accommodate exactly eight rational surfaces of the lowest finite toroidal mode number by choosing the integer j = 8, motivated by the numerical convergence study + reported in Appendix A. Local nonlinear simulations for STEP flat-top operating points on the above grids typically require 5 × 10 5 CPU-hours on the ARCHER2 high performance computing system; i.e. they are computationally expensive.

Flow shear considerations
Sensitivity of the linear growth rate to the ballooning parameter θ 0 = k x /(ŝk y ) is a reliable indicator of a mode's susceptibility to stabilisation by equilibrium flow shear.Linear analysis in [12] shows that the growth rate of the hybrid-KBM instability is highly sensitive to θ 0 and that the mode is stable for all θ 0 above a very small value.Hybrid-KBM turbulence should therefore be very sensitive to flow shear.Lacking any external momentum source from neutral beam injection in STEP, rotation is expected to be modest.Equilibrium flow shear will be present, however, at the diamagnetic level.We will assess the sensitivity of hybrid-KBM turbulence to the equilibrium flow shearing rate γ E , and estimate the diamagnetic flow shearing rate using (see e.g.[18]): where η i = L n /L T i is the ratio of the gradient scale lengths associated with density and ion temperature (which is equal for deuterium and tritium), and p i is the total ion pressure, and n i = n D + n T .From equation (1) the value of γ dia E at Ψ n = 0.49 is approximately γ dia E ≃ 0.05 c s /a.Details on the equilibrium flow shear implementation in CGYRO are reported in [19].Our main purpose is to assess challenges and uncertainties impacting on gyrokinetic calculations of turbulent transport at the STEP reference operating point; one crucial element is clearly to understand the sensitivity of turbulent fluxes to equilibrium flow shear.¶ CGYRO discretises velocity space using Laguerre and Legendre spectral methods in v and ξ, respectively (see [13] for details), where ξ = v ∥ /v with v ∥ the velocity component parallel to the equilibrium magnetic field direction.+ Small values of ∆k x were also required in local GK simulations of ITG turbulence in MAST in the presence of equilibrium flow shear.While the radial domain of the flux-tube was larger than MAST minor radius, the simulated turbulence properties were nevertheless found to be in remarkably close agreement with fluctuation measurements from Beam Emission Spectroscopy [17].

Hybrid-KBM turbulence results
We have performed a set of nonlinear simulations using the default grid, for γ E ∈ {0, 0.05, 0.1, 0.2} c s /a, where equilibrium flow shear is turned on well into the linear growing phase.
Previous works (see e.g.[20,21]) show that simulations of electromagnetic gyrokinetic turbulence do not always reach a saturated state.In finite simulations, however, it can be difficult to distinguish between turbulent fluxes that will never saturate and turbulent fluxes that might eventually saturate potentially at a very large value.In order to more clearly distinguish between simulations with runaway fluxes and those that saturate within the time simulated, this work introduces a rigorous definition of a saturated state using the augmented Dickey-Fuller test [22]; this tests the turbulent fluxes for statistical stationarity.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 typically chosen between 0.05 and 0.1.In this work, we apply the augmented Dickey-Fuller test to the time trace of the total heat flux and we define the time series to be stationary (i.e. the simulation reaches a robustly steady saturated state) when the p-value is below 0.1.Every simulation in this paper achieves a saturated state according to this criterion, and wherever time-averaged statistics are reported they are averaged over this saturated state.
Figure 3 shows the time averaged heat and particle flux values from nonlinear simulations with various flow shearing rates, γ E .Heat and particle fluxes both decrease with γ E .We note that the electromagnetic electron heat flux, Q e,em , dominates the total heat flux at all shearing rates considered.The second largest contribution comes from the electrostatic ion heat flux Q i,es , which is more than a factor of five smaller than Q e,em .The electrostatic and electromagnetic particle flux values are similar.At γ E = 0.05 c s /a ∼ γ dia E , heat and the particle fluxes exceed the maximum steadystate fluxes that could be provided from the sources assumed to be available in STEP (shown as a horizontal green line in Figure 3).Increasing the equilibrium flow shear to γ E = 0.1 c s /a brings the heat flux below the available source, though the particle flux remains above.We note that, unlike the heat source, fuelling can be more easily increased without compromising fusion gain.If the equilibrium flow shear is increased further the particle flux drops below the maximum available source.The same qualitative and quantitative trends of heat and particle fluxes with γ E are also observed in GENE simulations (see Appendix B for a comparison with GENE results).

Turbulence with γ E = 0
It is important to first discuss simulations without equilibrium flow shear, before shifting our focus to more physically relevant simulations with γ E > 0. The horizontal dashed lines indicate maximum fluxes that could be consistent with the reference STEP equilibrium, and are taken as the total heat and particle sources, respectively, divided by the area of the Ψ n = 0.49 surface.Open markers are used at γ E = 0 where the turbulent flux bursts (see Figure 4) influence the time averaged values.turbulent transport in these simulations) from the nonlinear simulation with γ E = 0.The lowest k y modes have the largest amplitudes in eδϕ/(ρ * T e ) and δA ∥ /(ρ * ρ s B 0 ), despite these low k y modes being only marginally unstable with the lowest k y > 0 mode linearly stable.After the initial linear growth phase, the heat flux contributions from various modes exceed 10 Q gB .Between t ≃ 150 a/c s and t ≃ 300 a/c s , both fields and fluxes at k y ρ s > 0.05 appear to saturate, though the amplitudes of the k y ρ s ≤ 0.05 modes increase weakly.The amplitude of δA ∥ for low k y modes exceeds the zonal δA ∥ amplitude at t ≃ 350 a/c s (see Figure 4), which coincides with a first heat flux burst.After the heat flux burst, the heat flux returns to values similar to the one observed in the time window between t ≃ 150 a/c s and t ≃ 300 a/c s .We note however that δϕ and δA ∥ at low k y modes (k y ρ s < 0.05) keep increasing slowly in this phase.At t ≃ 620 a/c s , the δϕ and δA ∥ amplitude of the non-zonal low k y modes exceeds the zonal amplitude, and a second heat flux burst occurs.Apart from the heat flux bursts, the heat flux value remains constant over a quite long time window for t > 150 a/c s .The augmented Dickey-Fuller test applied to the heat flux over t > 150 a/c s return a p-value of 0.01, thus supporting the stationarity of the heat flux time trace.We note that continuing this simulation further is extremely expensive due to the very small time step and it does not change the main outcome: turbulent fluxes are three orders of magnitude larger than values compatible with the available STEP sources in the absence of sheared flows.(This conclusion is also supported by similar results from GENE simulations that are reported in Appendix B). Figure 4 (d) shows the k y spectrum of the heat flux averaged over time from t = 150 a/c s to the end of the simulation.The heat flux spectrum peaks at k y ρ s ≃ 0.02 and decays at large k y .Figure 5 shows snapshot contour plots of δϕ and δA ∥ at the outboard mid-plane (θ = 0), taken from the last time step.The fluctuation amplitudes of eδϕ/(ρ * T e ) and δA ∥ /(ρ * ρ s B 0 ) are comparable, demonstrating the strongly electromagnetic nature of the turbulence * .The turbulent structures are highly elongated radially, extending to ∼ 100 ρ s , which corresponds here to approximately 50 cm.Similar large flux states have been found using various radial box sizes and radial resolutions as well as with different GK codes (see Appendix A for details), suggesting that for this local equilibrium large fluxes are a robust prediction of local gyrokinetic simulations.Given the large radial extent of the turbulent structures, however, global effects not captured in the local approach may also affect saturation; e.g.background equilibrium profile variation effects, or another recently proposed global mechanism whereby global zonal flows are enhanced by the radial transfer of entropy [23].Global gyrokinetic simulations that include δB ∥ may more accurately predict turbulent fluxes in this regime, and constitute a possible future research direction.We note, however, that local nonlinear simulations with modest flow shear produce less radially extended structures (e.g.see snapshots of δϕ and δA ∥ in Figure 7), and that further actuators and mechanisms (see §4) also result in turbulent states with lower transport that are more amenable to local gyrokinetics.

Turbulence with finite γ E
We now focus on the simulations with finite equilibrium flow shear.Figure 6 shows contributions to the heat flux from the lowest 16 k y modes as functions of time for the simulations with γ E = 0.05 c s /a and γ E = 0.1 c s /a, where equilibrium flow shear is turned on at t = 0 with the initial condition taken from after the linear growth phase in a simulation with γ E = 0.The heat flux decreases in the presence of equilibrium flow shear, with a sharp initial reduction followed by a phase of weaker decay.The two simulations take until approximately t ≃ 1000 a/c s to reach a saturated state.The k y spectrum of the time averaged heat flux in the saturated phase is shown in Figure 6 (c) and (d).At γ E = 0.05 c s /a, the total heat flux peaks around k y ρ s = 0.03 and decays for k y ρ s > 0.1 to a small tail contribution to the flux from high k y ρ s .At γ E = 0.1 c s /a, the spectrum is broader, with the main peak shifting to higher k y and a second distinct lower amplitude peak at k y ρ s ≃ 0.45.In both cases, low k y modes make a relatively large contribution to the heat flux, although the peak of the spectrum is broader and at larger k y compared to Figure 4 (d) for γ E = 0. Equilibrium flow shear clearly acts to suppress the formation of large amplitude radially extended turbulent structures.
The simulations with finite γ E are characterised by relatively strong zonal flows ( kx δϕ(k x , k y = 0)) and zonal fields ( kx δA ∥ (k x , k y = 0)), which reduce considerably  the radial extent of the turbulent structures, as shown in Figure 7 and 8.The amplitudes of eδϕ/(ρ * T e ) and δA ∥ /(ρ * ρ s B 0 ) are again comparable, as expected from electromagnetic turbulence.We highlight that the δA ∥ /(ρ * ρ s B 0 ) amplitude is smaller at γ E ≃ 0.1 c s /a than at γ E ≃ 0.05 c s /a.In addition, the radial extents of structures in δϕ and δA ∥ at γ E ≃ 0.1 c s /a are significantly smaller at ∼ 30ρ s (more easily visible if the zonal component is removed).The analysis presented in this section suggests that, while equilibrium flow shear reduces the turbulent fluxes, at the diamagnetic level it is insufficient to reduce the hybrid-KBM heat and particle fluxes to levels compatible with the assumed sources on the chosen surface from the STEP operating point.This motivates an assessment of the sensitivity of hybrid-KBM turbulent transport to key local equilibrium parameters, to  seek routes for optimising transport at the STEP flat-top operating point.

Sensitivity studies and mitigation mechanisms
This section is divided in two parts.In the first part, we assess the sensitivity of turbulent fluxes to the pressure gradient that drives the hybrid-KBM instability in nonlinear simulations with equilibrium flow shear.In the second part, we explore additional possible mechanisms other than equilibrium flow shear to mitigate the impact of hybrid-KBM.This includes some sensitivity scans to β e and q, with the aim of identifying more favorable regimes for a STEP flat-top operating point.
4.1.Sensitivity on local pressure gradient in simulations with γ E = 0.05 c s /a Section §3.3 shows that turbulent fluxes can exceed the target value for this STEP flattop operating point even in the presence of equilibrium flow shear.Here we explore the sensitivity of hybrid-KBM fluxes to the pressure gradient drive to assess whether target transport fluxes can be achieved through a modest change to the pressure gradient.
We perform nonlinear simulations with three different values of the equilibrium pressure gradient, L p,ref /L p ∈ {0.8, 1, 1.2}, where a/L p,ref denotes the reference value of the local pressure gradient.In this scan the pressure gradient is varied by scaling both the electron and ion density and temperature gradients simultaneously by the same factor.All other local equilibrium parameters are kept constant in this scan, except for β ′ ≡ β∂ ln p/∂ρ, which is scaled consistently with the pressure gradient to retain the important effect of β ′ on the linear hybrid-KBM instability [12].All simulations include equilibrium flow shear with a constant value of γ E = 0.05 c s /a.The numerical resolution of the simulation at a/L p = 1.2(a/L p ) ref is identical to that used in the nominal case, while the simulation at a/L p = 0.8(a/L p ) ref requires a higher k y,max ρ s = 1.28 to resolve the linear spectrum (discussed later).
Figure 9 shows the saturated heat and particle fluxes for the different values of a/L p .We observe a strong dependence of turbulent fluxes on the pressure gradient.The saturated heat flux at L p,ref /L p = 1.2 drops below the target flux.Whilst the electron electromagnetic heat flux dominates at all pressure gradients, the relative contribution from the electrostatic channel is significantly higher at L p,ref /L p = 1.2 compared to the other two cases.This sensitivity study to pressure gradient (and β ′ ) highlights the beneficial stabilizing effect of increasing β ′ on hybrid-KBM turbulent transport.Furthermore we note that the beneficial effect on transport at higher a/L p would be further reinforced if equilibrium flow shear were scaled proportionally with the ion pressure gradient as would be expected for diamagnetic flows, but γ E was held fixed in this scan.Decreasing the pressure gradient has the opposite effect and leads to much higher turbulent fluxes.Local equilibrium scans, like the a/L p scan presented here, reveal valuable insights, but in a truly consistent scan the whole global MHD equilibrium would evolve such that many local equilibrium parameters, like q, ŝ and β (and potentially the turbulent transport regime) might also change significantly.
Figure 9 is consistent with the strong, and somewhat counter-intuitive, dependence of the hybrid-KBM linear growth rate on a/L p (and β ′ ) that is illustrated in Figure 10.The linear analysis of [12] shows that the hybrid-KBM growth rate increases with density and temperature gradients when all the other parameters are held constant.On the other hand, we find here that the boost to the linear drive from the pressure gradient increase is more than compensated by increased stabilisation from the impact on local equilibrium geometry of higher β ′ (e.g. through the modification of magnetic drifts [24,25]).The growth rates at L p,ref /L p = 0.8 are larger than in the nominal case, except at low k y where they are comparable, while the growth rates in the L p,ref /L p = 1.2 case are smaller than in the nominal case for k y ρ s < 0.3, and they are comparable for 0.8 0.9 1.0

Possible mechanisms to mitigate hybrid-KBMs
We explore here the impact of local equilibrium parameters, β e and q, on hybrid-KBM transport fluxes.Unless it is explicitly stated otherwise, only the scan parameter is varied and other local parameters are held constant to isolate the effect of each parameter 0.000 0.025 0.050 0.075 0.100 0.125 0.150 0.175  and shed light on physical mechanisms that may mitigate the turbulent transport.In all the scans in this section we set γ E = 0.

Sensitivity to β e
We focus first on the effect of β.Nonlinear simulations at different values of β e are performed varying β ′ = β e /L p consistently (with the other local parameters kept constant).Numerical resolutions are chosen to properly resolve the linear spectrum of the dominant instability♯.
Figure 11 shows the total heat and particle fluxes from nonlinear simulations at different values of β e .At β e < 0.025, the particle and heat transport is mainly electrostatic and is driven by an ITG/TEM instability.As pointed out by [12], the hybrid-KBM couples to an ITG/TEM instability at low β.The electrostatic heat and particle fluxes remain above 100 MW/m 2 and 10 22 particles/(m 2 s) even at β e = β e,ref /4, and presumably reduced β ′ stabilisation plays a significant role in maintaining such large turbulent fluxes.In contrast to the low k y modes at nominal β, the amplitude of δA ∥ /(ρ * ρ s B 0 ) is smaller than the eδϕ/(ρ * T e ) one at lower β, which is consistent with a smaller electromagnetic drive.In particular, Figure 12 shows that the heat flux contribution from k y ρ s ≤ 0.1 is significantly smaller than the contribution from k y ρ s > 0.1 at lower β e , especially at β e = 0.005.
At β e = 0.16, above the nominal value, electromagnetic contributions dominate transport and the total heat flux drops considerably from its nominal value to Q tot ≃ 2 MW/m 2 , which is close to the assumed transport flux for this STEP flat-top operating ♯ In particular, the simulations with β e = 0.02 and β e = 0.005 are carried out with n ky = 128 and n kx = 64, while the simulation with β e = 0.16 is carried out with n ky = 32 and n kx = 256.point.This is a consequence of increased β ′ stabilisation of the hybrid-KBM, which affects all the unstable k y modes.The turbulent transport is dominated by MTMs, which emerge as the dominant linear instability at large β e , as detailed in [12].This local equilibrium at β e = 0.16 has considerably more favourable transport than the reference case at β e = 0.09, if it could be accessed through a route passing through lower β that avoided prohibitive transport losses from hybrid-KBMs.However, as noted in [12], a global equilibrium with this local β e would likely exceed the limiting β for effective control of Resistive Wall Modes.
Figure 13 suggests that the mitigation of hybrid-KBMs as β e may be related to the relative distance between the local equilibrium point on an ŝ − α diagram and the ideal n = ∞ ballooning boundary (where α = −Rq 2 β ′ ).The equilibrium points at β e = 0.02 and β e = 0.16 are further from the stability boundary than the reference equilibrium at β e = 0.09, potentially indicating a reduced electromagnetic contribution to the hybrid-KBM drive, and a subsequent reduction of heat and particle fluxes driven by the this component of the turbulence.

Sensitivity to q
The hypothesis that the n = ∞ ideal ballooning boundary can be used an indicator for the strength of the electromagnetic drive for the hybrid-KBM motivates further study of sensitivity to the safety factor.Figure 14 shows the heat and particle fluxes from a set of nonlinear simulations with different values of q.The heat and particle fluxes drop considerably when q is increased and approach the target values at q = 4.5.In all cases the heat flux is dominated by the electromagnetic electron contribution.We highlight that the simulations with q = 4.0 and q = 4.5 show a robust saturation of fluxes as well as turbulent structures that are considerably less radially extended than in the nominal case.For example, Figure 15 shows the time trace of the electrostatic  and electromagnetic heat fluxes from the simulation with q = 4.5 and snapshot contour plots of δϕ and δA ∥ from the outboard mid-plane at the end of the simulation.We also note in the simulation with q = 4.5 that the amplitude of δA ∥ /(ρ * ρ s B 0 ) is smaller than the amplitude of eδϕ/(ρ * T e ), thus suggesting a somewhat reduced importance of electromagnetic turbulence at higher q.
The nonlinear suppression of turbulent fluxes at higher q illustrated of Figure 14 is consistent with the strong reduction in linear growth rate values shown in Figure 16.The linear growth rate spectrum shows that a clear separation between modes at k y ρ s < 0.2 and k y ρ s > 0.2 appears as q is increased, with a frequency sign change in the low k y region with respect to the nominal case.This may suggest that the coupling with the TEM evolves with q, as also mentioned in [12].The mode frequency at k y ρ s > 0.2 is largely unaffected by the q increase.We highlight that, despite the strong reduction of the hybrid-KBM growth rate values at higher q, hybrid-KBM remains the dominant instability.In fact, the dominant mode at high q has a twisting parity, the particle flux saturates at a non negligible value, and the eδϕ/(ρ * T e ) fluctuation amplitude is larger than the δA ∥ /(ρ * ρ s B 0 ) one, suggesting a the presence of an additional electrostatic turbulence drive (see Fig. 15).We note that although hybrid-KBMs are linearly very weakly unstable at q = 4.5, they can still drive non-negligible turbulent transport nonlinearly (in the absence of equilibrium flow shear).Similarly to the β e stabilisation, the beneficial effect of the safety factor may be related to an increasing distance at higher q between the equilibrium point and the ideal ballooning stability boundary in the ŝ−α plane, as shown in Figure 13 (b).In this q scan the increase in this distance at higher q has a significant contribution from the movement of the boundary.Before concluding this section, we highlight that our analysis suggests a possible strategy to mitigate transport from hybrid-KBMs when designing STEP operating points: increasing the distance between the ideal ballooning stability boundary and the equilibrium point in the ŝ − α diagram.In all cases analysed here we find that turbulent transport is considerably reduced when the equilibrium is pushed further from the ideal ballooning limit in the direction that suppresses hybrid-KBMs in favour of other instabilities, such as TEMs or MTMs.In an important future work, we will perform a broader study of the saturation mechanisms at play in hybrid-KBM turbulence, including for local equilibria further from the ideal ballooning boundary.

The role of the subdominant MTM instability
The subdominant MTM instability, linearly characterised in [12], may become important when the hybrid-KBM is suppressed and it was seen to play an interesting role as β e was increased in §4.2.1.It is therefore important to evaluate turbulent transport from the subdominant MTM instability, and this is computationally tractable if we can isolate unstable MTMs from the dominant hybrid-KBM.As shown in [12], the hybrid-KBM instability can be artificially suppressed without any impact on MTMs linearly, by simply removing δB ∥ from the GK system of equations (so that only δϕ and δA ∥ are evolved); this change leaves MTMs as the dominant and only instability in this reduced system.Artificial suppression of the hybrid-KBM by removing δB ∥ therefore provides a convenient way to study isolated MTM turbulence in the STEP operating point.While it has been demonstrated that exclusion of δB ∥ has a negligible impact on the linear physics of MTMs in these STEP plasmas, the calculations of isolated MTM turbulence presented here also exclude potential nonlinear impacts associated with δB ∥ that may or may not be significant.
We have adopted this approach to study MTM-driven turbulence on the reference surface in STEP-EC-HD.Whilst it is unphysical to neglect δB ∥ , this study is nevertheless of considerable interest because (i) MTMs have been demonstrated to be unstable and even dominant over an extended range of binormal scales in various conceptual designs of fusion power plants based on the ST [3,6,26], (ii) MTMs may be important where the hybrid-KBMs is suppressed (e.g. at high β ′ as considered in §4.2.1), or on other radial surfaces in STEP-EC-HD, and (iii) MTMs have been demonstrated to have significant impacts on transport in other devices [27][28][29][30][31][32] and are consistent with pedestal magnetic fluctuation measurements in several experiments [33][34][35][36][37][38][39][40][41].Modelling low k y MTMs with conventional gyrokinetic codes is challenging due to the multiscale nature of the mode resulting in large grids and corresponding large computational cost: δϕ is highly extended in the field-line-following coordinate θ, which implies high radial wavenumbers, while δA ∥ is localised in θ and radially extended.Due to the very high computational cost, there are only a limited number of studies of MTM-driven turbulence in the literature (see e.g.[28,29,32]).
In our local GK simulations of MTM turbulence in STEP we find that the following numerical resolutions are sufficient: n θ = 32; n kx = 512; n ξ = 32; n v = 8; ∆k y ρ s = 0.02; ∆k x = 0.04; and 16 toroidal modes as MTMs are unstable only in a sub-region of the k y interval where the hybrid-KBM is unstable.Aside from the suppression of δB ∥ and the refined numerical grid, the simulation parameters are identical to those given in Table 1.The MTM simulation is performed with γ E = 0.
Figure 17 shows the time trace of the electrostatic and electromagnetic heat and particle fluxes from this nonlinear simulation, which ran for approximately fifty growth times of the most unstable k y (we note that the linear growth rate of the MTM is typically an order of magnitude smaller than that of the hybrid-KBM as shown in Figure 2).A robustly-steady saturated state is reached after t ≃ 3000 a/c s , where both the heat and particle fluxes are found to saturate at negligible levels †. Figure 17 also shows the time trace of the k y spectrum of δϕ and δA ∥ .It can be seen that the zonal mode in both δϕ and δA ∥ dominates over the non-zonal modes.In these simulations zonal fields [21,32], and local temperature flattening [42], are found to play a role in saturating this MTM instability at negligible heat flux values (see Appendix C).We mention that heat and particle flux predictions from the CGYRO simulation without † An appropriate parallel dissipation scheme was found to be essential to avoid runaway MTM fluxes [3] due to the onset of unphysical numerical instabilities with grid-scale oscillations in the parallel direction.This finding prompted an improvement to the parallel dissipation scheme in CGYRO, implemented at commit 399deb4c, that was used in all of the CGYRO simulations in this paper.δB ∥ are in good agreement with analogous simulations carried out with GENE and GS2 (see Appendix B).The dominance of zonal flows and fields is clearly observed in Figure 18, which shows contour plots of δϕ and δA ∥ at the outboard mid-plane taken from the last time point in the simulation.Turbulent eddies are strongly suppressed due to the presence of zonal flows and, consequently, turbulent transport is extremely modest.
It is interesting to note that although the linear growth rates of the subdominant MTM instability (see Figure 2) are similar to the linear growth rates of the hybrid-KBM instability at higher safety factor (see Figure 16), the turbulent fluxes driven by hybrid-KBMs are much larger than those driven by MTMs.This further highlights the importance of strategies to optimise the STEP equilibrium to suppress the hybrid-KBM in favour of other microinstabilities, perhaps using proximity to the ideal ballooning limit as a guide, as discussed in §4.2.

Conclusions
This paper presents the first nonlinear gyrokinetic simulations of turbulent transport in the core plasma of the conceptual STEP reference flat-top operating point, STEP-EC-HD, based on local nonlinear simulations at the q = 3.5 surface at Ψ n = 0.49.The linear microstability analysis reported in [12] has revealed that all unstable modes are strongly electromagnetic and at ion binormal scales (i.e.k y ρ s = O(1)), the fastest growing modes are hybrid-KBMs, and MTMs are subdominant.The turbulence is found to be strongly electromagnetic and dominated by hybrid-KBMs, which can drive large heat and particle fluxes that are very sensitive to equilibrium flow shear.Simulations with γ E = 0 are characterised by radially elongated turbulent structures and very large turbulent fluxes.Accurately predicting turbulent transport in this regime with γ E = 0 might require global gyrokinetic simulations that can capture equilibrium profile variation.Although substantial progress is being made in this direction (see e.g.[43]), most publicly available global codes do not, at the time of writing, retain magnetic compressional fluctuations (δB ∥ ), which are demonstrated in local gyrokinetic simulations to be essential for the hybrid-KBM to be unstable [12].When the hybrid-KBM is stabilised at long-wavelength by varying local equilibrium parameters, turbulent fluxes saturate at much more reasonable levels and turbulence appears to be well-described by the local gyrokinetic model (see e.g. Figure 15).
Local gyrokinetic simulations that include a modest level of equilibrium flow shear saturate more robustly at considerably lower turbulent fluxes (see §3).The saturated states are characterised by strong zonal flows and fields, and lower amplitude turbulent structures that are less extended radially; this turbulence appears to be better described by local gyrokinetics.Even in devices such as STEP, where rotation is expected to be modest due to the lack of external momentum injection, diamagnetic levels of flow shear can play an important role in saturating hybrid-KBMs.
Robustly saturated simulations at lower fluxes can also be obtained at γ E = 0 when the local equilibrium parameters act to suppress the linear drive for hybrid-KBMs.A study of the sensitivity of turbulent fluxes to the pressure gradient in §4 reveals a crucial competition between the linear drive of hybrid-KBMs and the stabilising impact of β ′ when β ′ is varied consistently.Strikingly, in the region of this local reference equilibrium, β ′ stabilisation overpowers the drive, and turbulent fluxes reduce or increase substantially when the pressure gradient is increased or reduced, respectively.This suggests that a favourable confinement regime exists at high β ′ , but accessing this regime would require crossing, or better evading, higher transport states at lower β ′ .Further sensitivity studies reveal that either increasing β e consistently with β ′ , or increasing q, leads to a substantial reduction in the hybrid-KBM turbulent fluxes, and expose potential transport mitigation mechanisms to supplement flow shear.These reductions in transport are found to correlate with distance between the equilibrium point and the ideal ballooning boundary in the ŝ − α diagram: at larger distances the transport fluxes and hybrid-KBM linear growth rates are reduced.This offers a simple approach to guide the optimisation of STEP scenarios, though it cannot capture important details of the turbulent transport ‡.Interestingly, we also note that heat and particle fluxes become more compatible with the available STEP sources when the local equilibrium is pushed to a regime where hybrid-KBMs are suppressed in favour of MTMs.Exploiting insights like those developed here will be essential for establishing routes to more optimal transport regimes.
Isolated turbulence from the linearly sub-dominant MTMs is also modelled in §5, by artificially suppressing hybrid-KBMs through the neglect of δB ∥ fluctuations.The MTM-driven heat flux saturates at negligible transport levels on the chosen surface in STEP-EC-HD, suggesting we should seek regimes where the hybrid-KBM is more stable, even if this comes at the expense of increasing the drive for MTMs.
In conclusion, first nonlinear local gyrokinetic simulations for a mid-radius surface in STEP-EC-HD predict that the plasma sits in a novel and poorly understood turbulent transport regime dominated by strongly electromagnetic hybrid-KBM turbulence.In the absence of equilibrium flow shear, simulations predict transport fluxes that are well above values compatible with STEP sources.Nevertheless, turbulent fluxes can be significantly reduced by including perpendicular flow shear (even at the diamagnetic level), and by other changes to the local equilibrium that increase distance from the n = ∞ ideal ballooning boundary (e.g. higher β ′ and higher q).Local gyrokinetics appears to be more suitable for describing the turbulence in these mitigated states.These first calculations will play an important role in informing the future optimisation of STEP.
The analysis performed here exposes a clear need to account for hybrid-KBM driven turbulent transport in future efforts to optimise STEP plasmas, which will likely require ‡ The low and high β e points in Figure 13 (a) which are at similar distances from the ideal ballooning boundary, but are in very different turbulent transport regimes as the low β e case is dominated by electrostatic turbulence while the high β e case is dominated by MTM turbulence several avenues to be pursued in parallel.Firstly, there is a clear need to improve the fidelity of hybrid-KBM turbulence simulations through the inclusion of physics likely to impact on the turbulent fluxes that has been neglected here, such as global physics, impurities, and fast particles.Secondly, it is urgent to hone the best available model(s) of hybrid-KBM turbulence and exploit them to predict the transport steady state; e.g. through computationally challenging coupled transport-turbulence simulations using local GK to model the hybrid-KBM turbulence or using higher fidelity reduced transport model able to capture hybrid-KBMs.A third priority is to improve our understanding of the basic physics associated with the onset of large fluxes as well as of the hybrid-KBM saturation mechanisms.Finally, it will be important to validate calculations of hybrid-KBM turbulence in detail against data from experiments in STEP-relevant regimes, i.e. in high β, low collisionality, low torque plasmas with substantial electron heating.We note that MAST-U will install additional NBI power and an Electron Bernstein Wave heating and current drive system [44], which should produce ST plasmas in regimes more relevant to STEP than it has been achieved previously.While these activities are important for advancing the ST route to a fusion power plant, they may also be useful for other devices seeking to operate in similar regimes.with long-lived turbulent structures that are radially extended.In this appendix, we show that this state of very large fluxes is not a numerical feature due to numerical resolution.A convergence study with respect to the radial box size is also presented for simulations with shearing rate of γ E = 0.1 c s /a.
Appendix A.1.Slightly reduced base grid in k y for resolution studies In order to reduce the computational cost of these resolution test simulations, the value of k y,max is decreased from k y,max ρ s = 0.95 used in the main text (e.g.see Figure 4) to k y,max ρ s = 0.63 for the simulations in this appendix.This value of k y,max is the largest k y value at which hybrid-KBMs are unstable (see Figure 2).In addition, as shown in Figure 4, modes with k y ρ s > 0.6 give negligible contributions to the total heat flux.Figure A1 compares the time trace of the total heat flux of the nominal simulation (see §3) and of a simulation with k y,max ρ s = 0.63.The time traces shows good quantitative agreement before t ≃ 300 a/c s , both reaching similarly large heat flux values in this timeframe.Around t ≃ 300 a/c s , both simulations show heat flux bursts that are qualitatively similar but differ in magnitude.After the heat flux burst, from t ≳ 400 a/c s there is better agreement and the heat flux values remain large at 200 − 300 Q gB .Therefore, apart from the heat flux burst, the two simulations agree quite well.Figure A1 compares also the k y spectra of the total heat flux averaged over the time interval t ∈ [400, 500] a/c s , and both spectra are similar at low k y .Some quantitative difference is observed at k y ρ s ≳ 0.2, where the heat flux values are slightly larger in the case of k y,max ρ s = 0.63, while the total contribution from k y ρ s ≳ 0.2 is very similar in the two cases.In order to reduce the computational cost of these test simulations, we use k y,max ρ s = 0.63 in this appendix, keeping the same value of k y,min to properly resolve the (most critical) low k y modes.
the turbulence.We note that an acceptable agreement on the linear growth rate values is observed among the three codes (see Figure 19 and Figure 20 of [12]).and n ξ = 24 (in GS2).The Sugama collision operator is used in the GENE simulation, while the linearised Fokker-Planck collision model of [50] is used in the GS2 simulation.
Figure B3 shows the time trace of the heat and particle fluxes from the simulations.The saturated level of the total heat flux from the three codes agrees very well, although the initial phase is different because of different initial conditions.The saturated heat flux level is approximately 10 −3 MW/m 2 in all the three simulations.The particle flux is negligible.All codes agree on the most important point that the MTM transport fluxes are essentially negligible.A higher level of quantitative agreement between codes is very difficult to achieve since this case is very close to marginal stability, with very low turbulent fluxes.

Appendix C. Saturation mechanisms of the MTM instability
Previous theories and simulations of microtearing turbulence have reported various saturation processes through energy transfer to long and short wavelengths [28,51], background shear flow [29], zonal fields [31,32], electron temperature flattening [42], and cross-scale interaction [30].In the STEP MTM simulations considered here, we find that electron temperature flattening and zonal fields play an important role in the saturation mechanism.As the electrons move swiftly along the perturbed magnetic field lines associated with magnetic islands at the rational surfaces, they undergo a periodic radial excursion, leading to a flattening of the electron temperature.Given that the electron temperature gradient provides the drive for the MTM instability, this flattening can locally stabilise the mode.We test whether this occurs in our simulations by removing the zonal electron temperature perturbations that cause the local temperature flattening, i.e. by redefining the zonal component of the electron distribution function as ⟨δf mod e ⟩ y,θ = ⟨δf e ⟩ y,θ − K(r)[m e v 2 /(2T e ) − 1.5]⟨F e0 ⟩ y,θ , where ⟨•⟩ y,θ denotes the fluxsurface average, F e0 is the electron background Maxwellian distribution function and K(r) is a function of the radial coordinate only, which is set at each time step such that ⟨δT e ⟩ y,θ = 0 [42].The resulting simulation, shown in Figure C1, returns a much larger heat flux than the reference case.Zonal fields can also provide a strong saturation mechanism for MTM turbulence by reducing the amplitude of non-zonal δA ∥ modes and the subsequent magnetic stochasticity.We also perform a nonlinear test simulation where the zonal fields are removed.As shown in Figure C1, this test simulation also returns much larger heat flux than the nominal simulation.Interestingly, the time trace of the heat flux from the two test simulations agrees very well, thus suggesting a competition (or a link) between zonal fields and local temperature flattening as saturation mechanisms.An additional simulation test also in Fig. C1 shows that zonal flows are not relevant for MTM turbulence saturation, at least in the case considered here.

Figure 1 :
Figure 1: Profiles from the STEP-EC-HD flat-top equilibrium: (a) contour plot of the poloidal magnetic flux, with the red contour denoting the q = 3.5 flux surface corresponding to normalised poloidal flux Ψ n = 0.49; (b) density and temperature profiles of electrons, deuterium and tritium as functions of Ψ n .

Figure 3 :
Figure 3: (a) Electrostatic electron and ion heat fluxes, Q e,es and Q i,es , as well as electromagnetic electron and ion heat fluxes, Q e,em and Q i,em , evaluated from nonlinear simulations at different values of γ E .(b) Electrostatic and electromagnetic particle fluxes, Γ es and Γ em , evaluated from nonlinear simulations at different values of γ E .The horizontal dashed lines indicate maximum fluxes that could be consistent with the reference STEP equilibrium, and are taken as the total heat and particle sources, respectively, divided by the area of the Ψ n = 0.49 surface.Open markers are used at γ E = 0 where the turbulent flux bursts (see Figure4) influence the time averaged values.

Figure 4 :
Figure 4: Time trace of the lowest 16 k y values of δϕ (a), δA ∥ (b) and Q tot (c) from the nonlinear simulation without flow shear.(d) k y spectrum of the total heat flux averaged over time between t = 150 a/c s and the end of the simulation time.

Figure 5 :Figure 6 :
Figure 5: Snapshot contour plots on the outboard midplane (at θ = 0) from the final time step of the simulation without equilibrium flow shear: (a) δϕ and (b) δA ∥ .

Figure 7 :Figure 8 :
Figure 7: Snapshot contour plots on the outboard midplane (at θ = 0) from the final time step of the simulation with γ E = 0.05 c s /a: (a) δϕ, and (b) δA ∥ .

Figure 9 :Figure 10 :
Figure 9: Electrostatic and electromagnetic electron and ion heat (a) and particle (b) fluxes from nonlinear simulations with different pressure gradient values where β ′ is varied consistently with L p,ref /L p , where a/L p,ref denotes the reference value of the total equilibrium pressure gradient.The horizontal dashed line indicates the target value.

k
y ρ s ≥ 0.3.The mode frequency values of the nominal and larger pressure gradient cases are comparable.

Figure 11 :
Figure 11: Total heat (a) and particle (b) fluxes from CGYRO nonlinear simulations at various values of β e with β ′ varying consistently.The dashed vertical line denotes the reference value of β e .Open markers are used for the reference simulation where the flux bursts influence the time averaged values.

Figure 12 :
Figure 12: Total heat flux k y spectrum from the nonlinear simulation with β e = 0.005 (a) and β e = 0.02 (b).

Figure 13 :
Figure 13: Ideal ballooning stability boundary in the ŝ − α plane of STEP-EC-HD.Panel (a) shows the position in the ŝ − α plane of the reference equilibrium and two cases with different β e .Panel (b)shows the ideal ballooning boundary for q = 3.5, q = 4 and q = 4.5 cases.

Figure 14 :
Figure 14: Electrostatic and electromagnetic electron and ion heat (a) and particle (b) fluxes from nonlinear simulations with different safety factor values.The horizontal dashed line denotes the value of the heat and particle fluxes at evaluated from the available assumed STEP sources.Open markers are used for the reference simulation where the flux bursts influence the time averaged values.

Figure 15 :
Figure 15: Top row: time trace of the electrostatic and electromagnetic heat (a) and particle (b) fluxes from the nonlinear simulation with q = 4.5.Bottom row: snapshot of δϕ (c) and δA ∥ (d) taken in the outboard midplane at the last time step.

Figure 16 :
Figure 16: Growth rate (a) and mode frequency (b) values from linear simulations at q = 3.5 (nominal value), q = 4.0 and q = 4.5.Stable modes are shown with open markers.

Figure 17 :
Figure 17: Top row: time trace of the electrostatic and electromagnetic heat (a) and particle (b) fluxes from a nonlinear simulation with δB ∥ = 0 (MTM instability).Bottom row: time trace of the k y spectrum of δϕ (c) and δA ∥ (d) from the same simulation.

Figure 18 :
Figure 18: Snapshot of δϕ (a) and δA ∥ (b) taken in the outboard midplane at the last time step from the simulation without δB ∥ .

Figure A1 :
Figure A1: (a) Time trace of the total heat flux from nonlinear simulations with k y,max ρ s = 0.95 and k y,max ρ s = 0.63.(b) k y spectrum of the total heat flux from the simulations averaged over the time interval t ∈ [400, 500] a/c s .

Figure B1 :Figure B2 :
Figure B1: Time trace of the total heat (a) and particle (b) fluxes from the CGYRO and GENE simulation with γ E = 0.

Figure B3 :
Figure B3: Time trace of the total heat (a) and particle (b) flux from CGYRO, GENE and GS2 nonlinear simulations without δB ∥ (subdominant MTM instability).

Figure C1 :
Figure C1: Time trace of the total heat flux from the nominal simulation (blue line) and three test simulations where the zonal flows (orange line), zonal fields (green line) or zonal temperature (red line) are removed from the system.