2.1. Z-pinch geometry

The set-up of the Z-pinch can be expressed in the usual cylindrical coordinates $(r,\theta,z)$ with basis vectors $(\hat {\boldsymbol {r}},\hat {\boldsymbol {\theta }},\hat {\boldsymbol {z}})$. A strong current runs in the longitudinal $\hat {\boldsymbol {z}}$-direction. It gives rise to a azimuthal magnetic field confining the plasma, assumed to have sufficiently low plasma $\beta$ so as not to not alter the longitudinally uniform magnetic vacuum field

(2.1)\begin{equation} \boldsymbol{B}=B(r)\hat{\boldsymbol{\theta}} \quad \mathrm{where}\quad \frac{1}{B(r)}\frac{{\rm d}B(r)}{{\rm d}r}\approx{-}\frac{1}{r}. \end{equation}

Because the plasma is free to rapidly equilibrate along magnetic field lines, the large-scale plasma density $n(r)$ and ion/electron temperatures $T_{i}(r)$ and $T_{e}(r)$ only depend on the radial coordinate $r$. For simplicity, we will also assume these satisfy $T_i=T_e=T(r)$ since a scaling argument reveals this to be the most unstable scenario (Ricci et al. Reference Ricci, Rogers, Dorland and Barnes2006). In the typical gyrokinetic fashion we are then interested in the small-scale fluctuations deriving their energy from the large-scale background plasma gradients. To encapsulate these it is advantageous to employ flux tube simulations (Candy, Waltz & Dorland Reference Candy, Waltz and Dorland2004), so we proceed to the local picture, expanding around the point $r=R$ and setting $B=B(R)$. Then it is advantageous to introduce the local gradient scale lengths

(2.2a,b)\begin{equation} \frac{1}{L_n} ={-} \left.\frac{1}{n(r)}\frac{{\rm d}n(r)}{{\rm d}r}\right|_{r=R}\quad \mathrm{and} \quad \frac{1}{L_T} ={-} \left.\frac{1}{T(r)}\frac{{\rm d}T(r)}{{\rm d}r}\right|_{r=R}, \end{equation}

