The Problem of Marginality in Model Reductions of Turbulence

Fusion plasmas display the property, common in dynamical systems that upon surpassing a critical threshold, an instability may promptly build up, inducing large fluxes which deplete the driving gradients and inhibit the instability. Background gradients thus hover in the vicinity of nonlinear near marginal thresholds [1]. Many strategies have been devised in modelling to mimic natural processes. All are not equivalent and different choices may critically affect the nature of computed statistical equilibria.

Forcing-dependent steady-states have indeed long been observed in a variety of systems. Systems with long-range interactions, either controlled in energy or in temperature display different equilibria [2, 3]. Swirling flows controlled through either imposed torque or imposed velocity display distinct steady states, as well as different dynamical regimes [4]. In magnetised fusion plasmas, auxiliary heating and current drive force the system out of equilibrium through a constant flux. Mimicking nature, ”flux-driven” (FD) forcing adiabatically imparts a volumetric flux to the system whose gradients self-consistently adapt, leading to the observation of large-scale transport events such as avalanches [5, 6, 7, 8, 9, 10, 11], nonlinear structures such as zonal mean flows [12, 13] or secondary patterns such as staircases [14, 15, 16, 17, 18, 19]. Such observations are absent or impaired when the system is driven through a body force, which amounts to imposing fixed mean gradients. This ”gradient-driven” (GD) strategy is widely used in direct numerical computations for it is computationally efficient. It indeed enforces, with respect to the FD approach, spatial and temporal scale separations between equilibrium and fluctuations and solves for the fluctuations only.

Known differences between FD and GD frameworks have been documented [20, 17, 18]. Whether these are of practical incidence in fusion-relevant configurations is non trivial. The matter is important for there are increasing requirements for fast, reduced, yet reliable models to explore the vast parameter space of magnetised plasma turbulence, interpret experimental results and forecast future large experiments such as ITER. Currently, the more advanced reduced models are based on quasilinear theory (QLT) [21, 22, 23]. With the advent of machine learning techniques, the ubiquitous closure problem of QLT is approached through data-driven techniques that use large-scale databases of first-principles GD computations. Systematic shortcomings, if any, within reference GD strategies are thus likely to be carried over to the reduced models. Comprehensive understanding of discrepancies between first-principles FD, GD and quasilinear approaches is thus important and timely, and the topic of this Letter.

To this end, we confront reference results from nonlinear FD gyrokinetics using the Gysela framework [24] to state-of-the-art GD nonlinear gyrokinetics and GD quasilinear calculations, using respectively Gkw [25] and QuaLiKiz [26, 27]. We further complement the study with confrontation with the QL flux-driven QuaLiKiz–Jetto framework. Several instances of QuaLiKiz locally compute at various radii flux-surface-averaged transport coefficients which are passed on to the Jetto integrated modelling suite [28] and used to evolve profiles through flux-driven transport equations. After a transport timescale, new profiles are fed to QuaLiKiz and the process loops.

Main results are: (i) steady-state predictions of fluxes moderately depend on the nature of forcing well above nonlinear threshold; (ii) near marginality however, nonlinear and quasilinear GD models sizeably underpredict turbulent heat transport. Under a driving flux, profiles display ’stiffness’, i.e. hover in the vicinity of their near-threshold flux-matching values. Large, hot devices such as ITER are expected to be stiff due to the temperature dependence of the gyroBohm heat flux scaling, making near marginality a regime which models must confront. Proximity to nonlinear thresholds implies additional complexity as it favours secondary pattern formation and mesoscale organisation. Despite this additional complexity, (iii) the underlying assumptions of QLT hold well across nonlinear regimes. We show that transport underprediction rather stems from the choice of closure, i.e. the nonlinear saturation rule. This work stresses the relevance of QLT for model reductions of turbulence whilst providing guidelines whereby reduced models can be improved. The nonlinear and QL approaches tested here are the current workhorse for estimating transport and confinement in turbulent fusion plasmas. This work thus has implications for present-day experimental data analysis and scenario extrapolation for fusion production. Novel saturation rules should strive to incorporate near marginal flux-driven specificities, often dubbed turbulence spreading [29, 11, 30], transport nonlocality or staircase organisation.

