Electromagnetic gyrokinetic instabilities in STEP

We present herein the results of a linear gyrokinetic analysis of electromagnetic microinstabilites in the conceptual high −β, reactor-scale, tight-aspect-ratio tokamak Spherical Tokamak for Energy Production, https://step.ukaea.uk. We examine a range of flux surfaces between the deep core and the pedestal top for two candidate flat-top operating points of the prototype device. Local linear gyrokinetic analysis is performed to determine the type of microinstabilities that arise under these reactor-relevant conditions. We find that the equilibria are dominated at ion binormal scales by a hybrid version of the kinetic ballooning mode (KBM) instability that has significant linear drive contributions from the ion temperature gradient and from trapped electrons, while collisional microtearing modes (MTMs) are sub-dominantly also unstable at similar binormal scales. The hybrid-KBM and MTM exhibit very different radial scales. We study the sensitivity of these instabilities to physics parameters, and discuss potential mechanisms for mitigating them. The results of this investigation are compared to a small set of similar conceptual reactor designs in the literature. A detailed benchmark of the linear results is performed using three gyrokinetic codes; alongside extensive resolution testing and sensitivity to numerical parameters providing confidence in the results of our calculations, and paving the way for detailed nonlinear studies in a companion article.


Introduction
Magnetically confined fusion is promising as a future power source.However, the viability of fusion power plants is strongly influenced by how well the thermal energy can be confined in the plasma.Often, the dominant process governing confinement is microinstability-driven plasma turbulence.The beneficial impacts of equilibrium geometry on the microstability properties of spherical tokamaks (STs) [1] were uncovered in early studies motivated by START, MAST and NSTX [2][3][4][5]; favorable magnetic drifts, allied with higher radial pressure gradients in STs, were found capable of suppressing some of the drift-wave instabilities that drive anomalous transport in other devices [6].Furthermore, in experiments with tangential NBI, the compact nature of the ST leads to high toroidal flows (see [7]) that can act to suppress turbulence, especially at ion Larmor scales (see recent review of transport and confinement in STs [8] and references therein), though it is anticipated that an ST power plant will have minimal momentum input and only modest externally driven flow.On the other hand, the higher trapping fraction in STs contributes to an increased drive for trapped electron modes (TEMs) at high density gradients, though this drive is mitigated if the magnetic drifts are favourable [3].In addition, the high β (the ratio of thermal pressure to magnetic pressure) accessible in STs (e.g., [9,10]) can drive electromagnetic modes unstable, which can significantly increase the core turbulent transport.
It is precisely these aforementioned electromagnetic instabilities which we expect to dominate transport in high-β, reactor-scale, tight-aspect-ratio tokamaks such as STEP (Spherical Tokamak for Energy Production) [11].To be economically competitive, ST power plant designs such as STEP require a high β, which further necessitates a high β ′ (the radial gradient of pressure) in a compact device.As a result, we also require sufficiently low turbulent transport in order to sustain these steep gradients and thus to maximise the self-driven bootstrap current and reduce the need for external current drive in a steady state device.In plasmas such as STEP where β and β ′ are sufficiently high, the curvature of the confining magnetic field and the plasma kinetic gradients can excite electromagnetic instabilities such as kinetic ballooning modes (KBMs) and microtearing modes (MTMs).The KBM is driven by electrons and ions at binormalscales approaching the ion Larmor radius (k y ρ i ≲ 1), propagates in the ion diamagnetic direction, and is closely related to the ideal ballooning mode of magnetohydrodynamics (MHD) [12].MTMs excite radially localised current layers on rational surfaces, are primarily driven unstable by the electron temperature gradient, and propagate in the electron diamagnetic direction.They generate magnetic islands on rational surfaces that tear the confining equilibrium flux surfaces and enhance electron heat transport through magnetic field line stochasticisation [13][14][15][16][17].In devices where β exceeds a certain threshold value, electromagnetic instabilities can become the fastest growing instabilities in the system and the dominant sources of transport in the plasma core [18].
These two instabilities, and the nonlinear interactions between them, will likely play a crucial role in setting the transport levels in the core of devices such as STEP, and dictate the confinement times attainable in next-generation STs such as STEP.Fully understanding the transport impacts of these modes is one of the major physics questions which must be answered to build confidence in the feasibility of designs of future ST power plants.
In this work, the first of two related papers, our contributions are: (a) to report on the main results of the gyrokinetic linear analysis of two candidate STEP equilibria, STEP-EC-HD-v5 ‡ (hereinafter STEP-EC-HD) and STEP-EB-HD-v4 § (hereinafter STEP-EB-HD), at various surfaces between the core and the pedestal top; (b) to identify the dominant and sub-dominant instabilities and elucidate the nature of these modes; and (c) to explore the resolution requirements for nonlinear simulations, ‡ SimDB UUID: 2bb77572-d832-11ec-b2e3-679f5f37cafe, Alias: smars/jetto/step/88888/apr2922/seq-1 § SimDB UUID: 8ea23452-dc00-11ec-9736-332ed6514f8d, Alias: twilson/jetto/step/88888/may2422/seq-1 thus paving the way for the companion work [19] (hereinafter referred to as Paper (II)), in which we will present the first local nonlinear turbulence simulations for a STEP conceptual design.We begin in §2 by introducing the STEP equilibria∥ and the associated plasma parameters, providing some motivation of the design choices and a brief discussion of how the equilibria compare to similar ST design points [18,20].In §3, we present the main results of the gyrokinetic linear analysis of the STEP-EC-HD equilibrium at four flux-surfaces between the core and the pedestal top.This analysis reveals the importance of two particular electromagnetic instabilities, a hybrid -KBM and a collisional MTM; in §4, §5, and §6, we explore the salient features of these modes (primarily focusing on one mid-radius surface in STEP-EC-HD and parameter scans around this flux surface).In §7, we present the results of a three-code microstability comparison for two surfaces in STEP-EC-HD and one surface in STEP-EB-HD.Finally, we present our conclusions and outlook in §8.

