Quasilinear modeling of heat flux from microtearing turbulence

The ability to obtain reliable predictions of turbulent transport in magnetically confined fusion plasmas is essential for efficient operation of present and the design of future devices. Gyrokinetic¹ simulations provide a high-fidelity foundation that incorporates the essential turbulence physics, but this approach is computationally intensive, commonly requiring 10⁴ or more CPU hours for a single nonlinear simulation. Toward the other end of the spectrum lie quasilinear (QL) solvers, where typically, the linear gyrofluid² or reduced linear gyrokinetic³ equations are solved, allowing for the fast (CPU minutes) predictions of how turbulent fluxes scale with varying input parameters. In modern variants, this technique is combined with neural networks to reduce computation time to far below a CPU second,^4–7 close to the speed necessary for real-time control and simultaneously offering a new path toward optimization of configurations for reduced turbulent transport.^8,9

At present, quasilinear (QL) transport models tend to focus on electrostatic, i.e., $E×B$ advective flux driven by electrostatic microinstabilities such as the ion-temperature-gradient (ITG) instability or trapped-electron modes (TEMs), which may be unstable at zero normalized electron pressure $β=8πne0Te0/B02$⁠. Here, $ne0$ and $Te0$ are the background electron density and temperature, respectively, and B₀ is the magnetic guide field. Inherently, magnetic instabilities, in particular, kinetic ballooning modes (KBMs) and the microtearing (MT) instability, require self-consistent treatment of shear-magnetic fluctuations $A∥$ at finite β. While some QL treatments account for the electrostatic fluxes produced by KBM turbulence, no similar proven model exists for MT turbulence.

Microtearing has been shown to arise in a variety of physical scenarios: at low toroidal wavenumbers k_y, it is occasionally observed in tokamak core plasma simulations,¹⁰ whereas it plays a more prominent role in the H-mode pedestal, where it often accounts for electron heat losses essential to pedestal evolution.^11–13 In DIII-D discharges with a high bootstrap-current fraction, MT turbulence was found to contribute.¹⁴ MT modes are routinely unstable in spherical tokamaks,^15–18 and reversed-field pinches may similarly see MT destabilized.^19,20 A practically inevitable form of MT activity, beyond pure turbulence regimes dominated by this instability, is the nonlinear excitation of a linearly stable MT mode,²¹ which can produce magnetic flutter transport in otherwise electrostatic turbulence regimes such as ITG.

Outside of these situations, experimental relevance of MT-driven turbulence and transport is still the subject of ongoing research. In particular, the core region of tokamak plasmas commonly includes fast-ion species, either from beam injection or, particularly in future devices, fusion alphas. Conceivably, the MT instability may be affected by their presence, as tearing physics undergoes stabilization.²² 

A rich trove exists in the literature on analytical theory of microtearing, beginning with its first appearance in 1975.²³ A variety of mode branches have thus been identified, in particular both (semi-) collisional^24–30 and collisionless MT,^20,31–33 with branches that do or do not rely on electrostatic fluctuations $Φ$⁠,³⁴ and MT driven by curvature^19,35,36 or the time-dependent thermal force.³⁷ In all cases, however, the MT instability requires finite values of β and of the electron temperature gradient. MT mode structures exhibit tearing parity, where $A∥$ is an even function along the guide field and $Φ$ is odd. This characteristic alone, however, is insufficient, as other instabilities such as ITG and TEM exist in tearing-parity states.^38,39

Theories that use QL techniques to gauge the efficacy of MT transport have been proposed (see, e.g., Refs. 40 and 41). So far, however, evidence of predictiveness of these models has been elusive. The present effort aims to design a QL MT model and demonstrate its predictive capability by comparing its parametric dependencies with those of gyrokinetic simulations of MT turbulence.

This paper is structured as follows: In Sec. II, a heuristic QL model for MT transport is proposed. Subsequently, this model is tested against gyrokinetic results in Sec. III by evaluating its terms using linear simulations and then comparing its transport predictions to those from nonlinear runs. The paper then concludes by summarizing key findings.

In strongly magnetized plasmas, the Larmor motion of particles about field lines resembles closed circles that flow along the guide field and drift across it. This description is the basis for the gyrokinetic framework, which averages over the gyromotion, thereby eliminating the gyrophase coordinate and ordering out fast processes on the order of the gyrofrequency.¹ The corresponding computational efficiency gain reaches many orders of magnitude.

