On the impact of temperature gradient flattening and system size on heat transport in microtearing turbulence

Microtearing instability is one of the major sources of turbulent transport in high-β tokamaks. These modes lead to very localized transport at low-order rational magnetic field lines, and we show that flattening of the local electron temperature gradient at these rational surfaces plays an important role in setting the saturated flux level in microtearing turbulence. This process depends crucially on the density of rational surfaces, and thus the system-size, and gives rise to a worse-than-gyro-Bohm transport scaling for system-sizes typical of existing tokamaks and simulations.


Introduction
Confinement in tokamaks is enabled by magnetic field lines that trace out nested toroidal surfaces. On rational surfaces, the field lines connect back to themselves after integer numbers of poloidal and toroidal turns; certain electromagnetic plasma instabilities, such as the microtearing modes of interest here, are localized near these rational surfaces, and break the nested topology by forming magnetic islands [1][2][3]. Microtearing modes may have a significant impact on confinement, especially in high-β spherical tokamaks [4] where electromagnetic instabilities tend to be more unstable [5][6][7], and hence * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. understanding microtearing transport is crucial for designing large spherical tokamak reactors such as STEP [8]. They also are believed to be important in the pedestal in conventional tokamaks [9].
Microtearing modes are characterized by radially-narrow parallel electron current layers that are driven resonantly at the rational surface and associated magnetic islands. While various branches of microtearing modes have been identified, including those driven by the time-dependent thermal force [10] or by curvature [11,12], and are present in various collisionality regimes [13][14][15][16], the electron temperature gradient remains a necessary condition for instability in all cases. Previous works have reported various microtearing saturation mechanisms: via energy transfer to long wavelengths [2] or to short wavelengths [17], via cross-scale coupling with electron temperature gradient modes [18], by background shear flow [6] or zonal fields [19], etc. However, despite these advances, predicting saturation levels remains a challenging task. We investigate a previously studied [17,20,21] microtearing turbulence case (known to be numerically tractable) in simple geometry, resolving only the ion scales, in the absence of background flow shear, and demonstrate a saturation mechanism where the magnetic islands associated with the resonant current layers flatten the electron temperature gradient, thereby reducing the linear drive at the rational surfaces.
The radial width of the resonant region at the rational surfaces is generally of the order of few ion Larmor radii and is set by the parallel correlation length in linear theory which scales with the square root of the mass ratio between ions and electrons [22]. Nonlinearly, the flux associated with these modes, although slightly broadened, is still localized near the rational surfaces. This is already known to be important in setting global flux levels in the pedestal [9]. In the most extreme case, if turbulent diffusivity is sufficiently large and localized near low-order rational surfaces, the system will remove the local driving gradients, increase the gradients away from the low-order rational surfaces, and saturate in a zero-flux state. In this work, we find a less extreme version of this process occurring in a standard microtearing regime. We make a scaling argument to quantify this effect and suggest that future reactor-size devices subject to microtearing turbulence may perform worse that expected.
In simulations of microtearing turbulence, runaway states quite often develop with extreme flux levels. In this work, we do not directly address how and where this process occurs, although we see this occur in certain of our simulations.
We proceed by demonstrating the strong electrontemperature-gradient flattening at low-order rationals in gyrokinetic simulations, and showing that this allows saturation by reducing mode drive. We test the impact of the various saturation mechanisms by suppressing zonal modulations and show the dominance of the temperature corrugations. Lastly, we consider system-size scaling and explain the origin of a non-gyro-Bohm scaling.