The STEP-EC-HD and STEP-EB-HD equilibria
STEP is a UK programme that aims to demonstrate the ability to generate net electricity from fusion.STEP is planned to be a compact prototype power plant (based on the ST concept) designed to deliver net electric power P > 100 MW to the national grid [21,22].The first phase of this ambitious programme is to develop a conceptual design of a STEP Prototype Plant (SPP) and STEP Plasma Reference (SPR) equilibria for preferred flattop operating points.The ST concept maximises fusion power P fus ∝ (κβ N B t ) 4 /A [23] and bootstrap current fraction f BS = I BS /I p in a compact device at relatively low toroidal field by allowing operation at high normalised pressure β N ≃ 4 − 5 and high elongation κ > 2.8.However, alongside these advantages, the ST concept also poses unique challenges, not only in terms of plasma microstability (the focus of this work) but also in terms of the engineering constraints; the compactness restricts significantly the available space for a solenoid, so the required plasma current of I p ≃ 20 MA has to be driven, ramped up and ramped down non-inductively.
STEP plasma concepts (with the global parameters given in Table 1) have mainly been designed [11] using the integrated modelling suite JINTRAC [24], to model transport and sources self-consistently in the core plasma with prescribed boundary conditions: simplified models are used for the pedestal boundary conditions, pellet fuelling, heating and current drive, and core transport uses an empirical Bohm-gyro-Bohm (BgB) model which has been tuned both to give dominant electron heat transport as observed experimentally in MAST, and also to give a desired β N .In the present design of the operating points, plasma confinement is largely assumed, with the confinement enhancement (or H-factor) over an energy confinement scaling law indicating the level of confinement required to achieve a particular non-inductive operating point satisfying a prescribed set of additional constraints.The primary drivers that constrain the confinement needed for a viable operating point include a specified fusion gain Q > 11 (a proxy for net electricity generation), a specified fusion power P fus > 1.5 GW, current drive efficiency validated against full wave modelling of either Electron Cyclotron (EC) or Electron Bernstein Wave (EBW) systems, P aux < 160 MW, a current profile consistent with MHD, vertical stability and divertor shaping constraints (see [25] for further details).Importantly, the STEP parameter regime is outside the range of validity of the most advanced reduced core transport models available, typically developed for present-day conventional tokamaks, which often do not capture the electromagnetic (EM) transport expected to prevail in STEP, as such, it is important to test the assumptions of BgB transport using linear and nonlinear GK simulations; a key thrust of this current work.
Here, we focus on two steady-state, non-inductive flat-top operating points, STEP-EC-HD and STEP-EB-HD, both of which are designed to deliver a fusion power P fus ∼ 1.8 GW.These two designs both use RF heating instead of neutral beams to generate the current drive, in order to maximise the wall area available for Tritium breeding and minimise the recirculating power fraction [11].There are modest differences between these equilibria because they use different RF current drive schemes: • STEP-EC-HD utilises only Electron Cyclotron Current Drive (ECCD) heating.
• STEP-EB-HD utilises a mixture of ECCD and Electron Bernstein Wave (EBW) heating.
Key global parameters of the preferred flat top operating points are shown in the highlighted columns of Table 1, and a contour plot of the magnetic flux surfaces in these two design points is shown in Figure 1, alongside the corresponding electron density and electron temperature radial profiles as functions of the normalised poloidal flux Ψ n .In Section §3 we perform linear microstability analysis on the surfaces Ψ n = 0.36, 0.49, 0.58, 0.71 in STEP-EC-HD.Our primary focus in Sections §4 to §6 will be on the q = 3.5 surface (Ψ n = 0.49) of STEP-EC-HD.Two surfaces from STEP-EC-HD and one surface from STEP-EB-HD are used in the three-code microstability comparisons reported in Section §7.
Table 1 also provides key global equilibrium parameters for two other recently developed conceptual burning ST plasma equilibria: the TDoTP high q 0 case [20] and an earlier concept BurST [18].
For each flux surface considered, a Miller parameterisation [26] was used to model the local plasma equilibrium.Miller parameters were fitted to the surface using pyrokinetics [27], a python library aiming to standardise gyrokinetic analysis between different GK codes and conventions.Pyrokinetics was also used throughout to facilitate the conversion of input files between the different GK codes used in this work (see the three code comparison reported in §7).Table 2 reports the local value of the normalised poloidal magnetic flux Ψ n , magnetic shear, ŝ = (ρ/q)dq/dρ, radial position, ρ = r/a, elongation and its radial derivative, κ and κ ′ , triangularity and its radial STEP-EC-HD STEP-EB-HD TDoTP-high-q 0 [20] 1: Basic global plasma parameters including the tokamak major radius, R geo , the aspect-ratio, A, the toroidal magnetic field at the tokamak magnetic axis, B T , the plasma current, I p , the electron density and temperature values at the magnetic axis, n e0 and T e0 , the fusion power, P fus , the total heating power, P tot , the ECCD heating power, P ECCD , the EBW heating power, P EBW , the radiated power, P rad , the fusion gain, Q, the normalised β, β N = βaB T /I p , the energy confinement times normalized to the multi-machine-based ITER-H98(y,2) scaling law, H98 (with radiation) and H98* (without radiation), the bootstrap fraction, f BS , and the technique used for heating and current drive, HCD, for the two baseline operating points examined in this work (STEP-EC-HD and STEP-EB-HD) are shown alongside parameters for comparable ST design points; the TDoTP high−q 0 baseline [20]; and an earlier prototype burning ST design BurST [18].Also shown are the plasma elongation κ and triangularity δ at the last closed flux surface.The red line denotes the flux surface corresponding to q = 3.5, while the blue line refers to the flux surface at q = 3.0 (considered only in STEP-EC-HD).(c) Electron density and electron temperature radial profiles of STEP-EC-HD and STEP-EB-HD.The coloured STEP-EC-HD surfaces in (a) and the vertical dashed lines in (c) denote the q=3 and q=3.5 surfaces in STEP-EC-HD, used for the bulk of our analysis in Sections §4 to §6.These surfaces together with the coloured STEP-EB-HD surface in (b) are used in three-code comparisons in Section §7.
derivative (the symbol ′ denotes derivative with respect to ρ), δ and δ ′ , the radial derivative of the Shafranov shift, ∆ ′ , the electron β, β e = 2µ 0 n e T e /B 2 T , and electron and deuterium density and temperature gradients at different flux surfaces corresponding to low order rational values of the safety factor q.For each surface, we report the value of the binormal wavenumber k y ρ s corresponding to the toroidal mode number n = 1, with ρ s = c sD /Ω D where c sD = T e /m D , Ω D = eB/m D and m D the deuterium mass.
Nominally, fives species (electron, deuterium, tritium, thermalised helium ash, and a heavy impurity species) are included in the integrated modelling of the STEP-EC-HD and STEP-EB-HD equilibria considered in this paper.The simulations performed in this paper are carried out with 3 kinetic species (electron, deuterium, and tritium) unless explicitly stated otherwise and a future work will explore the influence of fast α particles which are completely neglected in this analysis.

Overview of simulation results for various surfaces in STEP-EC-HD
We begin by finding the dominant linearly unstable modes at an initial ballooning angle of θ 0 = 0, i.e., those modes centered on the outboard midplane.For now, we focus our attention on STEP-EC-HD, with the results for STEP-EB-HD reported in §7.The linear simulations presented here are carried out with the gyrokinetic code Local parameters of all flux surfaces considered in this work, which include the surface at Ψ n = 0.36 (q = 3.0), Ψ n = 0.49 (q = 3.5), Ψ n = 0.58 (q = 4.0) and Ψ n = 0.71 (q = 5.0) of STEP-EC-HD and at Ψ n = 0.36 (q = 3.5) of STEP-EB-HD.Included also is the binormal wavenumber k n=1 y ρ s corresponding to the toroidal mode number n = 1.GS2 [28] (commit 675f0870).Later (in §7) we will verify the fidelity by comparing the main results obtained with GS2 against CGYRO [29] (commit 399deb4c) and GENE [30] (commit de99981) in a detailed three-code benchmark §7.Table 3 indicates grid parameters used in each code (see highlighted columns for GS2) for calculations that include (f B = 1) or neglect (f B = 0) the compressional magnetic perturbation δB ∥ .We find that neglecting δB ∥ is sufficient to suppress the dominant instability in our simulations, see §6.Otherwise, the physics included in f B = 0 is the same as that in f B = 1 simulations, evolving three kinetic species (electrons, deuterium and tritium).The linearized Fokker-Planck collision model of [31] is used to model collisions in the system.In this Section, we will report primarily on simulations of the dominant instability (using parameters in the f B = 1 column), and a thorough discussion of simulations of the subdominant instability (using parameters in the f B = 0 column) will be deferred to §6.
We begin by performing linear initial value calculations to find the dominant unstable modes (i.e., the fastest growing unstable mode) across a range of different of STEP-EC-HD.In CGYRO, n ξ is the number of Legendre pseudospectral meshpoints in the pitch-angle space and n ϵ is the number of generalized-Laguerre pseudospectral meshpoints.In GENE, n v ∥ and n µ are the number of grid points in the v ∥ and µ direction, respectively.In GS2, n ϵ is the number of energy grid points and n λ is the number of pitch-angles.