The present effort is based on simulations with the gyrokinetic turbulence code Gene;⁴² for the implemented equations and normalization, see Ref. 43. Here, the code is used to solve the Vlasov and field equations in $ŝ$-α flux-tube geometry assuming concentric flux surfaces (α = 0), where information of both pressure and the geometric profile is reduced to zeroth- and first-order terms: The background ion (electron) temperature is $Ti0$ (⁠$Te0$⁠), and its normalized radial gradient is $ωTi$ (⁠$ωTe$⁠), defined as the ratio of major radius R₀ to the gradient scale length $LTi$ (⁠$LTe$⁠). The radial variation of the safety factor q₀ is given by the normalized magnetic shear $ŝ=(r0/q0)dq/dr$⁠, where r₀ denotes the minor radius of the flux surface of interest. The local coordinate system is given by the radial, toroidal, and parallel (to the guide field) coordinates x = r, y, and z, respectively; in Gene, the first two are treated spectrally.

In a flux tube, the heat flux, normalized to $ne0Te0csρs2/R02$⁠, of species j is defined as

Commonly, the heat flux is split into an electrostatic—sometimes referred to as $E×B$—component (denoted by superscript es), corresponding to the $Φ$ term, and a magnetic—sometimes referred to as flutter—component (denoted by superscript em) corresponding to the $A∥$ term. Overbars denote gyroaverages that, for electrostatic and shear-magnetic fluctuations, take the form of J₀ Bessel functions, $k⊥$ is the perpendicular wavenumber, and the velocity space is spanned by the parallel velocity $v∥$ and the magnetic moment μ. $cs$ and $ρs$ refer to the ion sound speed and gyroradius, respectively. Angular brackets $⟨…⟩z$ indicate an average along the field line from $−π$ to π.

The heat flux may be related to the heat diffusivity, electrostatic or electromagnetic, via

where χ is normalized to $csρs2/R0$⁠. Using a mixing-length estimate $χ∝γ/k⊥2$⁠, the heat diffusivity can be related to the linear growth rate γ and a characteristic length scale—this is the most basic form of a quasilinear transport model. A more general approach leads to the expression

where η is an eigenmode label that covers all unstable modes at a given wavenumber and $C(ky)$ is a model constant, which is obtained by matching the nonlinear flux spectrum at one point in parameter space—other approaches use spectral shape functions instead. $wes(η,ky)=Qjes,lin(η)/|Φ(η)|2$ is the quasilinear weight, relating the volume-averaged heat flux and volume-averaged electrostatic potential of the linear eigenmode. The weighted flux-surface average is defined through integrals in the ballooning angle θ as

where integrals span ballooning angles from $−3π$ to $3π$ to avoid undue influence from artificially fine structures around the resonant surface.⁴⁴ 

The model given in Eq. (3) and similar variants have been used successfully for many physical situations, including $ωTe$ scalings in tokamak core TEM turbulence,^45,46 $ωTi$ scalings in ITG turbulence in low-shear stellarators,⁴⁷ and strongly driven turbulence.³⁹ Certain cases require additional physics models, such as low-shear tokamak geometry.⁴⁸ An equivalent approach is possible to compute quasilinear particle fluxes.^49–52 

Integrated modeling of fusion plasma discharges requires rapid calculation of fluxes; to this end, sophisticated solvers have been designed which use reduced models to obtain linear eigenmodes.^53–56 They have contributed key results to the literature, including comparing simulations with experimental impurity transport on Tore Supra,⁵⁷ predicting transport in ITER,^58,59 connecting velocity reversals to a transition between turbulence regimes,⁶⁰ comparing transport in JET and JT-60U,⁶¹ and evaluating steady-state feasibility on DIII-D.⁶² 

The above applies to modeling of electrostatic flux $Qes$⁠. Adjustments must be made in order to construct a model that recovers the dominant magnetic electron heat flux component of MT turbulence, where depending on the MT flavor in question, $Φ$ may or may not be affecting the instability. Thus, one may heuristically change the definition of the averaged perpendicular wavenumber given in Eq. (4) to

