Predicting the Dimits shift through reduced mode tertiary instability analysis in a strongly driven gyrokinetic fluid limit

Abstract

The tertiary instability is believed to be important for governing magnetised plasma turbulence under conditions of strong zonal flow generation, near marginal stability. In this work, we investigate its role for a collisionless strongly driven fluid model, self-consistently derived as a limit of gyrokinetics. It is found that a region of absolute stability above the linear threshold exists, beyond which significant nonlinear transport rapidly develops. Characteristically, this range exhibits a complex pattern of transient zonal evolution before a stable profile can arise. Nevertheless, the Dimits transition itself is found to coincide with a tertiary instability threshold, so long as linear effects are included. Through a simple and readily extendable procedure, tracing its origin to St-Onge (J. Plasma Phys., vol. 83, issue 05, 2017, 905830504), the stabilising effect of the typical zonal profile can be approximated, and the accompanying reduced mode estimate is found to be in good agreement with nonlinear simulations.

1. Introduction

Experimental fusion devices exhibit significantly higher transport than neoclassical predictions. The additional anomalous transport arises as a result of gyroscale microturbulence driven by various instabilities (Liewer 1985), such as the ion temperature gradient (ITG) mode (Choi & Horton 1980; Horton, Choi & Tang 1981; Conner & Wilson 1994) or the trapped electron mode (Kadomtsev & Pogutse 1970; Nordman, Weiland & Jarmén 1990). Moreover, the turbulent transport associated with these instabilities is very stiff. Thus, once instability is present, even a small increase in the plasma gradients will drastically increase transport levels, effectively freezing the gradients in place and restricting device performance (Ryter et al. 2011). This picture is expected to continue to hold for future fusion devices, and so being able to predict when this transport threshold is reached becomes of key importance to predicting overall behaviour and performance. This is of obvious importance for the understanding and design of experiments, possibly being particularly useful for optimisation. Here, especially stellarator devices spring to mind, since they possess a large degree of freedom in their magnetic geometry (Mynick 2006).

Naively, one might expect that the transport threshold should coincide with the linear instability threshold, since, fundamentally, these extract free energy from the plasma gradients to drive the turbulent transport. However, instead, it is found that finite transport actually commences at significantly steeper gradients. This apparent discrepancy traces its origin to self-generated poloidal zonal flows (Lin 1998; Diamond et al. 2005). Once the primary drift waves reach sufficient magnitude, such flows naturally arise through nonlinear interactions in what is known as a secondary instability (Rogers, Dorland & Kotschenreuther 2000). As the zonal flows become strong enough, they can then, in turn, nonlinearly stabilise the primary instability by shearing drift waves and decreasing their correlation length (Biglari, Diamond & Terry 1990). Because the zonal flows have a Landau-undamped component (Rosenbluth & Hinton 1998) they can, close to marginal stability and in the absence of collisions, persist for such a long time that the effective transport nearly vanishes. This is known as the Dimits regime, and the effective upshift of the critical gradient, i.e. the difference between the linear critical gradient and the observed critical gradient for the onset of turbulence, is known as the Dimits shift, both after their discoverer (Dimits et al. 2000).

Despite the qualitative picture of the Dimits shift as just outlined being somewhat firmly established, there are still some key features which are poorly understood. Thus a general quantitative prediction of the Dimits shift has proven elusive. To describe the ITG turbulence typically observed in experiments, it is necessary to employ full gyrokinetics to retain all relevant physics (Catto 1978; Frieman 1982; Abel et al. 2012). This, however, is a highly complex kinetic system, and attempting to thoroughly account for all the possibly relevant features necessary for a full description of the Dimits shift has proved a daunting task. Instead, much research has been undertaken for simpler systems which are analytically tractable, typically of the Hasegawa–Mima–Wakatani family (Hasegawa & Mima 1978; Hasegawa & Wakatani 1983), in order to gain the insight necessary to parse key features which could render the gyrokinetic problem solvable.

Many different features have been observed which could prove to be of relevance for the full problem. These include, but are not limited to, coupling to subdominant modes at unstable scales (Makwana et al. 2014; Pueschel, Li & Terry 2021), time-coherent localised soliton structures known as ferdinons (van Wyk et al. 2016, 2017; Ivanov et al. 2020), zonal-drift predator–prey-type interactions (Berionni & Gürcan 2011; Kobayashi & Rogers 2012) or the ability of a turbulent momentum flux to tear down or build up a decaying zonal profile (Kim & Diamond 2002; Ivanov et al. 2020). One feature which, however, repeatedly crops up in these studies is that instability causing drift waves to arise from an initially zonally dominated state, known as the tertiary instability (Rogers et al. 2000).

Despite seemingly being a natural candidate to explain the observed Dimits shift, based on findings from simpler systems, the importance of the tertiary instability for the Dimits shift has nevertheless been a topic of debate within the literature. St-Onge (2017) and Zhu, Zhou & Dodin (2020a), for example, based accurate predictions upon it, while Li & Diamond (2018) and Ivanov et al. (2020) on the contrary reported finding it unimportant. To help rectify this confusion, in this paper we will thus attempt to shed some light on the tertiary mode in the Dimits regime, investigating its relevance for the Dimits transition in a strongly driven fluid system directly derived from gyrokinetics.

In our investigations we will find that, just as Zhu et al. (2020a) stressed, in order to properly capture the behaviour of the tertiary instability in the marginally stable regime, the linear drive cannot be neglected. The tertiary instability should not be treated as a purely shear Kelvin–Helmholtz-like (KH) instability, but instead as a modified primary instability that includes such terms. Then, the tertiary instability alone seems sufficient to encapsulate the Dimits transition for the system under consideration. This is despite the fact that this system is ostensibly similar to the one recently studied by Ivanov et al. (2020), where the opposite case was found to hold, a discrepancy arising from the present absence of collisional zonal flow damping. Finally, we will see that a reduced mode scheme to approximate the tertiary instability can yield a simple but effective prediction (within 15 %–30 %). Furthermore this scheme seems readily extendable to more complete collisionless systems, including gyrokinetics itself, which will be the subject of an upcoming publication.

This paper is outlined as follows. The strongly driven gyrofluid system will first be introduced in § 2 and its key features will then be presented in 3. Next, we will in turn describe each of the present instabilities of the primary–secondary–tertiary paradigm (see Kim & Diamond 2002), noting their effects on the system as a whole. Guided by direct simulations presented in § 4, we will then home in further on the tertiary instability in § 5. There we will show that it can be employed to arrive at a very simple Dimits shift estimate, related to the one of St-Onge (2017), which could prove to be broadly applicable for other non-collisional systems as well. Finally, we will conclude with a brief summary and discussion in § 6.

2. Basic model

The Dimits shift was originally observed in, and is of most experimental relevance for, fully gyrokinetic simulations of tokamaks (Dimits et al. 2000). However, the intrinsic kinetic nature of this system makes analytical treatment of even just the tertiary instability intractable. Investigations have therefore focused on simplified problems (see e.g. Kolesnikov & Krommes 2005; Numata, Ball & Dewar 2007), hoping to find insights which can be extrapolated to the more complete problem. Naturally, these models all fail to capture much of the physics of the full gyrokinetic system because of their simplicity, possibly raising concerns about how valid such extrapolation will be. Therefore, we will here present another self-consistently closed gyrofluid system in two spatial dimensions, in the hope that it may prove yet another useful stepping stone to solidify and clarify the emerging picture of the Dimits shift when proceeding towards the full gyrokinetic problem.

2.1. Gyrokinetics and conventions

To arrive at the system of interest one starts from the usual electrostatic collisionless gyrokinetic equation in Fourier space, which we, in the vein of Plunk et al. (2014), express in non-dimensional form as

(2.1)\begin{equation} \left( \frac{\partial}{\partial t} + \textrm{i} \tilde{ \omega }_d \right) h_{\boldsymbol{k}} + \left\{ {\varPhi} , {h} \right\}_{\boldsymbol{k}} = \left( \frac{\partial}{\partial t} + \textrm{i} \tilde{ \omega }_* \right) \varPhi_{\boldsymbol{k}} f_0. \end{equation}

Here, $f_0$ is the ion Maxwellian distribution with charge q, mass m, temperature T, and mean thermal velocity $v_T = \sqrt { 2 T / m}$, $h$ is the non-adiabatic part of the ion fluctuations $\delta f_i$. Meanwhile, the gyroaverage in Fourier space is encapsulated by the Bessel function of the first kind $\textrm {J}_0 = \textrm {J}_0 (\sqrt {2} k_\bot w_\bot )$, where the normalised velocity $w=v/v_T$ and wavenumber $k$ are split into their parallel and perpendicular components $w_\parallel , k_\parallel , w_\perp , k_\perp$ with respect to the magnetic field. The gyroaverage enters (2.1) through the gyro-averaged electrostatic potential $\varPhi _{\boldsymbol {k}}=\textrm {J}_0\varphi _{\boldsymbol {k}}$, while the Fourier space Poisson bracket, in turn, is given by

(2.2)\begin{equation} \left\{ {a} , {b} \right\}_{\boldsymbol{k}}=\sum_{\boldsymbol{k}_1,\boldsymbol{k}_2}\left(k_{1y}k_{2x}-k_{1x}k_{2y}\right)a_{\boldsymbol{k}_1}b_{\boldsymbol{k}_2}\delta_{\boldsymbol{k},\boldsymbol{k}_1+\boldsymbol{k}_2}, \end{equation}

where $\delta _{\boldsymbol {k},\boldsymbol {k}_1+\boldsymbol {k}_2}$ is the Kronecker delta and the $x$- and $y$-coordinates are the radial and poloidal coordinates, respectively. After the introduction of a reference length scale $L_{\mathrm {ref}}$, the spatial and temporal dimensions are normalised to the typical ion gyroradius $\rho$ and the streaming time $v_T / L_{\mathrm {ref}}$ respectively, so that $\varphi = q \phi L_{\mathrm {ref}}/ T \rho$ is the dimensionless electrostatic potential. Furthermore the plasma $\beta$ is assumed small so that the magnetic field $\boldsymbol {B}=B\boldsymbol {b}$ satisfies $\boldsymbol {\nabla } \ln B \approx \boldsymbol {b} \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {b}$, which enables the velocity-dependent diamagnetic and magnetic drift frequencies

