Nonlinear microtearing modes in MAST and their stochastic layer formation

First nonlinear gyrokinetic simulations of microtearing modes in the core of a MAST case are performed on two surfaces of the high-collisionality discharge used in Valovi\v{c} et al. Nucl. Fusion 51.7 (2011) to obtain the favorable energy confinement scaling with collisionality, $\tau_E\propto\,\nu_*^{-1}$. On the considered surfaces microtearing modes dominate linearly at binormal length scales of the order of the ion Larmor radius. While the effect of electron collision frequency is moderate in linear simulations, a strong dependence on this parameter is found in nonlinear simulations at $r/a=0.5$, where $r$ and $a$ are the surface and tokamak minor radius, respectively. The dynamics of magnetic islands generated by microtearing modes is analysed, showing that the radial extent of the stochastic region caused by islands overlapping plays an important role in determining the saturation level of the microtearing mode driven heat flux. Local nonlinear gyrokinetic simulations show that the microtearing mode driven heat flux, $Q_e^\mathrm{MTM}$, is largely dominated by magnetic flutter and depends strongly on the magnetic shear, $\hat{s}$. Comparing two surfaces, $r/a=0.5$ and $r/a=0.6$, reveals that $Q_e^\mathrm{MTM}$ is negligible at $r/a=0.5$ ($\hat{s}=0.34$), with the electron temperature gradient driven heat flux, $Q_e^\mathrm{ETG}$, comparable to the experimental electron heat flux, $Q_e^\mathrm{exp}$, while $Q_e^\mathrm{MTM}$ is significantly larger and comparable to $Q_e^\mathrm{ETG}$ and $Q_e^\mathrm{exp}$ at $r/a=0.6$ ($\hat{s}=1.1$). Microtearing modes cause more experimentally significant transport in higher $\hat{s}$ regions and may influence (together with electron temperature gradient modes) the observed scaling of energy confinement time with collisionality (Valovi\v{c} et al. Nucl. Fusion 51.7 (2011)).