Electron Larmor radius scale modes
We begin by noting that there is no purely electron scale instability in the system at any of the core flux surfaces considered (see Fig. 2), a result which is due to the large β compared to conventional tokamaks.We also remark here that we see a distinct absence of the collisionless MTM which tends to dominate the instability spectrum at intermediate scales k y ρ s = O(1) in similar ST equilibria [18,32]; the absence of the collisionless MTM is due to the larger value of the density gradient owing to pellet fuelling (which strongly stabilises the collisionless MTM).

Ion Larmor radius scale modes
Approaching the binormal Deuterium Larmor radius scale k y ρ s ≲ 1, the stability landscape is somewhat more complicated.For clarity, we group the flux surfaces by spectra structure.
3.2.1.q = 3.5 (Ψ n = 0.49) From Figure 2, we note that the maximum growth rate occurs approximately at k y ρ s ≃ 0.4 for the surface at Ψ n = 0.49.Interestingly, the growth rate as a function of k y ρ s has two local maxima, similar to that seen in similar ST designs [18,32].The mode frequency is positive (i.e., the most unstable mode is propagating in the ion diamagnetic direction) for all unstable k y ρ s modes (i.e., those modes where γ > 0).The mode is stable as we approach the sub-ion Larmor radius scale (k y ρ s > 0.65, n > 130) but is unstable down to very long wavelengths (k y ρ s = 0.023, n = 5).
3.2.2.q = 4.0 (Ψ n = 0.58) and q = 4.5 (Ψ n = 0.71) Similarly to Ψ n = 0.49, the maximum growth rate occurs on both surfaces approximately at k y ρ s ≃ 0.4.For Ψ n = 0.71, the growth rate spectrum once again has two local maxima.For Ψ n = 0.58, there is a single local maxima but we observe a similar plateau structure in the growth rate spectrum.One key difference with respect to the q = 3.5 (Ψ n = 0.49) surface is that the longest wavelength unstable modes have weakly negative mode frequency (i.e., the mode is propagating in the electron diamagnetic direction).We note however, that the growth rate and frequency vary smoothly as k y ρ s decreases (c.f., the abrupt change of sign in the real frequency between the unstable modes k y ρ s < 0.65 and the stable modes k y ρ s > 0.65 on the Ψ n = 0.49 surface) suggesting that this is perhaps not a discrete mode transition but instead is a change in the nature of the dominant instability (see §4 for further discussion).
The maximum growth rate moves to slightly longer wavelengths k y ρ s ≃ 0.2 at Ψ n = 0.36, though it once again possess two local maxima.Again, we observe a slightly different dependence of the mode frequency on k y ρ s at Ψ n ≃ 0.36 compared to the q = 3.5 (Ψ n = 0.49) surface: the frequency increases at low k y , reaching a maximum around k y ρ s ≃ 0.3, the frequency then decreases and changes sign at k y ρ s ≃ 0.4.Once again, we note that this change in frequency occurs smoothly.The remainder of this manuscript is largely devoted to studying the linear instabilities identified in the STEP-EC-HD equilibrium, focusing in particular on the q = 3.5 flux surface (Ψ n = 0.49), unless otherwise explicitly indicated.

Hybrid-KBM instability
Based on previous results (see e.g., [8] and references therein), we might expect in these high-β plasmas that some of the instabilities in Figure 2 propagating in the ion diamagnetic direction are electromagnetic KBMs, especially where the local equilibrium profiles do not access second stability [33,34] This section is dedicated to studying the physics of the dominant instability identified in Figure 2.

Is the mode electromagnetic or electrostatic?
A sensible first step towards classifying and understanding this instability is to examine whether the mode is predominantly electrostatic or electromagnetic, this can be done by examining the eigenfunctions of the dominant instability identified in §3.In Figures 3-4, we plot the δϕ and δA ∥ eigenmode structures (both normalised to the maximum value of δϕ) as functions of ballooning angle θ, at k y ρ s ≃ 0.2 (Figure 3) and at k y ρ s ≃ 0.4 (Figure 4) for the flux surfaces with Ψ n = 0.36 and Ψ n = 0.49 respectively.We note that the amplitudes of δϕ and δA ∥ are comparable, thus suggesting that P.1 the mode is predominantly electromagnetic.
Electrostatic instabilities are typically characterised by |δA ∥ | ≪ |δϕ|.At both radial locations, the mode is strongly peaked around θ = 0, with even parity in ϕ and odd parity in A ∥ .Conventionally, even parity ϕ modes are called 'twisting parity' and odd parity ϕ modes are called 'tearing parity'.Therefore, P.2 the mode has twisting parity.
Based on properties P.1 and P.2, and in agreement with previous results in [18,32] with a similar parameter regime, the dominant instability may be associated with a KBM.We can investigate this further by examining whether the mode is indeed active where the equilibrium profiles do not access second stability.

Is the mode a KBM?
A general description of KBMs was presented by [35] and [36], in which the linear electromagnetic gyrokinetic equation is solved for the gyrokinetic distribution function Ψ n = 0.49, k y ρ s = 0.19 Ψ n = 0.49,∥k y ρ s = 0.19 Ψ n = 0.49,∥k y ρ s = 0.19 Ψ n = 0.36, k y ρ s = 0.39 Ψ n = 0.36,∥k y ρ s = 0.39 Ψ n = 0.36,∥k y ρ s = 0.39 Ψ n = 0.49,∥k y ρ s = 0.38 Ψ n = 0.49,∥k y ρ s = 0.38 in terms of the perturbed fields δϕ, δA ∥ , and δB ∥ , and the expression for the gyrokinetic distribution function is inserted into the field equations.This results in three coupled, linear, integro-differential equations which may then be solved in certain limits.In theory, one could analyze this system of equations to determine which design choices (e.g., shaping) are beneficial for KBM stability.However, the complexity of these equations make it difficult to assess whether kinetic effects have a net stabilising or net destabilising effect beyond simple limits [37,38].In a complex physical system such as an ST, accurately describing KBMs thus typically requires gyrokinetic simulations to explore the sensitivity of these modes.As a first step, we investigate the stability with respect to the ideal ballooning boundary, which is often used a simple proxy for KBM stability.