(2.3a,b)\begin{equation} \tilde{ \omega }_* = \omega_* \left( 1 + \eta \left( w^{2}-\frac{3}{2} \right) \right)\quad \mathrm{and} \quad \tilde{ \omega }_d = \omega_d \left( w_\parallel^{2} + \frac{ w_\perp^{2} }{ 2} \right) \end{equation}

to be succinctly expressed entirely in terms of the four parameters

(2.4a–d)\begin{align} \omega_* = \frac{k_yL_{\mathrm{ref}}}{\sqrt{2}L_n} = \omega_{*0} k_y,\quad \omega_d = \frac{\sqrt{2}k_yL_{\mathrm{ref}}}{R} = \omega_{d0} k_y,\quad \eta = \frac{L_n}{L_T},\quad \mathrm{and}\quad \tau=\frac{T_i}{T_e}, \end{align}

from the electron/ion temperatures $T_{e/i}$ and the characteristic density, temperature and magnetic curvature lengths

(2.5a–c)\begin{equation} L_n=\left(\frac{\textrm{d}\ln n}{\textrm{d} x}\right)^{{-}1},\quad L_T=\left(\frac{\textrm{d}\ln T}{\textrm{d} x}\right)^{{-}1}\quad \mathrm{and} \quad R=\left(\frac{\textrm{d}\ln B}{\textrm{d} x}\right)^{{-}1}, \end{equation}

all of which are negative by our convention.

To couple the potential $\varphi$ to the ion gyrocentre distribution $h$ and close the system, the electrons are taken to follow a modified adiabatic response (Dorland & Hammett 1993; Hammett et al. 1993) such that the quasineutrality condition becomes

(2.6)\begin{equation} \int \textrm{d}^{3}\,\boldsymbol{v} \textrm{J}_0 h_{\boldsymbol{k}} = n (1 + \tau \hat{\alpha}) \varphi_{\boldsymbol{k}}, \end{equation}

where $\hat {\alpha }$ is the operator

(2.7)\begin{equation} \hat{\alpha}a=a(x,y)-\frac{1}{L_y}\int_0^{L_y}a(x,y)\,{\textrm{d} y}, \end{equation}

i.e. an operator that is zero when acting on purely zonal $\boldsymbol {E}\times \boldsymbol {B}$ modes with $k_y=0$, and unity otherwise.

To serve our purpose of studying the Dimits shift, the gyrokinetic equation in the form of (2.1) clearly neglects both parallel variations and collisions. The former omission constitutes a considerable simplification from a spatially three-dimensional to a spatially two-dimensional system, but necessarily excludes the ITG slab mode. Instead the focus becomes a local description of the well-known bad-curvature-driven toroidal ITG instability (Beer 1995), which seems to be of most relevance for the Dimits transition (Dimits et al. 2000). The second omission is similarly made because, should collisions be included, their presence significantly muddies the waters. This is because a wide range of zonal flow behaviour then manifests, including bursty patterns (see Berionni & Gürcan 2011) or non-quasistatic flows (see Kobayashi & Rogers 2012), so that it can become somewhat difficult to identify a clear Dimits transition or even reliably define the Dimits shift. However, in their absence, Landau-undamped Rosenbluth–Hinton states (Rosenbluth & Hinton 1998) can produce static zonal flow states with zero transport, in principle (only limited by the finite simulation time available to find such a state) providing a clear cut distinction between systems within and outside the Dimits regime.

2.2. The strongly driven gyrokinetic fluid limit

Employing a subsidiary ordering such that

(2.8)\begin{equation} \frac{\partial_t}{\eta\omega_*} \sim \frac{\omega_d}{\omega_*} \sim k_\bot^{2} \sim \frac{\tilde{\varphi}^{2}}{\bar{\varphi}^{2}} \sim \frac{1}{\eta} \ll 1, \end{equation}

where the gyrophase-independent response and potential have been split into their zonal and non-zonal components like

(2.9a–d)\begin{equation} \bar{h} = \left(1-\hat{\alpha}\right) h, \quad \tilde{h} = \hat{\alpha} h, \quad \bar{\varphi} = \left(1-\hat{\alpha}\right) \varphi, \quad \tilde{\varphi} = \hat{\alpha} \varphi, \end{equation}

one finds (see Appendix A) that the gyrokinetic moment hierarchy self-consistently closes at second order, resulting in the renormalised equation system

(2.10)\begin{gather} \frac{\partial \tau \tilde{\varphi}_{\boldsymbol{k}}}{\partial t} + \left\{ {\bar{\varphi}} , {\tau \tilde{\varphi}} \right\}_{\boldsymbol{k}} + \left\{ {\boldsymbol{\nabla}_\bot^{2} \bar{\varphi}} , {\tilde{T}_\perp} \right\}_{\boldsymbol{k}} - k_\bot^{2} \left\{ {\bar{\varphi}} , {\tilde{T}_\perp} \right\}_{\boldsymbol{k}} - \left\{ {\bar{\varphi}} , {\boldsymbol{\nabla}_{\bot}^{2}\tilde{T}_\bot} \right\}_{\boldsymbol{k}} \nonumber\\ \quad = \textrm{i} \omega_* \left(1 - \eta k_\bot^{2} \right) \tilde{\varphi}_{\boldsymbol{k}} - \textrm{i} \omega_d \left( \frac{\tilde{T}_{{\parallel}\boldsymbol{k}}}{2} + \frac{\tilde{T}_{{\perp}\boldsymbol{k}}}{4} \right) - D_{\boldsymbol{k}} \tau \tilde{\varphi}_{\boldsymbol{k}}, \end{gather}

(2.11)\begin{gather} \frac{\partial k_x^{2} \bar{\varphi}_{\boldsymbol{k}}}{\partial t} + \left\{ {\boldsymbol{\nabla}_\bot^{2} \tilde{\varphi}} , {\tilde{T}_\perp} \right\}_{\boldsymbol{k}} - k_x^{2} \left\{ {\tilde{\varphi}} , {\tilde{T}_\perp} \right\}_{\boldsymbol{k}} - \left\{ {\tilde{\varphi}} , {\boldsymbol{\nabla}_\bot^{2} \tilde{T}_\bot} \right\}_{\boldsymbol{k}} = 0, \end{gather}

(2.12)\begin{gather}\frac{\partial \tilde{T}_{{\perp}\boldsymbol{k}}}{\partial t} + \left\{ {\bar{\varphi}} , {\tilde{T}_\perp} \right\}_{\boldsymbol{k}} = \frac{\textrm{i} \omega_* \eta \tilde{\varphi}_{\boldsymbol{k}}}{2} - D_{\boldsymbol{k}} \tilde{T}_{{\perp}\boldsymbol{k}}, \end{gather}

(2.13)\begin{gather}\frac{\partial \tilde{T}_{{\parallel}\boldsymbol{k}}}{\partial t} + \left\{ {\bar{\varphi}} , {\tilde{T}_\parallel} \right\}_{\boldsymbol{k}} = \frac{\textrm{i} \omega_* \eta \tilde{\varphi}_{\boldsymbol{k}}}{4} - D_{\boldsymbol{k}} \tilde{T}_{{\parallel}\boldsymbol{k}}. \end{gather}

Here, an ad hoc damping operator $D_{\boldsymbol {k}}$ acting on the non-zonal components, to be further discussed in § 3.2, has been added. This is to compensate for the loss of collisionless damping (Landau 1946) that occurs upon taking moments of the gyrokinetic equation. Note that the zonal components of the temperature do not enter, the system only consists of one zonal field $\bar {\varphi }$ and three non-zonal fields $\tilde {\varphi }$, $\tilde {T}_\perp$, $\tilde {T}_\parallel$, which, as a consequence of (2.8), differ by order like

(2.14)\begin{equation} \frac{\tilde{\varphi}}{\tilde{T}_\perp} \sim \frac{\tilde{\varphi}}{\tilde{T}_\parallel} \sim \frac{1}{\eta}. \end{equation}

However, combining (2.12) and (2.13) it is clear that the volume average of $\delta T=\tilde {T}_{\perp }-2\tilde {T}_{\parallel }$ transiently decays to zero under the action of $D_{\boldsymbol {k}}$. Nevertheless, we include this component in our simulations for completeness.

Some comments about (2.8) and its resulting system are now in order. First, apart from the additional separation of zonal and non-zonal components in the ordering scheme, this corresponds to a strongly driven limit with a high temperature gradient feeding a strong ITG instability and causing long wavelength turbulence to be dominant, previously studied separately in its linear (Plunk et al. 2014) and nonlinear (Plunk, Tatsuno & Dorland 2012) limits. Note that, although we call the limit ‘strongly driven’ since the drive term is large compared with the particle drift, stable modes still do exist, so one might alternatively call this limit non-resonant (in the linear fluid sense). As to the specific additional zonal/non-zonal separation within the ordering scheme, it is necessary for a consistent closure which includes both linear and nonlinear interactions. Beyond this, it also encapsulates the fact that only the former are so-called modes of minimal inertia (Diamond et al. 2005), being easily excited due to the density shielding of the adiabatic electron response. Furthermore, being Landau undamped, they can persist for long times, and so they are observed to be comparatively strong.

Secondly, all the present nonlinear terms affecting the drift waves involve zonal flows. Farrell & Ioannou (2009) have already shown that, beyond the Dimits regime, simple systems (specifically Hasegawa–Wakatani) can exhibit all relevant physics despite lacking drift wave self-interactions. Thus, it may be unsurprising that we here will find that the same can be true inside the Dimits regime. Beyond this, we note that the full nonlinear interaction is asymmetric between the different fields. While the governing equations for the non-zonal fields all include the typical $\boldsymbol {E}\times \boldsymbol {B}$-advection nonlinear $\left \{ {\bar {\varphi }} , {\cdot } \right \}$-term, both the zonal and non-zonal potentials $\bar {\varphi }$ and $\tilde {\varphi }$ are affected by an additional set of nonlinear diamagnetic drift finite Larmor radius terms coupling them to $\tilde {T}_\bot$. It should also be noted that, by ordering, there is no Reynolds stress, i.e. a term of the form $\left \{ {\varphi } , {\nabla ^{2}\varphi } \right \}_{\boldsymbol {k}}$, present. It has been pointed out that such a term greatly facilitates the construction of strong zonal flows (Diamond & Kim 1991), but as a consequence of zonal flows being unaffected by $D_{\boldsymbol {k}}$, zonally dominated states will here arise even though the Reynolds stress is absent.

