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. (Reference Plunk, Helander, Xanthopoulos and Connor2014), 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 Reference Dorland and Hammett1993; Hammett et al. Reference Hammett, Beer, Dorland, Cowley and Smith1993) 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 Reference Beer1995), which seems to be of most relevance for the Dimits transition (Dimits et al. Reference Dimits2000). 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 Reference Berionni and Gürcan2011) or non-quasistatic flows (see Kobayashi & Rogers Reference Kobayashi and Rogers2012), 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 Reference Rosenbluth and Hinton1998) 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 Reference Landau1946) 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. Reference Plunk, Helander, Xanthopoulos and Connor2014) and nonlinear (Plunk, Tatsuno & Dorland Reference Plunk, Tatsuno and Dorland2012) 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. Reference Diamond, Itoh, Itoh and Hahm2005), 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 (Reference Farrell and Ioannou2009) 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 Reference Diamond and Kim1991), 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. (Reference Plunk, Tatsuno and Dorland2012). 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. (Reference Ivanov, Schekochihin, Dorland, Field and Parra2020). 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.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. Reference Plunk, Helander, Xanthopoulos and Connor2014). 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 Reference Beer1995)

(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 Reference Kraichnan1967; Qian Reference Qian1986; Terry Reference Terry2004). In three-dimensional turbulence this is prevented by a scale balance of parallel streaming and turbulence, known as critical balance (Barnes, Parra & Schekochihin Reference Barnes, Parra and Schekochihin2011), 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 Reference Sugama1999). 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 Reference Sugama1999; Mishchenko, Plunk & Helander Reference Mishchenko, Plunk and Helander2018). Unfortunately, it is very hard to accurately reproduce the behaviour of these modes in a non-exotic way for a spectral fluid model (Sugama Reference Sugama1999), 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 Reference Kim and Diamond2002). 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. Reference Rogers, Dorland and Kotschenreuther2000). 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 Reference Kobayashi and Rogers2012; Kim, Min & An Reference Kim, Min and An2018, Reference Kim, Min and An2019). 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. Reference Zhu, Zhou and Dodin2020a). 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 Reference Kim and Diamond2002; Zhu, Zhou & Dodin Reference Zhu, Zhou and Dodin2018), and the latter a tertiary KH-like instability, analogous to the secondary instability (3.15) (Rogers et al. Reference Rogers, Dorland and Kotschenreuther2000; Ivanov et al. Reference Ivanov, Schekochihin, Dorland, Field and Parra2020).

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. Reference McMillan, Hill, Bottino, Jolliet, Vernay and Villard2011), and has already been noted in previous tertiary instability studies of simple systems (see e.g. Zhu et al. Reference Zhu, Zhou and Dodin2020a). 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. (Reference Ivanov, Schekochihin, Dorland, Field and Parra2020), 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. (Reference Ivanov, Schekochihin, Dorland, Field and Parra2020). 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. (Reference Rogers, Dorland and Kotschenreuther2000), nor ‘the Dimits regime should end when the zonal amplitude becomes too small’ as Zhu et al. (Reference Zhu, Zhou and Dodin2020a) 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.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 Reference Malkov, Diamond and Rosenbluth2001; Kobayashi & Rogers Reference Kobayashi and Rogers2012). 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 Reference Waltz, Dewar and Garbet1998; Kinsey, Waltz & Candy Reference Kinsey, Waltz and Candy2005; Kobayashi & Rogers Reference Kobayashi and Rogers2012) known as the quench rule.

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 (Reference St-Onge2017), 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.