The Ideal Ballooning Mode
The KBM instability is often associated with the MHD ideal ballooning mode (IBM) in the limit of n → ∞, which is derived from ideal MHD and thereby neglecting kinetic effects such as the finite Larmor radius and the effect of trapped particles.Despite making considerable simplifications to the physics, the IBM still describes the basic physics of the pure KBM instability; a competition between the stabilising effect of magnetic field line bending and the destabilising effect of a plasma pressure gradient combined with "bad" magnetic curvature.Moreover, IBM stability is much more easily assessed for a given plasma, and is sometimes used as a proxy for KBM stability in models such as the predictive pedestal model EPED (see e.g., [39]) and a good correlation is generally found in the pedestal of conventional tokamaks between the region where KBMs dominate and the region that is unstable to n → ∞ ideal ballooning modes (see e.g., [12,34] and discussion therein).
An approach pioneered by [40] allows one to calculate stability quickly and easily by integrating a one-dimensional differential equation for a given field line.This has been numerically implemented in GS2's module ideal ball.Moreover, for some fixed set of geometric parameters, ideal ball can be used to scan the normalised pressure gradient α = −Rq 2 dβ/dr, and magnetic shear ŝ ≡ ∂q/∂ψ to evaluate IBM stability for a given flux surface as a function of (ŝ, α).We therefore investigate where STEP-EC-HD and STEP-EB-HD are located with respect to the region of IBM stability and whether this is consistent with our suspicion that these equilibria are KBM dominated.
The results of these calculations are shown in Figure 5, where we see that all of the flux surfaces we have considered are well outside the unstable region.As such, we expect STEP plasma operating in these regions of parameter space to be stable to IBMs, a sensible proxy for KBM stability.This is an important piece of information about the dominant mode.

P.3
The dominant mode can be unstable in the region where the IBM is stable.
It is well known that kinetic effects can make the KBM unstable in the IBM-stable region, e.g. a finite ion temperature gradient was found to make KBM unstable below the beta threshold of the IBM [41].Therefore, we emphasise that P.3 is insufficient to  exclude the dominant mode from being labelled as a KBM, though it supports the need to explore broader mode properties, which will be a major focus of this paper.

Mode fingerprinting
Statements P.1 and P.2 indicate that the dominant mode has clear features consistent with the KBM instability, while P.3 seems to suggest a mode with different instability characteristics, i.e. this mode might be KBM-like, but it may also be coupling to some other modes in our system in order to be driven unstable.We can examine whether this might be the case by fingerprinting the mode.Mode 'fingerprints' [17] to identify the instabilities that cause transport losses in modern experiments from among widely posited candidates such as the KBM and others.The key idea underpinning mode fingerprinting is that analysis of both the gyrokinetic-Maxwell equations and gyrokinetic simulations of experiments reveals that each mode type produces characteristic ratios of transport in the density and heat channels.Thus, by examining the electron to ion heat and particle flux ratios, we might shed light on the nature of our instability.The important quantities for fingerprinting analysis are the particle and heat diffusivity, D α = Γ α /(dn α /dr), and , where α is the species label and Γ and Q denote the particle and heat flux respectively.The mode fingerprints identified in [17] are reported in Table 4 for the k y ρ s = 0.2 mode of q = 3.5 (Ψ n = 0.49).
Comparing our simulation results with the fingerprint identifiers given in Table 1 of [17] we see that our dominant instability does indeed have features in common with MHD-like modes (including the KBM) which cause very comparable diffusivities in all channels, and are characterised by |δE ∥ | ≪ |δϕ| (see Figures 3 and 4).However, we remark that Mode fingerprint χ e /χ i 0.86 D e /χ e 0.68 Table 4: Electron to ion heat diffusion coefficient ratio, χ e /χ i , and electron particle to heat diffusion coefficient ratio, D e /χ e , for the k y ρ s = 0.2 mode of q = 3.5(Ψ n = 0.49).
The heat and particle diffusion coefficients are computed as χ s = Q s /(∂p s /∂r) and D s = Γ s /(∂n s /∂r), respectively, where Q s and Γ s are the heat and particle fluxes of species s.
this fingerprint may also be also compatible with the ion-temperature gradient mode (ITG) and trapped electron mode (TEM).

P.4
The mode can be fingerprinted as a KBM or ITG/TEM.
Observation P.4 provides us with a way to reconcile P.3 with P.1 and P.2, the dominant mode is likely a hybrid instability.We now wish to study this hybrid instability.

Sensitivity to different physics parameters and hybridisation of the KBM
Based on our fingerprinting analysis, we have deduced that the dominant instability is likely a hybrid mode which shares the features of the KBM and of something reminiscent of an ITG or TEM.In the following, we attempt to further characterise the main instability by analysing the dependence on the inclusion of collisions, the numbers of species, the local gradients, the magnetic shear and safety factor, and β e .4.4.1.Pressure gradient When the local geometry is held fixed, the growth rate of a pure KBM increases with the total pressure gradient and β e .Thus, assessing the sensitivity of the dominant mode to these parameters allows us to test to what extent the mode is KBM-like.In Figure 6 we explore the dependence of the mode on the electron and ion temperature gradients, a/L T e and a/L T i , as well as on the density gradient ¶ a/L n , whilst all other parameters are held constant.Here, we note that both ion species temperatures are changed together.Figure 6 reveals another important characteristic of the hybrid mode, P.5 the growth rate is sensitive to changes in the pressure gradient (a KBM-like feature).
From comparing Figures 6 (a) and (b) we note that the growth rate appears to depend more strongly on the electron temperature gradient than on the ion temperature gradient for both of the binormal wavenumbers considered.In Figure 7, we show a comparison of the linear spectrum for two linear simulations with the same total kinetic pressure ¶ Note that we do not vary the electron and ion density gradients independently since quasineutrality requires n e = n i which in turn demands a/L ne = a/L ni globally).gradient but different electron temperature gradients; one with the nominal electron temperature gradient (blue markers) and one with zero electron temperature gradient (orange markers).In the second simulation, the total pressure gradient is kept fixed by putting the electron temperature gradient contribution into the ion temperature gradient.

P.6
The growth rate is sensitive to how the pressure gradient is varied -i.e. the mode is sensitive to the partitioning of the pressure gradient into electron and ion contributions.
If the instability was the MHD-like KBM, the growth rate would be the same and the two curves would be coincident + .However, what we observe is that the growth rate is much smaller in the case with zero electron temperature gradient across all scales.When understood alongside Figure 6 we thus deduce that this mode is indeed much more sensitive to changes in the electron temperature gradient than the ion temperature gradient.This asymmetry will have an impact on the transport properties associated with the mode, and suggests the KBM is hybridising with an electron instability (see e.g., [42]).Combining P.4, P.5 and P.6 suggests that while the mode is KBM-like, it deviates significantly from the pure KBM through kinetic effects.

Collisions
The impact of collisions on the growth rate of the hybrid -KBM is analysed in Figure 8 which shows the growth rate and mode frequency from GS2 linear simulations with and without collisions.We note that the growth rate values in the collisional case are lower than in the collisionless case, while the mode frequency is largely unaffected.That is, at the level of collisionality in this local equilibrium, collisions have a weakly stabilising effect on this dominant mode.
+ Strictly speaking, this is only true if T i = T e .In our simulations T i /T e = 1.03 and thus any deviation from MHD-like KBM behaviour might be expected to have more sensitivity to the ions -the opposite of what we observe.4.4.3.Species One can also explore the role of the different kinetic species in the simulation.Figure 9 shows the growth rate and mode frequency from GS2 simulations with two (electron and deuterium), three (electron, deuterium, and tritium) and five species (electron, deuterium, tritium, thermalised helium ash, and a heavy impurity).
In each instance, we ensure that the quasineutrality constraint is satisfied by making small adjustments to the value of the electron density gradient.Although there are some small quantitative differences between the three simulations, there is no large change to the linear properties of the hybrid -KBM observed when varying the species number * .In this paper we primarily consider linear simulations with three species (electrons, deuterium and tritium).The results of Figure 9 support this choice and have motivated using only two species in the nonlinear calculations reported in [19].