Thirdly, barring the splitting of the temperature moment into its separate parallel and perpendicular components, the nonlinear interaction is the same as Plunk et al. (2012). By a trivial modification of the results therein, the electrostatic energy conserved by the nonlinear interactions in (2.10)–(2.13) is therefore readily found to be given by

(2.15)\begin{equation} E_\varphi = E_{\bar{\varphi}} + E_{\tilde{\varphi}} = \frac{1}{2L_xL_y}\int |\boldsymbol{\nabla}\bar{\varphi}|^{2}\,\textrm{d}x\,\textrm{d}y + \frac{1}{2L_xL_y}\int \tau\tilde{\varphi}^{2} \,\textrm{d}x\,\textrm{d}y. \end{equation}

Finally, this strongly driven system seems formally far from the usual marginally unstable Dimits regime by virtue of its ordering, and one might question its relevance when investigating the Dimits shift. However, with a sufficiently large $D_{\boldsymbol {k}}$, marginal stability can be reinstated and a clear Dimits regime emerges. This system may thus act as a stepping stone, since its self-consistent closure means that the nonlinear interaction should closely resemble that of full gyrokinetics, at least in its range of validity. Indeed, it bears much resemblance to another, but highly collisional, gyrokinetic fluid limit recently studied by Ivanov et al. (2020). Beyond being collisionless, it differs in mainly three ways: (i) zonal flows are not subject to collisional dissipation, (ii) the nonlinear drift-wave self-interaction and Reynolds stress become too small to be of relevance and (iii) the zonal temperature perturbations cease to be dynamically relevant.

3. Key features

3.1. Primary instability

To arrive at a linear dispersion relation for the primary modes of the system (2.10)–(2.13), plane wave solutions proportional to $\exp (\lambda ^{p}_{\boldsymbol {k}}t + \textrm {i} \boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {r})$ are postulated, where $\lambda ^{p}_{\boldsymbol {k}}$ can be split into the growth rate $\gamma ^{p}_{\boldsymbol {k}}$ and frequency $\omega ^{p}_{\boldsymbol {k}}$ according to $\lambda ^{p}_{\boldsymbol {k}} = \gamma ^{p}_{\boldsymbol {k}} - \textrm {i} \omega ^{p}_{\boldsymbol {k}}$. A straightforward linear instability calculation then reveals the presence of a pure temperature mode (where $\tilde {\varphi }_{\boldsymbol {k}}=0$, but $\tilde {T}_\perp$ and $\tilde {T}_\parallel$ are non-zero) which is strictly damped,

(3.1)\begin{equation} \gamma^{pT}_{\boldsymbol{k}}={-}D_{\boldsymbol{k}}, \end{equation}

and two modes with the expected dispersion forms

(3.2)\begin{equation} \lambda_{\boldsymbol{k}}^{p\pm} = \gamma_{\boldsymbol{k}}^{p\pm} - \textrm{i} \omega_{\boldsymbol{k}}^{p\pm} ={-} D_{\boldsymbol{k}} + \textrm{i} \omega_* \frac{1 - \eta k_\bot^{2}}{2 \tau} \pm \frac{1}{2 \tau} \sqrt{ \eta \tau \omega_* \omega_d - \omega_*^{2} \left( 1 - \eta k_\bot^{2} \right)^{2} } \end{equation}

of the toroidal ITG mode inherent to the ordering (2.8) (Plunk et al. 2014). Note that, here and elsewhere, we reserve the $p$, $s$, $t$ superscripts for primary, secondary and tertiary quantities, and use $\pm$-superscripts to indicate the most/least unstable modes of each kind.

Because $D_{\boldsymbol {k}}$ generally introduces only a $k$-dependent shift in $\gamma$ towards lower values, it is useful to first consider $D_{\boldsymbol {k}} = 0$. Remembering that the definitions of $\omega _*$ and $\omega _d$ include a factor $k_y$, several features are readily apparent. The most unstable mode is as expected the purely radial streamer with $k_x=q=0$ satisfying

(3.3)\begin{equation} k_y^{2} = p^{2} = \frac{1}{3\eta}\left(2+\sqrt{1 + \frac{3 \eta \tau \omega_d}{\omega_*}}\right). \end{equation}

Note here the introduction of $q$ and $p$ which will henceforth be used for poloidal and radial wavenumbers, respectively. Now, when (3.3) is inserted into (3.2), it gives the expected bad-curvature ITG instability scaling (see Beer 1995)

(3.4)\begin{equation} \gamma_{\boldsymbol{k}}^{p+} \propto \frac{k_y}{\sqrt{\tau R L_T}}, \end{equation}

when the correction term under the root is taken to be small. When this term, on the other hand, is sufficiently large, the growth rate passes through zero. Thus, we find that only the wavenumbers within the annulus

(3.5)\begin{equation} \frac{1}{\eta}\left(1 - \sqrt{\frac{\eta \tau \omega_d}{\omega_*}}\right) < k_\bot^{2} < \frac{1}{\eta}\left(1 + \sqrt{\frac{\eta \tau \omega_d}{\omega_*}}\right) \end{equation}

can be unstable. Here, we see that $\eta$ pushes the instability to larger scales, while $\tau \omega _d/\omega _*$ controls the narrowness of the instability annulus and whether large scales are damped or not (indeed as $\eta \tau \omega _d/\omega _*$ exceeds 1, the annulus becomes a disk), see figure 1. Clearly (3.3) and (3.5) set the energy injection scale to be $\sqrt {(1/\eta )}$ in accordance with the subsidiary ordering (2.8), which is therefore justified a posteriori.

Figure 1. Normalised primary drift-wave growth rate $-\omega _{*0}^{-1}\gamma _{\boldsymbol {k}}^{p+}$ from (3.2) as a function of the radial and poloidal wavenumbers $q$ and $p$ for $\eta \tau \omega _d/\omega _*=0.2$, showing clearly the instability annulus (3.5). The instability boundary for primary unstable modes, in the presence of that $D_{\boldsymbol {k}}=D$ for which this configuration constitutes the Dimits threshold, is also shown.

We can arrive at a linear instability threshold for the temperature gradient when $D_{\boldsymbol {k}}=D$ is constant. Remembering that $\omega _*$ and $\omega _d$ both include a factor $k_y$, we find,upon inserting (3.3) into (3.2) and setting the result to $0$, that it becomes a condition on $\eta (D)$ for the most unstable mode to be marginally stable. Its $\eta >0$-solution is denoted by

(3.6)\begin{equation} \eta^{0} = \frac{8 \tau D^{2}}{9\omega_{d0}\omega_{*0}} \end{equation}

and it is found that, for $\eta$ larger than this value, $\gamma ^{p+}$ increases monotonically with $\eta$. When $D_{\boldsymbol {k}}$ is allowed to vary with respect to $p$, a correction to (3.6) appears, but monotonicity continues to hold. Therefore, for some $\eta ^{0}$ (typically close to (3.6) with $D=D_{\boldsymbol {p}}$) we have the necessary condition for instability

(3.7)\begin{equation} \eta > \eta^{0}. \end{equation}

Returning again to the unstable mode it naturally includes both potential and temperature perturbations, so upon inserting (3.2) into (2.10)–(2.13) the ratio between the two can be calculated to be

(3.8)\begin{equation} \frac{\tilde{\varphi}_{\boldsymbol{k}}}{\tilde{T}_{{\perp}\boldsymbol{k}}} = \frac{\tilde{\varphi}_{\boldsymbol{k}}}{2\tilde{T}_{{\parallel}\boldsymbol{k}}} = \frac{1}{\eta\tau}\left(1-\eta k_\bot^{2} - \sqrt{(1-\eta k_\bot^{2})^{2} - \frac{\eta\tau\omega_d}{\omega_*}}\right). \end{equation}

Using (3.5) to parametrise all the unstable modes with an angle $0\leqslant \theta \leqslant {\rm \pi}$ like

(3.9)\begin{equation} k_\bot^{2} = \frac{1}{\eta}\left(1 - \cos\theta\sqrt{\frac{\eta \tau \omega_d}{\omega_*}}\right), \end{equation}

one finds upon inserting this into the right-hand side of (3.8) that it reduces to the simple expression

(3.10)\begin{equation} \frac{\tilde{\varphi}_{\boldsymbol{k}}}{\tilde{T}_{{\perp} \boldsymbol{k}}} = \frac{\tilde{\varphi}_{\boldsymbol{k}}}{2\tilde{T}_{{\parallel} \boldsymbol{k}}} = \sqrt{\frac{\omega_d}{\eta\tau\omega_*}}\exp({\textrm{i}k_y\theta/|k_y|}), \end{equation}

since by our convention $\omega _{*0}<0$. Now, since the radial heat flux is given by

(3.11)\begin{equation} Q=\frac{1}{L_xL_y}\int \hat{x}\boldsymbol{\cdot} v_{\boldsymbol{E}} (\tilde{T}_\bot + \tilde{T}_\parallel) \,{\textrm{d} y} \,{\textrm{d} x} = \sum_{\boldsymbol{k}}\operatorname{Re}{\left(\textrm{i}k_y\tilde{\varphi}^{*}_{\boldsymbol{k}}(\tilde{T}_{{\perp} \boldsymbol{k}}+\tilde{T}_{{\parallel} \boldsymbol{k}})\right)}, \end{equation}

(3.10) implies that each mode provides a positive contribution to the total heat flux, since

(3.12)\begin{equation} \operatorname{Re}{\left(\textrm{i}k_y\tilde{\varphi}^{*}_{\boldsymbol{k}}(\tilde{T}_{{\perp} \boldsymbol{k}}+\tilde{T}_{{\parallel} \boldsymbol{k}})\right)}=3\sqrt{\frac{\omega_d}{\eta\tau\omega_*}}|\tilde{\varphi}_{\boldsymbol{k}}|^{2}k_y\sin{\left(\frac{k_y\theta}{|k_y|}\right)} \geqslant 0. \end{equation}

Additionally, (3.10) makes it clear that the potential component of the mode decreases compared with its temperature as $\eta$ increases. This linear result can be approached intuitively, since a strong temperature gradient causes $v_{\boldsymbol {E}}\boldsymbol {\cdot }\boldsymbol {\nabla } f_0$, the free energy source term in the gyrokinetic equation (where $v_{\boldsymbol {E}}$ is the $\boldsymbol {E}\times \boldsymbol {B}$-drift), to possess a larger temperature moment than density moment, the latter of which is most important for $\varphi$ through the quasineutrality condition (2.6).