where the ballooning angle is still integrated from $−3π$ to $3π$ unless stated otherwise. The quasilinear weight definition is altered to read

Here, n and $u∥$ denote the perturbed density and perturbed parallel flow, respectively. Thus, the complete model on which the results of this paper are based reads

As this expression effectively models the magnetic flutter transport described in Ref. 63, it is helpful to compare with their approach (compare Refs. 64–66),

Here, $Bx=ikyA∥$ is the radial magnetic field perturbation and $T∥e$ is the fluctuation of the parallel electron temperature. The last term, proportional to $ωTeBx2$⁠, tends to dominate away from the electrostatic limit.⁶⁷ Thus, it is expected that the fundamental gradient and β scalings are consistent between this model and Eq. (7), as the quasilinear approximation assumes γ to be a proxy for the squared amplitude, in Eq. (8) the magnetic fluctuation $Bx2$⁠.

A more fundamental electromagnetic quasilinear theory will require a separate effort and may need to address the different contributions to flutter transport from particles with different energies. Another possible improvement to the model in Eq. (7), which is beyond the scope of the present effort, is the inclusion of a triplet correlation time⁶⁸ that captures changes in saturation efficiency; this technique has been used successfully to capture finite-β transport trends.^69,70

To test the QL MT transport model, the MT scenario from Ref. 10 is used as a baseline. Physical input parameters are listed in Table I. Note that when converted to the units of the linear growth rate and frequency, $cs/R0$⁠, the collision frequency $νc=2.3031×10−5 ln Λ ne0R0/Te02=0.005$ becomes $νei=0.857$⁠. Here, $ln Λ$ is the Coulomb logarithm. In this expression for $νc$⁠, and only there, length scales are assumed to be in $m$ and temperatures in $keV$⁠.

With these parameters, there is generally one unstable eigenmode of MT type, and thus, the η sum in Eq. (7) can be omitted. Mode structures are of tearing parity, with $Φ$ extending much farther along the ballooning angle—i.e., to higher $|kx|$—than $A∥$⁠, which is mostly localized to $−π<θ<π$⁠. While at these parameters, the instability survives in the collisionless limit $νc=0$⁠, at the values studied here, collisions contribute to MT drive.

Figure 1 shows the linear growth rate and frequency spectra for different collisionalities. Little shift of the k_y location of the peak γ occurs with increasing $νc$⁠, while the spectrum broadens. The real frequency ω is largely unaffected by $νc$⁠. At the highest value of $νc=0.007$⁠, one begins to approach the fully collisional regime, where collisions no longer add additional drive and eventually act in a stabilizing fashion.

The equivalent data for variations in electron temperature gradient $ωTe$ are shown in Fig. 2. While $ωTe$ destabilized the MT mode in a manner similar to $νc$⁠, a moderate shift toward high k_y with increasing $ωTe$ is visible, and ω scales with the gradient. While the lowest value $ωTe=2.5$ yields finite growth, some degree of proximity to the critical gradient is evident.

As a third parameter essential to MT physics, the scaling with normalized pressure β, is presented in Fig. 3. Again, increasing β drives the mode and produces a spectral shift of γ to higher k_y, with no visible impact on ω.

Quasilinear fluxes can now be constructed on the basis of these linear simulations, matching nonlinear fluxes—by choosing appropriate $C(ky)$ as model constants—for each scan at the base case defined in Table I. More specifically, $C(ky)$ is scaled to match the quasilinear and nonlinear flux spectra at each wavenumber, resulting in identical nonlinear and quasilinear integrated fluxes at the base case for each scan.

For nonlinear simulations, typical numerical resolutions are 192 radial modes, 24 toroidal modes, 16 parallel, 32 parallel velocity, and eight magnetic moment grid points, with the perpendicular box spanning $150ρs$ radially and $314.159ρs$ toroidally. Fourth-order hyperdiffusion settings⁷¹ are D_z = 5 for parallel and $Dv=0.2$ for parallel velocity hyperdiffusion coefficients.

Well-behaved MT turbulence can be obtained with these settings. Fig. 4 shows the corresponding heat flux spectra: flux peaks at a rather low $ky=0.06$⁠. The nonlinear cross phases of $A∥×q∥e$ and $A∥×T∥e$ closely match those in linear simulations throughout the transport-relevant k_y range, as required for the validity of QL models. Here, $q∥e$ denotes the parallel electron heat flux.