Introduction
Previous theoretical and numerical works have shown the presence of a linear tearing instability at high mode numbers driven by an electron temperature gradient and denoted as microtearing mode (MTM) [1]. Recent nonlinear simulations have shown that the MTM instability can significantly contribute to the electron heat flux in the edge of H-mode plasmas as well as in the core of spherical tokamaks (see, e.g., Refs. [2][3][4][5][6]). The first theoretical description of the MTM instability has been proposed in Ref. [1] and extended later in Refs. [7][8][9][10][11][12]. In particular, Ref. [7] shows that the MTM instability separates into three regimes depending on the electron collision frequency, ν e = 4πn e e 4 ln λ/[(2T e ) 3/2 m 1/2 e ] (n e is the electron density, λ is the Coulomb logarithm, T e is the electron temperature and m e is the electron mass): a collisionless regime with ν e ≪ ω, a semi-collisional regime with ω ∼ (k ∥ v th,e ) 2 /ν e < ν e (with v th,e = 2T e /m e the electron thermal velocity and k ∥ the parallel wave vector), and a collisional regime with ν e ≫ ω, where ω is the MTM frequency. The work in Ref. [7] neglects the effect of the electrostatic potential in the collisionless and semi-collisional regimes. A following numerical analysis has extended this work by including the effect of electrostatic potential fluctuations, showing that these provide a strong destabilising effect [8], also confirmed in recent linear gyrokinetic simulations [13]. While first linear studies show that the mechanism driving the collisional MTM instability requires a velocity dependent collision frequency [9], unstable (collisional) MTMs have been found also when a velocity independent collision operator is considered [14], highlighting the presence of various driving mechanisms.
A magnetic perturbation δB mn associated with MTMs resonates at the rational surface with q = m/n, where q is the safety factor, m and n are the poloidal and toroidal magnetic perturbation mode number, respectively. Resonant modes can reconnect and form magnetic islands. An estimate of the island width is derived in Ref. [15], where B 0 is the unperturbed magnetic field, r the tokamak minor radius, R the tokamak major radius andŝ = (r/q)dq/dr is the magnetic shear. The distance between two rational surfaces with consecutive m and same n, corresponding to q(r m ) = m/n and q(r m+1 ) = (m + 1)/n, is approximated by ∆r ≃ 1/(nq ′ ) = r/(nqŝ) , where q ′ = dq/dr. In a typical flux-tube calculation resolving toroidal mode numbers {n 0 , 2n 0 , . . . , n = N n 0 } the minimum spacing between rational surfaces is given by [16] δr ≃ n 0 r n 2 qŝ .
If the magnetic island width is larger than the distance between two adjacent rational surfaces, a region of stochastic magnetic field lines can form [17]. Magnetic field stochasticity can provide a strong transport mechanism, as described in Ref. [18], thus accounting for the significant electron heat flux observed in nonlinear gyrokinetic MTM simulations [2]. In addition, electron heat transport consistent with the island overlap criterion has been observed in NSTX experiments [19]. While linear gyrokinetic simulations have been extensively carried out and show the presence of MTMs in many experimentally relevant scenarios [4,14,[20][21][22][23][24][25] and particularly in spherical tokamaks [26], the evaluation of the electron heat flux driven by MTMs requires one to perform nonlinear gyrokinetic simulations, which remain very challenging because of the high numerical requirements [3]. In recent years, significant effort has been devoted to understand the mechanisms behind the saturation of the MTM driven electron heat flux. For example, Ref. [2] shows that the MTM driven heat flux can be significantly reduced by equilibrium flow shear. Zonal fields [6] and local temperature flattening [27] have also been linked to the saturation of MTM turbulence. In Ref. [5], ion-scale MTMs are suppressed by electron-scale turbulence via cross-scale nonlinear interactions.
Refs. [2,28,29] shows that MTMs can drive significant electron heat transport in NSTX and links the collisionality dependence of the energy confinement time observed in NSTX to the collisionality dependence of MTM-driven heat flux. Analogously, Ref. [22] suggests a similar role played by MTM turbulence in MAST, in particular noting that, if MTMs were to dominate heat transport, lowering the electron collisionality would reduce the electron heat flux from MTMs and would be consistent with the observed energy confinement scaling τ e ∝ 1/ν * . While gyrokinetic nonlinear simulations of MTMs have been performed in NSTX to support this hypothesis [29], only gyrokinetic linear simulations have been carried out in Ref. [22].
In this work, we extend the linear study presented in Ref. [22] by performing a set of nonlinear gyrokinetic simulations of experimentally relevant cases built from the MAST discharge #22769 [22]: the simulations reported here are the first converged nonlinear simulations of MTM turbulence in MAST. All the cases considered here are characterised by a dominant ion scale collisional MTM instability, whose dependence on various parameters and, in particular, on the electron collision frequency is investigated. An electron temperature gradient (ETG) instability is found at electron scale. We also find that the MTM-driven heat flux Q MTM e is sensitive to magnetic shear and to electron collision frequency, ν e . Local gyrokinetic simulations on two neighbouring surfaces in this MAST plasma at r/a = 0.5 and r/a = 0.6 reveal that Q MTM e is negligible at r/a = 0.5 (lower magnetic shear,ŝ = 0.34), with Q ETG e comparable to Q exp e , while Q MTM e is substantial and comparable to both Q exp e and Q ETG e at r/a = 0.6 (higher magnetic shear,ŝ = 1.1).
In addition, this MAST equilibrium provides a useful reference for a scientific study of nonlinear saturation of MTMs in numerically tractable conditions. Saturation mechanisms are investigated, and indicate a significant contribution from zonal fields in the saturation process, in agreement with Ref. [6]. The level of saturated heat flux is found to strongly depend on the radial extent of the stochastic region due to magnetic island overlapping. This paper is organised as follows. The MAST reference case is introduced in Sec. 2, where the dominant microinstabilities are identified by means of linear gyrokinetic simulations. In Sec. 3, the main MTM instability is characterised and results from linear simulation scans in electron temperature gradient, density gradient and electron collision frequency are presented. Results of nonlinear simulations at r/a = 0.5 are reported in Sec. 4, where the main saturation mechanism is identified. The effect of the stochastic layer formation on heat flux is analysed in Sec. 5. The results of linear and nonlinear simulations carried out on the nearby r/a = 0.6 flux surface with higher magnetic shear are discussed in Sec. 6. Conclusions follow in Sec. 7.

The MAST reference case
The reference case is based on the MAST discharge #22769, which corresponds to the high collisionality discharge from a two point collisionality scan in MAST, where the energy confinement time was found to scale approximately as τ e ∝ ν −0.8 * [22], where ν * = ν e qR/(ϵ 3/2 v th,e ). The equilibrium and the kinetic profiles used here have been obtained by running TRANSP. Equilibrium magnetic flux surfaces for this discharge are shown in Fig. 1 at t = 0.2 s. A linear gyrokinetic analysis carried out in Ref. [22] shows that MTMs are the dominant linear microinstability at k y ρ s ≲ 1 (ρ s = c s /Ω i is the ion sound Larmor radius with c s = T e /m D the ion sound speed, Ω i = eB/m D the ion cyclotron frequency and m D the deuterium mass), thus suggesting a possible important role played by MTM turbulence in this discharge (only linear simulations were performed in Ref. [22]). Therefore this experimental case is of particular interest to characterise and analyse MTM turbulence and transport and is considered as a baseline case for the following analysis. Details on the equilibrium and profiles are reported in Ref. [22].
In this work, we perform local linear and nonlinear gyrokinetic simulations at two radial surfaces located in the tokamak core at r/a = 0.5 (depicted as a red line in Fig. 1) and r/a = 0.6. The safety factor and magnetic shear profiles in a region around r/a = 0.5 are also shown in Fig. 1. A Miller parameterisation [30], obtained by fitting the two radial surfaces using the pyrokinetics python library [31], is considered in the following. Local parameters at r/a = 0.5 and r/a = 0.6 are reported in table 1. In the following, we consider the surface at r/a = 0.5, while the surface at r/a = 0.6 is discussed in Sec. 6. Fig. 2 shows the growth rate and mode frequency of the dominant mode as a function of k y ρ s from linear simulations carried out by using the gyrokinetic code GS2 [32,33], with k y the binormal wave vector and ρ s = c s /Ω i the ion sound Larmor radius. We note in Fig. 2 the presence of two different instabilities at ion and electron Larmor radius scale. The sign of the mode frequency of both instabilities is negative (negative sign is used here for mode phase velocity in the electron diamagnetic direction). The maximum growth rate of the electron scale instability is two orders of magnitude larger than the one of the ion scale instability. The numerical resolution used in GS2 linear simulations is reported in table 2. Results of convergence tests are shown in Appendix A.
The real and imaginary components of δϕ and δA ∥ are shown in Fig. 3 as a function of the ballooning angle θ at the two different values of k y corresponding to the maximum growth rate of the ion and electron scale instabilities. The ion scale instability is characterised by δϕ extended along θ, while δA ∥ is very narrow around θ = 0. The mode has a tearing parity, i.e. δϕ has odd parity and δA ∥ even parity with respect to θ = 0. These are common features of MTMs [14]. The mode at electron scale has a twisting parity (δϕ even and δA ∥ odd) and both δϕ and δA ∥ are localised in the region of θ = 0. The electron scale modes are driven unstable by an ETG instability. This agrees with previous gyrokinetic linear simulations that have pointed out the presence   [22] at the radial surfaces corresponding to r/a = 0.5 (reference surface) and r/a = 0.6 *this surface is discussed in Sec. 6). The parameters δ, κ, ∆ ′ , ν e , a/L n and a/L T denote the plasma triangularity, the elongation, the Shafranov shift, the electron collision frequency, and the normalised inverse gradient lengths for density and temperature, respectively, with ρ * = ρ s /a.  Table 2: Numerical resolution used in GS2 and CGYRO linear simulations, with n θ and n r the number of grid points in the parallel and radial directions, respectively, and n ϵ the number of the energy grid points. In GS2 n λ is the number of pitch-angles, while in CGYRO n ξ is the number of Legendre pseudospectral meshpoints in the pitch-angle space. Results of CGYRO linear simulations are presented in Sec. 3. of the ETG instability in various MAST scenarios [34]. We note that the amplitudes of δϕ and δA ∥ , both normalised to max(δϕ), are comparable at k y ρ s = 0.5, while δϕ ≫ δA ∥ at k y ρ s = 20. This is consistent with the ion (electron) scale instability being electromagnetic (electrostatic).