3.2. The damping operator

Having determined the linear properties of (2.10)–(2.13) we are now in a position to discuss $D_{\boldsymbol {k}}$ in greater detail. Examining the linear growth rate (3.2) it is clear that, as long as there exists bad magnetic curvature providing finite $\omega _d$, and in the absence of artificial dissipation $D_{\boldsymbol {k}}$, the primary instability is present at $\eta k_\bot ^{2}=1$ given any arbitrarily small density and temperature gradients. Furthermore, all arbitrarily small scales are completely undamped and so can act as a reservoir of energy. Numerically, this means that even though every unstable $\lambda _{\boldsymbol {k}}^{p+}$-mode in the injection range is accompanied by a damped $\lambda _{\boldsymbol {k}}^{p-}$-mode, without $D_{\boldsymbol {k}}$ the system could nonlinearly diverge while exhibiting large-scale energy pileup typical of two-dimensional turbulence with its inverse energy cascade (Kraichnan 1967; Qian 1986; Terry 2004). In three-dimensional turbulence this is prevented by a scale balance of parallel streaming and turbulence, known as critical balance (Barnes, Parra & Schekochihin 2011), but no such mechanism is available here.

Given what was just outlined, in order to prevent non-physical absolute instability and the excitation of arbitrarily small or large scales, the necessity of including some kind of $D_{\boldsymbol {k}}$ is apparent. Physically, this is meant to represent the Landau-type damping present in weakly collisional toroidal ITG but which was lost upon only considering its moments to arrive at (2.10)–(2.13) (Sugama 1999). Although our ordering (2.8) implies that the kinetic damping is small, it is nevertheless non-zero and so dynamically relevant, particularly for the marginally stable small scales it firmly stabilises. Its inclusion is further justified since the ordering (2.8), although ‘strongly driven’, nevertheless allows the primary instability to also be weak, so that $D_{\boldsymbol {k}}$ can stabilise the system so it exhibits a Dimits regime.

As to the specific form of $D_{\boldsymbol {k}}$ which we will employ in this paper, we will always include a constant component $D$, present for all non-zonal modes. The reason for this is that, at least within and close to the Dimits regime, we have found its inclusion to be sufficient to prevent a large-scale energy pileup. This form has some physical justification, in that a rigorous linear analysis of the full kinetic mode reveals, beyond the normal mode whose Landau damping can be approximated as viscous dissipation, the presence of an algebraically decaying continuum mode. In the marginally stable regime of interest, the continuum modes of the sidebands should thus be dominant, since they decay much more slowly (Sugama 1999; Mishchenko, Plunk & Helander 2018). Unfortunately, it is very hard to accurately reproduce the behaviour of these modes in a non-exotic way for a spectral fluid model (Sugama 1999), and so, in the absence of better alternatives, a flat decay can be used to model this. Beyond this component it is natural to include some kind of hyperviscosity $k_\bot ^{\alpha }$ in $D_{\boldsymbol {k}}$. However, this seems to have little effect on the key results of this paper, presumably because of how sharply peaked in $k$-space the linear instability is, and so we typically do not include it. For generality we will nevertheless allow it in our instability calculations.

3.3. Secondary instability

Within the primary–secondary–tertiary hierarchy, the secondary instability develops once the primary drift waves have grown sufficiently for slight flow inhomogeneities to amplify through a shearing interaction to magnify small zonal perturbations (Kim & Diamond 2002). It is analytically best treated via a Galerkin truncation of (2.10)–(2.13) into the 4-mode system consisting of the most unstable purely radial primary mode, its two sidebands and a single zonal mode

(3.13a–c)\begin{align} \boldsymbol{p}=(0,p),\quad \boldsymbol{r}^{{\pm}}=({\pm} q,p),\quad\boldsymbol{q}=(q,0). \end{align}

The potential and temperature amplitudes of the primary mode are then fixed and taken to be much larger than other variables so that the linear terms of the sidebands can be ignored compared with the nonlinear interaction with the primary mode, which now become linearised. Then it is straightforward to obtain the KH-like secondary dispersion relation

(3.14)\begin{equation} {\lambda^{s}}^{2} + 4 p^{2} q^{2} \left| \tilde{T}_{\bot \boldsymbol{p}} \right|^{2} \left( \frac{2 q^{2} }{\tau} - \operatorname{Re}\left(\frac{\tilde{\varphi}_{\boldsymbol{p}}}{\tilde{T}_{\bot \boldsymbol{p}}}\right)\right) = 0. \end{equation}

Inserting the potential/temperature ratio of the unstable mode (3.8) into (3.14) is now natural since it is this mode which initiates the secondary instability in the primary–secondary paradigm, and this results in

(3.15)\begin{equation} {\lambda^{s}}^{2} + \frac{4 p^{2} q^{2} \left| \tilde{T}_{\bot \boldsymbol{p}} \right|^{2}}{\eta\tau} \left( 2\eta q^{2} + \eta p^{2} - 1 + \operatorname{Re}\sqrt{(1-\eta p^{2})^{2}-\frac{\eta\tau\omega_d}{\omega_*}} \right) = 0. \end{equation}

The last, stabilising term in (3.15) is similar to the opposite of the destabilising term of $\gamma _{\boldsymbol {k}}^{p+}$ in (3.2), but gives rise to a discontinuity instead of a bifurcation at

(3.16)\begin{equation} \eta p^{2} = 1\pm\sqrt{\frac{\eta\tau\omega_d}{\omega_*}}, \end{equation}

a feature which can clearly be seen in figure 2 where the secondary growth rate $\gamma ^{s+}$ is plotted. As $p$ approaches the lower threshold from below, the radical in (3.15) approaches $0$, causing $\gamma ^{s+}$ to rapidly increase until its derivative with respect to $p$ discontinuously flattens.

Figure 2. Normalised secondary growth rate $\sqrt {\tau \eta }|\tilde {T}_{\perp \boldsymbol {p}}|^{-1}\gamma ^{s+}$ as a function of $q$ and $p$ for $\eta \tau \omega _d/\omega _*=0.7$. The bifurcation of (3.16) is clearly visible.

Now the absolute requirement on the zonal wavenumber in order for an unstable primary mode to be secondary unstable is established by (3.15) to be

(3.17)\begin{equation} q^{2} < \frac{1-\eta p^{2}}{2\eta}. \end{equation}

This means that modes with $p^{2}>\eta$, and in particular the most unstable mode satisfying (3.3), are completely stable to the secondary instability and so can continue to grow unabated until another channel for zonal flow generation is established.

One might suspect that the zonal flow would be initiated by those less unstable primary modes with smaller $p$ which are able to satisfy (3.17) once they have grown to a sufficient amplitude. However, in simulations it is instead observed that, since the primary growth rate is rather sharply peaked around (3.3), these modes do not grow fast enough to be dynamically relevant at this stage. Instead, it seems that, since the small $q$-sidebands of the most unstable primary mode grow at nearly the same rate, it is their mutual sideband–sideband interaction which jump starts the zonal growth. This is evidenced by the fact that the initial zonal growth rate remains mostly unchanged even when all modes but the most unstable primary mode and its sidebands are set to 0. Thus we conclude that the secondary instability of this form is, in fact, presently irrelevant in the Dimits regime.

3.4. Local tertiary instability

Turning now to the final stage of the primary–secondary–tertiary hierarchy, once the zonal flow has grown enough to quench the drift waves the tertiary instability is that instability which allows the drift waves to re-emerge from a zonally dominated state (Rogers et al. 2000). In analysing this instability we will consider two separate limits, one localised and one de-localised (i.e. localised in $k$-space).

It is a well-known feature of tertiary modes that they localise to regions of zero zonal shear rate, $\partial _x^{2}\bar {\varphi } = 0$ (Kobayashi & Rogers 2012; Kim, Min & An 2018, 2019). Therefore, we consider a poloidal band of modes ($k_y=p$) subject to a large amplitude zonal flow localised around such a point. Taking $D_{\boldsymbol {k}}=D$ again, in real $x$-space equations (2.10)–(2.13) then become

(3.18)\begin{gather} \left(\partial_t+\textrm{i}p\partial_x\bar{\varphi}\right)\tilde{T}_\bot - \frac{\textrm{i}\eta\omega_*\tilde{\varphi}}{2} ={-}D\tilde{T}_\perp, \end{gather}

(3.19)\begin{gather}\left(\partial_t+\textrm{i}p\partial_x\bar{\varphi}\right)\tilde{T}_\parallel{-} \frac{\textrm{i}\eta\omega_*\tilde{\varphi}}{4} ={-}D\tilde{T}_\parallel, \end{gather}

(3.20)\begin{gather} \tau\left(\partial_t+\textrm{i}p\partial_x\bar{\varphi}\right)\tilde{\varphi} - \textrm{i}\omega_*\left(1-\eta(p^{2}-\partial_x^{2})\right)\tilde{\varphi} +\textrm{i}\left(\frac{\omega_d}{4}+2p\partial_x^{3}\bar{\varphi}+2p\partial_x^{2}\bar{\varphi}\partial_x\right)\tilde{T}_\bot\nonumber\\ \quad + \textrm{i}\frac{\omega_d}{2}\tilde{T}_\parallel={-}D\tilde{\varphi}. \end{gather}

If we consider a narrow region in which $\partial _x\bar {\varphi }$ and $\partial _x^{3}\bar {\varphi }$ are approximately constant and allow ourselves to consider the mode to also be localised around the most unstable primary mode with $k_x\approx 0$, we therefore find the local tertiary dispersion form

(3.21)\begin{equation} \lambda^{t\pm}_L ={-}D+\textrm{i}\left(\frac{\omega_*-\eta\omega_*p^{2}}{2\tau}-p\partial_x\bar{\varphi}\right) \mp \frac{\omega_*}{2\tau}\sqrt{\frac{\eta\tau}{\omega_*}(\omega_d+8p\partial_x^{3}\bar{\varphi})-\left(1-\eta p^{2}\right)^{2}}. \end{equation}