which are assumed to be greater than the ion gyroradius,  $L_n\sim L_r\gg \rho _i$. We also adopt a more suitable local flux tube coordinate system $(x',y',z')$

(2.3a–c)\begin{equation} \hat{\boldsymbol{x}}'=\hat{\boldsymbol{r}}, \quad\hat{\boldsymbol{y}}'={-}\hat{\boldsymbol{z}}, \quad \hat{\boldsymbol{z}}'=\hat{\boldsymbol{\theta}}, \end{equation}

in which $x'\in (-L_x/2,L_x/2)$ is the radial coordinate, $y'\in (-L_y/2,L_y/2)$ the bilinear ‘poloidal’ coordinate and $z'$ the parallel coordinate, where $L_x,L_y\ll R$ are the flux tube box widths. We will, however, drop the prime notation for future convenience and simply write $(x,y,z)$. Note that, as mentioned previously, parallel variations vanish in the regime of interest and so all quantities will only depend on the perpendicular coordinates $x$ and $y$.

2.2. Gyrokinetics

The collisionless gyrokinetic equation (Catto Reference Catto1978; Frieman Reference Frieman1982; Abel et al. Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013) for the ion-scale dynamics can in Fourier space and dimensionless form here be expressed as (see Plunk et al. Reference Plunk, Helander, Xanthopoulos and Connor2014)

(2.4)\begin{equation} \left( \frac{\partial}{\partial t} + {\rm i} \omega_{ds\boldsymbol{k}} \right) g_{s\boldsymbol{k}} + \left\{ {\varPhi_s} , {g_s} \right\}_{\boldsymbol{k}} = {\rm i} (\omega_{*s\boldsymbol{k}} - \omega_{ds\boldsymbol{k}}) Z_s \varPhi_{s\boldsymbol{k}} f_{0s} \end{equation}

for the component with wave vector $\boldsymbol {k}$. Here, the macroscopic and microscopic spatial dimensions are normalised to $R$ and the ion gyroradius $\rho _i=m_iv_{Ti}/\sqrt {2}eB$, respectively, while the temporal dimension is normalised to the ion streaming time $\sqrt {2}R / v_{Ti}$. The subscript $s=i/e$ denotes the ion/electron species with charge $Z_s e=\pm e$, mass $m_s$ and macroscopic Maxwellian distribution $f_{0s}$ with mean thermal velocity $v_{Ts} = \sqrt { 2 T / m_s}$. The corresponding gyrocentre distribution is given by $g_{s\boldsymbol {k}}={\rm J}_{0s\boldsymbol {k}}\delta f_{s\boldsymbol {k}}R/\rho _i$, where $\delta f_{s\boldsymbol {k}}$ is the fluctuating distribution. Here, the Bessel function of the first kind

(2.5)\begin{equation} {\rm J}_{0s\boldsymbol{k}} = {\rm J}_{0} (\sqrt{2m_s/m_i} k_\bot v_{s\bot}), \end{equation}

where $v_s=v/v_{Ts}$ is the normalised velocity, encapsulates the Fourier space gyroaverage. It also enters through $\varPhi _{s\boldsymbol {k}}={\rm J}_{0s\boldsymbol {k}}\varphi _{\boldsymbol {k}}$, where the dimensionless electrostatic potential $\varphi$ is obtained from the electric potential $\phi$ as $\varphi = e \phi R / T \rho _i$. Naturally, the normalised velocity and wavenumber are both split into their parallel and perpendicular components $v_{s\parallel }, k_\parallel, v_{s\perp }, k_\perp$ with respect to the magnetic field.

Next, in (2.4), the Fourier space Poisson bracket, representing the $E\times B$-nonlinearity, is given by

(2.6)\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. Finally, $\boldsymbol {\nabla } \mathrm {ln} B = \boldsymbol {b}\boldsymbol {\cdot }\boldsymbol {\nabla }\boldsymbol {b}=\hat {\boldsymbol {x}}/R$ since magnetic $\beta$ is small, so the velocity-dependent diamagnetic and magnetic drift frequencies can be expressed as

(2.7a,b)\begin{equation} \omega_{*s\boldsymbol{k}} ={-}\frac{k_yR}{\sqrt{2}L_n} \left( 1 + \eta \left( v_s^{2} -\frac{3}{2} \right) \right)\quad \mathrm{and} \quad \omega_{ds\boldsymbol{k}} ={-} \sqrt{2}k_y \left( v_{s\parallel}^{2} + \frac{v_{s\perp}^{2}}{2} \right), \end{equation}

where $\eta =L_n/L_T$ and the local gradients $L_n$ and $L_T$ are given by (2.2a,b).

To complement the gyrokinetic equation and provide closure, the gyrocentre distribution $g_{s\boldsymbol {k}}$ is related to the potential $\varphi _{\boldsymbol {k}}$ via the quasineutrality condition

(2.8)\begin{equation} \sum_{s=i,e}Z_s\int\,{\rm d}^{3}v_s {\rm J}_{0s\boldsymbol{k}} g_{s\boldsymbol{k}} = \sum_{s=i,e} n \left( 1 - \varGamma_{0s\boldsymbol{k}} \right) \varphi_{\boldsymbol{k}}, \end{equation}

where $\varGamma _{0s\boldsymbol {k}}$ can be expressed in terms of the modified Bessel function ${\rm I}_0$ through

(2.9)\begin{equation} \varGamma_{0s\boldsymbol{k}}=\int \,{\rm d}^{3}v_s{\rm J}_{0s\boldsymbol{k}}^{2}f_{0s}={\rm I}_0\left(\frac{m_sk_\bot^{2}}{m_i}\right) \exp\left({-\frac{m_sk_\bot^{2}}{m_i}}\right). \end{equation}

As a final note, it is plain that the parallel streaming term has deliberately been neglected in (2.4), constituting a significant reduction from spatial three dimensions to spatial two dimensions. This is because the entropy mode, to be described in § 2.3, is uniform in the parallel direction. Although ostensibly a significant simplification, in the regime of interest the entropy mode is massively dominant. In nonlinear simulations where parallel streaming is included slab-mode-driven parallel structures are thus observed to only attain amplitudes of at most $0.5\,\%$ of their uniform counterparts, justifying this simplification.

2.3. The entropy mode

Before we delve further into the details of the tertiary instability we must be familiar with how the primary linear instability is determined. We thus look for solutions with an $\exp (\lambda ^{p}_{\boldsymbol {k}} t)$-time dependence, neglecting the nonlinear mode coupling through the Poisson bracket. The partial time derivative $\partial _t$ in the gyrokinetic equation (2.4) is symbolically replaced by the primary complex frequency $\lambda ^{p}_{\boldsymbol {k}}=\gamma ^{p}_{\boldsymbol {k}}-{\rm i}\omega ^{p}_{\boldsymbol {k}}$, which is then multiplied by ${\rm J}_{0s\boldsymbol {k}}$. The resulting equation is then integrated over velocity space and inserted into the quasineutrality condition (2.8), resulting in the primary dispersion relation

(2.10)\begin{equation} D^{p}(\lambda^{p}_{\boldsymbol{k}},\boldsymbol{k}) = \sum_{s=i,e}\left( \frac{1}{n} \int \,{\rm d}^{3}v_s \frac{{\rm i}(\omega_{*s} - \omega_{ds}) {\rm J}_{0s\boldsymbol{k}}^{2} f_{0s}}{\lambda^{p}_{\boldsymbol{k}}+{\rm i}\omega_{ds}} - \left( 1 - \varGamma_{0s\boldsymbol{k}} \right) \right) = 0, \end{equation}

where the $\boldsymbol {k}$-dependence enters through $\omega _{*s}, \omega _{ds}, {\rm J}_{0s\boldsymbol {k}}$ and $\varGamma _{0s\boldsymbol {k}}$.

Based on the definition (2.7a,b) for $\omega _{*s}$ and $\omega _{ds}$, it is plain that two independent parameters, $R/L_n$ and $\eta$, span the Z-pinch configuration space to uniquely determine $\lambda ^{p}_{\boldsymbol {k}}$.Footnote ¹ $\omega _{*s}$ and $\omega _{ds}$, being proportional to $k_y$, furthermore makes (2.10) break down for modes with $k_y=0$ unless $\lambda ^{p}_{\boldsymbol {k}}=0$. But these modes make up the zonal potential

(2.11)\begin{equation} \bar{\varphi}=\frac{1}{L_y}\int_0^{L_y}\varphi \,{{\rm d}y}, \end{equation}

and the zonal $E\times B$-flow $\partial _x\bar {\varphi }$. Thus zonal flows are, as previously mentioned, indeed linearly conserved.

Returning to (2.10), a key feature is the inclusion of kinetic electrons in lieu of a simplified adiabatic response as is frequently employed for efficiency in gyrokinetic toroidal ITG simulations (Dorland & Hammett Reference Dorland and Hammett1993) where the Dimits shift has conventionally been studied (Dimits et al. Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000). This is necessary because the dominant electrostatic instability at small gradients, originally considered by Kadomtsev (Reference Kadomtsev1960), is the so-called entropy mode. This mode, like the magnetohydrodynamic (MHD) interchange mode (Rosenbluth & Longmire Reference Rosenbluth and Longmire1957), exchanges $n$ and $T$ to modify the entropy $S=\ln (T^{5/2}/p)$ while leaving the pressure $p=nT$ unchanged. It is the plasma analogue of a fluid thermal convective instability, enabled by oppositely directed displacements of ions/electrons, and so cannot be encapsulated when the electron dynamics is excluded (Ware Reference Ware1962). Despite existing even in the absence of a temperature gradient, i.e. $\eta =0$, this mode is occasionally also referred to as the drift-temperature-gradient mode (see e.g. Kesner Reference Kesner2000), since it shares the same drive term as the ITG mode. Thus, it is unsurprising that its associated turbulence, also subject to the same $E\times B$-nonlinearity, can exhibit that shear stabilisation which characterises the Dimits shift (Kobayashi & Rogers Reference Kobayashi and Rogers2012).

As to the solution of (2.10), since it is a kinetic integral equation that cannot be solved analytically, it must be treated numerically. Although one can derive a comparable dispersion relation in a gyrofluid model, both Ricci et al. (Reference Ricci, Rogers, Dorland and Barnes2006) and Kobayashi & Rogers (Reference Kobayashi and Rogers2012) note that it significantly overestimates the stability threshold. We must therefore solve (2.10), but instead of explicitly doing so it is convenient to initialise a random state and let it evolve according to the gyrokinetic equation (2.4). The largest growth rate component will then eventually dominate, allowing for easy determination of $\gamma ^{p}$ by fitting an exponential to the system's long-term evolution as outlined in § 4.

The result of such a calculation, using the Z-pinch implementation (Bañón Navarro, Teaca & Jenko Reference Bañón Navarro, Teaca and Jenko2016) of GENE (Jenko et al. Reference Jenko, Dorland, Kotschenreuther and Rogers2000), can be seen in figure 1, where (a) shows the range of gyrokinetic instability compared with the gyrofluid result of Ricci et al. (Reference Ricci, Rogers, Dorland and Barnes2006). It confirms that the critical gradients are indeed lower than the gyrofluid value, and indicates that $R/L_n$ is a ‘good’ instability parameter, i.e. one that uniformly destabilises the plasma, unlike $\eta$ or $R/L_T$. In (b,c), on the other hand, the entropy-mode growth rate is shown as a function of the wavenumber $\boldsymbol {k}$ for $\eta =0.25$ when $R/L_n=1.4$ and $R/L_n=1.8$, which, as we will see in § 3.1, are within the Dimits regime and at the Dimits threshold, respectively. The result is typical of the Dimits regime, with the mode exhibiting a wide unstable range. However, the growth rate also exhibits a sharp peak around some radial streamer, here $\boldsymbol {k}\approx (0,0.6)$, while other modes constitute a broader, lower growth rate shoulder.

Figure 1. (a) Region of collisionless gyrofluid entropy-mode instability (blue) and collisionless gyrokinetic entropy-mode instability (green). The latter is seen to become unstable at significantly lower gradients. Shown in (b,c) is the $\eta =0.25$ entropy-mode (primary) growth rate $\gamma ^{p} R/c_s$ as a function of the radial and poloidal wavelengths $k_x,k_y$ for (b) $R/L_n=1.4$ and (c) $R/L_n=1.8$. On increasing $R/L_n$, the unstable range is seen to widen significantly, but its maximum remains mostly stationary, peaked around $k_y=0.6$.

2.4. The interchange mode

So far we have been focusing on the gyrokinetic equation and the resulting entropy mode, but the full picture is more complex. At larger gradients there are several MHD instabilities present, the first of which to appear as gradients are increased is the ideal interchange mode (Rosenbluth & Longmire Reference Rosenbluth and Longmire1957). Of course, we are technically concerned with the completely collisionless case where an MHD treatment is not justified, but the collisionless analogue of this instability can be found by employing the Chew–Goldberger–Low model (Chew, Goldberger & Low Reference Chew, Goldberger and Low1956). It possesses a very similar dispersion relation with a simple instability criterion given by

(2.12)\begin{equation} \frac{R}{L_n}\geq \frac{7}{2(1+\eta)}, \end{equation}

when electrons/ions have the same temperature (Ricci et al. Reference Ricci, Rogers, Dorland and Barnes2006).

For our purposes, the interchange mode with its characteristic $\boldsymbol {k}$-independent growth rate eliminates the possibility of self-generated small-scale and small-amplitude zonal flow stabilisation at larger gradients. Simulations in this range are therefore ill-behaved, exhibiting secular growth at large scales with non-convergent fluxes, and so we dare not attempt to extract any information from these simulations. Still, Z-pinch MHD-like instabilities like this could, in principle, also be shear-flow stabilised. For example, it was predicted that the $m=1$ kink mode (Newcomb & Kaufman Reference Newcomb and Kaufman1961) could be stabilised (Shumlak & Hartman Reference Shumlak and Hartman1995), which was also experimentally confirmed (Shumlak et al. Reference Shumlak, Golingo, Nelson and Den Hartog2001, Reference Shumlak, Adams, Blakely, Chan, Golingo, Knecht, Nelson, Oberto, Sybouts and Vogman2009). These flows are, however, of a different type, being large scale and externally imposed instead of small scale and self-generated. As such, they affect the background equilibrium, and so we will not be considering the range given by (2.12).

3.1. The Dimits shift

Turbulent bursts like those of § 3 have continually plagued investigations of the Dimits shift by rendering time-resolved statistics of transport frustratingly imprecise and possibly misleading (Pueschel & Jenko Reference Pueschel and Jenko2010). In essence, statistics averaged over different, similar sections of a burst cycle yield very different results. This can make it particularly difficult to determine whether a given configuration belongs to the Dimits regime or not (for example, compare St-Onge Reference St-Onge2017; Zhu et al. Reference Zhu, Zhou and Dodin2020b), which unfortunately also holds true here. Although variance will remain high, the only remedy for a more unambiguous determination of the Dimits shift is to average over much longer time intervals that encapsulate multiple complete burst cycles. The particle flux $\varGamma$ thus obtained can be seen in figure 3 when $\eta =0.25$. $\varGamma$ is seen to increase exponentially with increasing $R/L_n$. However, a discontinuity from a negligible level of ${\sim } 10^{-2}\varGamma _{{\rm GB}}$ to the gyro-Bohm value ${\sim }\varGamma _{{\rm GB}}$ is clear at $R/L_n\approx 1.8$, which we interpret to correspond to the Dimits threshold $R/L_D$.

Figure 3. In black: particle flux $\varGamma$ as a function of the density gradient $R/L_n$ for $\eta =0.25$. In blue: the fraction $\varTheta _{0.1}$ of the simulation time after the initial transport peak where the transport level is higher than 10 % of the initial transport peak, i.e. the fraction of time spent outside of quiescence. The transport is well described by two exponential fits (dashed red) of slopes $4\pm 0.3$ and $4.9\pm 0.4$, discontinuously jumping between the two at $R/L_n\sim 1.8$. This corresponds to that same point at which $\varTheta _{0.1}$ becomes greater than 0 and will be identified as the Dimits threshold.

To provide a quantitative measure of turbulent bursts we introduce

(3.2)\begin{equation} \varTheta_\kappa=\lim_{t\rightarrow\infty}\frac{1}{t}\int_{t_i}^{t} \,{\rm d}t' H\left(\varGamma-\kappa\varGamma(t_i)\right), \end{equation}

where $H$ is the Heaviside function and $t_i$ is the time when the initial, linear growth phase ceases, i.e. $\varGamma (t_i)$ is taken to be a representative level of turbulent transport. $\varTheta _\kappa$, measuring the fraction of time the transport level exceeds $\kappa \varGamma (t_i)$, is also plotted in figure 3 when $\kappa =0.1$. It is seen that this fraction begins to increase from $0$ at $R/L_D$, gradually increasing until full turbulence ($\varTheta _{0.1}\sim 1$) develops well above. The picture is the same for $0.01\lesssim \kappa \lesssim 0.3$, the choice of which only affects how rapidly $\varTheta _\kappa$ increases above $R/L_D$. Although initially rare, the emergence of bursts causing appreciable transport of up to ${\sim }0.3\varGamma (t_i)$ thus characterises the Dimits threshold. This clear signal was therefore used to unambiguously determine $R/L_D$.

Still analysing figure 3, the observed exponential growth of $\varGamma$ with respect to $R/L_n$ is slightly different above and below $R/L_D$, $\partial (\log \varGamma /\varGamma _{{\rm GB}})/\partial (R/L_n)=4\pm 0.3$ and $4.9\pm 0.4$, respectively. This is because turbulent bursts with $\kappa \gtrsim 0.01$ dominate $\varGamma$ above $R/L_D$. Below this point, as we will describe in § 4, $\varGamma$ is driven by drift modes that couple together quasilinearly into marginally (un)stable tertiary modes. These can be rendered stable or unstable by even a small perturbation, like the dynamically evolving $\nu _x, \nu _y$. That this is enough to sustain small transport levels has been confirmed by running simulations where $\nu _x, \nu _y$ were kept constant at comparable values to those obtained through the gyroLES procedure. Then $\varGamma =0$ was eventually obtained for $R/L_n\lesssim b$, where the value $1.2\lesssim b\lesssim 1.6$ depended on the choice of $\nu _x, \nu _y$.

One might be worried that $\nu _x, \nu _y$ similarly affects $\varGamma$ at and above the Dimits threshold $R/L_D$. Driven by turbulent bursts, it was, however, found to only be significantly affected when $\nu _x, \nu _y$ were increased much further. This $\varGamma$-insensitivity at $R/L_D$ echoes collisional findings (Weikl et al. Reference Weikl, Peeters, Rath, Grosshauser, Buchholz, Hornsby, Seiferling and Strintzi2017), but there is still cause for alarm since bursts resemble the continuous turbulence at larger $R/L_n$-values. Unfortunately, this turbulence is prone to an unbalanced inverse energy cascade known to plague two-dimensional turbulence (Kraichnan Reference Kraichnan1967; Qian Reference Qian1986; Terry Reference Terry2004). The resulting pile up, as energy tries to transfer towards unresolved large scales, means that turbulent transport levels cannot be trusted. Thankfully, the collisionless Dimits regime is characterised by marginally unstable ‘quasilinearly’ evolving states, which are resolved. That this includes $R/L_D$, i.e. the point where these states fail to remain uniformly quiescent, is confirmed in figure 4. It clearly demonstrates that the discontinuous $\varGamma$-increase at $R/L_n$ remains persistent, and consistent, as the resolution is increased.

Figure 4. Nonlinear particle flux as a function of the density gradient around the observed Dimits threshold $R/L_D\sim 1.8$ of figure 3 for different numerical implementations: standard runs as used in this article (black squares), runs with twice as large box length (blue triangles) and runs with twice the $\boldsymbol {k}$-space resolution (red circles) (the latter two are visually offset for clarity). Here, the uncertainty is calculated by splitting the full runs into 10 smaller pieces. Since the observed Dimits threshold is seen to be consistent across configurations, the standard configuration employed is seen to be sufficient to determine the Dimits threshold.

3.2. Zonal shear

A common quantitative measure of the stabilising effect of zonal flows for turbulence suppression is the box-averaged zonal shear magnitude $\langle |\partial _x^{2}\bar {\varphi }|^{2}\rangle ^{1/2}$ (Waltz et al. Reference Waltz, Kerbel and Milovich1994; Diamond et al. Reference Diamond, Itoh, Itoh and Hahm2005). For example, it has been noted that when turbulence is strongly suppressed the quench rule $\langle |\partial _x^{2}\bar {\varphi }|^{2}\rangle ^{1/2}\sim \gamma ^{p}$ holds (Waltz, Dewar & Garbet Reference Waltz, Dewar and Garbet1998; Kinsey, Waltz & Candy Reference Kinsey, Waltz and Candy2005). However, as can be seen in figure 5, the nonlinear shear significantly exceeds $\gamma ^{p}$, attaining a value of ${\sim } 1$ in the Dimits regime. It is possible that part of this discrepancy can be traced to zonal dissipation, which would push the flow towards quench rule levels, where turbulence, as found by Ivanov et al. (Reference Ivanov, Schekochihin, Dorland, Field and Parra2020), acts to reinforce the zonal flow in the Dimits regime. Indeed, Kobayashi & Rogers (Reference Kobayashi and Rogers2012) also observed that zonal shear levels decreased with increasing collisionality in the Z-pinch Dimits regime until quench rule levels heralded continuous turbulence.

Figure 5. Time- and box-averaged nonlinear zonal shear $\langle |\partial _x^{2}\bar {\varphi }|^{2}\rangle ^{1/2}$ after the transient period as a function of the density gradient $R/L_n$ for $\eta =0.25$ (blue). Also plotted is the entropy-mode growth rate $\gamma ^{p}$ (green), which is significantly lower. The zonal shear range where a sinusoidal zonal profile can be stable (approximately obtained by using a strongly stabilising $k_x=0.4$-mode as discussed in § 4) is furthermore shown in red. It is seen that the nonlinear shear rate remains relatively constant throughout the Dimits regime, at a level around the lower boundary of sinusoidal profile stability, before increasing thereafter.

Another heuristic estimate involving zonal shear gives rise to simple scalings of the kind $\varGamma \propto 1-\epsilon \langle |\partial _x^{2}\bar {\varphi }|^{2}\rangle ^{1/2}/\gamma ^{p}$ with some $\epsilon$ values that are sometimes employed to model zonal shear quenching (Waltz et al. Reference Waltz, Dewar and Garbet1998; Kinsey et al. Reference Kinsey, Waltz and Candy2005). Such estimates imply that only when $\langle |\partial _x^{2}\bar {\varphi }|^{2}\rangle ^{1/2}\gtrsim \gamma ^{p}$ will zonal flows have a significant effect on turbulent saturation (Xanthopoulos et al. Reference Xanthopoulos, Merz, Görler and Jenko2007; Pueschel & Jenko Reference Pueschel and Jenko2010). As is clear in figure 5, this presently holds, but $\langle |\partial _x^{2}\bar {\varphi }|^{2}\rangle ^{1/2}$ continues to increase faster than $\gamma ^{p}$ above the Dimits threshold. This behaviour is inconsistent with a a scaling of $\varGamma \propto 1-\epsilon \langle |\partial _x^{2}\bar {\varphi }|^{2}\rangle ^{1/2}/\gamma ^{p}$, which would predict that increasing $R/L_n$ will eventually eliminate transport. Even accounting, in the same way as Hahm et al. (Reference Hahm, Beer, Lin, Hammett, Lee and Tang1999), for the fact that the effective shearing decorrelation rate decreases when the zonal flow varies rapidly in time does not change this observation.

Now an increase of $\langle |\partial _x^{2}\bar {\varphi }|^{2}\rangle ^{1/2}$ as the driving gradients are increased like the one here has also been observed in ITG simulations by Terry et al. (Reference Terry, Li, Pueschel and Whelan2021). This was, however, noted as being inconsistent with the tertiary picture of the Dimits shift, seemingly under the assumption that greater zonal shear should be more stabilising. But that this need not be the case has already been observed (see e.g. Kinsey et al. Reference Kinsey, Waltz and Candy2005; Kobayashi & Rogers Reference Kobayashi and Rogers2012). To presently reaffirm this, figure 5 also depicts the approximate range of zonal shear amplitudes where a sinusoidal profile can be completely stable (determined as outlined in §§ 4 and 5). A clear upper boundary where stable profiles cease to exist is observed. Consistent with the tertiary picture, this boundary intersects with the lower boundary at approximately the Dimits threshold, leaving no stable profiles at all.

3.3. Localisation

A persistent feature of marginal drift waves in the presence of a zonal flow is their localisation to regions of zero zonal shear, i.e. where $\partial ^{2}_x\bar {\varphi }=0$. This is in accordance with the intuitive picture that the zonal flow is unable to decorrelate radial streamers at such points, effectively only transporting them in the poloidal $y$-direction (Ivanov et al. Reference Ivanov, Schekochihin, Dorland, Field and Parra2020). However, this necessary condition is not sufficient for drift waves to congregate, since they respond asymmetrically to minima $\partial _x^{3}\bar {\varphi }>0$ and maxima $\partial _x^{3}\bar {\varphi }<0$; significant localisation is only observed at minima (McMillan et al. Reference McMillan, Hill, Bottino, Jolliet, Vernay and Villard2011; Zhu et al. Reference Zhu, Zhou and Dodin2020a). This picture is presently reaffirmed in figure 6, showing that the zonal-averaged drift wave density $\langle \tilde {n}^{2}\rangle ^{1/2}$ is clustered at precisely these points. We note that that there are only a few of them since the zonal flows form staircase states (Dif-Pradalier et al. Reference Dif-Pradalier, Diamond, Grandgirard, Sarazin, Abiteboul, Garbet, Ghendrih, Strugarek, Ku and Chang2010; Peeters et al. Reference Peeters, Rath, Buchholz, Camenen, Candy, Casson, Grosshauser, Hornsby, Strintzi and Weikl2016), i.e. broader rectangular shapes of nearly constant shear with sharp transitions from positive to negative values. These clearly become more developed and broad as $R/L_n$ approaches the Dimits threshold from below. We also note for later that the ion/electron temperature moments predominantly self-organise in such a way to predominantly have opposite sign, something we will clarify further in §§ 4.3–4.4.

Figure 6. Mean drift wave density perturbation $\tilde {n}$ as a function of the radial coordinate $x$ in the presence of the quasistationary zonal flow $\partial _x\bar {\varphi }$ and zonal fluctuating temperature $\bar {T}$ for $\eta =0.25$ and (a) $R/L_n=1.2$, (b) $R/L_n=1.5$ and (c) $R/L_n=1.8$. As can be seen, the ion temperature (blue) and the electron temperature (red, dotted) tend to display opposite signs over a majority of the region. Meanwhile, it can be seen that as the gradient is increased towards the Dimits threshold the stable zonal flow by necessity increasingly resembles a staircase state with broadening ‘steps’ of nearly constant zonal shear $\partial ^{2}_x\bar {\varphi }$. Meanwhile, the drift waves increasingly localise around zonal flow minima, i.e. where the zonal shear $\partial ^{2}_x\bar {\varphi }$ vanishes and $\partial ^{3}_x\bar {\varphi }>0$.

It was recently conjectured by Zhang & Krasheninnikov (Reference Zhang and Krasheninnikov2020) that localising well-defined staircases like these may provide effective transport barriers of drift waves and drifting coherent structures like ferdinons (Ivanov et al. Reference Ivanov, Schekochihin, Dorland, Field and Parra2020), reflecting and trapping them in the vicinity of zonal extrema. Although coherent structures are not observed, figure 7 shows that, during quiescent periods of zonal profile cycling, this seems to be in effect. However, it can also be seen that when transport flares up drift waves encompass the full volume. Meanwhile, tertiary streamer modes remain localised, only marginally broadening. This raises doubts about the present validity of the further claim of Zhang & Krasheninnikov (Reference Zhang and Krasheninnikov2020) that localisation may be more important than the zonal shearing-reduced streamer growth rate for transport quenching. Furthermore, closer to marginal stability turbulence essentially dies away, leaving only very low amplitude tertiary streamers like those depicted in figure 7 with minimal self-interactions. In § 4, we will see that this must be explained by the fact that their drive then is insufficient to overcome zonal shear.

Figure 7. A snapshot from a simulation with $\eta =0.25$ and $R/L_n=1.8$ during a (a) quiescent and a (b) turbulent period. Plotted is the drift wave potential $\tilde {\varphi }$ together with $\partial _x^{2}\bar {\varphi }$ (solid black) and $\partial _x^{3}\bar {\varphi }$ (dashed red) of the zonal potential $\bar {\varphi }$, with the corresponding most unstable tertiary drift wave modes in the lower-most panel. For the zonal profile of (a) the drift waves are localised around the same points as the tertiary modes that feed them while for that of (b) they fill the space more uniformly while the tertiary modes remain localised.

4.1. Tertiary instability simulations

When (4.10) is expanded, the result involves intractable products of integrals like (2.10) with ever more complex integrands. Since (2.10) must already be solved numerically, the same holds true here, but it becomes exceedingly cumbersome to do so directly as $q_{G}$ is increased. To remedy this, the alternate procedure of instead evolving the gyrokinetic equation directly as mentioned in § 2.3 is instead employed.

A zonal profile, i.e. the set of $2(2j+1)$ zonal distributions $\bar {g}_{s\boldsymbol {k}}$ with potential modes $\bar {\varphi }_{\boldsymbol {k}}$, is either extracted from nonlinear simulations or initialised from scratch into some desired configuration. Then a single poloidal wavenumber $p$ is chosen and a new nonlinear simulation, with a Fourier decomposition of $2j+1$ $k_x$-modes and $2$ $k_y$-modes ($k_y=0$ and $k_y=p$), is initialised with this zonal profile and small-amplitude drift waves. The system is then allowed to evolve nonlinearly, but with the zonal profile frozen in place. Since drift wave self-interactions are not present by the triplet condition (4.1), the linear tertiary growth rate eventually converges. For this to be stable, some small $\nu _x$ must, however, be included. Thankfully, the specific value is generally unimportant for present purposes unless specific comparison with nonlinear observations must be made.

To begin investigating the present tertiary instability using the aforementioned scheme, in figure 8 the tertiary growth rate $\gamma ^{t}$ of different zonal profiles, extracted from nonlinear simulations, are compared with the primary growth rate $\gamma ^{p}$. The profiles are obtained just above the $\eta =0.25$ Dimits threshold $R/L_D\approx 1.8$, where turbulent bursts first manifest. It is seen that $\gamma ^{t}$ is much greater for zonal profiles of turbulent periods compared with their marginally unstable quasistationary counterparts with $\gamma ^{t}\sim 0.01\gamma ^{p}$. However, in both cases, the growth rate is seen to peak around the same value $p\approx 0.6$ where $\gamma ^{p}$ attains its maximum.

Figure 8. Tertiary growth rates of different, randomly selected zonal profiles, obtained from nonlinear simulations with $\eta =0.25$ and $R/L_n=1.9$, as a function of the poloidal wavenumber $p$. Blue lines correspond to zonal profiles taken during turbulent periods, while red (inlaid) lines correspond to profiles taken during quiescent periods. Additionally, the primary growth rate $\gamma ^{p}$ is plotted in green for comparison. As can be seen, the quiescent profiles are only unstable for a few $p$-values clustered around ${\sim }$0.6, the most primary unstable point, with growth rates severely reduced to approximately 1 % of the primary growth rates. In comparison, the turbulent profiles are much more unstable to a broader range of $p$-values.

Having observed that quasistationary zonal flows are very nearly tertiary stable, the question becomes how robust this stability is. For the fluid model under consideration in the previous work of Hallenbert & Plunk (Reference Hallenbert and Plunk2021) this concept was key. There, frequent turbulent bursts sustained drift waves at sufficient amplitude to affect stable zonal profiles, markedly altering them over time. Therefore, to prevent subsequent turbulent bursts, the zonal profile had to not only be individually tertiary stable, but also belong to a family of similar profiles that also were stable. Here, the necessity of such stringent stability is smaller. As observed in figure 2, the Z-pinch exhibits very little drift wave activity over extended periods in the Dimits regime, during which the zonal flow remains very nearly unchanged. As such, individual profile stability is therefore re-emphasised.

To verify that individual profile stability is sufficient, figure 9 depicts $\gamma ^{t}$ of different zonal flow profiles when a rescaling factor $a$ is used to (artificially) adjust the zonal flow amplitude. Here, $\gamma ^{t}$ is observed to increase relatively swiftly as $a$ is moved away from $1$ for quasistationary profiles. On the other hand, this process may stabilise the turbulent profiles somewhat, although not to comparable levels of the quasistationary profiles. We interpret this as meaning that quasistationary flows, in being established when turbulent bursts end, are delicately pushed towards configurations of greater tertiary stability. Note again that increasing the zonal shear, i.e. increasing $a$, does not necessarily imply that $\gamma ^{t}$ decreases. This fact was already noted by Kobayashi & Rogers (Reference Kobayashi and Rogers2012), who, however, misattributed this fact to KH-type zonal flow breakup instead of their reduced efficacy in stabilising drift waves, something we will expand upon in § 5.1.

Figure 9. Tertiary growth rate $\gamma ^{t}$, normalised to the primary growth rate $\gamma ^{p}$, of different zonal flow profiles obtained from nonlinear simulations with $\eta =0.25$ and $R/L_n=2.0$ as their amplitude is multiplied by a factor $a$. Quiescent zonal profiles (red triangles) are seen to typically be of very nearly that amplitude which is most stabilising, i.e. when $a=1$. The same is not true for profiles from turbulent periods (blue squares), which, beyond generally being less stabilising, may be made more stabilising if their amplitude is altered.

4.2. Tertiary instability of a single zonal mode

To gain some basic understanding of the tertiary instability, we now turn to the minimal case of a sinusoidal zonal profile, i.e. a single zonal mode $\bar {\varphi }_{\boldsymbol {q}}$ which we take to be real. As is apparent from the presence of $\bar {g}_{s\boldsymbol {q}}$ in the definition (4.9a,b) of $\boldsymbol{\mathsf{C}}_2$, the zonal flow $\partial _x\bar {\varphi }$ alone is insufficient to determine $\gamma ^{t}$; sideband coupling is modified by $\bar {g}_{s\boldsymbol {q}}$. One reason for this is because $\bar {g}_{s\boldsymbol {q}}$ modifies the background gradients, but, in principle, any kinetic modification of $\bar {g}_{s\boldsymbol {q}}$ also affects the coupling. This fact, necessarily not included in previous studies of fluid systems like Zhu et al. (Reference Zhu, Zhou and Dodin2020a), Ivanov et al. (Reference Ivanov, Schekochihin, Dorland, Field and Parra2020) or Hallenbert & Plunk (Reference Hallenbert and Plunk2021), has perhaps not been adequately appreciated. To investigate this effect in a controlled way, we will in the next sections vary the zonal distribution $\bar {g}_{s\boldsymbol {q}}$, holding $\bar {\varphi }_{\boldsymbol {q}}$ fixed by employing a multiplicative factor.

Before exploring this in detail, we can immediately deduce one way in which $\bar {g}_{s\boldsymbol {q}}$ does not affect $\gamma ^{t}$. In Appendix A we show that the tertiary growth rate of a single zonal mode $\bar {\varphi }_{\boldsymbol {q}}$ is unchanged under the transformation $\bar {g}_{s\boldsymbol {q}} \rightarrow \bar {g}_{s\boldsymbol {q}}^{*}\bar {\varphi }_{\boldsymbol {q}}^{2}/|\bar {\varphi }_{\boldsymbol {q}}|^{2}$, which also leaves $\bar {\varphi }_{\boldsymbol {q}}$ unchanged. Expressing $\int \,\textrm {d}^{3}v_s\textrm {J}_{0s\boldsymbol {q}}\bar {g}_{s\boldsymbol {q}}$ as $\varrho _s\textrm {e}^{\textrm {i}\vartheta _s}\bar {\varphi }_{\boldsymbol {q}}$, this transformation acts to transform the phase like $\vartheta _s\rightarrow -\vartheta _s$. However, quasineutrality (2.8) makes it clear that the distribution parts $\bar {g}_{s\mathrm {out}}$ that are ${\rm \pi} /2$ out of phase with $\bar {\varphi }_{\boldsymbol {q}}$ are opposite, i.e. $\operatorname {Im}{(\varrho _i\textrm {e}^{\textrm {i}\vartheta _i})}=-\operatorname {Im}{(\varrho _e\textrm {e}^{i\vartheta _e})}$. One species distribution therefore trails the potential while the other leads it, and we conclude that $\gamma ^{t}$ does not depend upon which is which. This result can be interpreted in light of the tertiary localisation highlighted in § 3.3; $\bar {g}_{s\mathrm {out}}$ satisfies $\partial _x\bar {g}_{s\mathrm {out}}=0$ at the points $\partial ^{2}_x\bar {\varphi }=0$ that are tertiary unstable, and so are not able to affect the tertiary instability by modifying the background gradients there.

The preceding discussion did not depend upon the Galerkin truncation. However, before we further investigate the role of $\bar {g}_{s\boldsymbol {q}}$ in § 4.3 for $\gamma ^{t}$, it is prudent to investigate its effect. Thus figure 10 depicts how $\gamma ^{t}$, obtained as outlined in § 4.1, changes as $q_{G}$ is varied for some example distributions $\bar {g}_{s\boldsymbol {q}}$ with $\bar {\varphi }_{\boldsymbol {q}}=3.5$ (a typical nonlinear value) and $\bar {\varphi }_{\boldsymbol {q}}=35$. Also plotted for comparison is the primary sideband growth rate $\gamma ^{p}_{\boldsymbol {p}+\boldsymbol {q}_{G}}$. For small to ‘normal’ zonal amplitudes $\gamma ^{t}>\gamma ^{p}_{\boldsymbol {p}+\boldsymbol {q}_{G}}$ seems to uniformly hold in accordance with the view that the tertiary instability mixes primary modes (Hallenbert & Plunk Reference Hallenbert and Plunk2021). When significantly increased, $\gamma ^{t}<\gamma ^{p}_{\boldsymbol {p}+\boldsymbol {q}_{G}}$ is occasionally found to hold when $\boldsymbol {q}_{G}$ is small. Presumably this is explained by immediate subdominant mode coupling (see e.g. Terry, Baver & Gupta Reference Terry, Baver and Gupta2006; Hatch et al. Reference Hatch, Terry, Jenko, Merz and Nevins2011), the matching of which becomes less resonant when more modes are included, i.e. when $\boldsymbol {q}_{G}$ is increased.

Figure 10. Single zonal mode tertiary growth rate $\gamma ^{t}$, normalised $\gamma ^{p}$, of some example distributions when $q_{G}$ is varied. Here, $\eta =0.25$, $R/L_n=1.8$, $q=0.4$ and $p=0.6$, while the distributions are of the form (4.11a,b) with (a) $\bar {\varphi }_{\boldsymbol {q}}=3.5$, $\alpha _1+\alpha _2 \textrm {i}=(0.5, 2+\textrm {i}, 0, -1)$, and (b) $\bar {\varphi }_{\boldsymbol {q}}=35$, $\alpha _1+\alpha _2 \textrm {i}=(1+\textrm {i}, 0.4+\textrm {i}, 0.1+0.25\textrm {i}, -1+\textrm {i})$. Also plotted in green is the primary growth rate $\gamma ^{p}_{\boldsymbol {p}+\boldsymbol {q}_{G}}$ of the sideband mode with $k_x=q_{G}$. The configurations are chosen to give a wide spread of converged $\gamma ^{t}$, which occurs around $q_{G}/q\sim 10$ and $q_{G}/q\sim 30$ for $\bar {\varphi }_{\boldsymbol {q}}=3.5$ and $\bar {\varphi }_{\boldsymbol {q}}=35$, respectively.

Continuing, we note that the point $q_{G}$ at which $\gamma ^{t}$ converges effectively constitutes a measure of how many sidebands are coupled together into a combined tertiary mode. Intuitively, it increases with increasing $\bar {\varphi }_{\boldsymbol {q}}$, since in the opposite weak coupling limit $\bar {\varphi }_{\boldsymbol {q}}\rightarrow 0$ different wavenumbers decouple. For the typical zonal amplitude of figure 10(a) this point occurs at $k_x\approx 12q=4.8<6.4$, which as stated in § 3 is the smallest scale included in our simulations, offering some further evidence that our numerical implementation was sufficient to capture the dynamics of the Dimits regime.

4.3. Role of ion/electron response phase

The question now is how the choice of $\bar {g}_{s\boldsymbol {q}}$, as uncovered in § 4.2, affects $\gamma ^{t}$. As a first facet, we choose to investigate the relative phase of the ion/electron response. For this, both distributions are taken to be simple representative Maxwellians, i.e.

(4.11a,b)\begin{equation} g_{e\boldsymbol{q}} = \frac{M(\alpha_1+\alpha_2 {\rm i})}{({\rm \pi}\alpha_3)^{3/2}} {\rm e}^{{-}v_e^{2}/\alpha_3} \quad \mathrm{and} \quad g_{i\boldsymbol{q}} = \frac{M}{({\rm \pi}\alpha_4)^{3/2}} {\rm e}^{{-}v_i^{2}/\alpha_4}, \end{equation}

where $M$ is a multiplicative factor that ensures $\bar {\varphi }_{\boldsymbol {q}}$ attains its desired value. It can be calculated using quasineutrality (2.8) as

(4.12)\begin{equation} M = \frac{\bar{\varphi}_{\boldsymbol{q}}(2-\varGamma_{0i\boldsymbol{q}}-\varGamma_{0e\boldsymbol{q}})}{\int\, {\rm d}^{3} v_i{\rm J}_{0i\boldsymbol{q}}g_{i\boldsymbol{q}}/M - \int \,{\rm d}^{3} v_e{\rm J}_{0e\boldsymbol{q}}g_{e\boldsymbol{q}}/M}. \end{equation}

Here

(4.13)\begin{equation} \frac{n_{ge\boldsymbol{q}} }{ n_{gi\boldsymbol{q}}} = \frac{\int \,{\rm d}^{3}v_e g_{e\boldsymbol{q}} }{ \int \,{\rm d}^{3}v_i g_{i\boldsymbol{q}}} = \alpha_1+\alpha_2{\rm i} \end{equation}

is the cross-species phase whose influence we want to investigate. On the other hand, $\alpha _3=\alpha _4=0.2$ are ‘effective temperatures’, the choice of which will be explained in § 5.

Figure 11 shows how $\gamma ^{t}$, of the poloidal band $p=0.6$ at the Dimits threshold $\eta =0.25$ and $R/L_n=1.8$, changes as $\alpha _1/\alpha _2$ are varied. Four different cases, corresponding to $\varphi _{\boldsymbol {q}}=0.35, 3.5, 35$ and $350$, are shown, each exhibiting a radically different appearance (note the usage of the result of Appendix A for a reduction to the half-plane). Two features, however, persist. The first is an extreme instability $\gamma ^{t}\gg \gamma ^{p}$, centred at the point $(\alpha _1,\alpha _2)=(1,0)$ where the ion/electron responses are completely in phase that spreads with increasing $\bar {\varphi }_{\boldsymbol {q}}$. The second feature is the stabilisation $\gamma ^{t}<\gamma ^{p}$ of opposing responses with $\alpha _1<0$ so long as $\alpha _2$ is not too large.

Figure 11. Tertiary growth rate $\gamma ^{t}$ for $\eta =0.25$, $R/L_n=1.8$ and $p=0.6$ in the presence of a sinusoidal zonal profile of $q=0.4$ and amplitude (a) $\bar {\varphi }_{\boldsymbol {q}}=0.35$, (b) $\bar {\varphi }_{\boldsymbol {q}}=3.5$, (c) $\bar {\varphi }_{\boldsymbol {q}}=35$ and (d) $\bar {\varphi }_{\boldsymbol {q}}=350$, as a function of the relative phase $\alpha _1+\alpha _2\textrm {i}=\int \,\textrm {d}^{3}v_e g_{e\boldsymbol {q}} / \int \,\textrm {d}^{3}v_i g_{i\boldsymbol {q}}$ of the zonal ion/electron responses of (4.11a,b). Although the stabilisation changes significantly as the zonal amplitude is varied, two noteworthy features persist: maximum instability occurs for $(\alpha _1,\alpha _2)=(1,0)$ and the stabilising region extends out from $(\alpha _1,\alpha _2)=(-1,0)$, which is consistently (although not necessarily most) stabilising. Note that $\gamma ^{t}$-values larger than $2\gamma ^{p}$ are not shown and, as outlined in § 4.2, $\gamma ^{t}$ is unchanged when $\alpha _2\rightarrow -\alpha _2$.

Based on the result above, quiescent zonal profiles should be expected to exhibit similar opposing ion/electron responses. Indeed, the opposite sign of the ion and electron zonal temperature profiles of figure 6 already hinted at this. To firmly establish that this is the case, multiple simulations scattered closely around the point $\eta =0.25$ and $R/L_n=1.7$ were conducted to extract long-term stable zonal distributions. The result can be seen in figure 12, where the distribution of (a) the zonal shear contribution $2q^{2}|\bar {\varphi }_{\boldsymbol {q}}|$ and (b) $\mathrm {arg}(n_{ge\boldsymbol {q}}/n_{ge\boldsymbol {q}})$ for the first $15$ zonal modes are plotted. It is striking that $q=0.3$ and $q=0.4$-modes (whose shear contribution typically is the largest) indeed develop the expected phase shift $\mathrm {arg}(n_{ge\boldsymbol {q}}/n_{ge\boldsymbol {q}})\approx {\rm \pi}$ while the unstable range $\mathrm {arg}(n_{ge\boldsymbol {q}}/n_{ge\boldsymbol {q}})\approx 0$ constitutes an ‘avoidance zone’ they never attain. Confusingly, other modes are seen to exhibit similar ‘avoidance zones’ centred elsewhere or appear nearly uniformly distributed.

Figure 12. (a) Box-averaged zonal shear mode contribution $2q^{2}|\bar {\varphi }_{\boldsymbol {q}}|$ of different quiescent zonal profiles obtained from nonlinear simulations with $(\eta,R/L_n)=(0.25,1.7)$, with three coloured red, green and blue in particular. (b) The argument of the corresponding relative phase $n_{ge\boldsymbol {q}}/n_{gi\boldsymbol {q}}$. (c) The corresponding tertiary growth rates for $p=0.6$ of the modified profiles $\bar {\varphi }_{< q_\mathrm {cutoff}}$ and $\bar {\varphi }_{>q_\mathrm {cutoff}}$ that are obtained via the removal of modes not satisfying $q< q_\mathrm {cutoff}$ (triangles) and $q>q_\mathrm {cutoff}$ (squares). Modes with $q=0.3$ and $q=0.4$ typically develop a phase shift of ${\sim }{\rm \pi}$ and stabilise $\gamma ^{t}$ inordinately compared with their shearing. This can be seen by how $\gamma ^{t}$ greatly increases around $q_\mathrm {cutoff}\sim 0.3 - 0.4$ when large-scale modes are excluded. Similarly, when small-scale modes are excluded, $\gamma ^{t}$ instead decreases there.

The picture is clarified when we selectively filter out Fourier components to produce some high- and low-pass-filtered zonal fields $\bar {\varphi }_{>q_\mathrm {cutoff}}$ and $\bar {\varphi }_{< q_\mathrm {cutoff}}$, the tertiary growth rate $\gamma ^{t}$ of which are shown in figure 12(c). Here, $\bar {\varphi }_{>q_\mathrm {cutoff}}$ is seen to exhibit $\gamma ^{t}$ which rapidly increases from ${\sim }0$ to above $\gamma ^{p}$ as $q_\mathrm {cutoff}$ increases to exclude the $q=0.3$ and $q=0.4$-modes. A similar trend is observed for $\bar {\varphi }_{< q_\mathrm {cutoff}}$ when $q_\mathrm {cutoff}$ traverses in the opposite direction, although here, $\gamma ^{t}$ also displays large fluctuations in the range of $0.5\gamma ^{p} - 1.2\gamma ^{p}$. Combined, it therefore seems that the phase-shifted modes contribute much stronger to the total stabilisation than other modes, based on what one would suspect solely upon their shear. For the purposes of stabilisation, these modes are the most important component and can be thought of as a foundation to which modes at smaller-scale modes can contribute, either constructively or destructively.

4.4. Velocity space dependence

We noted in § 4.2 that $\gamma ^{t}$ is affected by velocity space structures. Therefore, we now want some understanding of how $\gamma ^{t}$ depends upon

(4.14)\begin{equation} \int\,{\rm d}^{3}v_sv_s^{2n}\bar{g}_{s\boldsymbol{q}}, \quad n\in\{1,2,3,\ldots\}, \end{equation}

for some representative $\bar {g}_{s\boldsymbol {q}}$-distributions. The reason why, as we noted in § 4.1, is that the relevant velocity moments for the tertiary theory entering (4.10) have the form

(4.15)\begin{equation} \int\,{\rm d}^{3}v_s\frac{v_s^{2n}\bar{g}_{s\boldsymbol{q}}}{L_s}, \quad n\in\{1,2,3,\ldots\}, \end{equation}

where the resonant denominator $L_s$ originates from $\boldsymbol{\mathsf{A}}_s^{-1}$, see (4.8a,b)–(4.10). In general, $L_s$ has a complex dependency on both $\lambda ^{t}$ and $v_s^{2}$, which is why we instead focus on the simpler moments (4.14) to obtain a rough assessment of how velocity space structures affect $\gamma ^{t}$.

As a representative case we initially lift the most stable distributions from figure 11(b,c), having the form (4.11a,b), and look for even more stable configurations by adding an additional part $\bar {g}_{s\boldsymbol {q}}^{T}=(a_{s1}+a_{s2}i)p_1(v_{s\bot }^{2})\exp (-v_s^{2})$, where $p_1$ is the second degree polynomial satisfying

(4.16a,b)\begin{align} \int\,{\rm d}^{3}v_sv_{s\bot}^{2}p_1(v_{s\bot}^{2}) {\rm e}^{{-}v_s^{2}}=1,\quad \int\,{\rm d}^{3}v_sp_1(v_{s\bot}^{2}) {\rm e}^{{-}v_s^{2}}=\int\,{\rm d}^{3}v_sv_{s\bot}^{4}p_1(v_{s\bot}^{2}) {\rm e}^{{-}v_s^{2}}=0, \end{align}

or $\bar {g}_{s\boldsymbol {q}}^{\chi }=(a_{s3}+a_{s4}\textrm {i})p_2(v_{s\bot }^{2})\exp (-v_s^{2})$, where $p_2$ is the second degree polynomial satisfying

(4.17a,b)\begin{align} \int\,{\rm d}^{3}v_sv_{s\bot}^{4}p_2(v_{s\bot}^{2}) {\rm e}^{{-}v_s^{2}}=1,\quad \int\,{\rm d}^{3}v_sp_2(v_{s\bot}^{2}) {\rm e}^{{-}v_s^{2}}=\int\,{\rm d}^{3}v_sv_{s\bot}^{2}p_2(v_{s\bot}^{2}){\rm e}^{{-}v_s^{2}}=0. \end{align}

These additions should be thought of as a ‘pure temperature’ perturbation and a ‘purely higher moment’ perturbation respectively, whose stabilisation effect we are interested in.

The resulting modification $\Delta \gamma ^{t}$ to $\gamma ^{t}$ when $\bar {g}_{s\boldsymbol {q}}$ is modified with the inclusion of $\bar {g}_{s\boldsymbol {q}}^{T}$ and $\bar {g}_{s\boldsymbol {q}}^{\chi }$, keeping $\bar {\varphi }_{\boldsymbol {q}}$ fixed using (4.12), can be seen in figures 13 and 14,respectively. It is very clear that the inclusion of either term could only marginally reduce $\gamma ^{t}$ by some $0.01\gamma ^{p}$. At greater amplitudes its inclusion even becomes strongly destabilising. Of course we cannot guarantee that it is impossible to employ some trial distribution with a different velocity space structure to significantly stabilise $\gamma ^{t}$. Nevertheless, a range of different exploratory investigations with different trial distributions hints that the velocity space stabilisation cannot move $\gamma ^{t}$ much at all below its minimum value among the configurations of figure 11. We will apply this evidence to simplify the critical gradient model developed later.

Figure 13. Tertiary growth rate modification $\Delta \gamma ^{t}$ when an additional part $\bar {g}_{s\boldsymbol {q}}^{T}=(a_{s1}+a_{s2}\textrm {i})p_1(v_{s\bot }^{2})\exp (-v_s^{2})$, with $p_1(v_{s\bot }^{2})$ being the second degree $v_{s\bot }^{2}$-polynomial satisfying (4.16a,b), is added to the initial zonal response of figure 11. In (a,b) $(\alpha _1,\alpha _2,\bar {\varphi }_{\boldsymbol {q}})=(-1,0,3.5)$, corresponding to the most stable point of figure 11(b), while in (c,d) $(\alpha _1,\alpha _2,\bar {\varphi }_{\boldsymbol {q}})=(0,0.5,35)$, corresponding to the most stable point of figure 11(c). The inclusion of the ‘pure temperature perturbation’ $\bar {g}_{s\boldsymbol {q}}^{T}$ is seen to predominantly destabilise ($\Delta \gamma ^{t}>0$) the tertiary instability, at most stabilising ($\Delta \gamma ^{t}<0$) some ${\sim } 2\,\%$ of the primary growth rate $\gamma ^{p}$.

Figure 14. The equivalent result to figure 13, but where instead $\bar {g}_{s\boldsymbol {q}}^{\chi }=(a_{s3}+a_{s4}\textrm {i})p_2(v_{s\bot }^{2})\exp (-v_s^{2})$, with $p_2(v_{s\bot }^{2})$ being the second degree $v_{s\bot }^{2}$-polynomial satisfying (4.17a,b), is added. Once again, the inclusion of a higher moment perturbation is seen to predominantly destabilise ($\Delta \gamma ^{t}>0$) the tertiary instability, with marginal stabilisation ($\Delta \gamma ^{t}<0$) not exceeding ${\sim }2\,\%$.

5.1. Adapting the prediction

The astute reader will already have identified the problem with (5.2) as written: it does not prescribe how one should choose $\bar {g}_{s\boldsymbol {q}}$, being formulated for a fluid model. This can, in principle, be fixed by extending the optimisation to find maximally stabilising distributions, but the vast possibility space makes this massively computationally inefficient. We must therefore reduce our attention to a minimal representative set of trial distributions

(5.5)\begin{equation} \bar{g}_{s\boldsymbol{q}}=h(\alpha_{s1},\alpha_{s2},\ldots,\alpha_{sj_s},v_{s\bot}^{2},v_{s\parallel}^{2}), \end{equation}

parametrised by some $j_e+j_i$ variable parameters $\alpha _{sl}$ with $l\in 1,2,\ldots,j_s$. To complete our system, these can, analogously to $\bar {\varphi }_{\boldsymbol {q}}$ and $q$ in (5.2)–(5.3), then be fixed as

(5.6)\begin{equation} \frac{\partial\gamma^{t}}{\partial\alpha_{sl}}=0, \end{equation}

which is added to (5.1)–(5.4) to find the most stable zonal mode within the space of trial distributions.

The question now is what trial distributions to substitute into the right-hand side of (5.5). Noting that Maxwellian velocity dependence affords a predictable control over density and temperature moments, we took our trial distributions to be those of (4.11a,b) with $j_e+j_i=4$. Naturally,this form is quite restrictive. Thankfully, § 4.4 strongly hints that, when minimising $\gamma ^{t}$, the inclusion of finer velocity space structures in our functional dependence can safely be neglected at an accuracy cost of perhaps only some $0.02\gamma ^{p}$. Their propensity to cause instability makes them, perhaps counter-intuitively, unimportant for the Dimits shift, since in the Dimits regime, where zonal flows evolve towards stability, they are unfavoured to emerge. This is massively favourable and signifies that our choice of trial distributions is sufficient.

As a further stroke of luck, exploratory investigations revealed that $\alpha _3\approx \alpha _4\approx 0.2$ consistently yielded the lowest tertiary growth rates for typical values $q\lesssim p$ and $\bar {\varphi }_{\boldsymbol {q}}\sim 1-10$. With minimal impact $\alpha _3=\alpha _4=0.2$ could thus be set, bypassing the need to dynamically solve (5.6) for these parameters. We could therefore restrict our attention from the full kinetic distribution possibility space to just a cross-species phase in the form of $\alpha _1$ and $\alpha _2$, which, however, massively affects the tertiary growth rate as § 4.3 and figure 11 highlights. In the end, our modified Dimits shift prediction is finally obtained as the solution $R/L_n$ to the equation system

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

(5.8)\begin{gather} \frac{\partial\gamma^{t}}{\partial \bar{\varphi}_{\boldsymbol{q}}}=0, \end{gather}

(5.9)\begin{gather} \frac{\partial \gamma^{t}}{\partial q}=0, \end{gather}

(5.10)\begin{gather} \frac{\partial\gamma^{t}}{\partial\alpha_{1}}=0, \end{gather}

(5.11)\begin{gather} \frac{\partial\gamma^{t}}{\partial\alpha_{2}}=0, \end{gather}

(5.12)\begin{gather} \gamma^{t}\left(\frac{R}{L_n}\right) = 0. \end{gather}

5.2. Prediction comparison

Finally, proceeding to put the prediction into action, the $\gamma ^{t}$-minimum of a given $R/L_n$-value, with poloidal wavenumber $p$ given by (5.7), was obtained through a steepest-descent walk of the parameters $\varphi _{\boldsymbol {q}}$, $q$, $\alpha _1$ and $\alpha _2$ that simultaneously solved (5.8)–(5.11).

To enhance our understanding of the zonal flows explored by this procedure we turn to figure 15. There, we can see the 32M $\gamma ^{t}$ as a function of $\bar {\varphi }_{\boldsymbol {q}}$ and $q$ with $\alpha _1$ and $\alpha _2$ determined by (5.10) and (5.11). Strikingly, as $q>p$ all distributions $\bar {g}_{s\boldsymbol {q}}$ satisfy $\gamma ^{t}>\gamma ^{p}$ except when $\bar {\varphi }_{\boldsymbol {q}}$ is small. This is the behaviour of KH-like modes (Zhu, Zhou & Dodin Reference Zhu, Zhou and Dodin2018), with the generalised KH instability criterion that the radial wavenumber must exceed the poloidal of gyrokinetics, $k_x^{2}>k_y^{2}$, already well known (Kim & Diamond Reference Kim and Diamond2002; Diamond et al. Reference Diamond, Itoh, Itoh and Hahm2005). Thus we can easily identify the unstable upper-right region of figure 15 as KH dominated, with an associated $k_xk_y|\bar {\varphi }|$-dependent growth rate.

Figure 15. Growth rate $\gamma ^{t}$, normalised to $\gamma ^{p}$, at the Dimits threshold $\eta =0.25$ and $R/L_n=1.8$ of the tertiary mode with $p=0.6$. Here, the sinusoidal zonal profile $\bar {\varphi }_{\boldsymbol {q}}$ with radial wavenumber $q$ has $\bar {g}_{s\boldsymbol {q}}$ correctly out of phase for maximum stability, as seen in figure 11. Strong instability commences only when $q>p$ and, unless $q$ is large, only for large zonal flows (although not shown, $\gamma ^{t}$ explodes above $2\gamma ^{p}$ in this region). Below this threshold the zonal flow is uniformly stabilising, to varying degrees. Note the peculiar feature that $\gamma ^{t}$ exhibits two separate local minima with respect to $\bar {\varphi }_{\boldsymbol {q}}$ for $q\gtrsim 1.3$.

In contrast to the KH range, stabilising distributions with $\gamma ^{t}<\gamma ^{p}$ exist when $q< p$ and in a limited region of small $\bar {\varphi }_{\boldsymbol {q}}$ when $q>p$. Initiating the aforementioned steepest-descent walk at a point with $q< p$, we eventually reach the most stable configuration at $q\approx 0.55$ and $\bar {\varphi }_{\boldsymbol {q}}\approx 5$ where $\gamma ^{t}\approx 0.01\gamma ^{p}$. Were we to employ the 4M-reduction to calculate $\gamma ^{t}$ instead of the full 32M of figure 15, this point would be moved upwards to $\bar {\varphi }_{\boldsymbol {q}}\approx 35$. Studying figure 10(b), this hints that the direct coupling to conjugate mirror modes of Pueschel et al. (Reference Pueschel, Li and Terry2021), captured by the 4M-description, may easily be spoiled by the inclusion of further sidebands. This view is further strengthened because nonlinearly $\bar {\varphi }_{\boldsymbol {q}}$ also typically reaches values where a larger $q_{{G}}$ is required for $\gamma ^{t}$ to converge.

Knowing how to produce a single, minimal $\gamma ^{t}(R/L_n)$-value, we arrive at our prediction by solving (5.12). Because of the presence of multiple undriven modes, the convergence of $\gamma ^{t}$ is exceptionally slow around the Dimits threshold. Therefore, the algorithm is commenced at $R/L_n$ well into the expected turbulent range where $\gamma ^{t}(R/L_n)$ and $\partial \gamma ^{t}/\partial (R/L_n)$ can more quickly be calculated. Then Newton's algorithm is deployed for $\gamma ^{t}(R/L_n)$ to find the solution to (5.12), ameliorating the problem of slow convergence.

How the result of this entire procedure, both 4M and 32M, compares with the result of nonlinear simulations, run for up to $t=8000$ and categorised through the measure $\varTheta _{0.1}$ of (3.2), can be seen in figure 16. The 4M-prediction is found to match the observed Dimits threshold remarkably well. However, as mentioned in § 5.1, the 4M solution of (5.7)–(5.12) produces high zonal amplitudes that would cause coupling to additional sidebands if these were included. The 32M-solution does not suffer from this problem, and also does not differ much from the observed Dimits threshold. Curiously, the 32M critical gradient of the latter is actually lower than the 4M one, a feature that will be discussed in § 6. In either case, the key result is that a reduced-mode tertiary instability prediction is able to accurately capture the Dimits transition.

Figure 16. Characterisation of the qualitative long-term behaviour of nonlinear simulations, with (possibly intermittent) unstable or continuously quiescent zonal flows, as a function of the system configuration $R/L_n$ and $\eta =L_n/L_T$. The estimated Dimits threshold, obtained as the solution to (5.7)–(5.12) using both a reduced 4M and a full 32M (i.e. $q_{G}=15q$) calculation, is also plotted, exhibiting remarkable agreement.