MAST
In this work, we focus only on the MTM instability, even though, as demonstrated in Appendix B, ETG modes drive most of the turbulent transport at r/a = 0.5 in this particular MAST equilibrium. We highlight that the aim of this paper is mainly to investigate the saturation mechanism of this MTM instability and the role of the stochastic layer in MAST rather than to provide an accurate prediction of the heat flux in this MAST case.
Since ν e ≃ ω, this MTM instability sits between the collisionless and semi-collisional regimes of Ref. [7]. This regime has been numerically studied in Ref. [8] and analytically addressed in a recent work reported in Ref. [35]. Both works show that the growth rate of MTMs strongly depends on the electron collision frequency. Since the collisionality in the MAST reference case is ν * ≃ 0.1, trapped electron effects may also provide an additional drive for the MTM instability, as described in Ref. [10]. The MTM is a tearing instability that leads to the formation of magnetic islands. Fig. 4 shows a Poincaré map of the magnetic field at three different amplitudes of δA ∥ at k y ρ s = 0.5, which corresponds to the most unstable MTM. We highlight that the amplitude of δA ∥ used in Fig. 4 is chosen only for representative purposes, as the actual value of δA ∥ for a linear unstable mode grows exponentially until nonlinear effects cause saturation. The magnetic island forms at a rational surface located at x = 0. We note from Fig. 4 that the magnetic island width increases from w island ≃ 0.5 ρ s to w island ≃ 2 ρ s as δA ∥ is increased from 5 × 10 −3 ρ * ρ s B 0 to 8 × 10 −2 ρ * ρ s B 0 , which is in agreement with the prediction of Eq. (1). The distance between adjacent rational surfaces at k y ρ s = 0.5 is given by Eq. (2) and it is ∆r = 1/(sk y ) ≃ 6 ρ s , which is a factor of three larger than the island width at max |δA ∥ | = 0.08 ρ * ρ s B 0 . If δA ∥ is further increased, the magnetic islands generated by this mode at different rational surfaces overlap partially, therefore generating a region of stochastic magnetic field. The formation of a stochastic layer can strongly enhance heat transport, as shown later in Sec. 5. We note that multiple toroidal modes are evolved in nonlinear simulations, thus reducing the distance between adjacent resonant surfaces as compared to the distance between rational surfaces for fixed n.

Linear characterisation of the MTM instability
We explore here the sensitivity of the MTM instability to various parameters and, in particular, to the electron collision frequency. Linear simulations are carried out with the gyrokinetic codes CGYRO [36] and GS2, using a similar numerical resolution in the two codes, as reported in table 2. A benchmark of the reference case is shown in Fig. 5, where linear simulations with and without δB ∥ are considered. A good agreement between CGYRO and GS2 growth rate and mode frequency values is observed in the region 0.2 < k y ρ s < 0.7 where MTMs are unstable. We note that the MTM instability is weakly affected by δB ∥ , in agreement with Ref. [14]. The modes at k y ρ s ≤ 0.2 are stable when δB ∥ ̸ = 0 and unstable when δB ∥ = 0. Parallel magnetic fluctuations are therefore important to suppress this low k y ion temperature gradient (ITG) instability, which has a positive mode frequency sign (phase velocity in the ion diamagnetic direction) and is also observed in the electrostatic limit. Parallel magnetic fluctuations are retained in all the following linear and nonlinear simulations, and therefore the ITG mode at low k y is stable. The effect of the electron temperature gradient is investigated in Fig. 6, where growth rate and mode frequency values are shown at different values of a/L Te . The growth rate depends on the electron temperature gradient, as predicted by early  analytical works [1,9]. We note that the dependence on a/L Te is non-monotonic and the local maximum growth rate is reached at the reference value of a/L Te . A non-monotonic dependence has also been observed in previous MAST linear simulations, as shown in Ref. [14], where a resonance mechanism occurring at ν e ≃ ω is proposed as a possible explanation, as well as in ASDEX Upgrade and JET linear gyrokinetic simulations [24].
In addition, we also show in Fig. 6 the growth rate value from linear simulations with adiabatic passing electrons and kinetic trapped electrons at two different values of a/L Te . We note that these modes are stable when adiabatic passing electrons are considered, hence pointing out a minor role played by trapped electrons. Fig. 6 shows that MTMs are stable at intermediate a/L Te values, while another instability appears at large a/L Te , associated with a transition in the mode frequency, which remains in the electron diamagnetic direction. The eigenfunctions corresponding to the mode at k y ρ s = 0.5 and a/L Te = 4.2 are shown in Fig. 7. The electrostatic potential is very elongated in the ballooning angle, similarly to the MTM instability. On the other hand, the mode has a twisting parity (δϕ is even and δA ∥ odd), it is unstable also in the electrostatic limit and it is driven unstable by kinetic passing electrons, as shown in Fig. 6. The instability appearing at large a/L Te is an ETG mode characterised by k y ρ s ∼ 1 and k x ρ s > 1, which is similar to the long wavelength ETG instability described in Ref. [37] (see Appendix C for further details on this instability).
The results of a density gradient scan are presented in Fig. 8, where the growth rate and mode frequency at different k y values are shown as a function of a/L n . The growth rate decreases as the density gradient increases, similar to recent findings in linear simulations of MTMs for a high-β spherical tokamak conceptual power plant design [25].
The effect of the electron collision frequency is investigated in Fig. 9, where growth rate and frequency values are shown at various values of ν e and k y . The MTM instability is suppressed at low collisionality in favour of ITG (at low k y ) and ETG (at high k y ), thus confirming the collisional nature of this MTM instability, which is stable in the collisionless limit. The MTM growth rate increases with ν e , until it reaches its maximum value around ν e = 0.42 c s /a. At this value of collision frequency, which is approximately half the value of the collision frequency in the reference case, the growth rate is a factor of two larger than in the reference case. The growth rate value decreases when the collisionality is further increased. Since the maximum growth rate occurs at ν e ≃ ω * e , with ω * e the diamagnetic electron frequency, this may suggest a resonance mechanism similar to the one observed in the electron temperature scan. We note that the ITG instability is suppressed when ν e > 0.4 c s /a, while the onset of the ETG instability shifts at higher k y values when ν e is increased (see Appendix C for further details). We also note that future magnetic confinement fusion devices, based on the high-β spherical tokamak concept, are expected to have a lower collisionality than the case considered here and, therefore, the corresponding scenario may be located in the phase space region where the collisional MTM growth rate increases with the electron collision frequency and where the MTM instability connects to other instabilities, such as ITG or ETG modes, although the higher β values expected in high-β reactor-scale spherical tokamak scenarios might change the location and/or the presence of a maximum of γ(ν e ). In the following, we consider the most unstable case at r/a = 0.5 occurring when ν e = 0.42 c s /a. A comparison between CGYRO and GS2 linear simulations at ν e = 0.42 c s /a is shown in Fig. 10. A good agreement is observed over the entire k y range both in the growth rate and mode frequency values. The driving mechanism of the MTM instability of Ref.
[1] requires a velocity dependent collision frequency [9]. In Fig. 11, we compare the results of linear simulations with and without a velocity dependence in the collisional operator. The MTM instability is retrieved only when the velocity dependence is retained. Modes at k y ρ s > 0.1 are stable when ν e (v) = ν e (v th,e ). The mode at k y ρ s = 0.1 is stable in the reference case and unstable when an energy independent collision operator is considered. As an aside, we note, however, that earlier linear gyrokinetic studies of MAST in Ref. [14] found that MTMs remained unstable when the energy dependence of the collision operator was similarly removed. Previous analytical works suggest also an important stabilising effect on the MTM instability from the ion dynamics [11]. Results of linear simulations with adiabatic ions are shown in Fig. 11. Growth rate values increase slightly when considering adiabatic ions, so the stabilising effect from kinetic ions is very weak in the case considered here.
Finally, the dependence on the inclusion of the electrostatic potential fluctuations is investigated. Fig. 11 shows that MTMs are stable in the simulation with δϕ = 0, i.e. electrostatic potential fluctuations are essential for the drift-tearing mode to be unstable in this case. This is in agreement with the numerical calculations of Ref. [8]. We note that inclusion of δϕ was also found to be destabilising in previous MAST MTMs linear simulations [14]. This is expected from theoretical calculations in Refs. [7,35], where in Ref. [7] it is noted that the electrostatic potential plays an increasing role for MTMs in more collisional limits.