As is easily seen, this expression is precisely the linear dispersion of the primary mode (3.2) Doppler shifted by $p\partial _x\bar {\varphi }$ and with a zonal shear modified magnetic curvature. We note that the real part of this expression vanishes when the driving gradients are removed, meaning that no tertiary instability exists at all in their absence. We are thus dealing here with only a modified primary, extracting energy from the background gradients, rather than from the zonal flow like the KH instability (see Zhu et al. 2020a). This is because the fundamental ordering (2.8) eliminates both the Reynolds stress and the zonal temperature. If present, the former would give rise to a true tertiary KH instability (Kim & Diamond 2002; Zhu, Zhou & Dodin 2018), and the latter a tertiary KH-like instability, analogous to the secondary instability (3.15) (Rogers et al. 2000; Ivanov et al. 2020).

Returning to the specific expression (3.21), we see that the tertiary instability is asymmetric with respect to zonal flow velocity minima, $\partial ^{3}_x\bar {\varphi }>0$, and maxima, $\partial ^{3}_x\bar {\varphi }<0$; the former are destabilising while the latter are stabilising. This asymmetry matches gyrokinetic observations that zonal flow minima are significantly more prone to turbulent transport (McMillan et al. 2011), and has already been noted in previous tertiary instability studies of simple systems (see e.g. Zhu et al. 2020a). As we will see in § 4, the same holds true here: turbulence consistently localises around the points where $\partial _x^{2}\bar {\varphi }=0$ and $\partial _x^{3}\bar {\varphi }>0$ in accordance with (3.21). Nevertheless, results like (3.21) cannot be taken at face value. In the closely related system of Ivanov et al. (2020), where the presence of zonal temperature perturbations considerably complicates the picture, the equivalent expression for the growth rate fails to match what is observed in simulations.

3.5. Four-mode tertiary instability

We now proceed to study the tertiary instability of a sinusoidal profile corresponding to the mode $\bar {\varphi }_{\boldsymbol {q}}$. In order to gain further insight we employ, for simplification, the same 4-mode (4M) Galerkin truncation (3.13a–c) as we did for the secondary instability, even though in general it is less justifiable here. Naturally with so few modes this analysis will fail to capture any intricate localisation effects, but we will nevertheless be able to discern some important features of the tertiary instability. After all, (3.21) employed strong approximations (e.g. fully neglecting $k$-space coupling), which may in general not be satisfied.

Assuming without loss of generality that $\bar {\varphi }_{\boldsymbol {q}}$ is real (all results herein will only depend upon its magnitude), i.e.

(3.22)\begin{equation} \bar{\varphi}=2\bar{\varphi}_{\boldsymbol{q}}\cos{qx}, \end{equation}

the 4M tertiary dispersion relation can after some algebra, analogous to that of § 3.3, be expressed as the product of three polynomials like

(3.23)\begin{equation} \mathcal{D}_{{P}\boldsymbol{r}}\left(\lambda^{t}_{4M}\right)\mathcal{D}_{T}\left(\lambda^{t}_{4M}\right)\mathcal{D}_{\mathrm{MP}}\left(\lambda^{t}_{4M}\right)=0. \end{equation}

This factorised form of the dispersion relation separates the linear tertiary modes into different groups, corresponding to zeros of each of the three factors. The equation $\mathcal {D}_{{P}_{\boldsymbol {r}}}=0$ is the unmodified primary dispersion relation of the sidebands $\boldsymbol {r}^{\pm }$ with solutions (3.1) and (3.2), corresponding to a solution of (2.10)–(2.13) where the primary mode is absent (i.e. the $\boldsymbol {p}$-mode is 0) and the sidebands are of equal amplitude. Next

(3.24)\begin{equation} \mathcal{D}_{T}=(\lambda^{t}_{4M}+D_{\boldsymbol{p}})(\lambda^{t}_{4M}+D_{\boldsymbol{r}})+2p^{2}q^{2}\bar{\varphi}_{\boldsymbol{q}}^{2}=0 \end{equation}

is the dispersion relation of two stable pure temperature modes affected by the zonal flow, and finally

(3.25)\begin{align} \mathcal{D}_{\mathrm{MP}} = (\lambda^{t}_{4M}-\lambda^{p+}_{\boldsymbol{p}})(\lambda^{t}_{4M}-\lambda^{p-}_{\boldsymbol{p}}) (\lambda^{t}_{4M}-\lambda^{p+}_{\boldsymbol{r}})(\lambda^{t}_{4M}-\lambda^{p-}_{\boldsymbol{r}})+4p^{4}q^{4} \bar{\varphi}_{\boldsymbol{q}}^{4}+Ap^{2}q^{2}\bar{\varphi}_{\boldsymbol{q}}^{2}=0, \end{align}

where

(3.26)\begin{align} A(\lambda^{t}_{4M}) &= 4(\lambda^{t}_{4M}+D_{\boldsymbol{p}})(\lambda^{t}_{4M}+D_{\boldsymbol{r}}) +\frac{\eta\omega_*\omega_d}{\tau}-2\frac{\omega_*^{2}}{\tau^{2}}(1-\eta p^{2})^{2}\nonumber\\ &\quad +2\textrm{i}\frac{\omega_*}{\tau}\left[(\eta p^{2}-1)(2\lambda^{t}_{4M}+D_{\boldsymbol{p}}+D_{\boldsymbol{r}})-\eta q^{2}(\lambda^{t}_{4M}+D_{\boldsymbol{r}})\right], \end{align}

is the dispersion relation of the 4M zonal flow modified primary mode. Henceforth, we focus on this latter equation, since the modified primary will prove to be the most unstable tertiary mode.

Let us consider the dispersion relation of the modified primary in the large zonal flow limit. Expanding (3.25) in orders of $\bar {\varphi }_{\boldsymbol {q}}$ like

(3.27)\begin{equation} \lambda^{t}_{4M}=\lambda_1^{t}+\lambda_{1/2}^{t}+\textit{O}(1), \end{equation}

we find, after collecting terms up to order $\sqrt {\bar {\varphi }_{\boldsymbol {q}}}$ and using (3.2), that (3.25) can be reduced to

(3.28)\begin{equation} \left(\left(\lambda^{t}_{4M}\right)^{2}+2p^{2}q^{2}\bar{\varphi}_{\boldsymbol{q}}^{2}\right)^{2}+B \left(\lambda^{t}_{4M}\right)^{3}+\left(2B-4\textrm{i}\frac{\omega_*}{\tau}\eta q^{2}\right)p^{2}q^{2}\bar{\varphi}_{\boldsymbol{q}}^{2}{\lambda^{t}}_{4M}=0, \end{equation}

where

(3.29)\begin{equation} B={-}(\lambda^{p+}_{\boldsymbol{p}}+\lambda^{p-}_{\boldsymbol{p}}+\lambda^{p+}_{\boldsymbol{r}}+\lambda^{p-}_{\boldsymbol{r}})=2(D_{\boldsymbol{p}}+D_{\boldsymbol{r}})+\frac{\textrm{i}\omega_*}{\tau}\left(\eta(q^{2}+2p^{2})-2\right). \end{equation}

At leading $\textit {O}(\bar {\varphi }_{\boldsymbol {q}}^{2})$-order (3.28) yields the purely oscillating solutions

(3.30)\begin{equation} \lambda_1^{t} ={\pm} \sqrt{2}\textrm{i}pq\bar{\varphi}_{\boldsymbol{q}}, \end{equation}

similar to the modes of (3.24). In order to find the real part of these modes we then have to proceed to order $\textit {O}(\bar {\varphi }_{\boldsymbol {q}})$ since the $\textit {O}(\bar {\varphi }_{\boldsymbol {q}}^{3/2})$-part identically vanishes. At that order (3.28) yields

(3.31)\begin{equation} \lambda_{1/2}^{t}={\pm}\sqrt{\pm\frac{\eta\omega_* q^{2}}{\sqrt{2}\tau}}\sqrt{pq\bar{\varphi}_{\boldsymbol{q}}}. \end{equation}

Combining these results, in the large $\bar {\varphi }_{\boldsymbol {q}}$-limit we therefore have the four solutions

(3.32)\begin{gather} \lambda_{4M}^{t\pm}=\sqrt{2}\textrm{i}pq\bar{\varphi}_{\boldsymbol{q}}\pm \sqrt{-\frac{\eta\omega_* q^{2}}{\sqrt{2}\tau}}\sqrt{pq\bar{\varphi}_{\boldsymbol{q}}}, \end{gather}

(3.33)\begin{gather}\lambda^{t}_{4M}={-}\sqrt{2}ipq\bar{\varphi}_{\boldsymbol{q}}\pm \sqrt{\frac{\eta\omega_* q^{2}}{\sqrt{2}\tau}}\sqrt{pq\bar{\varphi}_{\boldsymbol{q}}}, \end{gather}

of which only the first is unstable since $\omega _*<0$. Do keep in mind that we will continue to use the $+$-superscript for the most unstable mode, regardless of whether $\bar {\varphi }_{\boldsymbol {q}}$ is large or not.

Converting the solutions corresponding to (3.32) from $k$-space to real space using the zonal profile (3.22) we find that they take the form

(3.34)\begin{equation} \tilde{\varphi}=\left(1\pm\sqrt{2}\sin{qx}\right)\operatorname{Re}{\left(\exp({\textrm{i}py+\lambda^{{\pm}}_{4M}t})\right)}. \end{equation}

It is apparent that the $x$-envelope, being given by the first factor, predominantly localises the unstable mode around minima of the zonal flow velocity and the stable mode around maxima, entirely in accordance with the picture that these points are tertiary (de-)stabilising as outlined in § 3.4. Furthermore, it is seen that, despite not being sufficiently localised for the treatment of § 3.4 to be justified, this result nevertheless agrees with the large $\bar {\varphi }$-limit of (3.21) up to numerical constants.

Now, turning to the opposite small $\bar {\varphi }_{\boldsymbol {q}}$-limit, we are interested in how the unstable primary mode is modified by the presence of a small zonal flow. Taylor expanding (3.25) around $\lambda ^{t}_{4M}=\lambda ^{p+}_{\boldsymbol {p}}$ one straightforwardly obtains the solution

(3.35)\begin{equation} \lambda^{t+}_{4M}=\lambda^{p+}_{\boldsymbol{p}}-Cp^{2}q^{2}\bar{\varphi}_{\boldsymbol{q}}^{2}, \end{equation}

where