Gysela resolves ion Larmor radius scale turbulence and collisional transport in global tokamak geometry with no scale separation between equilibrium and fluctuations. A centrally peaked heat source drives a deuterium plasma out of equilibrium, which converges towards a steady temperature profile $T$ on energy confinement times. Two paradigmatic simulation regimes are considered in the electrostatic regime with Boltzmann electrons. Both cases are run in the so-called “local limit” [31], at $\rho_{\ast}=\rho_{i}/a\leq 1/250$, where comparison to local Gkw is fair. ”Local” simulations exploit to numerical advantage the assumption of a scale separation between equilibrium and fluctuations and solve a limited subset of the whole plasma volume which twists around the torus due to the magnetic shear of the background magnetic equilibrium. In this approximation, local fluxes bijectively relate to local gradients. Here $\rho_{i}$ and $a$ are respectively the local Larmor radius and minor radius of the torus. The time averaged (over several tens of turbulence correlation times in each case) radial profiles of normalised temperature gradients $R/L_{T}=-R\partial_{r}T/T$ and zonal flow shear are plotted in Fig.1 for the ‘Above Threshold’ (top) and ‘Near Marginal’ (bottom) cases. The shaded areas represent temporal standard deviation. The large deviation from the mean shearing rate in the near marginal case (bottom, right axis) results from the meandering of staircases, which have already been reported to play an important role in this regime of parameters [18]. The linear (black hourglass symbols) and nonlinear (red squares) thresholds are discussed in next section. In order to broadly span parameter regimes, main plasma parameters vary significantly between cases, as illustrated in Table 1.

For each set of local values of the Gysela parameters, two thresholds are to be distinguished: linear $R/L_{T}^{\text{lin}}$ represents the normalised temperature gradient above which an unstable mode grows, at vanishing $\mathbf{E}\times\mathbf{B}$ shear. The inclusion of self-generated flow shear –nonlinear in essence– introduces a second threshold $R/L_{Tc}$, indicative of the nonlinear saturation of turbulence [32], in particular by zonal flows. Practically, $R/L_{Tc}$ is computed in local GD frameworks as the minimum local temperature gradient which provides a non-vanishing nonlinear heat flux. It is estimated with Gkw through a series of increasing $R/L_{T}$ nonlinear computations, whilst imposing local $\mathbf{E}\times\mathbf{B}$ shear values from flux-driven Gysela. We call “above threshold” the regime such that $R/L_{T}^{\text{lin}}<R/L_{Tc}\ll R/L_{T}$, characterised by a static, smooth zonal mean shear $\gamma_{E}=\partial_{r}^{2}\langle\phi\rangle/B$ and subdominant zonal fluctuations. The flux-surface average of the electric potential is $\langle\phi\rangle$, $\tilde{\phi}$ are its fluctuations and $R$ is the tokamak major radius. We call “near marginal” regimes such that $R/L_{T}^{\text{lin}}<R/L_{T}\lessapprox R/L_{Tc}$. Proximity from below to $R/L_{Tc}$ is a hallmark of near marginality, featuring a well-defined staircase pattern of flows and associated temperature corrugations which meanders within the time interval –hence the large shear deviation– and significant avalanching activity [18].

QLT is valid [33] in the low Kubo number limit $K=\tau_{\text{jump}}/\tau_{\text{int}}<1$, ratio of a jumping time $\tau_{\text{jump}}$ of particles from one turbulent eddy to the next over a nonlinear eddy-particle interaction time $\tau_{\text{int}}$. Kubo numbers are estimated from the Gysela flux-driven data in various ways, summarised in Tab.2. With analogy to incompressible fluids, particles trapped in turbulent convective cells explore eddies in a typical turn-over time given by the local vorticity $\tau_{\text{int}}^{\text{trap}}\sim B/\langle|\nabla_{\bot}^{2}\tilde{\phi}|^{2}\rangle^{1/2}$. Whilst they undergo this vortical motion they also drift in about $\tau_{\text{jump}}^{\star}\sim L_{\theta}(eB/\nabla T)$ from one turbulent structure to the next at the typical speed of the local diamagnetic velocity. Alternatively, the slower evolution of the potential field provides a relevant correlation time $\tau_{\text{jump}}^{{\text{corr}}}$ for turbulent fluctuations, trade-off between unstable growth and nonlinear saturation. It must be computed as a Lagrangian correlation time, in the co-moving frame of the eddies to correct for the Doppler shift induced by their phase velocity. We estimate it through image registration, following the toroidal shift between 3-dimensional turbulence snapshots of $\tilde{\phi}$. It is compared to turbulence-driven stochastic transport times of particles which, assuming a diffusive ansatz for $\mathbf{E}\times\mathbf{B}$ fluctuations, drift across eddies in about $\tau_{\text{int}}^{{\text{diff}},x}=L_{x}/\langle|\tilde{v}_{E,x}|^{2}\rangle^{1/2}$, with $x=\{r,\theta\}$ and $L_{x}$ the transverse correlation lengths computed from Gysela outputs. Interestingly, employing a Eulerian correlation time would result in severe Kubo number underestimation, locally up to factors of 25.