Nonlinear simulation results and saturation mechanism
We present here the results of the first nonlinear simulations of MTM turbulence carried out in a MAST case. The simulations are performed using CGYRO at the radial surface located at r/a = 0.5. Nonlinear simulations are carried out without equilibrium E × B flow shear, whose effect is investigated in Ref. [38], which shows that it has a relatively weak effect on the saturated heat flux value in this MAST reference case. We note that these nonlinear simulations are computationally quite expensive since a high numerical resolution is required to properly resolve the MTM instability, especially in the radial direction (see table 3). The simulation with the highest numerical resolution considered in this work required approximately 10 5 CPU-hours on the ARCHER2 supercomputer (Edinburgh, United Kingdom).  Table 3: Numerical resolution used in CGYRO nonlinear simulations. The quantities n ky , k y,min and L x represent the number of evolved k y modes, the minimum evolved finite k y value and the radial extent of the flux tube domain, respectively. The simulation "higherŝ" is discussed in Sec. 5. The linear scan presented in the previous section shows that the maximum MTM growth rate is achieved at approximately ν e ≃ 0.42 c s /a. We note that the heat flux driven by MTMs is negligible for ν e > 0.6 c s /a, despite MTMs being linearly unstable (see Fig. 9). Therefore, there is no contribution to the heat flux from the MTM instability in the reference MAST case at r/a = 0.5 (most of the turbulent heat flux is driven by the ETG instability as shown in Appendix B) despite MTMs being unstable. This shows that the presence of linearly unstable MTMs is not a sufficient condition to drive significant heat flux. Fig. 12 (b) shows the electromagnetic and electrostatic contribution to the saturated heat flux at different values of electron collision frequency. When ν e decreases from 0.63 c s /a to 0.42 c s /a, the heat flux increases by an order of magnitude, while the maximum growth rate of the linear MTM instability increases by less than a factor two. At ν e = 0.42 c s /a and ν e = 0.21 c s /a, the heat flux saturates approximately at Q tot ≃ 0.02 Q gB , where Q gB = ρ 2 * n e T e c s is the gyro-Bohm heat flux, corresponding to Q tot ≃ 0.002 MW/m 2 . We note that the heat flux at ν e < 0.6 c s /a saturates at a value that is more than a factor of two smaller than the saturated heat flux driven by the ETG instability in the reference case (see Appendix B). Fig. 12 (b) shows that the electromagnetic heat flux largely dominates over the electrostatic contribution, in agreement with previous nonlinear gyrokinetic simulations of MTMs [2,3].