Regarding the prominence of zonal features in saturation, both zonal flows⁷² and zonal fields—i.e., zonal $A∥$—form intermittently and transiently throughout the quasi-stationary state at the default parameter settings. This suggests that a difference in the critical gradient between linear and nonlinear simulations may exist; compare Ref. 73 in the case of zonal-flow-mediated ITG turbulence, which will be discussed below, which would be expected to cause QL predictions to break down below the nonlinear critical gradient.

Figure 5 contains both the nonlinear and quasilinear fluxes over a scan of $νc$⁠, matched at the default 0.005. Clearly, good agreement is observed between the turbulent fluxes and the model predictions obtained from Eq. (7). The only data points that are overpredicted by the QL model lie at the lowest collisionality, where fluxes approach marginality, in line with the expectation of model limitations.

A similar picture emerges when varying the temperature gradient $ωTe$ (see Fig. 6). Again, QL fluxes track nonlinear fluxes, while at the lowest gradient $ωTe=2.5$⁠, the QL model overpredicts the turbulent flux, which has become marginal. As discussed above, this is due to simulations approaching the nonlinear critical gradient, at which point a standard QL model such as the one described here can no longer predict fluxes.

Note that some variability of the nonlinear data in Fig. 6, compared with the more monotonic behavior in, e.g., Fig. 5, is caused by shorter simulation run times: for the $ωTe$ scan, nonlinear data was taken directly from Ref. 10, whereas the above collisionality scan is based on new simulations with better temporal statistics.

Finally, Fig. 7 shows nonlinear and quasilinear data for a scan in β. In the above cases, where $νc$ or $ωTe$ was scanned, the results were essentially unaffected by the ballooning angle integration range in Eq. (5). This changes for the β scan, where a separate curve indicates the results obtained with $⟨k⊥2⟩$ integration from $−13π$ to $13π$⁠, where the nonlinear scaling is matched more accurately than for the default $−3π$ to $3π$⁠. The discussion whether layers around the resonant surface matter and, consequently, which integration range should be used, will be deferred to future work. Regardless, both curves capture the comparatively weak dependence of heat flux on β, demonstrating that the important scaling of magnetic transport with β is included consistently in the model, accurately capturing the physics discussed in Ref. 63.

Below, the findings of this paper are summarized.

A quasilinear transport model for microtearing turbulence has been proposed which relies on the magnetic potential eigenfunction and a mixed quasilinear weight with both electrostatic and electromagnetic proxies. This QL model, written out in Eqs. (7) and (5), was then tested against nonlinear simulations of MT turbulence in an idealized tokamak core plasma scenario. As long as turbulence remains sufficiently above the nonlinear critical gradient, flux dependencies on the collision frequency, the electron temperature gradient, and the normalized electron pressure are all recovered quantitatively. These results lay the foundation for efficient computation of MT fluxes, with applications including the inter-ELM (edge-localized modes) evolution in the tokamak pedestal.

Various improvements on this approach can be envisioned, chief among them a fundamental derivation of an electromagnetic quasilinear theory that accounts for MT. Additionally, changes in saturation efficiency in zonal-flow- or zonal-field-mediated MT turbulence may be tracked by including a triplet correlation time factor in the model. Furthermore, future research will need to determine whether other flavors of MT instability produce flux scalings that are captured similarly well by the present model. Beyond pure MT turbulence, mixed regimes where other ion-scale instabilities affect the dynamics⁷⁴ will need to be considered, and reduced models be designed, which account for mode-interaction. Similarly, interaction of microtearing turbulence with electron-scale turbulence will require attention.^75,76

Once an encompassing QL model for MT turbulence has emerged, a fast solver will be required which allows its deployment in integrated modeling and real-time control.

The authors wish to thank M. Hamed and J. Citrin for helpful discussions and H. Doerk for making simulation data available for analysis. This work was supported by the U.S. Department of Energy, Office of Science, under Grant No. DE-FG02–04ER-54742, and the Office of Fusion Energy Sciences Scientific Discovery through Advanced Computing (SciDAC0) program, under Award No. DE-SC0018429.

The data that support the findings of this study are available from the corresponding author upon reasonable request.