Three Kubo numbers, combinations of the above nonlinear times are plotted in Fig.2 for both ‘Above Threshold’ and ‘Near Marginal’ regimes. The various definitions for $K$, coherent, provide the following picture: (i) on the basis of order unity Kubo numbers, QLT should be marginally valid. Yet, as shown below, key assumptions at the heart of the QLT reduction remain valid throughout nonlinear evolution, which strengthens the case for quasi-linear integrated modelling. Interestingly also, (ii) consistently larger $K$ values near marginality stress the more percolative nature of transport there. Avalanching emerges as a key theme to distinguish between ‘Above Threshold’ and ‘Near Marginal’ regimes as they likely underpin the larger $K$ values computed near marginality. It is a likely indication that incorporating avalanching and its zonal mean flow regulation may significantly alter transport predictions and improve model behaviour near marginality.

Heat fluxes are computed with QuaLiKiz and Gkw from the Gysela time-averaged steady-state profiles plotted in Fig.1 with the same degree of approximations (electrostatic & Boltzmann electrons). All codes have different normalisations, with Gysela e.g. being normalized to $\hat{Q}^{\text{Gys}}=Q/(n_{0}T_{0}v_{Ti0})$ with $v_{Ti0}=(T_{0}/m_{i})^{1/2}$ the ion thermal velocity. In Fig.3–(a) and (b), fluxes from all codes are in the units of $\hat{Q}=\hat{Q}^{\text{Gys}}/\hat{Q}^{\text{gB}}$, with $\hat{Q}^{\text{gB}}=\rho_{\star,50}^{2}\,a/R$ and values of $\rho_{\star,50}$ and $R/a$ given in Tab.1. Consistent rescaling factors have been applied to QuaLiKiz and Gkw fluxes. The stabilising effect of zonal flow shear is accounted for, locally imposing the reference Gysela  shear values. Non-axisymmetric (turbulent) contributions to heat fluxes are shown in Fig.3–(a) and (b). Turbulence spreading and profile corrugations, inherent to flux-driven complexity are absent or hindered in QL or GD approaches. A fair comparison thus requests that the Gysela reference data be significantly coarse-grained before being handed over to Gkw and QuaLiKiz as input profiles (including $\mathbf{E}\times\mathbf{B}$ shear profiles) and after computations, when comparing e.g. fluxes. Practically, Gysela observables are time averaged well over a typical linear growth time ($\geq 30,000\,\Omega_{ci}$), at steady-state and a radial smoothing is performed over $60\,\rho_{i}$ to smear out visible FD specificities in the profiles —FD turbulence correlation lengths and profile corrugations are indeed about $10\rho_{i}$ [18]. Here, $\Omega_{ci}$ denotes the ion cyclotron frequency. Sensitivity to gradient fluctuations, inherent to FD approaches is further assessed by additional scans in $R/L_{T}$ and $\mathbf{E}\times\mathbf{B}$ shear within one temporal standard deviation.