(3.36)\begin{equation} C=\frac{A(\lambda^{p+}_{\boldsymbol{p}})(\lambda^{p+}_{\boldsymbol{p}}-\lambda_{\boldsymbol{r}}^{p+})^{{\dagger}}(\lambda^{p+}_{\boldsymbol{p}}-\lambda_{\boldsymbol{r}}^{p-})^{{\dagger}}}{(\lambda^{p+}_{\boldsymbol{p}}-\lambda^{p-}_{\boldsymbol{p}})|(\lambda^{p+}_{\boldsymbol{p}}-\lambda^{p+}_{\boldsymbol{r}})(\lambda^{p+}_{\boldsymbol{p}}-\lambda^{p-}_{\boldsymbol{r}})|^{2}}. \end{equation}

Let us now employ a subsidiary ordering in $q$ to see how $\operatorname {Re}(C)$ behaves in the two limits $q^{2}\ll 1$ (keeping $q^{2}\gg \bar {\varphi }_{\boldsymbol {q}}^{2}$) and $q^{2}\gg 1$ (keeping $q^{2}\bar {\varphi }_{\boldsymbol {q}}^{2}\ll 1$) to see whether the most unstable mode is initially stabilised or destabilised by the presence of zonal flows at these scales. Because it is apparent that the denominator of (3.36) is positive, and because the numerator of $C$ becomes

(3.37)\begin{equation} -2\frac{\omega_*^{2}\eta^{2}q^{4}}{\tau^{2}}(|\lambda_{\boldsymbol{p}}^{p+}|^{2}+D_{\boldsymbol{r}}\gamma^{p+}_{\boldsymbol{p}}), \end{equation}

for large $q$-values, it is apparent that $\operatorname {Re}({C})$ is negative for small-scale zonal flows, which therefore are destabilising already at small amplitude. As $\bar {\varphi }_{\boldsymbol {q}}$ is further increased we furthermore know that the mode under consideration transitions into (3.32), and so we can conclude that small-scale zonal flows are always destabilising. In fact it is numerically found that the transition to instability with increasing $q$ occurs much before this limit, already as $q\sim p$, as can be seen in figure 3. Note that this point is still much below that where the KH-like $(qp\bar {\varphi }_{\boldsymbol {q}})^{\alpha }$-scaling of (3.32) develops, and thus the mode is still ostensibly ‘more primary’ in character.

Figure 3. Four-mode tertiary growth rate $\gamma ^{t+}_{4M}$ of the most unstable poloidal band satisfying (3.3) for $(\omega _{d0},\omega _{*0},\tau ,D,\eta )=(-0.8,-1.02,1,0.5,1.8)$ normalised by the unmodified primary growth rate $\gamma ^{p+}_{\boldsymbol {p}}$, as a function of the zonal wavenumber $q$ for $\bar {\varphi }_{\boldsymbol {q}}=0.1$, 0.3 and 0.6, respectively

If $q$ on the other hand is small, the numerator of $C$ becomes

(3.38)\begin{equation} 16\left(\gamma_{\boldsymbol{p}}^{p+}\right)^{2}\left(\textrm{i}\frac{\omega_*}{2\tau}-\frac{\omega_*^{2}(1-\eta p^{2})}{4\tau^{2}\gamma^{p+}_{\boldsymbol{p}}}\right)\eta q^{2}. \end{equation}

This has a positive real part for those modes with $\eta p^{2}>1$, including the most unstable primary mode satisfying (3.3), and so large-scale zonal flows are initially stabilising for these modes. It should be noted that, despite (3.38) reversing sign for modes of smaller $p$, the zonal flow does initially not destabilise large-scale drift waves. As can be seen from the linear growth rate (3.2), for these values of $p$ the small $q$ sidebands are in fact more unstable than the pure drift mode. Thus, upon repeating the calculation above, but instead expanding around $\lambda ^{t}_{4M}=\lambda _{\boldsymbol {r}}^{p+}$, one finds precisely the opposite stabilisation effect on the sidebands. These results can be summarised by the observation that, at small zonal amplitude, the tertiary mode constitutes a sort of weighted average of its constituent primary modes.

With the asymptotic behaviour of (3.32) and (3.35) in hand, it is clear that the initial stabilisation (3.35) of small-amplitude zonal flows must reverse as the amplitude is increased, and there necessarily exists some zonal flow amplitude which is most stable. Precisely this can be seen in figure 3, where the most unstable 4M tertiary growth rate $\gamma _{4M}^{t+}$ for the example system $(\omega _{d0},\omega _{s0},\tau ,D,\eta )=(-0.8,-1.02,1,0.5,1.8)$ with linear instability threshold $\eta ^{0}=1.64$ (which will be the focus of the remainder of this paper), is plotted. In accordance with (3.37) and (3.38), as $\bar {\varphi }_{\boldsymbol {q}}$ begins to be increased the tertiary mode initially stabilises for $q\lesssim p$, and destabilises for $q$ of greater magnitude.

3.6. Role of the tertiary instability for the Dimits transition

Extrapolating the consequences of the findings above to the dynamics of the zonally dominated states typical of the Dimits regime, some conclusions can be drawn. If zonal profile conditions are not ideal, the said profile will fail to suppress the tertiary instability. Then drift waves will grow in amplitude to eventually affect the zonal profile. While such conditions prevail, the zonal profile will evolve through different configurations in a process we will refer to as zonal profile cycling. In this process, energy will continue to be injected into the drift waves at a faster rate for more tertiary unstable zonal profiles. Thus the profile should be observed with higher probability in a state of low tertiary growth. Indeed, it is expected that a state of absolute tertiary stability could be sustained indefinitely. In conclusion, we argue that the tertiary instability therefore preferentially selects a set of zonal profiles which will predominantly appear as the system evolves.

Because the zonal flows by construction are linearly undamped, a tertiary stable zonal profile can emerge that sustains the system in a state of suppressed turbulence, so long as the decaying residual drift-wave activity, in turn, does not affect it too much. We will refer to such profiles as robustly stable. Now, from the 4M result above, we can extrapolate that the tertiary instability for our system exhibits only a finite ability to be stabilised. Naturally, this means that the number of robustly stable zonal profiles should decrease as the driving gradient $\eta$ is increased. At some point, none remain, so turbulence and transport must arise. If this point indeed corresponds to the Dimits transition, then the only the relevant features needed to explain the Dimits shift should be the tertiary instability and the ability of the zonal flows to cycle through stabilising profiles.

Of course, it is possible that, even in the absence of collisional zonal damping, a stable zonal state cannot be attained. Another possibility is that some nonlinear mechanism continues to reduce transport above the tertiary instability threshold, making the Dimits threshold of appreciable transport not coincide with the tertiary threshold. An example of such a feature, already observed and explored in other systems, would be e.g. the ferdinons of Ivanov et al. (2020). For our system, however, this does not seem to be the case, and, as we will see, the tertiary instability alone seems sufficient to explain the Dimits shift. That is, below a certain point $\eta =\eta ^{\mathrm {NL}}$, tertiary stable zonal profiles are always able to form and completely quench transport, while above it they cease to manifest and the time-averaged transport levels rapidly increase with $\eta$.

In conclusion, it should be noted that precisely what profiles are robust is a delicate and highly non-trivial question, which nevertheless ultimately decides when the Dimits transition occurs. Thus the tertiary instability should not enter the Dimits shift picture via so simple a rule as ‘the Dimits regime should end when the zonal amplitude becomes too large’ as envisioned by Rogers et al. (2000), nor ‘the Dimits regime should end when the zonal amplitude becomes too small’ as Zhu et al. (2020a) stated. Although the latter may hold when collisional damping limits the zonal amplitude, in general it is the much more nebulous question of ‘can a robust zonal profile be reached and sustained during the subsequent transient period of decay’ which must be answered.

4. Nonlinear simulation results

In order to thoroughly investigate the strongly driven system, it was first simulated pseudospectrally for several configurations on a square grid using a sixth-order Runge–Kutta–Fehlberg method including $512\times 256$ modes with $0\leqslant k_x,k_y\leqslant 5p_m$ and with dealiasing using the 3/2-rule. Sensitivity scans with regards to the number of modes and the minimum wavenumbers found this selection to be well beyond what was necessary for convergence within the Dimits regime, so long as the most unstable drift wave $p_m$ was included.

The nonlinear simulations usually employed a $\boldsymbol {k}$-independent $D_{\boldsymbol {k}}$, and were initiated with small Gaussian noise of (normalised) energy density $10^{-8}$. As can be seen in figure 4, the expected behaviour is then observed, where primary modes emerge until they are strong enough to nonlinearly engage the zonal modes. In accordance with the secondary mode analysis of 3.3, sideband–sideband interactions here play a vital role for initial zonal mode growth to occur, and so necessarily the second zonal harmonic, i.e. the mode with twice the smallest wavenumber $q_{\mathrm {min}}$, is primarily engaged. Following this, the largest-scale zonal modes also begin to grow appreciably, until they in turn reach sufficient amplitude for nonlinear interactions to quickly shuffle energy from the unstable modes to higher and higher $q$-sidebands. These, in turn, engage higher $q$ zonal modes to affect the primary growth, and the growth phase ceases. This typically occurs when both the drift-wave and zonal flow energy densities reach a comparable magnitude of around 1.

Figure 4. (a) Drift-wave and (b) zonal mode amplitudes during the initial growth phase for the system of figure 3. The blue line corresponds to the mode with smallest radial wavenumber $q_{\min }$, the red line to its second harmonic $2q_{\min }$ and other modes are denoted by black dashed lines. After an initial linear phase, the sideband–sideband interaction of $\tilde {\varphi }_{(q_{\min },p)}$ excites $\bar {\varphi }_{(q_{\min },0)}$ and $\bar {\varphi }_{(2q_{\min },0)}$, which thus grow at a rate proportional to $\sim |\tilde {\varphi }_{(q_{\min },p)}|^{2}$, which is plotted with a dotted blue line for comparison. Modes of higher and higher $q$ are then excited one by one, until the zonal flows reach a magnitude comparable to the drift waves, which are then suppressed.

It is important to note that, as a consequence of there being no direct coupling between drift waves of differing poloidal wavenumber, the system typically stratifies into separate $p$-layers only interacting with each other via their influence on the zonal flow. In the Dimits regime, where necessarily $D_{\boldsymbol {p}}\sim \gamma _{\boldsymbol {p}}^{p+}$, only those few modes with $p$ around $p_m$, satisfying (3.3), are linearly unstable. Consequently, the layer corresponding to the dominant primary mode becomes solely dynamically important, as is borne out in simulations. Although the zonal profile could excite other bands through the tertiary instability, since the primary band is the most tertiary unstable this is not borne out in practice. More layers become important only once they become primary unstable at larger $\eta$-values, occurring above the point at which the transition to continuous transport occurs.