Trapped particles
Thus far, our sensitivity study has shown once again that our dominant mode has many properties in common with the KBM, but also has some non-KBM-like properties.The fingerprinting analysis in §4.3 suggested that the non-KBM-like properties might be due to coupling to a ITG or TEM (P.4).A further investigation of the dominant mode is presented in Figure 10, which shows the growth rate and mode frequency as functions of k y ρ s from four GS2 linear simulations: (i) the nominal simulation, (ii) a simulation with hybrid electrons, where the passing electrons are treated adiabatically (i.e., the passing particles have a Maxwellian response to δϕ perturbations), while trapped electrons are treated kinetically, (iii) a simulation with adiabatic electrons, and (iv) a simulation with adiabatic ions.We note that the hybrid  and kinetic electron curves (i.e., (i) and (ii)) follow each other closely, although the hybrid electron curve is not suddenly stabilised at k y ρ s > 0.6 (note also there is no sudden change in mode frequency).We thus determine that; P.7 the dominant instability has a substantial drive from trapped electrons.
Furthermore, we note that the simulation with adiabatic electrons is marginally stable at all binormal scales, once again highlighting the importance of kinetic electrons.The ion dynamics also provide an important drive for this instability, this can be seen by noting that the growth rate is strongly reduced in the simulation with adiabatic ions.

Labelling the dominant mode
The careful analysis presented here has revealed that the dominant instability is KBM- Figure 10: Growth rate (a) and mode frequency (b) as functions of k y from GS2 linear simulations with kinetic ions and electrons (blue makers), and hybrid electrons (orange markers), adiabatic electrons (green markers) and adiabatic ions (red markers).Solid and open markers refer to unstable and stable modes, respectively.The simulation with adiabatic electrons is electrostatic and there are no unstable modes.
like, but also has properties that suggest this mode is hybridising with other modes.For example, sensitivity scans in β e and β ′ (discussed later in §5.2) show that this instability can also be tracked to the electrostatic limit where it connects to a ion temperature gradient (ITG) mode.Furthermore Figure 10 clearly highlights the importance of ion and trapped electron dynamics, indicating hybridisation of the KBM with ITG and TEM drive mechanisms.Henceforth, based on the properties uncovered above, we will choose to refer to this mode as a hybrid-KBM.
It is important to highlight that ultimately the name hybrid -KBM is just a convenient label to refer to the properties stated below.
• The mode generally propagates in the ion-diamagnetic direction (Figure 2).
• The mode eigenfunction is strongly peaked in ballooning space and the mode has twisting parity (P.2).
• The mode is driven by the pressure gradient (P.5) and by β (see discussion in §5.2).
• The mode is unstable in a regime well below the ideal n = ∞ MHD limit.
• The mode is sensitive to how the pressure gradient is varied e.g., it varies more strongly with the electron temperature gradient than with the ion temperature gradient (P.6).
• The mode is driven by trapped electrons (P.7).
• The mode couples smoothly to an electrostatic instability (see discussion in §5.2).
• The mode requires access to δB ∥ drive in order to be unstable (see discussion in §6).

Stabilising the hybrid KBM
It should be emphasized here that understanding the nature of this mode (rather than simply naming it) is of the utmost importance for studying the high performance phase in conceptual ST reactors similar to STEP.Earlier studies for similar high beta conceptual burning ST plasmas [18] found MTMs dominating over several ranges in k y ρ s , and inferred that MTMs could cause substantial transport.Here, however, for STEP STEP-EC-HD and STEP-EB-HD (see §7 for STEP-EB-HD results) we find that these hybrid -KBMs dominate at all scales across a range of equilibria at various surfaces between the deep core and pedestal top.Paper (II) shows that this hybrid -KBM instability is responsible for driving most of the heat and particle transport in the STEP plasmas considered here.It is thus important to understand the nature of this mode and find strategies to mitigate it.

Sensitivity to θ 0
Shear E × B flows are known to play a stabilising role in gyrokinetics; a result which has been established theoretically [43] and borne out experimentally [44].Sheared E × B flows in linear local gyrokinetic simulations can be modelled by introducing a time-dependence into the radial wavenumber k x , which corresponds to the ballooning parameter, θ 0 = k x,0 /k y ŝ.For a mode in ballooning space at a given k y , the dependence of growth rate on θ 0 is a useful indicator of the mode's susceptibility to flow shear stabilisation [45].If the mode is stable at some θ 0 , then when flow shear advects the mode it can be moved into a stabilising region, reducing its effective growth rate.
Figure 11 shows the growth rate and frequency of the dominant mode at k y ρ s = 0.2 and k y ρ s = 0.3 as functions of θ 0 for the surface at Ψ n = 0.49 of STEP-EC-HD.At k y ρ s = 0.3, the growth rate is strongly suppressed as θ 0 increases, with the hybrid -KBM instability being stable already at θ 0 ≥ π/8.The growth rate at k y ρ s = 0.2 decreases as θ 0 increases from 0 to π/4.At θ 0 ≃ π/4, the mode is stable and remains close to the marginal stability until θ 0 ≃ 3π/4.At k y ρ s = 0.2, a different instability propagating in the electron drift direction (actually the previously subdominant MTM which has no such simple dependence on θ 0 [46]) appears at θ 0 ≃ π.
The high sensitivity of the hybrid -KBM instability to θ 0 suggests a possible important effect of flow shear, a relationship which is explored further in Paper (II).We note that a strong dependence of a KBM-like instability on θ 0 was also observed in a similar STEP conceptual design [18], where it was noted that in the local equilibrium studied flow shear effects may largely suppress transport from KBMs♯.
♯ It is important to remark that the stiffness of the pure KBM counters this argument by suggesting that even a small increase in drive could compensate for any stabilisation.It is also worth mentioning that in the pedestal of conventional aspect ratio tokamaks, it has been noted that owing to the stiffness of KBM transport, the KBM may still play a role in limiting gradients close to the critical value even when the mode is marginally stable [17]  All simulations in this paper in this paper are performed at θ 0 = 0 unless explicitly stated otherwise.

Sensitivity to β e and β ′
Motivated by the electromagnetic nature of the dominant instability, we study the sensitivity of the mode with respect to β e .Since the mode has many features in common with the KBM, we would expect the mode to be stable below some finite value of β e , and for the growth rate to scale with β e above this threshold.

5.2.1.
Varying β e with β ′ fixed In Figure 12, we study the impact of varying β e whilst the other parameters (notably β ′ ) are held fixed (the nominal value is denoted by a red vertical line).We note that in this case the mode is stabilised as β e is dropped (and the growth rate increases when β e is increased) which indicates that the dominant mode is accessing the electromagnetic component of the drive terms (the electromagnetic component of the drive is reduced at lower β, while the stabilising effect of β ′ is retained).
Later, in §6, we will see that the hybrid -KBM necessitates the inclusion of parallel magnetic fluctuations δB ∥ in order to access the electromagnetic drive.This result is in line with earlier works [42,47] that find parallel magnetic fluctuations act to destabilise the KBM and an absence of δB ∥ effects lead to a decrease of the KBM growth rate up to a factor of 6 [37].However, this is even more severe in simulation of the hybrid -KBM, which we find to be everywhere stable when δB ∥ is neglected.