Without inclusion of $\mathbf{E}\times\mathbf{B}$ shear stabilisation [i] QuaLiKiz heat fluxes are overestimated with respect to Fig.3 by over an order of magnitude (not shown here). With the inclusion of shear, [ii] at locations of low or vanishing $\mathbf{E}\times\mathbf{B}$ shear and despite the $60\,\rho_{i}$ smoothing, gradient-driven (standalone) QuaLiKiz commonly displays [subplots (a) and (b)] variations by factors in heat fluxes from one radial position to the next whilst Gkw exhibits much less sensitivity to $\mathbf{E}\times\mathbf{B}$ shear stabilisation. Interestingly, this large sensitivity of QuaLiKiz to shear is mitigated when called within the integrated framework of Jetto [28] to allow for a flux-driven QL profile evolution, driven by a central source that mimics the one of Gysela [subplot (c)]. In the regime above threshold [iii] reasonable agreement in computed fluxes is found across fidelity hierarchy. Conversely near marginality, [iv] despite significant smoothing, heat flux discrepancies in Fig.3–(b) are well outside allowed gradient sensitivity and fluctuation ‘error’ bars. Secular growth of zonal flows in Gkw is responsible for the observed large flux underprediction. This echoes previous observations casting concern on near marginal gradient-driven predictions [17]. The soundness of separating fluctuations from mean is thus clearly questioned near marginality.

The conclusions above are further confirmed when comparing FD Gysela to the FD quasilinear framework of QuaLiKiz–Jetto, in the same two regimes. In subplot (c), profiles from QuaLiKiz–Jetto are evolved until heat fluxes match the nonlinear Gysela reference fluxes. The figure of merit now becomes how close quasilinear profiles are at flux equilibrium with those of Gysela. Remarkable [v] profile agreement above threshold echoes the agreement in fluxes, displayed in panel (a). Near marginality however, [vi] a large over-prediction of the temperature gradient is required for QuaLiKiz–Jetto to carry the same flux as Gysela. Differently stated, both GD nonlinear fluxes and QL fluxes (irrespective to the forcing framework) are underpredicted near marginality, confirming that a constitutive ingredient is missing in that regime, at least with Boltzmann electrons.

Flux discrepancies between Gysela and Gkw likely stems from disregarding the feedback of fluctuations on an assumed fixed ”equilibrium”. This enforces local and single-valued flux–gradient relations, underestimates turbulence spreading, avalanching and mesoscale organisation. All of which contribute to transport, especially near marginality. Discrepancies with QuaLiKiz may either come from violating assumptions central to QLT –linearity of fluctuations– or by inheriting shortcomings akin to those of Gkw, through the choice of saturation rules.

To disentangle these questions, we compute in Fig.4, for all three approaches, the complex argument $\alpha_{n}$ of the $n$-th Fourier component of the heat flux [panel (a)] –a proxy for the linear cross-phase between electric potential and pressure fluctuations. They depend on the toroidal wave number $n$ labelling each eigenmode which is related to the normalised poloidal wave vector $k_{\theta}\rho_{i}$ through: $k_{\theta}\rho_{i}=(a/r)\,n\,q\,\rho_{\star}$. In the present cases, the $n=25$ toroidal wave number corresponds to, at mid-radius $a/r=2$: $k_{\theta}\rho_{i}=25*1.4*2/250=0.28$ and $k_{\theta}\rho_{i}=25*1.7*2/350\approx 0.24$ for the Near Marginal and Above Threshold cases, respectively. In Gysela, $Q_{n}^{\text{Gys}}=-i|Q_{n}|\mathrm{e}^{i\alpha_{n}}=\int\hat{v}_{E,n}^{r}\frac{2}{3}\left(\hat{P}^{\ast}_{||,n}/2+\hat{P}^{\ast}_{\bot,n}\right)\mathrm{d}\theta$, where parallel and transverse components of the pressure and the radial component of the $\mathbf{E}\times\mathbf{B}$ drift are computed from 3D output data. The $n$-Fourier components of the heat flux $|Q_{n}|$ and of the squared potential $|\phi_{n}|^{2}$, i.e. the saturation rule are respectively shown in panels (c) and (d). Panel (b) displays their ratio, a proxy for the dispersion relation, sometimes called ‘QL flux integrals’.