Simulation set-up
Our numerical investigation uses flux-tube and global Gene gyrokinetic simulations [23,24]. In the flux-tube version [25], the background quantities and their gradients are assumed constant over its radial domain [26]. The global simulations, on the other hand, can accommodate more realistic background profile variations, but in general may be more computationally expensive, and furthermore, setting appropriate boundary conditions, sources etc can prove difficult [24,27]. By comparing these two simulation types, we are able to determine whether the more complicated global physics is playing an important role [28].
The simulations use a field-aligned coordinate system [25] where x is the radial coordinate, y the binormal coordinate and z the parallel coordinate. The binormal wavenumber k y = nq 0 /x 0 , where n is the toroidal mode number and q 0 is the safety factor at the radial position x 0 where the simulation is centered. Parallel velocity v ∥ and magnetic moment µ are the velocity space coordinates. We use the simple microtearing-dominated equilibrium used in [17,21], modelling the outer core region of a tokamak. Concentric circular flux-surface geometry [29] is considered with an inverse aspect ratio ϵ = x 0 /R = 0.15. An electron-ion mass ratio of m i /m e = 1836, temperature ratio of T i,0 /T e,0 = 1 and normalized electron pressure of β e = 0.4% are considered. To model collisions, the linearized Landau operator is used with an electron-ion collision frequency ν ei /(v th,e /R) = 0.02, where v th,s = (T s,0 /m s ) 1/2 is the thermal velocity of species s. δB ∥ fluctuations are not included.
In the flux-tube simulations, a safety factor of q 0 = 3 and magnetic shearŝ = 1 are considered. The inverse of the density, ion temperature and electron temperature background gradient scale lengths, normalized to the major radius R, are R/L n = 1, R/L T i = 0 and R/L Te = 4.5, respectively. The standard nonlinear flux-tube simulation considered in this work has a minimum binormal wavenumber of k y,min ρ i = 0.02 (L y = 314.2ρ i ), radial width of L x = 150ρ i and grid resolutions given by It is run until normalized time 2000R/v th,i (max lin. growth rate = 0.018v th,i /R) and saturates to give a gyro-Bohm normalized electron electromagnetic heat flux of Q e,em /Q GB = 7.9 (see figure 5).
The corresponding global simulations, centered at x 0 = 0.5a, span a radial width of either L x = 0.15a or 0.3a, where a is the tokamak minor radius. A quadratic q-profile of the form q(x) = 1.5 + 6(x/a) 2 is considered. The radial background temperature and density profiles are of the form A s = exp[−κ As ϵ ∆A s tanh((x − x 0 )/(a∆A s ))] where A s represents the temperature or density of species s; κ n = 1, κ Ti = 0, κ Te = 4.5 and ∆n = ∆T i = ∆T e = 0.3. The numerical resolutions for the simulation with ρ ⋆ = ρ i /a = 0.004 and L x = 0.3a are N x × N y × N z × N v ∥ × N µ = 128 × 36 × 16 × 36 × 16. Krook heat and particle sources (see [30] for details) are also employed with a source rate of γ h = γ p = 0.015v th,i /R, and with radial smoothing over 0.09 a applied so these operators maintain the global-scale profiles, but do not damp finer-scale corrugations. The result is a radially smooth heating profile, and thus a radially smooth quasi-steady state heat flux, as would be expected in experiment, and consistent with local simulations. Doubling or halving this source rate is found to have little effect on the time-averaged density and temperature profiles, and the heat flux-levels change only by 20% at most. This is unlike the global simulations of microtearing of [9] that had heat fluxes that were sensitive to the source level, and had very radially peaked heat fluxes near low-order rationals.

T e flattening at low-order rational surfaces
Modes at a specific toroidal mode number k y create magnetic islands around the resonantly driven current layers at their respective mode rational surfaces (MRSs). Note that the distance between MRSs for a given k y is 1/(ŝk y ) in flux-tube simulations. The MRS of all k y radially align at the lowest-order mode rational surface (LMRS), where the magnetic islands can persist even in the turbulent phase. For the standard nonlinear flux-tube simulation, this can be seen at the LMRSs  [31,32]. Each color denotes an individual field line. Away from the LMRSs, the MRSs of each k y are radially misaligned and the overlapping magnetic islands give rise to ergodic regions.
As the electrons move swiftly along the parallel direction following the perturbed magnetic field associated with the islands at the low-order MRSs, they also undergo periodic radial excursions. This leads to a short-circuit of the perturbed T e profile, leading to its flattening. This can be seen in figure 2(a), where the green curve denoting the time-averaged effective temperature gradient ω eff Te is plotted as a function of the radial coordinate for the standard nonlinear simulation. ω eff Te is defined as the sum of the contributions from the background temperature gradient and the time-averaged zonal perturbed temperature gradient, i.e.
The perturbed temperature is defined as where δf e is the perturbed electron distribution function, and the flux-surface average, denoted by ⟨·⟩ yz , extracts the zonal part. Time-averaged ω eff Te in global simulations (both L x = 0.15a and 0.3 a) too show similar flattenings at low-order rational surfaces, as shown in figure 2(b) for the runs with L x = 0.3a. Note that for the two low-ρ ⋆ global simulations having L x = 0.3a, the time-average is taken over the initial saturated state before they 'run away'; more on this in section 6.
One may also understand the temperature flattening as a consequence of turbulence self-interaction-a mechanism where modes that are significantly extended along the field line 'bite their tails' at the rational surfaces [33,34]. In the case of microtearing modes, the parallel electron heat current density q e,∥ =´v 3 ∥ δf e d 3 v which is extended along the field line, interacts with the A ∥ of the same eigenmode, to drive zonal parallel electron temperature perturbations ⟨δT e,∥ ⟩ yz , leading to its flattening at MRSs. Here, δT e,∥ is defined as δT e,∥ = (m e /n 0 )´v 2 ∥ δf e d 3 v for convenience. Taking the v 2 ∥ moment and the flux-surface average of the gyrokinetic Vlasov equation, one arrives at an equation for the time evolution of the zonal δT e,∥ . Since we are interested in the microtearing mode, only the electromagnetic (∝ A ∥ ) nonlinear term need to be considered. Furthermore, the gyro-average over A ∥ may be ignored given that it has a radially broad structure as shown in figure 3. These approximations helps one to obtain the relation where the constant C = B 0 /|∇x × ∇y|. The linear structures of q e,∥,ky andÂ ∥,ky for k y ρ i = 0.04 are plotted with dashed lines in figure 3. The product of the two, proportional to a linear heat flux contribution, drives a zonal δT e,∥ that leads to the flattening of the parallel electron temperature at each MRS. The same process repeats for the perpendicular electron temperature. However, note a significant broadening of the timeaveragedq e,∥,ky in nonlinear simulation, also shown in figure 3 with a solid red line. A detailed description of this nonlinear broadening mechanism is given in [34,35] and can be summarized as follows. The radially narrow linear eigenmode structures lead to extended tails in k x −Fourier space and in ballooning representation (called 'giant tails' [36]). However, in a nonlinear simulation, only the first few linearly coupled k x -Fourier modes starting from k x = 0 of the eigenmode are able to retain their linear characteristics, i.e. their high amplitudes and relative phase differences with the k x = 0 mode, whereas the Fourier modes further away in the tail undergo a significant reduction in their amplitudes as a result of nonlinear couplings, implying a broadening in real space. The width of the flattened electron temperature is therefore also broadened.