5.2.2.
Varying β e and β ′ consistently In Figure 13 we study the effect of varying β e whilst also varying β ′ consistently with the local equilibrium, following an approach similar to that outlined in [32].Figure 13   unstable even in electrostatic limit (β e = 0, β ′ = 0).The smooth variation of the growth rate (and the real frequency) with β e indicates that the mode is coupling to some electrostatic instability which prevents stabilisation.This is another feature of the hybrid -KBM which markedly distinguishes it from the simple KBM (which would be stable in an electrostatic simulation).We remark that the behaviour seen here is consistent with coupled KBM-ITG, and similar behaviour has been seen in both theory [32] and experiment [48].Figure 12 also shows that the growth rate is reduced at higher β e due to β ′ stabilisation, noting that a mode transition occurs as β e passes through β e ∼ 0.13 (note the abrupt change of sign of frequency and the further reduction of the growth rate).At these values of β e and β ′ , the hybrid -KBM is fully stabilised, revealing the underlying subdominant MTM instability (see §6).As discussed in Section 5.1, hybrid -KBMs are strongly ballooning with growth rates that peak strongly at θ 0 = 0, while MTM growth rates depend very differently on θ 0 , with a typically weak dependence at low k y [18].For the reference equilibrium the hybrid -KBM growth rate is much more unstable at θ 0 = 0 than the highest MTM growth rate at any θ 0 .Figure 12 indicates that increasing β e together with β ′ results in stabilisation of the dominant hybrid -KBM.
As a complement, Figure 14 shows the linear growth rate spectrum for three different values of β e with consistently varied β ′ -the nominal case (orange markers) and a lower and higher β e case (blue and green markers).In the case with higher β e , the hybrid -KBM instability vanishes and the MTM instability (see §6) becomes the most unstable mode in the system.In the case of lower β e , the ITG instability drives an unstable mode with a higher growth rate than the most unstable hybrid mode at the

Implications for the current ramp
We remark here that, as mentioned in §2, one of the major challenges for STEP is the need to generate the required plasma current of I p ≃ 20 MA.Although Figures 13 and 14 demonstrate that a high beta regime free of the hybrid-mode exists, it is less clear how the hybrid mode could be avoided on the approach to such a flat-top during the I p ramp.During the current ramp, the plasma equilibria will evolve continuously from a β e = 0 state, where it will be dominated by electrostatic instabilities, up to the reference β e where it is dominated by the hybrid -KBM.However, we have seen in Figures 13 and 14 that the hybrid -KBM becomes active at much smaller β e than that which we are aiming to achieve in the flat top.Thus, in getting to this equilibria one must first pass through a region where the hybrid -KBM is active, and this could shut down the evolution of the plasma e.g., if the turbulent transport were too large to sustain the profiles (see Paper (II) for further discussion).It is currently not clear whether it's possible to avoid the onset of the hybrid -KBM completely or how much heating and fueling would be required burn through it should it appear earlier in the current ramp.It is also worth remarking that the β e = 0.16 case will probably exceed the resistive wall mode control limit and thus is likely not viable for other reasons.

Safety factor and magnetic shear
The hybrid -KBM also shows some sensitivity to changes to the local equilibrium parameters such as the magnetic ŝ and the safety factor q, (see Figure 15).Studying the sensitivity of the hybrid -KBM to ŝ and q reveals that, although we have placed great emphasis on the properties of this mode that distinguish it from the IBM, some of the physical intuition we can develop from the ideal theory is still useful for these hybrid -KBMs.Panel (a) of Figure 15 shows how the dominant mode is destabilised by increasing ŝ, consistent with the ideal ballooning mode behaviour (i.e. a KBM-like behaviour).Interestingly, inspection of the mode frequency in this scan reveals an isolated mode transition occurring only at ŝ = 0, to a mode that propagates in the electron diamagnetic direction.The dependence of the mode on q (Panel (b) of Figure 15) is slightly more complicated but is also consistent with the behaviour of the ideal ballooning mode.As q increases the stability boundary of the ideal ballooning mode moves towards higher ŝ and lower β ′ .Therefore, increasing q will make access to the second stability region easier (which is why we see the mode stabilised with increasing q from the reference value).At sufficiently low q, the stability boundary gets pushed to higher β ′ enabling the equilibrium to lie in the first stability region (see [34] for a more careful discussion of first and second stability) which is consistent with the stabilisation of the KBM-like dominant mode.We note that the global MHD equilibrium is not varied consistently in these cases.These results in Section §5 motivate the need to analyse the subdominant instability.

Subdominant MTM instability
We now turn our attention to the subdominant MTM instability, which may play an important transport role, especially if the hybrid mode is effectively stabilised by flow shear.MTMs generate magnetic islands on rational surfaces that tear the confining flux surfaces and generate heat transport primarily through the electron channel [14,16].We note that local GK simulations have revealed MTMs as the dominant microinstabilities in the wavenumber range k y ρ s < 1 locally at mid-radius (where β e ∼ 5% − 10%) in several spherical tokamak plasmas (see [8] and references therein).GK studies for the high performance phase in conceptual ST reactors have also found MTMs dominant over an extended range of binormal scales and likely to have significant impacts on transport [18,32].The presence of the fastest growing tearing mode can be investigated by enforcing the tearing (odd) parity of the perturbed distribution function, exploiting the up-down symmetry of the Miller equilibrium.This test is carried out with GS2 at θ 0 = 0 and the results are shown in Figure 16 (see orange curve).We thus see that there are in fact unstable modes with tearing parity (e.g., MTMs), but on this surface in STEP-EC-HD these are always subdominant to the hybrid -KBM.
Another way to obtain the MTM from an initial value solver as the fastest growing mode in our system specifically, without forcing the parity of the eigenmode, is to simply switch off compressive magnetic perturbations i.e., we exclude the δB ∥ contribution to the GK equation.Figure 16 shows that the simulation neglecting δB ∥ (green) is equivalent to the nominal simulation with a tearing parity initial distribution function (orange).Essentially, removing δB ∥ fluctuations from the system stabilises the hybrid -KBM, whilst having no impact on the MTM and thus leaving the MTM as the dominant mode.
The eigenfunctions of ϕ and A ∥ corresponding to the MTM at k y ρ s = 0.14 are shown in Figure 17.We find that the A ∥ fluctuation is significantly larger than the electrostatic fluctuations close to the inboard midplane, as expected for MTMs.The electrostatic potential eigenfunction exhibits a clear multiscale structure (ion-scale in k y , electron scale in k x ), with a narrow oscillatory structure in θ overlaying a much broader oscillation.The A ∥ function is more strongly peaked about θ = 0, with subsequent peaks occurring along the field line at θ mod 2π = 0, the outboard midplane.Similar MTM eigenfunctions extended in ballooning angle have been seen in simulations of MAST [13] and NSTX [14] discharges and BurST [18].The extended nature of these modes, requiring a parallel domain θ ∈ [−70π, 70π], coupled with a very small growth rate, means that even linearly resolving the subdominant MTM can become very computationally expensive.Nonlinear simulations involving the MTMs in Figure 17 will be computationally challenging in these STEP plasmas, owing to the intrinsic multiscale character of the MTM in the radial direction (which is linked to its multiscale character in θ in ballooning space).Figure 18 illustrates how the MTM growth rates (for modes at k y ρ s = 0.1 and k y ρ s = 0.3) depend on θ 0 , showing that the MTM growth rate (particularly at k y ρ s = 0.1) is much less sensitive than the hybrid -KBM growth rate (see Figure 11 of §5.1) this therefore suggests that these MTMs should be much less susceptible than hybrid -KBM modes to flow shear stabilisation.We note that tokamak regimes exist where turbulent transport from MTMs is affected by flow shear stabilisation [14].Recent theoretical work has identified an important local equilibrium parameter that helps explain this [49], and the relevance of this parameter in experiments and numerical simulations is explored in [46].The insights gained here may be helpful in the future optimisation of STEP design points.