Clearly, linear cross-phases display reasonable agreement, within factors of 2 and across all regimes. They are not responsible for the flux discrepancies. Factors of disagreement –and avenues of improvement for QL modeling– are essentially twofold: [i] $90\%$ of the heat flux is carried by modes $n\in[\![5,50]\!]$ in the above threshold regime; all approaches provide similar conclusions. Near marginality however, flux is carried in Gysela and Gkw through $n\in[\![5,70]\!]$, twice the amount of active modes with respect to QuaLiKiz. Furthermore, [ii] saturation rules are clearly responsible for much of the observed flux discrepancies. In the above threshold regime, flux spectra agree well [panel (c)] –though this results from a surprising compensation: the potential, or saturation rule [panel (d)] is under-estimated, the dispersion relation [panel (b)] is over-estimated yet heat flux spectra in the above threshold regime are in reasonable agreement. No such compensation occurs near marginality; there severe under-estimation of the potential spectrum is clearly responsible for the under-prediction of fluxes.

The success of reduced models especially hinges on the reproduction of nonlinear gyrokinetic fluxes [27]. Flux under-prediction in the dynamically important regime of near marginal stability is thus a matter of importance. At the heart of this paper lies the fact that FD and GD models provide significantly different flux predictions with Bolzmann electrons close to marginal stability $R/L_{Tc}$, which underpins basic discrepancies in how nonlinear saturation of turbulence is modelled. This should foster renewed interest in ways to complete QL or GD models near marginality. The robustness of linear features [34] across fidelity hierarchy and across turbulent regimes is encouraging perspective and provides constructive directions whereby reduced models could be improved.

Discrepancies between Gysela and Gkw are strong indications that turbulence spreading and mesoscale patterning are in fact central to accurate transport predictions near marginality. Larger Kubo numbers near marginality, as shown in Fig.2 reinforce this point, which is also differently stressed by recent works in the plasma edge [11, 30]. The fact that QuaLiKiz behaves better near marginality than Gkw is likely because the QL closure does not include the long-lived zonal flows that quench transport near marginality in GD approaches [17]. In QuaLiKiz, stable modes are neglected and the turbulent intensity spectrum –chosen as a double power law in $k_{\theta}$ for unstable modes– is fitted [35, 27] onto databases of GD nonlinear computations, similar to those presented here with Gkw. With this procedure, QL models inherit the shortcomings of the primitive GD models onto which they are adjusted. The near marginal transport shortfall in QuaLiKiz indeed largely comes from issues with the QL closure, i.e. the choice of saturation rule.

This provides directions for improvement. In physical terms, near marginal regimes require description of transport below or at linear stability and coupling to modes presently predicted as stable in QuaLiKiz. It also requires to model the self-advection (spreading) of turbulent domains and the possibility of non-monotonic flux-gradient relations. Several routes can be explored. In the spirit of current frameworks, QL models could be trained on near marginal flux-driven databases such as provided by the likes of Gysela. This would likely lead to close QL models with regime-dependent turbulent intensity spectra. Alternatively to present closures, QL models could also be coupled to dynamic equations for the turbulence intensity, e.g. in the form of reaction–diffusion [36] or $k-\epsilon$ equations [37], enriching accessible nonlinear dynamics.

Ongoing studies are concerned with further characterising our “Above Threshold” and “Near Marginal” regimes [38] when kinetic features of electron dynamics are present. This is important to assess relevance for ITER extrapolations. Electron dynamics is known to locally modify turbulence organisation near low order rational surfaces [39]. Interestingly however, key features of near marginal turbulence with Boltzmann electrons (flow patterning, shear effectiveness and staircase organisation) which are central here to the near marginal regime robustly endure in kinetic electron regimes [40]. The present work provides a framework of understanding; it will be expanded to various parameter regimes.

The authors acknowledge stimulating discussions with P.H. Diamond and participants at the 2019 and 2021 ”Festival de Théorie” in Aix-en-Provence. This work has been carried out within the framework of the EUROfusion Consortium and was supported by the EUROfusion Theory and Advanced Simulation Coordination (E-TASC) initiative under the TSVV (Theory, Simulation, Verification and Validation) “L-H transition and pedestal physics” project (2019–2020) and TSVV “Plasma Particle/Heat Exhaust: Gyrokinetic/Kinetic Edge Codes” (2021–2025). It has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement No 633053. The authors gratefully acknowledge funding from the European Commission Horizon 2020 research and innovation program under Grant Agreement No. 824158 (EoCoE-II). The views and opinions expressed herein do not necessarily reflect those of the European Commission. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958. We acknowledge PRACE for awarding us access to Joliot-Curie at GENCI@CEA, France and MareNostrum at Barcelona Supercomputing Center (BSC), Spain.