Microtearing stability with corrugated background gradients
Now, we consider the linear stability of microtearing modes when the effective electron temperature gradient ω eff Te has local flattenings at LMRSs. This is equivalent to the tertiary instability analysis of zonal flows [37], except with a fixed temperature corrugation rather than a zonal flow pattern.
We perform the same nonlinear local simulation but with only two k y -the zonal k y = 0 mode and an unstable microtearing mode k y = 0.02ρ −1 i . The zonal mode is reinitialized at each time-step such that the effective gradient ω eff Te remains constant in time and has local flattenings at LMRSs as shown by the dashed magenta curves in figure 4(a), and resembling the time-averaged ω eff Te in the standard nonlinear simulation shown by the green curve. The k y ρ i = 0.02 mode is then initialized to a low seed-level amplitude and let to evolve in time atop the background of the fixed effective gradient. The simulation is then repeated for different values of ω eff Te at the MRS, and the growth rate of the microtearing mode is measured for each case.
Since the resonant current drive leading to the microtearing instability is also localized at the MRS, we expect the growth rate of the microtearing modes considered in these tertiary instability simulations to be set mostly by the effective gradient ω eff Te,MRS at the MRS, i.e. the temperature gradient away from MRS is of little significance. This is verified in figure 4(b) by the close match between the growth rate obtained from the tertiary instability simulations plotted as a function of ω eff  in standard nonlinear simulation is set by the critical gradient of the instability. For k y > k y,min , while the modes at LMRSs are made almost fully stable, those at MRSs away from the LMRSs, with lesser flattenings, are made less stable. In general, by reducing the local drive of microtearing modes at the rational surfaces, the system saturates to a state with lesser flux.

Removing zonal modulations
To further investigate the role of electron temperature flattening on saturation, a nonlinear flux-tube simulation is run while eliminating any local modifications to the temperature gradient. This is achieved by redefining the zonal component of the electron distribution function as where f M,e is the electron background distribution function-a homogeneous local Maxwellian. Note that K(x) is only a function of x and is set at each time-step such that ⟨δT e ⟩ yz = 0, and therefore ω eff Te = R/L Te throughout the simulation. The heat flux Q e,em /Q GB = 28.1 in this simulation is many times higher than Q e,em /Q GB = 7.9 in the original standard nonlinear simulation, as shown in figure 5, confirming that electron temperature flattening indeed plays a significant role in turbulence suppression.
Another way to reduce the electron temperature flattening is by weakening the self-interaction process by increasing the parallel length L z = 2π N pol of the simulation volume [33,34,38], where N pol indicates the number of times the fluxtube wraps around poloidally before connecting back to itself. Given that self-interaction essentially results from 'modes biting their tails', to correctly capture its effects, the domain length along a field line at a rational surface must be correctly captured by modelling the full flux-surface, and hence N pol and the minimum toroidal mode number n min both need to be set to 1 [33]. Therefore, by increasing N pol , we are unphysically weakening the self-interaction process and the resulting electron temperature flattening. Doubling N pol (which also doubles the density of LMRSs) is found to weaken the temperature flattening from ≈70% to ≈20% at LMRSs, leading to an increase in the flux level as shown in figure 5.
While these results confirm that the local flattening of electron temperature is crucial for correctly predicting the saturated turbulent state, the fact that these simulations, either with fully eliminated or weakened electron temperature flattenings, did saturate, indicates the presence of other, less dominant, saturation mechanism(s). Deleting the zonal electrostatic potential Φ or the zonal A ∥ in simulations changes the flux levels at most by 12%, implying that zonal flows and fields do not play a significant role in saturation in the case considered. For similar parameters, [17] shows that the free-energy flow to short wavelengths could be another saturation mechanism.