Code comparison on surfaces in STEP-EC-HD and STEP-EB-HD
Careful benchmarking is essential for ensuring the fidelity of GK simulations in nextgeneration reactor design, and to identify (and ideally rectify) issues that may arise in simulations using any single code (see e.g., the discussion of the numerical instability in paper II).Furthermore, this benchmarking also paves the way for the detailed nonlinear investigation of the companion article.In this section, we compare the results of CGYRO, GENE and GS2 linear GK simulations carried out at the radial surfaces corresponding to q = 3.0 (STEP-EC-HD and STEP-EB-HD) and q = 3.5 (STEP-EC-HD only) equilibria.As previously, we compare simulations and results (linear eigenvalues and eigenmodes) for both the hybrid -KBM instability and the subdominant MTM instability.As before, the numerical resolutions used in these simulations are listed in Table 3, where again we resort to different resolutions for simulations of the hybrid -KBM instability and simulations of the subdominant MTM instability.These simulations evolve three species (electron, deuterium and tritium) and include both perpendicular and parallel magnetic fluctuations, δA ∥ and δB ∥ for simulations of the hybrid -KBM whilst including only δϕ and δA ∥ for MTM simulations.In each code, we try to use the most advanced physics model available whilst also ensuring results are comparable by adopting as similar approaches as possible.We have therefore used the Sugama collision model [50] in both CGYRO and GENE.† † The linearized Fokker-Planck collision model of [31] is used in GS2.
Figure 19 compares the growth rate and the mode frequency at Ψ n = 0.49 of STEP-EC-HD, and a reasonable agreement is found between all three codes.We note that it is of no great surprise that there is some variation between the growth rates since; e.g., each code employs differing discretisations of the 5D space and schemes for parallel dissipation.
The comparison for the subdominant MTM instability is shown in Figure 20.Retrieving a good agreement here is much more challenging, since the growth rates are relatively small and therefore more sensitive to the different numerical implementations and dissipation employed in the three codes (see [28][29][30] for code-specific details).In addition, it was found that numerical convergence in these simulations required a very high pitch-angle resolution (CGYRO) and a very high θ resolution (GENE) in order to capture the parallel structure of these very extended modes (Appendix A).The resolutions used in these simulations are as listed in Table 3.We also note that the maximum growth rate differs by less than 20 % when a lower resolution is considered, thus motivating the lower numerical resolution used in some of the nonlinear simulations of Paper (II).As shown in Figure 21, we can see that all three codes show a good agreement on the eigenfunctions for both the dominant and subdominant modes.† † A more advanced exact Landau collision operator [51] available in GENE has not been used here.The three code comparison is also carried out on the q = 3.0 flux surface of STEP-EC-HD; the dominant instability is shown in Figure 22 and a comparison for the subdominant MTM instability is shown in Figure 23; and also for the q = 3.5 surface of STEP-EB-HD; the dominant instability is shown in Figure 24 and a comparison for the subdominant MTM instability is shown in Figure 25.
We conclude by noting that there is a good agreement between the three codes in all the considered cases.

Conclusions
In this paper, we have presented the results of local linear microinstability studies of the thermal plasma on a range of flux surfaces from the core to the pedestal top in the two preferred STEP flat-top operating points.We find that the linear spectra is dominated by a hybrid mode, sharing features of the KBM, ITG, and TEM instability, at the ion Larmor scale, with weakly-growing subdominant MTMs present at similar scales.The local equilibria examined here (q = 3.0, 3.5, 4.0, 5.0 in STEP-EC-HD and q = 3.5 in STEP-EB-HD) were found to be completely stable to electron scale modes.
A summary of the dominant microinstabilities from some of our simulations is given in Table 5 alongside the results from some similar conceptual designs of burning ST plasmas, namely the TDoTP high q 0 equilibrium taken from [20]; and an earlier prototype for a burning ST reactor BurST taken from [18].Shown here is a summary for only the q = 3.5 surface in each equilibrium.†  -The mode is sensitive to θ 0 .
-The mode growth rate is smaller at larger values of β and β ′ .
-The mode growth rate is less sensitive to β and β ′ at low β values.
-The mode requires access to δB ∥ drive in order to be unstable.
• There is no unstable branch of collisionless MTMs at intermediate scales (k y ρ s ∼ 4) unlike in [20] and [18] • A collisional MTM is unstable but is always subdominant to the KBM-like instability.
It is important to remark that for the equilibria examined in this work, the  confinement is assumed and this is a substantial caveat.In a natural extension to the linear work presented in this paper, a parallel companion article [19] explores whether or not these assumed plasma profiles can be sustained by the available heating and fuelling, i.e. whether the turbulence driven by the instabilities studied in this paper are compatible with the assumptions in the design of STEP-EC-HD and STEP-EB-HD.Linear analysis presented here has uncovered different routes to stabilising the hybrid -KBM which could be very useful in future exploration of optimised operating points.We have shown that the hybrid -KBM can be stabilised by increasing β ′ .Due to the strong sensitivity on θ 0 , the mode might be suppressed by E × B flow shear, and this is explored in [19] by means of nonlinear turbulent simulations.To obtain further information on the data and models underlying this paper please contact PublicationsManager@ukaea.uk.

Appendix A. Numerical resolution convergence
Here we discuss the numerical resolution convergence studies in GS2 linear simulations for the dominant hybrid -KBM instability at the q = 3.5 flux surface of STEP-EC-HD.Similar resolution convergence scans were performed also with CGYRO and GENE at all the radial surfaces considered in this work, and for both the dominant and subdominant modes.These detailed convergence tests are also used to inform the resolutions used in the nonlinear simulations which are the focus of Paper (II).
Figure A1 shows the growth rate and mode frequency at different parallel grid resolutions.Convergence is achieved for n θ ≥ 32 at low mode numbers, while a higher resolution (n θ ≥ 64) is required at k y ρ s > 0.4.
Convergence with respect to the grid extent in ballooning space (which is equivalent to the radial grid resolution in the flux tube) is controlled in GS2 by the parameter nperiod, and is investigated in Figure A2.Growth rate is only slightly affected by this parameter, as expected since for the hybrid -KBM both δϕ and δA ∥ are very localised around θ = 0 (see Figure 3).It should of course be noted once again that the MTM has much more stringent conditions on the radial and parallel grid resolutions (see Table 3).
Velocity space resolution convergence is tested in Figures A3 and A4, where the number of passing pitch-angles and the number of energy grid points are varied.Over these grid parameters there is a weak dependence of the growth rate and mode frequency on the velocity space resolution over the entire k y spectrum.