Now, initial saturation amplitudes of both zonal and drift waves exhibit a very slight dependence on the initialisation amplitude. This is because differing primary/sideband growth rates causes there to be more or less initial energy to distribute, depending on the amount of time the primary mode has had to grow before the sidebands trigger zonal growth. Nevertheless, in the Dimits regime the zonal amplitude usually quickly returns to a system-configuration-dependent typical amplitude, as the tertiary analysis of (3.5) suggests. Once there, a rectangular zonal shear profile resembling a less developed version of the staircase states observed by e.g. Dif-Pradalier et al. (2010), Kobayashi & Rogers (2012) or Peeters et al. (2016), quickly develops, which can then efficiently suppress the amplitude of drift waves, typically on the order of 2 magnitudes (but occasionally much more). All that then remains are the localised tertiary modes which will eventually die if stable or grow back if unstable. Observing the heat flux after a long time has passed, as in figure 5, the presence of a Dimits shift is thus revealed, since stable states only exist close to the linear threshold $\eta ^{0}$.

Figure 5. Long-time-averaged heat flux, given by (3.11), as a function of $\eta$ for $(\omega _{d0},\omega _{s0},\tau ,D)=(-0.8,-1.02,1,0.5)$. The linear instability threshold occurs at $\eta ^{0}\approx 1.66$, yet finite heat flux only commences beyond $\eta ^{\mathrm {NL}}\approx 1.9$, constituting a clear Dimits shift. Between these points, the system relaxes to completely stable purely zonal states.

4.1. Drift-wave bursts

As can be seen figure 6, when the instability parameter $\eta$ is increased away from the linear threshold $\eta ^{0}$ while other parameters are kept the same, the initial zonal profiles attained are commonly completely tertiary stable. However, as $\eta$ is further increased this usually ceases to be the case, since the primary/sideband coupling then fails to remain strong enough for the primary mode to decay together with the damped sidebands. Consequently, a spreading turbulence burst destroys the initial zonal state, cycling through zonal profiles until another stable state is reached. The rapidity by which such bursts occur, and the time until a stable zonal profile is attained, both rapidly increase with $\eta$, unless the cycling by happenstance quickly produces a stable state. At even larger $\eta$-values, no stable state is ever attained. Furthermore, the typical burst amplitude is also reduced as a result of less efficient drift-wave quenching.

Figure 6. Time evolution for the configuration of figure 5 of the zonal (black dashed) and non-zonal (red) energy densities $E_{\bar {\varphi }}$ and $E_{\tilde {\varphi }}$, given by (2.15), with increasing instability $\eta$ above the linear threshold $\eta ^{0}$, where $\eta ^{0}$ is the linear instability threshold and $\Delta \eta ^{\mathrm {PR}}$ is the predicted Dimits transition as introduced in § 5. As the system becomes more unstable it is observed to take longer to arrive at a completely stable zonal state, while simultaneously exhibiting more rapid bursty behaviour.

Note that this entire burst pattern is, on the surface, similar to the zonal/drift predator–prey interactions commonly observed in many systems as a result of zonal damping (see e.g. Malkov, Diamond & Rosenbluth 2001; Kobayashi & Rogers 2012). However, it differs fundamentally in that the turbulent bursts are not typically accompanied by large zonal amplitude swings. Instead, it traces its origin to tertiary mode localisation.

To see how this is the case, some snapshots of a typical burst are displayed in figure 7. A zonal profile exhibits tertiary modes predominantly localised at the points where the conditions $\partial _x^{2}\bar {\varphi }=0$, and $\partial _x^{3}\bar {\varphi }>0$ are satisfied, in accordance with (3.21). Eventually, the one at $x\approx 32$ grows enough to affect the zonal profile at this point, resulting in a central flattening of the zonal amplitude. While the zonal shearing rate $\partial _x^{2}\bar {\varphi }$ remains 0 in the process, $\partial _x^{3}\bar {\varphi }$ is reduced, except at the boundary of the full mode, causing a central tertiary instability reduction. The tertiary mode now becomes more unstable at its boundary, where the condition $\partial _x^{3}\bar {\varphi }>0$ is still maintained. As a result, the zonal flattening continues and the tertiary mode broadens behind a propagating zonal front, eventually encompassing much of the domain and destroying the zonal profile.

Figure 7. Snapshots for the system of figure 3 at (a) $\gamma ^{p+}_{\boldsymbol {p}}t=33$, (b) $\gamma ^{p+}_{\boldsymbol {p}}t=35.8$, (c) $\gamma ^{p+}_{\boldsymbol {p}}t=37.5$ and (d) $\gamma ^{p+}_{\boldsymbol {p}}t=47$ depicting drift waves on the left and the zonal potential (blue), the zonal flow shear $\partial _x^{2}\bar {\varphi }$ (red) and its derivative $\partial _x^{3}\bar {\varphi }$ (dotted red) on the right. These depict a turbulent burst originating as an unstable tertiary mode at $x\approx 32$ where $\partial _x^{2}\bar {\varphi }=0$ and $\partial _x^{3}\bar {\varphi }>0$ at $\gamma _{\boldsymbol {p}}^{p+} t=33$, broadening and growing in amplitude between two tertiary unstable propagating zonal fronts at $\gamma _{\boldsymbol {p}}^{p+} t=35.8$ until the drift waves encompass the whole volume at $\gamma _{\boldsymbol {p}}^{p+} t=37.5$, rapidly modifying the zonal profile until a new zonally dominated state can be reinstated at $\gamma _{\boldsymbol {p}}^{p+}t=47$, but which exhibits seeded tertiary modes at $x\approx 17$ and $x\approx 64$ which will eventually repeat this process.

After a period of zonal profile cycling, a stable zonal profile is eventually re-established, which again quenches the drift waves down to the original amplitude. Typically, the tertiary instability is now localised to different points, where seeded drift waves can eventually repeat the process. However, since these points were initially tertiary stable it takes a long time for a localised mode to fully develop. This explains the large swings in transport levels at marginally unstable $\eta$-values observed in figure 6.

As a final remark it is worth noting that, during an entire typical burst process, the box-averaged zonal shear magnitude $|\partial _x^{2}\bar {\varphi }|$ remains comparable to the primary growth rate $\gamma ^{p+}_{\boldsymbol {p}}$. This is a typical result also obtained in previous investigations of zonally dominated states (Waltz, Dewar & Garbet 1998; Kinsey, Waltz & Candy 2005; Kobayashi & Rogers 2012) known as the quench rule.

5. Reduced mode Dimits shift estimate

Having identified the importance of the tertiary instability for the Dimits transition in §§ 3.4, 3.5 and 4, we now turn our attention to the problem of attempting to predict the size of the Dimits shift using tertiary instability analysis. For the system under consideration it is clearly necessary that there exists tertiary stable zonal profiles for the system to be located within the Dimits regime. This has been exploited before by Zhu et al. (2020a) and Zhu, Zhou & Dodin (2020b) in other systems where full stability characterises the Dimits regime to match the Dimits shift threshold to a tertiary transition. However, this matching could only be done in hindsight by tuning a representative zonal curvature value by hand, and did thus not constitute a prediction.

To predict the Dimits shift, the major problem one encounters is, as we outlined in § 3.6, the full multitude of possible zonal profiles, and accounting for how these can be generated through nonlinear interactions with the drift waves. Some profiles may fail to form nonlinearly, while others fail to be robustly stable enough to persist while the residual drift-wave decay. Thus, it may not be sufficient that a tertiary stable zonal profile exists for the system to be in the Dimits regime. Thoroughly accounting for this seems a herculean task and instead some major simplifying assumptions have to be employed.

An example of a simple method doing just that is the heuristic prediction of St-Onge (2017), which relies on the same tertiary 4-mode Galerkin truncation as (3.13a–c) in lieu of accounting for the complex interplay of the many modes constituting the full zonal flow profile. In a modified Terry–Horton system it was postulated that a typical mode could be approximated by that maximally coupled tertiary 4M satisfying the condition $\lambda ^{t+}_{4M}=\lambda ^{t-}_{4M}$ (in our notation). Approximating the coupling of the primary mode to its sidebands through the 4M-interaction alone, and assuming that it is the most unstable primary mode that determines when the Dimits regime ends, the Dimits transition was then taken to occur when this cluster of maximally coupled modes became unstable.

St-Onge (2017) found this prediction to be in excellent agreement with the observed Dimits transition. However, the sensitivity of this transition to numerical dissipation, the small transport levels immediately beyond this point, and the slow evolution made it somewhat difficult to definitively classify a state as stable or turbulent close to this threshold in subsequent reproductions of this system by Zhu et al. (2020b). On the other hand, for our purposes, the major flaw of this prediction is the fact that, with the inclusion of 3 modes for each $\boldsymbol {k}$, there ceases to form anything resembling maximally coupled modes. Nevertheless, because of the simplicity of such a scheme we now look for a similar zonal profile reduction to arrive at a reduced mode prediction.

The first key feature of the present system when attempting to arrive at a simplified prediction is the aforementioned stratification observed in the nonlinear simulations. As mentioned, it is the poloidal band corresponding to the most unstable primary mode $\boldsymbol {p}$ which goes tertiary unstable first, and thus solely determines when the Dimits regime ends. Secondly, we recall the result of §§ 3.5 and 3.6 that the tertiary instability frequently acts in such a way as to push the zonal amplitude $\bar {\varphi }_{\boldsymbol {q}}$ towards its most stabilising value. At least as long as the 4M-interaction is dominant, the zonal profile should therefore repeatedly evolve into a state similar to this one.

For the final piece of the puzzle a much more reductive simplification is employed which relegates this method firmly into being an estimate. We choose to approximate a typical full zonal profile with a single mode, with wavelength $q$, which is the most 4M tertiary stabilising. Of course any real profile will include many more modes of varying amplitude. Nevertheless, some of these will act to destabilise and some to stabilise the tertiary modes, and so we assume that their cumulative effect can be approximated by a representative mode. Indeed the surprising qualitative similarity between the 4M tertiary growth rate in the high $\bar {\varphi }$-limit (3.32) and the local tertiary (3.21) hints that this approximation may be less far fetched than it seems.