CGYRO nonlinear simulations
Since the aim of the present work is to analyse the properties of MTM turbulence in cases built from a MAST scenario, from this point onwards we focus our studies on more strongly driven MTM turbulence at r/a = 0.5 that arises at ν e = 0.42 c s /a (which is half the nominal collisionality on this surface).
The sensitivity to some numerical parameters is investigated by carrying out a set of nonlinear simulations with different parallel grid resolution, different values of k y,min (the minimum finite k y mode evolved in the simulation) and different size of the radial flux-tube domain. In our nonlinear simulations k y,min is varied to resolve the linearly unstable MTMs, but k y,max is kept constant. The value of k y,max is chosen to exclude the range of k y where ETG modes dominate (see Appendix D for further discussion). The saturated heat flux level for each of these simulations is shown in Fig. 13. The simulation with half the number of points in the parallel direction predicts the same heat flux value within the error bar. Also the simulation with k y,min = k y,min,ref /2 predicts the same heat flux level within the error bar. On the other hand, the simulation with k y,min = 2k y,min,ref predicts a much lower heat flux. This is partially expected as modes at k y ρ s ≃ 0.2 are MTM unstable. The importance of low k y modes is highlighted in Fig. 14, where |δϕ(k x , k y )| 2 and |δA ∥ (k x , k y )| 2 spectra are shown as a function of k x and k y . The maximum value of |δϕ(k x , k y )| 2 and |δA ∥ (k x , k y )| 2 occurs at k y ρ s ≃ 0.15 and small k x values. In the simulation with k y,min = 2k y,min,ref , the |δϕ(k x , k y )| 2 and |δA ∥ (k x , k y )| 2 spectra are under resolved and the heat flux is therefore underestimated. We note that varying k y,min may potentially affect the formation of stochastic layers due to magnetic island overlapping, as discussed in Sec. 5.
The δϕ spectrum peaks at low k y values and is extended in k x [see Fig. 14 (a)], i.e. significant δϕ amplitude is observed at high k x . The electrostatic potential fluctuations are therefore very narrow radially and elongated in the binormal direction, as shown in Fig. 15 (a). On the opposite, the δA ∥ spectrum is more narrow in k x [see Fig. 14   thus resulting in much longer radial perturbations, as shown in Fig. 15 (b), which remain, however, smaller than the radial size of the flux tube domain. Given the presence of elongated δA ∥ structures, the effect of the radial extent of the flux tube domain is tested. Fig. 13 shows that no significant difference on the heat flux level is observed in the simulation with L x = L x,ref /2. This is also expected from Fig. 15 (b), where the radial extent of A ∥ fluctuations is considerably smaller than the radial size of the flux tube domain, which therefore could be reduced by a factor of two without affecting δA ∥ . We note that, because of the low magnetic shear and k y,min values, the L x and L y values used here are comparable to the size of MAST. This may question the applicability of the local approximation and global gyrokinetic simulations may be required to accurately predict the MTM driven heat flux, which is outside the scope of the present work in which we are interested in the saturation of MTMs within the local gyrokinetic framework. The |δϕ(k x , k y )| 2 and |δA ∥ (k x , k y )| 2 spectra in Fig. 14 reveal the presence of a very strong zonal component, which is more than an order of magnitude higher than the largest non zonal δϕ and δA ∥ modes. Zonal flows and/or zonal fields may therefore provide an important saturation mechanism here. In particular, the perturbed magnetic shear from zonal A ∥ perturbations,s = qR/B(dδB y /dx) ≃ 0.1, is comparable to the local magnetic shear,ŝ = 0.34, thus suggesting a potential important role of zonal fields. Saturation occurs via nonlinear interaction, where the nonlinear source term in the gyrokinetic equation can be written as [36] where x and y are the radial and binormal coordinates, h a is the non adiabatic perturbed distribution function of species a, and χ a is the generalised field potential, with R the guiding-center position, ρ = b × v/Ω ca and ⟨·⟩ R denoting the gyro-average (see Ref. [36] for details). Different simulation tests are carried out to investigate the effect of ⟨δϕ(k x , k y = 0)⟩ θ and ⟨δA ∥ (k x , k y = 0)⟩ θ in the nonlinear term of Eq. (4). The time trace of the total heat flux from these tests is shown in Fig. 16. Removing ⟨δϕ(k x , k y = 0)⟩ θ in Eq. (4) has little effect on the saturated heat flux level, thus excluding any effect of zonal flows on the saturation mechanism. On the other hand, removing ⟨δA ∥ (k x , k y = 0)⟩ θ leads to a substantial increase of the heat flux, pointing out the main role played by zonal fields, in agreement with Ref. [6].
The k x spectrum of the zonal fields, ⟨A ∥ (k x , k y = 0)⟩ θ , and their shear, ⟨k 2 x A ∥ (k x , k y = 0)⟩ θ , are shown in Fig. 17 centered around k x = 0. We note that the zonal A ∥ spectrum peaks approximately at |k x ρ s | ≃ 0.15 and decays quickly as |k x | increases. The shear of the zonal fields shows a broader k x spectrum, with its maximum value occurring in the region between |k x ρ s | ≃ 0.15 and |k x ρ s | ≃ 0.3. Therefore, we expect an important role played by the low k x zonal field modes on the saturation mechanism. The effect of the low k x zonal field modes on the heat flux is tested in Fig. 16 (b) by removing ⟨δA ∥ (k x > k x0 )⟩ θ in Eq. (4) with different values of k x0 . Lowering k x0 increases the damping of zonal fields and their associated shear, and the heat flux increases when k x0 ρ s is set below 0.25. We note that a transition to very large heat flux values is observed only when the zonal field modes with |k x0 ρ s | > 0.11 are removed from the nonlinear source term. We also mention that removing the zonal A ∥ modes with |k x ρ s | ≤ 0.11 while retaining the modes with |k x ρ s | > 0.11 causes a transition to large heat flux values, thus confirming the important of the low k x zonal A ∥ modes. Fig. 17 shows that the maximum of ⟨δA ∥ ⟩ θ and ⟨k 2 x δA ∥ ⟩ θ spectra is well resolved in simulations with L x = L x,ref and L x = L x,ref /2. On the other hand, the resolution worsens in a simulation with L x = L x,ref /4 = 1/(sk y,min ), which is the minimum L x value that can be considered keeping all other parameters fixed. This simulation is affected by heat flux oscillations and convergence difficulties. The k x resolution is therefore important here to correctly capture the saturation mechanism via zonal fields.

Magnetic islands interaction and local shear effect
We analyse here the effect of magnetic islands overlapping and subsequent formation of a stochastic layer. Fig. 18 shows a Poincaré map of the magnetic field on a radial section of the flux tube domain generated from the reference nonlinear simulation with ν e = 0.42 c s /a. Several magnetic islands with different mode numbers can be clearly distinguished. The largest island width corresponds to the mode at k y /k y,min = 2, in agreement with the δA ∥ spectrum shown in Fig. 14, where the largest non-zonal δA ∥ amplitude is achieved at k y /k y,min = 2. The island width is smaller at higher k y and the subsequent perturbation of the equilibrium magnetic field is weaker. Since the minimum radial separation between adjacent islands is proportional to 1/n 2 ∝ 1/k 2 y (see Eq. 3), the effect of these high-k y islands can nonetheless be important. For example, the radial separation of resonant surfaces for the k y /k y,min = 5 islands is ∆r = 1/(ŝk y ) ≃ 8.4 ρ s and these islands appear in Fig. 18 at x ≃ 4 ρ s , x ≃ 13 ρ s and x ≃ 21 ρ s . The k y /k y,min = 7 islands appear at x ≃ 2.2 ρ s , x ≃ 8 ρ s , x ≃ 14 ρ s and x ≃ 19.6 ρ s (not clearly visible in Fig. 18). The k y /k y,min = 3 and k y /k y,min = 9 islands also appear at a radial surface near x = 14 ρ s . The overlap of these magnetic islands generates a layer of stochastic magnetic field lines in proximity of x = 14 ρ s , as shown in Fig. 19 (a). We note that the stochastic layer around x = 14 ρ s is quite narrow in the radial direction and is surrounded by regions of weakly perturbed magnetic field. In fact, the overall magnetic field stochasticity is relatively low, with regions of well separated islands and almost unperturbed magnetic field. This is consistent with the MTM driven heat flux at ν e = 0.42 c s /a (where the MTM instability is linearly most unstable) saturating at a relatively small value.
The main role of zonal fields is to saturate the non zonal δA ∥ at low amplitude, thus reducing the width of the magnetic islands and therefore the stochastic layer size (see Sec. 4). A strong zonal A ∥ can directly affect the formation of a stochastic layer. This is shown in Fig. 19 (b), where the Poincaré map is generated without including the zonal A ∥ (the nonlinear simulation does include the zonal A ∥ ). By comparing Fig. 19 (a) and (b), we note that the magnetic islands at x ≃ 12 ρ s and x ≃ 14 ρ s merge together in the case without zonal A ∥ , thus extending the radial size of the stochastic layer. We also note that the zonal A ∥ slightly shifts the radial position of rational surfaces. For example, the resonant surface with k y /k y,min = 5 moves from x ≃ 12.5 ρ s to x ≃ 13.5 ρ s by adding the zonal A ∥ .
Heat flux transport caused by stochastic magnetic field lines depends on both the size of magnetic islands and the separation between adjacent resonant surfaces. The island width is estimated from δA ∥ using Eq. (1), while the minimum separation between adjacent resonant surfaces is given by Eq. (3). By following Refs. [2,16], we compare in Fig. 20 (a) w island and δr from the nonlinear simulation at ν e = 0.42 c s /a. We note that δr is larger than w island at all the k y modes evolved in the nonlinear simulations, in agreement with the low level of field lines stochasticity shown in Fig. 18.
The resonant radial surface spacing is inversely proportional to the magnetic shear. A larger stochastic layer is therefore expected to form at higher magnetic shear. We consider therefore an additional case at higher magnetic shear,ŝ = 2ŝ ref ≃ 0.7. The growth rate as a function of k y in this new case is shown in Fig. 21. The maximum growth rate of the MTM instability is higher than in the reference case withŝ = 0.34 and it occurs at a lower k y value. A nonlinear simulation at higherŝ is carried out with the numerical resolution listed in table 3. A lower value of k y,min is used here to account for the low k y unstable modes. Fig. 20 (b) shows the comparison between w island and δr in the simulation at higherŝ. The island width at different k y values is comparable to the one from the nonlinear simulation atŝ = 0.34. In fact, the amplitude of A ∥ fluctuations from the two simulations with different magnetic shear is comparable. This is in agreement with the nonlinear MTM theory developed in Ref. [39], which predicts |δB/B| ≃ ρ e /L Te ≃ 3.5 × 10 −5 (the temperature gradient is the same in the two simulations). This value is close to the one obtained from nonlinear simulations at ν e = 0.42 c s /a, i.e. |δB/B| ≃ 6 × 10 −5 . On the other hand, the resonant surface spacing is a factor two smaller in the simulation at higherŝ. This leads to an important qualitative change in Fig. 20 (b), where w island is larger than δr for k y ρ s > 0.3. Consequently, island overlapping is more effective and a stochastic layer is expected to extend over a wider region. This is clearly shown by the Poincaré map in Fig. 22. The radial extent of the stochastic layer is much larger than in the reference case (see Fig. 18).
Most of the magnetic islands visible in Fig. 18 in the low magnetic shear simulation are destroyed in the high magnetic shear simulation by the presence of a stochastic layer. Consequently, the heat flux increases in the higher magnetic shear simulation (see Fig. 23) and largely overcomes the heat flux driven by the ETG instability. We note that w island values in the low and high magnetic shear simulations are similar and the formation of a stochastic layer is mainly caused by a factor of four reduction of the minimum adjacent resonant surface spacing. Following Ref. [18], a magnetic diffusion coefficient is introduced, where r i is the radial position of a field line, l is the distance along the field line and N is the number of field lines considered in the average. The magnetic diffusion coefficient converges to a well-defined value at large N . The magnetic diffusivity is computed in all the nonlinear simulations by considering the full radial extent, N = 400 field lines and integrating along the perturbed field line for 2000 poloidal cycles. An estimate of the electron heat transport due to stochastic magnetic field lines can be derived from D m [16], where f p ≃ 1 − r/R is the fraction of passing particles (magnetically trapped particles do not contribute to the stochastic transport [18]). Fig. 23 compares the electron heat flux given by Eq. (7) to the electromagnetic electron heat flux calculated from nonlinear simulations. The trend is well reproduced and the predictions of Eq. (7) are in qualitative agreement with the heat flux predicted by nonlinear simulations. Importantly, the order of magnitude increase in the heat flux observed at higher magnetic shear is reproduced by Eq. (7). We note that the electromagnetic electron heat flux is entirely due to the stochastic magnetic diffusivity in the case ofŝ = 2ŝ ref ≃ 0.7, while a smaller contribution from the stochastic transport is observed in the simulations with the nominal magnetic shear value, though the stochastic contribution remains nonetheless important. The minimum spacing between resonant surfaces depends on n 0 and, therefore, on k y,min [see Eq. (3)]. As a consequence, different values of k y,min may affect the stochastic layer and the subsequent saturated heat flux level. Fig. 23 shows the stochastic heat flux in the case with k y,min /k y,min,ref = 2.0 and k y,min /k y,min,ref = 0.5. The value of Q e,stochastic is much smaller when k y,min /k y,min,ref = 2.0 than in the reference case and it is significantly lower than the heat flux predicted by the nonlinear simulation. As discussed in Sec. 4, a smaller heat flux is expected as the lowest k y unstable mode is not included in the simulation. The discrepancy between Q e,stochastic and the electron heat flux from the simulation suggests that the stochastic magnetic diffusivity is underestimated by the simulation with k y,min /k y,min,ref = 2, and the (very low) heat flux is driven by a different mechanism. In fact, the minimum resonant surface spacing increases with k y,min , thus reducing the magnetic island overlap. On the other hand, Q e,stochastic at k y,min /k y,min,ref = 0.5 is similar to Q e,stochastic at k y,min = k y,min,ref , despite the smaller resonant surface spacing. In fact, the amplitude of δA ∥ (k y ) decreases as the number of k y modes increases, thus preserving the condition w island < δr. This comparison highlights that convergence on k y,min should be carefully verified when performing nonlinear MTM simulations in order to avoid a possible underestimation of the stochastic heat flux, even when all the low k y unstable modes are included in the nonlinear simulation, especially when w island ≃ δr. The condition w island ≃ δr can also be used to indicate when substantial stochastic electron heat flux might arise from MTMs.
6. Gyrokinetic analysis at r/a = 0.6: experimentally significant Q MTM e at higher magnetic shear In this section we briefly discuss the results of linear and nonlinear simulations carried out at the radial surface corresponding to r/a = 0.6. The local parameters of this surface are shown in table 1. We note that this surface shares similar local parameter values with respect to the surface at r/a = 0.5, extensively analysed in the previous sections, except for a higher value of magnetic shear,ŝ = 1.1. Based on the results of Sec. 5, we expect a more important role played by MTMs turbulence at r/a = 0.6 than at r/a = 0.5 (although we note that the value of β, important for MTMs, is slightly smaller at r/a = 0.6 than at r/a = 0.5).
The growth rate values as a function of k y at r/a = 0.6 are shown in Fig. 24 along with the results at r/a = 0.5 for comparison. MTMs and ETG modes dominate at ion and electron scale, respectively. We note that the maximum growth rate of both the MTM and ETG instabilities is larger at r/a = 0.6 than at r/a = 0.5. The MTM instability range extends at lower k y values at r/a = 0.6. Despite these quantitative differences, the growth rate spectrum on the two surfaces is qualitatively similar.
Similarly to Sec. 3, we focus our analysis on the ion scale MTM instability. Fig. 25  shows the time trace of the total heat flux, which is largely dominated by the magnetic flutter, from a nonlinear gyrokinetic simulation at r/a = 0.6. For comparison purposes, Fig. 25 shows also the total heat flux at r/a = 0.5. We note that the total heat flux driven by MTMs at r/a = 0.6 saturates at approximately Q tot ≃ 0.2 Q gB , corresponding to Q tot ≃ 0.01 MW/m 2 , which is two orders of magnitude larger than the MTM-driven heat flux at r/a = 0.5. We also remark that Q MTM e at r/a = 0.6 is comparable to the ETG-driven heat flux (see Appendix B) and the total heat flux driven by MTMs and ETG modes is close to the experimental electron heat flux in this MAST discharge at r/a = 0.6 (see Fig. 2(e) of Ref. [22]).
The analysis at r/a = 0.6 supports a role for MTMs complementing ETG modes in generating experimentally significant transport in higherŝ regions in MAST: MTMs together with ETG might therefore be expected to influence the observed scaling of energy confinement time with collisionality, τ E ∝ 1/ν * [22], where candidate mechanisms to explain such a scaling have been proposed based on MTMs [28] and ETG modes [40]. We note, however, that our single scale computations of ETG and MTMs neglect multiscale interactions between modes that may be important in setting the overall turbulent transport [5]. However, computationally demanding multi-scale nonlinear simulations are outside the scope of the present work.

Conclusions
The MTM instability can provide significant electron heat flux transport in the highβ core of spherical tokamaks as well as in the edge of conventional aspect ratio tokamaks. An accurate prediction of MTM-driven heat flux requires one to perform expensive nonlinear gyrokinetic simulations, which are often very challenging because of their numerical requirements and convergence difficulties. This motivates improved understanding of the saturation and transport mechanisms to aid the development of cheaper reduced models. In this work, the results of the first gyrokinetic nonlinear simulations of MTM turbulence in a MAST scenario are presented, therefore extending the linear analysis of Ref. [22]. Linear simulations show that MTMs are the dominant linear instability at k y ρ s < 1 in all the cases considered here and that this MTM instability is sensitive to electron temperature gradient, density gradient, magnetic shear and electron collision frequency. At r/a = 0.5 the magnetic shear is low (ŝ = 0.34) and MTMs drive a negligible fraction of the total electron heat flux. On the other hand, slightly further out (r/a = 0.6) the magnetic shear is larger (ŝ = 1.1) and the MTM-driven heat flux is much larger and experimentally relevant. This MTM instability requires a velocity dependent electron collision frequency, it is weakly affected by ion dynamics or parallel magnetic fluctuations, and it is strongly destabilised by the inclusion of the electrostatic potential. A comparison between CGYRO and GS2 linear simulations is also carried out, showing an overall good agreement both in the reference case and in a case of stronger MTM drive. The heat flux driven by MTMs at r/a = 0.5 is negligible in the nonlinear simulation with the reference value of electron collision frequency as compared to the ETG driven heat flux, but reducing ν e by a factor of two results in Q MTM e rising by an order of magnitude, but still only to a level where Q MTM e < 0.5 Q ETG e . A strong zonal ϕ and A ∥ is observed in all the nonlinear simulations. While the effect of zonal ϕ on the saturated level of heat flux is weak, a much larger heat flux is obtained when the zonal A ∥ is removed from the nonlinear source term, thus pointing out the importance of zonal fields in the saturation mechanism of this MTM instability in the considered MAST case.
The MTM instability leads to the formation of magnetic islands at resonant surfaces. These magnetic islands can overlap if their width exceeds the radial separation between adjacent resonant surfaces, thus generating a layer of stochastic magnetic field lines. The effect of the stochastic layer formation and its radial extent on the heat flux is analysed in various nonlinear simulations. In the reference MAST case at r/a = 0.5, the island widths are smaller than the minimum spacing between resonant surfaces and this is consistent with the saturation at low heat flux level. On the other hand, the heat flux increases by more than an order of magnitude when the value of the magnetic shear at r/a = 0.5 is doubled. In this case, the island width exceeds the radial spacing between rational surfaces and a radially extended stochastic layer forms.
Heat transport caused by stochastic magnetic field lines is quantified through the magnetic diffusivity D m given in Eq. (6). A reasonable agreement is found between the electromagnetic heat flux predicted by nonlinear simulations and the stochastic heat flux at different values of electron collision frequency and magnetic shear, thus suggesting that the MTM driven heat flux in the cases considered here is mostly due to stochastic magnetic field diffusivity and confirming the important role played by the formation of stochastic layers. The criterion w island ≳ δr can therefore be used to indicate when substantial electron heat transport might be expected to arise from stochastic fields. The saturated level of heat flux is shown to depend on the radial extent of the stochastic layer, which in turn depends on both the magnetic island width, proportional to the magnetic fluctuation amplitude, and the radial separation between adjacent resonant surfaces. Nonlinear simulations with similar magnetic perturbation amplitude but different radial separation between adjacent resonant surfaces are shown to saturate at very different values. It is also interesting to note that reduced Rechester-Rosenbluth based models of electron heat transport from stochastic magnetic fields [18] have been shown to model transport in MAST discharges reasonably well [41].
Nonlinear simulations at r/a = 0.6 at higher magnetic shear find much larger Q MTM e that is similar in magnitude to Q ETG e and experimentally relevant, thus supporting a role for MTM turbulence complementing ETG modes in generating experimentally significant transport in higherŝ regions of MAST. Both MTMs and ETG modes might therefore be expected to influence the observed scaling of energy confinement time with collisionality, τ E ∝ 1/ν * [22]. We remark that our simulations do not account for the multi-scale interaction between ETG modes and MTMs. Multi-scale simulations are computationally very expensive and outside the scope of the present work. Future work will be performed in MAST scenarios to specifically address the multi-scale interaction between MTMs and ETG.
We also note that further work is required to derive a relation between magnetic fluctuation amplitude, adjacent resonant surfaces separation and saturated heat flux value, which will extend current quasi-linear theories (see, e.g., Ref. [42]) and lead towards reliable and accurate MTM driven heat flux predictions from reduced models. of the total heat flux from GS2 and CGYRO nonlinear simulations of the reference MAST case at r/a = 0.5 and r/a = 0.6. The heat flux saturates at Q tot /Q gB ≃ 0.05 on the r/a = 0.5 surface and at Q tot /Q gB ≃ 0.2 on the r/a = 0.6 surface, where Q gB is approximately 0.09 MW/m 2 at r/a = 0.5 and 0.05 MW/m 2 at r/a = 0.6. We note that the ETG-driven heat flux values are of the same order of magnitude of the heat flux computed by using TRANSP for this MAST discharge (see Fig. 2 of Ref. [22]). We highlight that the ETG heat flux is approximately two orders of magnitude larger than the MTM-driven heat flux at r/a = 0.5 while it is comparable to the MTM-driven heat flux at r/a = 0.6.
Appendix C. Binormal ion scale ETG instability Fig. 9 shows that, when the electron collision frequency value is decreased, two different instabilities appear: an instability with positive mode frequency at k y ρ s < 0.4, which is identified as an ITG instability, and an instability with negative mode frequency at k y ρ s > 0.4, which we show here to be the electron scale ETG instability that extends into the ion scale k y region at low electron collision frequency. Fig. C1 shows the growth rate and the mode frequency values of unstable modes with k y ρ s > 0.8 from GS2 linear simulations at the reference value of electron collision frequency and at a lower value. The growth rate of the ETG instability is higher at lower ν e than in the reference case and the ETG instability extends to modes at k y ρ s ≃ 1.
The real and imaginary parts of δϕ and δA ∥ at k y ρ s = 1.0 are shown in Fig. C2 for the simulation with ν E = 0.21 c s /a. Electrostatic potential mode structure is quite extended along θ and it is similar to the one in Fig. 7, which shows δϕ(θ) at k y ρ s = 0.5 in the case of large electron temperature gradient values. Therefore, while the MTM instability is suppressed at large a/L Te values, the ETG instability extends into ion scale k y region, similarly to what is observed when the value of ν e is decreased. These ETG modes, which are unstable at k y ρ s ≃ 1 and very extended along the magnetic field line, are similar to the ones characterised in Ref. [37]. We note that decreasing the collisionality decreases the electron detrapping frequency and therefore increases the drive at lower k y and ω values, as shown in appendix A of Ref. [34].
Appendix D. Heat flux spectrum and dependence on k y,max The value of k y,max in the nonlinear simulations of the present work is chosen such that the linear k y spectrum of MTMs is well resolved, while excluding the region in k y where ETG is the dominant mode. At ν e = 0.42, MTMs are unstable in the k y range 0.1 ≲ k y ρ s ≲ 0.7, while ETG modes are unstable at k y ρ s ≳ 0.9. This motivated our choice of k y,max ρ s ≃ 0.7. We briefly investigate here the effect of k y,max on the heat flux by comparing the results of four simulations with k y,max ρ s ∈ {0.53, 0.67, 0.81, 0.94}. We note that the simulation with k y,max ρ s = 0.94 includes a weakly unstable ETG mode. The saturated value of heat flux from these simulations is shown in Fig. D1 (a). The total heat flux depends weakly on k y,max at k y,max ρ s < 0.9. On the other hand, the total heat flux increases by a factor of two when k y,max ρ s = 0.94. In this last case, we also note an increase of the electrostatic heat flux, which is driven by the high k y unstable ETG. Fig. D1 (b) shows the total heat flux k y spectrum from the simulation with k y,max ρ s = 0.94. Although the heat flux peaks at low k y , significant contribution comes from the high k y region. In particular, we note that the heat flux decreases monotonically from k y ρ s ≃ 0.3 to k y ρ s ≃ 0.8 and increases at k y ρ s > 0.8 due to the ETG contribution. We highlight that only a very small fraction of the unstable ETG spectrum is included in this nonlinear simulation. Multi-scale nonlinear simulations, outside the scope of this work, are required to carefully account for MTM and ETG driven turbulent fluxes.