Appendix B. Dependence on the collision model
Here we briefly discuss the dependence of our simulation results on the collision model.In particular, we observe that the linear spectrum is more sensitive to the plasma composition (the number of species evolved) when the Sugama collision operator [50], a sophisticated approximation to the full linearized, gyro-averaged Fokker-Planck operator, is used instead of the full linearized Fokker-Planck collision operator which is used in all of our GS2 simulations.
We first compare the growth rate and mode frequency values from CGYRO and GS2 simulations carried out by evolving two species (see Figure B1) and three species (see Figure B2) and considering different collision operator models: CGYRO simulations are performed by using the Lorentz and the Sugama collision operators, while GS2 uses the linearized Fokker-Planck collision operator.In all of these simulations, we note that there is only a very weak dependence on the type of collision model employed.
However, Figure B3 shows that growth rate values are more sensitive to the collision model when five species are considered.Interestingly, although there a relatively good agreement observed between CGYRO and GS2 simulations results when the Lorentz operator is used in CGYRO, we find that the growth rate values are smaller (by approximately 20%) when the Sugama collision operator is used in CGYRO.It is important to note however that there is no qualitative difference in the main instability.We will also see in Paper (II) that reducing the linear growth rates by such a small amount makes no significant difference to the transport properties on the surface considered.

Figure 1 :
Figure1: Magnetic flux surfaces of the of the STEP-EC-HD (a) and STEP-EB-HD (b) equilibria.The red line denotes the flux surface corresponding to q = 3.5, while the blue line refers to the flux surface at q = 3.0 (considered only in STEP-EC-HD).(c) Electron density and electron temperature radial profiles of STEP-EC-HD and STEP-EB-HD.The coloured STEP-EC-HD surfaces in (a) and the vertical dashed lines in (c) denote the q=3 and q=3.5 surfaces in STEP-EC-HD, used for the bulk of our analysis in Sections §4 to §6.These surfaces together with the coloured STEP-EB-HD surface in (b) are used in three-code comparisons in Section §7.

Figure 2 :
Figure 2: Growth rate (a) and mode frequency (b) as functions of k y ρ s from GS2 linear simulations of STEP-EC-HD at various radial surfaces corresponding to low q rational surfaces.The considered k y values cover a range corresponding to toroidal mode numbers between n = 2 and n = 5000.Frequency values are shown only for unstable modes with n < 200.Solid and open markers refer to unstable and stable modes, respectively.The growth rate of stable modes is set to zero.No unstable modes are found below n = 5.

Figure 5 :
Figure 5: Ideal ballooning stability boundary in the ŝ − α plane of STEP-EC-HD at Ψ n = 0.49 (black line) and at Ψ n = 0.36 (blue line) and of STEP-EB-HD at Ψ n = 0.35 (red line).The makers denote the equilibrium value of ŝ and α for each surface.

Figure 6 :
Figure 6: Growth rate from GS2 linear simulations as a function of a/L Te (a), a/L T i (b) and a/L n (c) at k y ρ s = 0.2 (black line) and k y ρ s = 0.3 (blue line).The red dashed line represents the reference value.All the other parameters except a/L Te in (a), a/L T i in (b) and a/L n in (c) are kept fixed.Results at the surface Ψ n = 0.49 of STEP-EC-HD.

Figure 7 :
Figure 7: Growth rate (a) and mode frequency (b) as functions of k y ρ s from GS2 linear simulations with nominal (blue markers) and zero (orange markers) electron temperature gradient and same pressure gradient.Results at Ψ n = 0.49 of STEP-EC-HD.

Figure 8 :Figure 9 :
Figure 8: Growth rate (a) and mode frequency (b) as functions of k y ρ s with and without collisions.Results from GS2 simulations at Ψ n = 0.49 of STEP-EC-HD.

Figure 11 :
Figure 11: Growth rate (a) and mode frequency (b) as functions of θ 0 at k y ρ s = 0.2 (blue line) and k y ρ s = 0.3 (red line) for the surface at Ψ n = 0.49 of STEP-EC-HD.Unstable and stable modes are represented by filled and open markers, respectively.

Figure 12 :
Figure 12: Growth rate (a) and mode frequency (b) from GS2 linear simulations as functions of β e at k y ρ s = 0.2 (black line) and k y ρ s = 0.3 (blue line).The red dashed line represents the reference value.The value of β ′ is kept constant while varying β e .Filled and open markers refer to unstable and stable modes, respectively.Results at Ψ n = 0.49 of STEP-EC-HD.

Figure 13 :Figure 14 :
Figure 13: Growth rate (a) and mode frequency (b) from GS2 linear simulations as functions of β e at k y ρ s = 0.2 (black line) and k y ρ s = 0.3 (blue line).The red dashed line represents the reference value.The value of β ′ is consistently varied with β e .Results at Ψ n = 0.49 of STEP-EC-HD.

Figure 15 :
Figure 15: Growth rate from GS2 linear simulations as a function of ŝ (a) and q (b) at k y ρ s = 0.2 (black line) and k y ρ s = 0.3 (blue line).The red dashed line represents the reference value.Results at Ψ n = 0.49 of STEP-EC-HD.

Figure 16 :Figure 17 :
Figure 16: Growth rate (a) and mode frequency (b) as functions of k y ρ s at Ψ n = 0.49 of STEP-EC-HD from GS2 linear simulations.The blue line refers to the nominal simulation with δB ∥ fluctuations (f B = 1), the orange line to the simulation with tearing parity enforced in the distribution function and the green line to the simulation without δB ∥ fluctuations (f B = 0).Only unstable modes are shown.

Figure 18 :
Figure 18: Growth rate (a) and mode frequency (b) as functions of θ 0 at k y ρ s = 0.1 (blue line) and k y ρ s = 0.3 (red line) from GS2 linear simulations with δB ∥ = 0.

Figure A1 :
Figure A1: Growth rate (a) and mode frequency (b) as functions of k y from GS2 linear simulations at Ψ n = 0.49 of STEP-EC-HD with different parallel grid resolution.

Figure A2 :Figure A3 :
Figure A2: Growth rate (a) and mode frequency (b) as functions of k y from GS2 linear simulations at Ψ n = 0.49 of STEP-EC-HD with different values of nperiod.

Figure A4 :
Figure A4: Growth rate (a) and mode frequency (b) as functions of k y from GS2 linear simulations at Ψ n = 0.49 of STEP-EC-HD with different values of n ϵ .

Figure B1 :
Figure B1: Comparison of the growth rate (a) and mode frequency (b) as functions of k y from CGYRO simulations with two different collision operators (Lorentz and Sugama) and two species.Also shown is the growth rate and mode frequency from the GS2 simulation with two species.

Figure B2 :Figure B3 :
Figure B2: Comparison of the growth rate (a) and mode frequency (b) as functions of k y from CGYRO simulations with three different collision operators (Lorentz and Sugama) and three species.Also shown is the growth rate and mode frequency from the GS2 simulation with three species.

Table 3 :
Numerical resolution used in CGYRO, GENE and GS2 linear simulations of the full model (f B = 1) and the model without δB ∥ reveals that this hybrid -KBM mode remains