Combining these pieces, we conclude that the Dimits transition should occur at around that point $\eta ^{\mathrm {PR}}$ where our single zonal mode ceases to be able to stabilise the most unstable primary mode, leaving no robustly stable zonal state attainable. Thus we can express our Dimits shift prediction $\Delta \eta ^{\mathrm {PR}}$ as

(5.1)\begin{equation} \Delta\eta^{\mathrm{PR}}=\eta^{\mathrm{PR}}-\eta^{0}, \end{equation}

where $\eta ^{\mathrm {PR}}$ is the solution to the constrained optimisation problem

(5.2)\begin{gather} \frac{\partial \gamma^{p+}_{\boldsymbol{p}}}{\partial p}=0, \end{gather}

(5.3)\begin{gather}\frac{\partial\gamma^{t+}_{4M}}{\partial |\bar{\varphi}_{\boldsymbol{q}}|}=0, \end{gather}

(5.4)\begin{gather}\frac{\partial \gamma^{t+}_{4M}}{\partial q}=0, \end{gather}

(5.5)\begin{gather}\gamma^{t+}_{4M}(\eta^{\mathrm{PR}}) = 0. \end{gather}

In figure 8 this method can be seen in action for the configuration of § 4. Although we again stress that (5.1) is clearly a non-rigorous estimate of the Dimits transition, the broader accuracy of which has to be confirmed by comparison with nonlinear simulations we will perform in § 5.1, this predicted Dimits transition at $\eta ^{\mathrm {PR}}=1.86$, corresponding to a Dimits shift of $\Delta \eta ^{\mathrm {PR}}\approx 0.23$, is indeed close to the point $\eta ^{\mathrm {NL}}\approx 1.9$ of figure 5 where the drift waves observed in nonlinear simulations fail to vanish.

Figure 8. (a) The primary instability growth rate $\gamma ^{p+}_{\boldsymbol {r}}$, as given by (3.2), for the sideband modes $\boldsymbol {r}^{\pm }$ of the 4M system (3.13a–c) with poloidal wavenumber (3.3) (note the differing positive and negative scales), (b) the most stabilising zonal amplitude $\bar {\varphi }_{\boldsymbol {q}}$ as given by (5.3), (c) the corresponding 4M-tertiary growth rate $\gamma ^{t+}_{4M}$ as a function of $\eta$ and $q$, all for the configuration of figure 5. The resulting 4M-tertiary Dimits prediction $\eta ^{\mathrm {PR}}$ is indicated by a dashed line, constituting a predicted Dimits shift $\Delta \eta ^{\mathrm {PR}}=\eta ^{\mathrm {PR}}-\eta ^{0}$ of ${\sim }0.23$.

Now, in principle it should be possible to apply the estimation method as just outlined to other systems, up to and including gyrokinetics, so long as one is mindful of what zonal profiles are typically observed and whether some slight modification is necessary to account for these. However, it can only be expected to be useful so long as collisional damping is low enough for the Dimits regime to be characterised by sufficiently stable zonal states, so that the Dimits transition coincides with the point of tertiary destabilisation. Should this not be the case, some other method would have to be employed such as e.g. that of Ivanov et al. (2020), asymptotically accurate in the highly collisional limit, which investigates whether the effect of the turbulence upon the zonal flow either counteracts or reinforces its collisional dissipation.

5.1. Comparison of prediction and nonlinear results

The question now is to what extent the prediction as just outlined in § 5 is generally accurate. To investigate this question we are greatly aided by the observed poloidal stratification and how the dominant primary band is the most tertiary unstable. This means that we can restrict ourselves to only include zonal modes and the most unstable poloidal band in nonlinear simulations, which constitutes a tremendous speedup enabling us to investigate a very wide range of different configurations. This reduction was found to have no appreciable effect on the observed Dimits shift for any of the many disparate test cases, and seems to be uniformly valid for this investigation.

The specific way in which the Dimits shift for a configuration with a given $(\omega _{d0}, \omega _{*0}, \tau ,\eta ^{0})$ was determined can be described as follows. A set of simulations with increasing $\eta$ were performed and were allowed to run continuously until a fully stabilising profile was obtained and all drift-wave amplitudes died down below their original values. If this had not occurred within the time $t_{\mathrm {end}}=3000\gamma ^{-1}$, the simulation was stopped, and the Dimits transition taken to be the final $\eta$ where a stable state arose.

Across thousands of simulations and nearly a hundred configurations, only a handful of times did a simulation with a higher $\eta$ reach stability while one below it failed to do so, and in the few cases for which this occurred all were located right at the Dimits threshold. As outlined in § 3.6, this is because the space of robustly stable zonal configurations rapidly shrinks to 0 at the true Dimits threshold $\eta ^{\mathrm {Di}}$, beyond which no stable states exist. Thus, although it occasionally happens that a stable profile is quickly obtained, the average time $t_{\mathrm {avg}}$ to attain one of these stable states in nonlinear simulations diverges at $\eta ^{\mathrm {Di}}$,

(5.6)\begin{equation} \lim_{\eta \to \eta^{\mathrm{Di}}}t_{\mathrm{avg}}=\infty. \end{equation}

Therefore, the observed Dimits transition $\eta ^{\mathrm {NL}}$, obtained in nonlinear simulations, typically varies very little with respect to $t_\mathrm {end}$ once it is sufficiently large,

(5.7)\begin{equation} \lim_{t_{{\mathrm{end}}}\to \infty}\frac{\partial \eta^{\mathrm{NL}}}{\partial t_{\mathrm{end}}}=0. \end{equation}

With these observations we feel justified in claiming that this method of determining the Dimits shift is robust for our system.

A comparison of the predicted Dimits shift $\Delta \eta ^{\mathrm {PR}}$ and the nonlinearly obtained Dimits shift $\Delta \eta ^{\mathrm {NL}}$ for configurations within the wide range given by $\omega _{d0}\in [-10,-10^{-5}]$, $(\omega _{*0},\tau )\in [-10,-10^{-2}]$, $\eta ^{0}\in [10^{-0.5},10^{2.5}]$, where all these parameters have been simultaneously randomly chosen, can be seen in figure 9. The prediction is seen to generally underpredict the Dimits shift by approximately 15 % barring a few outliers varying by some 10 %. Nevertheless, this is a favourable result because of how very consistent it remains across such a wide range of different configurations, while ultimately being so simple. Indeed no trend for the deviation with respect to any parameter can be discerned.

Figure 9. Dimits shift mismatch between the reduced mode prediction $\Delta \eta ^{\mathrm {PR}}$ of § 5 and what is observed in nonlinear simulations $\Delta \eta ^{\mathrm {NL}}$ as a function of the configuration parameters $\omega _d$, $\omega _*$, $\tau$ and $\eta ^{0}$ for a set of configurations where all parameters were simultaneously randomly chosen. The prediction is generally seen to underpredict the actual shift by some 5 %–30 % but otherwise remain consistent across configurations.

Finally, as mentioned in § 3.2, the nonlinear findings presented in § 4, and which were used to check the accuracy of prediction (5.1), used a constant (with respect to $k_y$) dissipation $D_{\boldsymbol {k}}$. In principle, the validity of these findings could cease to hold for other types of damping, since its inclusion would modify both the prediction and nonlinear Dimits transition slightly. Nevertheless this seems to be an unimportant effect in so far as the matching between the two, just as for St-Onge (2017), appears to remain mostly intact for different reasonable hyperviscosities, as observed across multiple exploratory simulations. Noting this, we feel content that using a constant $D_{\boldsymbol {k}}$ is general enough for our purposes and choose not to investigate this in any further detail.

6. Discussion

What has been demonstrated in this paper is that, through an appropriate asymptotic expansion, the moment hierarchy of gyrokinetics can self-consistently be closed at second order, resulting in a strongly driven gyrofluid system with both linear drive and nonlinear drift/zonal interactions. With the ad hoc introduction of additional damping, meant to encapsulate the effect of Landau damping, this system can be made to exhibit a Dimits shift.

Within the system studied in this paper, the Dimits transition was found to correspond to that point at which the zonal profile can no longer form a robustly tertiary stable profile to eventually kill all drift waves. This is an entirely different mechanism to the one found by Ivanov et al. (2020) in a very similar system. The difference arises because, in that system, being fundamentally collisional, zonal flows decay and no self-stable zonal flows exist. Instead the Dimits regime is characterised by zonal staircases continually rebuilt by momentum influx from small drift-wave turbulence, yielding low mean transport rates. At high collisionality Ivanov et al. (2020) could accurately predict the Dimits transition as that point where momentum flux turned destructive, but this scheme failed at lower collisionality. Likewise, were we to include zonal dissipation in our model our scheme would increasingly fail with increasing collisionality, since it relies on stable zonal states. Thus collisionality is clearly both a qualitatively and quantitatively important parameter to consider when one attempts to study the Dimits regime.

Although similar, the transport bursts observed around the Dimits transition in this system are not strictly speaking the well-known zonal drift-wave predator–prey oscillations observed by e.g. Malkov et al. (2001) or Kobayashi, Gürcan & Diamond (2015), since the zonal energy typically only varies slightly during a drift-wave burst. Furthermore, due to the lack of drift-wave self-interactions it is clear that these bursts do not arise as a result of ferdinons (van Wyk et al. 2016). Instead, it seems that these bursts trace their origin to the movement of tertiary unstable points as the zonal profile is modified, combined with the fact that new localised tertiary modes take so long to emerge, being close to marginal stability within the Dimits regime.

Similarly to the prediction of St-Onge (2017), it is possible to employ a consistently and comparatively accurate estimate of when the Dimits transition should occur by approximating the full zonal profile with a single mode, if it is taken to be that mode which is most stabilising. Naturally, the validity of such a simplification may not in general be taken for granted; one has to be reasonably sure that the system under consideration has small enough collisional damping and that other kinds of nonlinear behaviour do not dominate during the Dimits transition so that the tertiary instability is still of prime importance. Since $\boldsymbol {E}\times \boldsymbol {B}$-shearing by strong zonal flows should remain the dominant nonlinear interaction, however, it seems likely that the Dimits transition continues to coincide with a tertiary stability threshold. Then one might, in more general collisionless systems, be able to approximate the typical zonal profiles with some reduced mode scheme, adapting (5.2)–(5.5) in some simple way, to maintain that computational simplicity that would make our theoretical method of Dimits shift estimation a practically useful and predictive tool. Because of this, future work aims to investigate whether this state of affairs holds for fully gyrokinetic systems.