Effect of system-size
Given that microtearing turbulence saturation via temperature flattening happens primarily at rational surfaces, the separation distance between them and the corresponding finite systemsize effect is crucial. To study this in detail, we first perform two system-size scans using global simulations, choosing radial width L x = 0.15a for one scan and L x = 0.3a for the other. As mentioned in the previous section, to correctly capture the effects of self-interaction including temperature flattening in these simulations, the minimum toroidal mode number n min is set to 1.
In general, the heat flux increases with increasing systemsize as shown in figure 5. In the later part of this section, we explain why this is consistent with a saturation mechanism that depends on zonal temperature flattening at MRSs. Note that two sets of global simulations have different flux levels, and both are also below the local simulation flux levels. This difference cannot be explained in terms of the zonal-temperaturegradient flattening mechanism explored here, so suggests there is an additional system-size effect dependent on how the radial domain is treated.
For the global simulation with L x = 0.3a at the two larger system sizes denoted by open squares, the fluxes saturate to a turbulent steady state for at least a duration of 400R/v th,i , after which they undergo a 'runaway' similar to what has been reported in [9]. The runaway persists even when the grid resolutions are increased and the time-step is decreased, which suggests the origin of these may be physical. Inspecting the explosive behavior in detail, we find that a very localized intense process occurs at the current sheet of the magnetic island at the LMRS. In both local simulations, and simulations with a narrower global domain, runaway is not seen, indicating that this process is somehow sensitive to boundary conditions or spatial background nonuniformity.
To isolate the cause of this system-size dependence, i.e. the increase in heat flux with increasing system-size, we use flux-tube simulations. Note that by setting N pol = 1 and by choosing k y,min ρ i = n min q 0 (a/x 0 )ρ ⋆ corresponding to the fundamental mode (n min = 1) of the tokamak, one can simulate the full flux-surface of interest in the tokamak using flux-tube. This helps to capture the correct parallel connection length along the magnetic field lines when the LMRSs correspond to an integer rational surface. Furthermore, the radial distance between rational surfaces is also correctly captured. See figures 6(b) and (c) for comparing the distance between (L)MRSs in flux-tube and global simulations respectively. Given that self-interaction is essentially 'modes biting their tails' at rational surfaces, the flux-tube framework offers the possibility to accurately capture the effect of self-interaction, while neglecting other finite ρ * effects such as profile shearing (see [33,34] for more details). Given that the separation distance between the LMRSs is 1/(ŝk y,min ) in flux-tube simulations, we can study this particular system-size effect through a scan in k y,min ρ i . Both L x /ρ i and k y,max ρ i are kept fixed in this scan.
As k y,min is decreased, the radial density of regions with flattened electron temperature (see figure 2(a)), and hence weaker linear drive, at low-order MRSs decreases. Concurrently, flux increases, as shown by the blue asterisks in figure 5. That is, the temperature flattening mechanism becomes less effective in large systems. For the k y,min ρ i = 0.02 case, the k y ρ i = 0.04 mode contributing most to the flux has six MRSs, three of which at LMRSs experience ∼70% flattening, and the other three at second-order MRSs experience ∼10% flattening. Whereas for k y,min ρ i = 0.04, the k y ρ i = 0.04 mode sees a ∼70% temperature flattening at every MRS, so mode stabilization is much more effective. When the electron temperature flattenings are eliminated, there is still some non-gyro-Bohm scaling (black markers in figure 5), but this is less consistent.
We suggest a crude model to understand the increase in flux with increasing system-size that persists in the local limit. In the turbulent steady state, when the electron heat flux Q e becomes radially constant, one defines the pointwise diffusivity via χ e ≡ Q e / (dT e /dx) and the radial average ⟨ dT e dx so ⟨1/χ e ⟩ −1 x is the effective average diffusivity. This analysis also holds in the source-free region of a tokamak or global simulation. In the local limit, boundary conditions impose zero average temperature fluctuation, thus Q e = dT e,0 /dx⟨1/χ e ⟩ −1 x . Microtearing modes are now modelled to lead to regions of high diffusivity near each MRS, which reinforce at LMRSs, resulting in the temperature-gradient corrugations seen in simulations.
The set of 'relevant' MRSs, i.e. associated with all the modes contributing significant flux, becomes more radially dense with decreasing k y,min , because the number of toroidal modes increases, and each mode has associated rationals separated by 1/ŝk y . This is illustrated in figure 6. In figure 6(a), k y −spectra of the electron electromagnetic heat flux for the flux-tube simulations are plotted for the three considered values of k y,min .Q e,em,ky is defined such that the total flux Q e,em = ∑ kyQ e,em,ky dk y , where dk y = k y,min ρ i . In figure 6(b), the positions of all MRSs associated with k y ρ i ⩽ 0.1 (i.e. contributing significant flux) are marked for each of the three simulations, clearly indicating how the MRSs become radially dense with decreasing k y,min . The opposite is true for LMRSs (the MRSs common to all k y s and separated by 1/ŝk y,min ), denoted by thicker markers in figure 6(b), which become more spaced apart with decreasing k y,min . For comparison, the corresponding plot for global simulations with L x = 0.3a is shown in figure 6(c).
As the harmonic mean of diffusivity sets flux levels, concentrating the diffusivity at widely-spaced MRSs (large k y,min , small system-size) leads to lower flux than distributing it more evenly at a larger number of closely-spaced MRSs (low k y,min , large system-size). In an extreme limit, microtearing creates infinite local diffusivity and completely flattens gradients near each MRS, but elsewhere the diffusivity is a small constant χ b . The effective average diffusivity, crudely assuming no overlap between flattened regions, is χ b /(1 − WN), where W is the proportion of the radius flattened by each toroidal mode and N is the number of toroidal modes. This leads to a scaling Q e ∝ 1/(1 − w/ρ ⋆ ) where w a dimensionless small parameter that depends on parameters other than ρ ⋆ ; note that the flux has a singularity at small enough ρ ⋆ , where the flux explodes. This might be related to the increase in flux with system size seen (figure 5) in our MT simulations. This is also analogous to the avalanche transport arising when transport windows caused by fast-particle-driven modes overlap across much of the tokamak [39].
The width W of the high-transport region was assumed fixed in this simple picture, but actually may increase with higher flux, and the larger overlap may cause a runaway situation; this may be tied to failure to reach saturation in certain microtearing simulations.
Apart from system-size scaling, another possible consequence of electron temperature flattenings and magnetic islands at low-order rational surfaces is the potential to seed the growth of NTMs [40]. The possibility for microturbulence to excite NTMs via nonlinear coupling has been demonstrated in the past [41]. Furthermore, to experimentally verify our results, one may measure the electron temperature and look for flattenings near rational surfaces, similar to previous investigations of ITG turbulence [42]. One may also be able to measure low toroidal mode number magnetic perturbations associated with the microtearing islands in external magnetic coils; for instance, the radial magnetic perturbation associated with the islands is (δB x /B 0 )/(ρ i /R) ≃ 0.14 in the standard nonlinear simulation.

Conclusions
In conclusion, the fast motion of electrons across the magnetic islands at the LMRSs short-circuit the electron temperature, resulting in local electron temperature flattening, which then decreases the local linear drive of microtearing modes and allows lower saturated transport levels than when this process is artificially suppressed. The spacing and width of the lowconfinement regions near low-order rationals are crucial, and this provides a pathway to understand microtearing saturation (or lack thereof); one direct consequence is that microtearing turbulent transport and its study are more important in larger future devices than previously thought.
Note that microtearing saturation mechanisms that depend on zonal corrugations of other quantities, such as zonal fields or flows at rationals, would also be expected to show the same non-gyro-Bohm dependence, due to the same densityof-rationals argument.
We note that there are various important unanswered questions about size-scaling and global effects for microtearing transport. For example, the transport level is also found to be sensitive to the nature (global versus local) and radial extent of the simulation domain. This is perhaps not surprising given that the field perturbations of microtearing modes are radially extended, even though they cause localized transport. Also, we do not know how the runaway process occurs, except that it is the result of a localized intense event in the current sheet around a magnetic island. Only some of the simulations here seem subject to this process, and this also provides some clues.