Predicting the Z-pinch Dimits shift through gyrokinetic tertiary instability analysis of the entropy mode

Abstract

The Dimits shift, an upshift in the onset of turbulence from the linear instability threshold, caused by self-generated zonal flows, can greatly enhance the performance of magnetic confinement plasma devices. Except in simple cases, using fluid approximations and model magnetic geometries, this phenomenon has proved difficult to understand and quantitatively predict. To bridge the large gap in complexity between simple models and realistic treatment in toroidal magnetic geometries (e.g. tokamaks or stellarators), the present work uses fully gyrokinetic simulations in a Z-pinch geometry to investigate the Dimits shift through the lens of tertiary instability analysis, which describes the emergence of drift waves from a zonally dominated state. Several features of the tertiary instability, previously observed in fluid models, are confirmed to remain. Most significantly, an efficient reduced-mode tertiary model, which previously proved successful in predicting the Dimits shift in a gyrofluid limit (Hallenbert & Plunk, J. Plasma Phys., vol. 87, issue 05, 2021, 905870508), is found to be accurate here, with only slight modifications to account for kinetic effects.

1. Introduction

Transport associated with small-scale plasma density/temperature-gradient-driven instabilities (Liewer 1985) has proven to be highly significant for nuclear fusion experiments, limiting performance by contributing significant, often stiff, transport (Horton 1999; Mantica et al. 2009; Ryter et al. 2011). Knowing the associated turbulent transport levels that dictate confinement properties is thus imperative, but an accurate prediction of these requires numerically expensive gyrokinetic simulations (Lin 1998; Dimits et al. 2000; Parker et al. 2004). Among other things, linear instability growth rates have therefore been used to expedite this analysis, forming the basis of simpler mixing length estimates (Waltz, Kerbel & Milovich 1994; Bourdelle et al. 2007; Pueschel et al. 2016).

However, near marginal linear stability where reactors may be expected to operate due to the aforementioned stiff transport, such a treatment is called into question by the existence of self-generated poloidal zonal flows (Diamond & Kim 1991; Lin 1998; Diamond et al. 2005). These efficiently suppress radial streamers by $E\times B$-flow shear decorrelation (Biglari, Diamond & Terry 1990; Lin 1998) or zonal-flow-catalysed energy transfer (Terry et al. 2018, 2021). This process can be so efficient that transport becomes strongly quenched, resulting in an effective upshift of the critical gradient, known as the Dimits shift (Dimits et al. 2000). Since it is well known that zonal shear quenching exhibits a complex dependence upon the magnetic geometry (Kinsey, Waltz & Candy 2007; Belli, Hammett & Dorland 2008), so should the size of the Dimits shift. This opens up the possibility of linearly unstable configurations that nevertheless are Dimits stable, which may prove valuable for reactor optimisation, particularly for stellarators with their large configuration space (Boozer 1998).

Unfortunately, a single general description of the Dimits upshift has proven elusive, with many competing explanations having been put forth (see e.g. Ivanov et al. 2020; Pueschel, Li & Terry 2021). This is because there exists a wealth of Dimits regime physics, even in the simpler fluid models usually under investigation. Although a full overview is beyond the present scope, a few examples are illustrative. Starting with non-local transport, Ivanov et al. (2020) studied ‘ferdinons’, which are radially propagating soliton-like coherent structures resembling those previously reported by van Wyk et al. (2016) in fully developed gyrokinetic turbulence. In a similar vein, Qi, Majda & Cerfon (2020) observed transient disordered radial avalanches similar to those of McMillan et al. (2009), which are known to coexist with coherent structures in fully developed turbulence (McMillan, Pringle & Teaca 2018). Despite providing a dominating transport contribution, it is unknown precisely how the two are related and what their effect upon the Dimits shift is.

Next, although the time-averaged transport by definition is small in the Dimits regime, this does not preclude the possibility of transient turbulent bursts (see Berionni & Gürcan 2011), of which the aforementioned avalanches are an example. Thus predator–prey-type energy oscillations between zonal flows, in the role of predator, and drift waves, as prey (Diamond et al. 1994), are prevalent. These oscillations are observed to vary qualitatively in e.g. amplitude, frequency or energy channels, even within a single system as the Dimits regime is traversed (Malkov, Diamond & Rosenbluth 2001; Berionni & Gürcan 2011; Kobayashi, Gürcan & Diamond 2015; Zhu, Zhou & Dodin 2020a). However, they eventually give way to continuous finite transport around the Dimits threshold, but it remains unknown to what extent a causal link exists.

Related to the previous features, it is well known that strong zonal flows cause drift waves and their associated turbulence to localise around those limited regions of zero $E\times B$ flow shear in the Dimits regime (Kobayashi & Rogers 2012; Kim, Min & An 2018, 2019; Zhang & Krasheninnikov 2020). In effect, this creates transport barriers with limited drift waves between these points, limiting overall transport. However, an associated nonlinear effect also seems to be that zonal flows eventually saturate into so called $E\times B$-staircases (Dif-Pradalier et al. 2010; Peeters et al. 2016; Garbet et al. 2021) that may extend well above the Dimits regime. This raises doubts that the breakdown of these barriers are sufficient to capture the Dimits shift. Nevertheless, although the formation and sustainability of these staircases are poorly understood, they are such an omnipresent feature that they surely must be accounted for.

Despite the complexity the examples above hint at, attempts continue to be made to encapsulate and predict the Dimits shift. Broadly, we can classify these into two categories. The first of these focuses on the turbulent dynamics above the Dimits threshold and attempts to extrapolate ‘downwards’ with the inclusion and emphasis of some particular piece of physics. The advantage of this approach is how it connects with already existing transport modelling. Unfortunately, efficient conventional quasilinear modelling fails in the Dimits regime, since mean zonal shear exceeds instability growth rates (Pueschel et al. 2016, 2021).

Thankfully, progress continues to be made on this front, with e.g. the combined work of Pueschel et al. (2021) and Terry et al. (2021) predicting a Dimits shift of the appropriate size. This was achieved by extending the turbulence saturation model of Terry et al. (2018) while including an observation of Hatch et al. (2016): nonlinear coupling is favoured between unstable modes and their conjugate side band mirror modes, even when these are not present linearly. The reason is that the three-wave correlation time is inversely proportional to the frequency mismatch, which is minimal for this coupling. Although this approach is promising, some piece of physics not captured by this scheme may still be vital. It remains to be seen if extrapolation candidates like this can prove generally apt.

The second family of Dimits shift descriptions, to which the present work belongs, instead focuses on the extreme quiescence very close to the linear instability threshold. Then the linearised behaviour of small-amplitude drift waves, in the presence of an essentially static zonal flow, is extrapolated ‘upwards’, i.e. towards the onset of turbulence. Of course, this may only prove truly justifiable in the collisionless regime where the zonal flow does not decay collisionally, but it has been observed that a significant collision rate is required to significantly alter the Dimits shift (Weikl et al. 2017). Unfortunately, it is not clear if the Dimits threshold coincides with a sharp linearised stability boundary. Should this not be the case, investigations of this second kind may nevertheless be useful to those of the first kind, clarifying properties of the Dimits regime.

As a result of these investigations, some success predicting the Dimits shift has been found in individual cases. However, these cases all involve some kind of simplified fluid model, predominantly of the Hasegawa–Mima–Wakatani kind (Hasegawa & Mima 1978; Hasegawa & Wakatani 1983). Frequently, although not always (see e.g. Ivanov et al. 2020), the analysis underlying these emphasises the so-called tertiary instability (St-Onge 2017; Zhu et al. 2020a; Zhu, Zhou & Dodin 2020b). The name of this instability arises from its place in the primary–secondary–tertiary hierarchy (Rogers, Dorland & Kotschenreuther 2000; Kim & Diamond 2002; Diamond et al. 2005). Here, the primary instability is the initial linear drift wave instability, the secondary instability amplifies drift wave inhomogeneities to then give rise to dominant zonal flows, which the tertiary instability finally affects to let drift waves re-emerge. Frequently, the tertiary instability, like the secondary instability, has been thought to essentially be a plasma nonlinear Kelvin–Helmholtz-type (KH) instability (see e.g. Kim & Diamond 2002; Numata, Ball & Dewar 2007; Kobayashi & Rogers 2012). Subsequent investigations have, however, revealed this assumption to typically be false under experimentally relevant conditions. Instead, the dominant tertiary instability is in fact a primary instability modified by zonal shear, and so it is vital that linear terms driving the primary instability are included for its treatment (Zhu et al. 2020a; Zhu & Dodin 2021).

Although eminently fruitful in elucidating key features, a fluid model cannot fully encapsulate the dynamics of gyrokinetics. Thus, a Dimits shift prediction obtained from these models must typically be modified and validated for the kinetic case, which may not be straightforward. Owing to the inherent complexity there has been scant work on this front, which the present article will attempt to begin rectifying. In previous work (Hallenbert & Plunk 2021), building upon that of St-Onge (2017), a reduced-mode tertiary analysis was employed to predict the Dimits shift of ion temperature-gradient (ITG) turbulence in a collisionless gyrofluid limit. Here, we proceed to extend that model to the kinetic case, avoiding the question of complex geometries by focusing on the simplest possible one: the Z-pinch.

Previous gyrokinetic simulations of Z-pinch entropy-mode-driven turbulence, while investigating and highlighting various features of near-marginal turbulence, have already established that the Z-pinch indeed exhibits a Dimits shift (Kobayashi, Rogers & Dorland 2010; Kobayashi & Gürcan 2015; Kobayashi et al. 2015). Here, we will again observe a similar dynamics in nonlinear simulations using GENE (Jenko et al. 2000), but also employ said code to perform a sophisticated kinetic tertiary instability analysis and arrive at a Dimits shift prediction that matches the nonlinear findings well.

This article is outlined as follows. The necessary fundamental physics of the subsequent analysis is provided in § 2, with the Z-pinch geometry in § 2.1, its collisionless gyrokinetics in § 2.2 and some relevant instabilities in §§ 2.3–2.4. A description of nonlinear simulations and some key results for this system then follow in § 3. Then, the tertiary instability is presented in § 4 and its numerical implementation discussed in § 4.1, with detailed analysis following in §§ 4.2–4.4. With this in hand, the reduced-mode Dimits prediction is presented, adapted for the present system, and finally tested in §§ 5–5.2, before a final brief summary and discussion in § 6.

2. Basic physics

The focus of the present article will be collisionless gyrokinetics in the Z-pinch geometry. This is advantageous for an exploratory investigation since the dominant electrostatic instability does not exhibit any parallel (along the magnetic field) component, i.e. with parallel wave number $k_\parallel =0$ (Simakov, Hastie & Catto 2002), which effectively reduces the system from spatially three-dimensional to two-dimensional. Furthermore, zonal flows are linearly conserved, simplifying and emphasising the tertiary picture of zonal flows. This is because of the lack of geodesic acoustic modes (Winsor, Johnson & Dawson 1968; Qiu, Chen & Zonga 2018) and their accompanying zonal/drift wave predator–prey-type intermittent oscillations (Miki et al. 2007; Kobayashi & Gürcan 2015), leaving the Rosenbluth–Hinton remnant fraction (Rosenbluth & Hinton 1998) identically one. Thus, this system is an excellent test bed for clarifying how key findings of simplified models (see e.g. Kolesnikov & Krommes 2005a; Numata et al. 2007; St-Onge 2017; Zhu et al. 2020a) are affected by the presence of kinetic effects.

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. 2006). 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 2004), 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 1978; Frieman 1982; Abel et al. 2013) for the ion-scale dynamics can in Fourier space and dimensionless form here be expressed as (see Plunk et al. 2014)

(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}}$.¹ $\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 1993) where the Dimits shift has conventionally been studied (Dimits et al. 2000). This is necessary because the dominant electrostatic instability at small gradients, originally considered by Kadomtsev (1960), is the so-called entropy mode. This mode, like the magnetohydrodynamic (MHD) interchange mode (Rosenbluth & Longmire 1957), 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 1962). 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 2000), 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 2012).

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. (2006) and Kobayashi & Rogers (2012) 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 2016) of GENE (Jenko et al. 2000), can be seen in figure 1, where (a) shows the range of gyrokinetic instability compared with the gyrofluid result of Ricci et al. (2006). 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 1957). 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 1956). 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. 2006).

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 1961) could be stabilised (Shumlak & Hartman 1995), which was also experimentally confirmed (Shumlak et al. 2001, 2009). 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. Nonlinear simulations

Before proceeding further to discuss the tertiary instability, it is worthwhile describing how the nonlinear simulations are performed and show some typical results. As previously mentioned, gyrokinetic simulations were performed using the GENE code (Jenko et al. 2000). Naturally, it is desirable to use a minimal numerical implementation for the multitude of runs needed to categorise a broad parameter range as belonging to the Dimits regime or not. Thus nonlinear simulations were performed with 32 points in $v_{\parallel }$-space and 12 points in $\mu$-space. These values were chosen to yield $\gamma ^{p}$-values within a margin of $0.01v_{\rm th}/R$ of the converged value close above the linear instability threshold and $0.0001v_{\rm th}/R$ around the Dimits threshold. As to the Fourier space decomposition, 128 $k_x$-modes and 16 $k_y$-modes with smallest wavenumber $k_x=k_y=0.1$ were found to be suitably minimal without compromising key results.

Although collisions were absent in the simulations, explicit numerical dissipation had to be included for stability. To this end, dynamically tuned gyrokinetic large eddy simulation (gyroLES) hyperviscous dissipation $\nu _x k_x^{4}+\nu _y k_y^{4}$ (Morel et al. 2011, 2012; Bañón Navarro et al. 2014) was employed to reduce the excitation of small-scale modes. A strength of this approach was its ability to capture of energy flows to smaller, unresolved scales in periods of turbulent flare up, while flexibly reducing numerical dissipation in the marginal Dimits regimes. Typical values for quiescent, zonally dominated states were $\nu _x\sim 10^{-4}$ and $\nu _y\sim 10^{-2}$, both of which flared up to ${\sim }0.05$ in periods of significant turbulence. Although schemes like this have already been noted as capable of allowing, and being consistent with, a Dimits shift (Morel et al. 2011), dissipation was also artificially made not to affect zonal modes. This eliminated the kind of predator–prey-type dynamics (see e.g. Kobayashi et al. 2015) that would complicate the tertiary picture.

Some example simulations of the kind just outlined can be observed in figure 2, which depicts the box-averaged particle flux²

(3.1)\begin{equation} \varGamma = \frac{1}{L_xL_y}\int_0^{L_x}\int_0^{L_y} \,{\rm d}x\,{\rm d}y \hat{\boldsymbol{x}} \boldsymbol{\cdot} v_{\boldsymbol{E}} \int \,{\rm d}^{3}v_s\delta f_s, \end{equation}

where $v_{\boldsymbol {E}}=\hat {\boldsymbol {z}}\times \boldsymbol {\nabla }\varphi$ is the $E\times B$-drift, and the zonal flow $\partial _x\bar {\varphi }$ over time. Three different configurations are shown: the first below the Dimits threshold, the second just above it and the last far above it. The dynamics is observed to vary significantly in a way consistent with previous findings (see e.g. Kobayashi & Gürcan 2015). When $R/L_n$ is small, a transient period of transport gives rise to an essentially stable zonal flow, which effectively quenches drift waves and their associated transport. On the other hand, when $R/L_n$ is large the zonal flow continually evolves and fails to quench turbulent transport. At intermediary values, the two types of dynamics coexist, with quasi-stationary zonal flows eventually succumbing in turbulent bursts before being re-established in a continually repeating pattern. Dubbed zonal flow cycling, this pattern was previously observed for the gyrofluid limit described in Hallenbert & Plunk (2021), differing significantly only in that the cycling here is much slower.

Figure 2. Time-average particle transport rate $\varGamma$, normalised to the gyro-Bohm level $\varGamma _{{\rm GB}}$, and the zonal flow $\partial _x\bar {\varphi }$ over time for $\eta =0.25$ and (a) $R/L_n=1.6$, (b) $R/L_n=2.0$ and (c) $R/L_n=2.5$. These values correspond to below, just above and well above the Dimits threshold. In case (a) the potential is virtually unchanged after it initially forms and the transport is very low, while in case (c) the potential is continuously modified and the transport is high. In the intermediate case (b) the system alternates between these two states over extended periods, resulting in bursty transport, whose mean magnitude is set by the non-quasistationary levels, and zonal profile cycling.

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 2010). 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 2017; Zhu et al. 2020b), 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. 2017), 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 1967; Qian 1986; Terry 2004). 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. 1994; Diamond et al. 2005). 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 1998; Kinsey, Waltz & Candy 2005). 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. (2020), acts to reinforce the zonal flow in the Dimits regime. Indeed, Kobayashi & Rogers (2012) 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. 1998; Kinsey et al. 2005). 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. 2007; Pueschel & Jenko 2010). 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. (1999), 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. (2021). 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. 2005; Kobayashi & Rogers 2012). 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. 2020). 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. 2011; Zhu et al. 2020a). 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. 2010; Peeters et al. 2016), 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 (2020) that localising well-defined staircases like these may provide effective transport barriers of drift waves and drifting coherent structures like ferdinons (Ivanov et al. 2020), 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 (2020) 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. The tertiary instability

Having described how nonlinear simulations were performed and some typical results, we now consider a static, high-amplitude zonal profile $\bar {\varphi }$ and proceed to extend the analysis of § 2.3 to capture its tertiary instability. Naturally, we can no longer neglect the nonlinear term, but (2.6) makes it clear that the $E\times B$ nonlinearity $\left \{ {\varPhi _s} , {g_s} \right \}_{\boldsymbol {k}}$ only couples modes together if they satisfy the triplet condition

(4.1)\begin{equation} \boldsymbol{k}=\boldsymbol{k}'+\boldsymbol{k}''. \end{equation}

Thus, zonal modes with $\boldsymbol {k}=(k_x,0)$ only couple modes sharing the same poloidal wavenumber, and for small-amplitude drift waves the nonlinear self-interaction between modes of different $k_y$ can be neglected. Because of this, the tertiary instability is characterised by a single poloidal wave vector $\boldsymbol {p}=(0,p)$. The corresponding discrete spectral decomposition into wave vectors

(4.2)\begin{equation} \boldsymbol{p}, \boldsymbol{p}\pm\boldsymbol{q}, \boldsymbol{p}\pm 2\boldsymbol{q} \cdots , \end{equation}

where $\boldsymbol {q}=(q,0)$ is purely radial, then gives rise to a combined tertiary mode

(4.3)\begin{equation} {\rm e}^{\gamma^{t} t+{\rm i}py}\sum_{s,n} a_{sn}{\rm e}^{{\rm i}lqx}g_{s\boldsymbol{p}+n\boldsymbol{q}}, \end{equation}

where $l$ runs over the integers and $a_{sl}$ has to be determined by insertion into the gyrokinetic equation.

Galerkin truncating (4.2) after some $\boldsymbol {q}_{G}=j\boldsymbol {q}$ we have a finite set of modes whose radial wavenumbers can be combined into a $(2j+1)$-dimensional vector

(4.4)\begin{equation} \boldsymbol{\mathsf{Q}} = \begin{bmatrix}-q_{G}\\ -q_{G}+q\\\cdots\\q_{G}-q\\q_{G}\end{bmatrix}, \end{equation}

whose elements we denote by $Q_l$, where $-j\leq l\leq j$. Then the zonal and non-zonal Fourier components of other variables can also be combined into corresponding vectors

(4.5)\begin{equation} \left. \begin{aligned} \bar{\boldsymbol{\mathsf{g}}}_s & =\begin{bmatrix}g_{s-\boldsymbol{q}_{G}}\\ g_{s-\boldsymbol{q}_{G}+\boldsymbol{q}}\\\cdots\\g_{s\boldsymbol{q}_{G}-\boldsymbol{q}}\\g_{s\boldsymbol{q}_{G}}\end{bmatrix}, \quad \tilde{\boldsymbol{\mathsf{g}}}_s=\begin{bmatrix}g_{s\boldsymbol{p}-\boldsymbol{q}_{G}}\\g_{s\boldsymbol{p}-\boldsymbol{q}_{G}+\boldsymbol{q}}\\\cdots\\g_{s\boldsymbol{p}+\boldsymbol{q}_{G}-\boldsymbol{q}}\\g_{s\boldsymbol{p}+\boldsymbol{q}_{G}}\end{bmatrix},\\ \boldsymbol{\bar{\varphi}} & =\begin{bmatrix}\varphi_{-\boldsymbol{q}_{G}}\\\varphi_{-\boldsymbol{q}_{G}+\boldsymbol{q}}\\\cdots\\\varphi_{\boldsymbol{q}_{G}-\boldsymbol{q}}\\\varphi_{\boldsymbol{q}_{G}}\end{bmatrix}, \quad \boldsymbol{\tilde{\varphi}}=\begin{bmatrix}\varphi_{\boldsymbol{p}-\boldsymbol{q}_{G}}\\\varphi_{\boldsymbol{p}-\boldsymbol{q}_{G}+\boldsymbol{q}}\\\cdots\\\varphi_{\boldsymbol{p}+\boldsymbol{q}_{G}-\boldsymbol{q}}\\\varphi_{\boldsymbol{p}+\boldsymbol{q}_{G}}\end{bmatrix}, \end{aligned} \right\} \end{equation}

with elements $\tilde {g}_{sl}, \bar {g}_{sl}, \tilde {\varphi }_{sl}$ and $\bar {\varphi }_{sl}$, for the potential and gyrocentre distributions. Letting $\boldsymbol{\mathsf{J}}_{0s}$, and $\boldsymbol {\varGamma }_{0s}$ denote the corresponding diagonal matrices, i.e. with $l$th entry $\textrm {J}_{0sll}=\textrm {J}_0(\sqrt {2(p^{2}+Q_l^{2})m_s/m_i}v_{s\bot })$ and $\varGamma _{0sll}=\textrm {I}_0((p^{2}+Q_l^{2})m_s/m_i) \exp ({-(p^{2}+Q_l^{2})m_s/m_i})$, respectively, we can thus write the set of coupled equations given by (2.4) and (2.8) for these modes as

(4.6)\begin{equation} \boldsymbol{\mathsf{A}}_s\tilde{\boldsymbol{\mathsf{g}}}_s = \boldsymbol{\mathsf{B}}_s\boldsymbol{\mathsf{J}}_{0s}\boldsymbol{\tilde{\varphi}}, \end{equation}

and

(4.7)\begin{equation} \sum_{s=i,e}Z_s\int \,{\rm d}^{3}v_s \boldsymbol{\mathsf{J}}_{0s} \tilde{\boldsymbol{\mathsf{g}}}_s = \sum_{s=i,e} n \left( \boldsymbol{\mathsf{I}} - \boldsymbol{\varGamma}_{0s} \right) \boldsymbol{\tilde{\varphi}}, \end{equation}

where $\boldsymbol{\mathsf{I}}$ is the identity matrix, the matrices $\boldsymbol {A}_s$ and $\boldsymbol {B}_s$ are given by

(4.8a,b)\begin{equation} \boldsymbol{A}_s=(\lambda+{\rm i}\omega_{ds\boldsymbol{p}})\boldsymbol{\mathsf{I}}+\boldsymbol{\mathsf{C}}_{1s}, \quad \boldsymbol{B}_s={\rm i}Z_sf_{0s}(\omega_{*s\boldsymbol{p}}-\omega_{ds\boldsymbol{p}})\boldsymbol{\mathsf{I}}+\boldsymbol{\mathsf{C}}_{2s}, \end{equation}

while the Poisson bracket is encapsulated by the matrices $\boldsymbol{\mathsf{C}}_{1s}$ and $\boldsymbol{\mathsf{C}}_{2s}$ with elements

(4.9a,b)\begin{equation} C_{1smn}={-}p\sum_{l={-}j}^{j}Q_l{\rm J}_{0sll}\bar{\varphi}_l\delta_{m(l+n)}, \quad C_{2smn}={-}p\sum_{l={-}j}^{j}Q_l\bar{g}_{sl}\delta_{m(l+n)}. \end{equation}

The procedure is now analogous to that of the primary instability. Thus we multiply (4.6) by $\boldsymbol{\mathsf{A}}^{-1}$ and substitute the result into (4.7) to find that, in order for the solution to possess non-zero $\boldsymbol {\tilde {\varphi }}$, the tertiary dispersion relation

(4.10)\begin{equation} \det\left[\sum_{s=i,e}\left(\frac{Z_s}{n}\int \,{\rm d}^{3}v_s\boldsymbol{\mathsf{J}}_{0s}\boldsymbol{\mathsf{A}}_s^{{-}1}\boldsymbol{\mathsf{B}}_s\boldsymbol{\mathsf{J}}_{0s}-(\boldsymbol{\mathsf{I}}-\boldsymbol{\varGamma}_{0s})\right)\right]=0 \end{equation}

must be satisfied. The exact solution can then be recovered in the continuum limit where $q_{G}\rightarrow \infty$, although we will only be concerned with the Galerkin-truncated system.

The similarity of (2.10) and (4.10) as presented here is deliberate: we consider the tertiary instability to be a linear confluence of primary modes within a poloidal band, coupled together by the zonal flow. Although this allows the description of a KH-like mode, as the tertiary mode has frequently been thought to be (see e.g. Kolesnikov & Krommes 2005b; Numata et al. 2007), it also allows modes of an entirely different kind. This is the modified primary instability, which differs in that the energy feeding the instability is supplied by the background gradients instead of the zonal flow, a feature elucidated and emphasised by Zhu et al. (2020a).

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 (2021) 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 (2012), 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. (2020a), Ivanov et al. (2020) or Hallenbert & Plunk (2021), 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 2021). 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 2006; Hatch et al. 2011), 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. A reduced-mode tertiary instability Dimits shift prediction

We now have all the tools in hand to proceed to what motivated the present work: an attempt to expand and apply the reduced-mode Dimits shift estimate outlined in Hallenbert & Plunk (2021), which proved successful in the gyrokinetic strongly driven fluid limit system. Expressed in the present notation, the predicted Dimits threshold was obtained from the critical solution of the equation system involving the four system parameters $p$, $q$, $\bar {\varphi }_{\boldsymbol {q}}$ and $R/L_n$ that is given by

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

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

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

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

Here, (5.4) is the ‘most important’, as its solution specifies a $R/L_n$-value that is taken to correspond to the Dimits threshold.

We will now go through (5.1)–(5.4) to make their meaning clear. Starting from the bottom, (5.4) succinctly expresses that the Dimits threshold corresponds to a tertiary stability threshold. Equation (5.1) next conveys the expectation that the tertiary mode of the poloidal band containing the most unstable primary mode also is the most unstable, and thus destabilises first. The effect of typical zonal flows is then approximated by that single maximally stable zonal mode which is set through (5.2) and (5.3), where the former fixes its amplitude and the latter its wavelength. This, as explained in Hallenbert & Plunk (2021), attempts to capture how the random nature of turbulence makes it plausible that the system will continue to explore the landscape of zonal flows until it finds a zonal profile that can possibly stabilise the system.

Rather than accepting the above description wholesale, it is prudent to review its applicability and validity. To this end we recall that, as indicated by figures 2, 8 and 9, $\gamma ^{t}$ is marginal for the quasistationary flows of the Dimits regime and increases substantially above it. Thus it seems the transition is indeed characterised by some condition like (5.4). Similarly, figure 8 makes it patently clear that 5.1 is a suitable and robust condition. That leaves just the reduction to a maximally stable single zonal mode of (5.2) and (5.3). But, here, figure 12 shows that a few modes give rise to most of the stabilisation, suggesting that a reduced-mode description should prove apt. The seeming evolution towards stable flows hinted at by figure 9, then finally affirms that its stabilisation should be maximised.

Before proceeding, it is important to note that in Hallenbert & Plunk (2021), the system (5.1)–(5.4) was solved with a minimal 4-mode³ (4M) $q_{G}=q$ Galerkin truncation with an eye towards computational efficiency for kinetic systems such as this one. Unfortunately, as will be described in § 5.1, there are some doubts about the present validity of such a severe reduction. Therefore, we will also make use of a less efficient, but more justifiable (see the convergence of figure 10), 32-mode (32M) reduction with $q_{G}=15q$) for the present prediction scheme.

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 2018), 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 2002; Diamond et al. 2005). 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. (2021), 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.

6. Discussion

In this article we have performed nonlinear flux tube simulations and employed extensive tertiary instability analysis for the entropy-mode-driven Z-pinch using GENE. This enabled us to extend the fluid model Dimits shift prediction of Hallenbert & Plunk (2021), which is essentially an efficient reduced zonal tertiary stability optimisation routine to mirror the observation that the zonal profile typically evolves towards more stable configurations. The inclusion of kinetic effects encouragingly alters the prediction minimally, and the result closely matches the observed Dimits threshold. Furthermore, the only apparent obstacle preventing this prediction from being adapted for toroidal geometries is a suitable means of accounting for the fact that only Rosenbluth–Hinton residual flows (Rosenbluth & Hinton 1998) can be stable. Hopefully,these can be captured by a similar trial distribution approach as in the present work, but more work on this front needs to be carried out.

Compared with the gyrokinetic fluid limit of Hallenbert & Plunk (2021), one key difference of Dimits regime dynamics is apparent. There, the Dimits range was characterised by large but transient turbulent bursts giving way to completely stable zonal profiles. Here, such stability can only occur at the very most stable part of the Dimits range. Instead, the Dimits range is characterised by a localised, low-amplitude drift wave with insufficient amplitude to significantly alter the zonal profile. Correspondingly, zonal profile evolution is much slower, and the need for ‘robustness’ lessened, see also § 4.1.

While large turbulent bursts do not occur within the Dimits regime of Z-pinch gyrokinetics as in the fluid limit, they emerge to dominate transport above the Dimits threshold but below higher $R/L_n$ where continuous turbulence eventually develops. These bursts seem to be of the same nature as those in Hallenbert & Plunk (2021), exhibiting little change in zonal energy, in contrast to conventional predator–prey-type zonal drift wave oscillations (Malkov et al. 2001; Kobayashi et al. 2015). Such oscillations arise as a result of zonal collisional damping. Higher collisionalities, although stabilising primary modes, more importantly also dampen zonal flows. This facilitates zonal/drift wave oscillations, which effectively induces higher transport and can even remove the Dimits shift entirely (Kobayashi & Rogers 2012; Weikl et al. 2017; Ivanov et al. 2020). Without such a mechanism here, as can be seen in figure 2(c), the present bursts instead coincide with periods of zonal profile cycling between profiles that are tertiary unstable at different points.

For computational efficiency, the present Dimits shift prediction can use a 4M reduction. This constitutes a simplification in which only immediate sideband coupling is considered that draws parallels to recent toroidal ITG work by Pueschel et al. (2021) and Terry et al. (2021), who found that energy transfer rapidly diminishes in the chain of sidebands. The reason why is that direct sideband coupling can engage stable conjugate sideband mirror modes, and by a nearly zero nonlinear frequency mismatch prolonging this three-wave correlation time, this is favoured (Hatch et al. 2016).

It is helpful to contrast the view above with the conventional shearing picture. Therein energy is continually shuffled, through shear, to smaller and smaller radial scale sidebands before being dissipated. Despite this energy transfer difference, both are consistent with the tertiary picture, with poloidal bands exhibiting collective, coherent behaviour. Nevertheless, Terry et al. (2021), based on findings of Whelan et al. (2019) that finite $\beta$-effects initially diminishes transport despite zonal flow shearing then becoming weakened, considers it likely that the dominant stabilising mechanism within the Dimits shift is direct mirror mode coupling instead of zonal shearing. For the collisionless Z-pinch, however, we have three pieces of evidence that indicate that this is not the case. First, as seen in figure 1, the Z-pinch possesses a very broad range of unstable sidebands even well inside the Dimits regime. Second, we noted that removing small-scale modes of an initially stable profile is sufficient to destabilise it. Finally, stabilisation via direct coupling in the Dimits range requires large-amplitude zonal modes, whose tertiary growth rates were observed in figure 10 to become unstable as more modes were included.

With this information in hand it is natural to question the validity of a 4M-prediction, which ostensibly captures direct coupling. Therefore, we completed the analogous 32M-prediction for comparison, which could also encapsulate zonal shear stabilisation. Since the maximally stabilising zonal profile found in each prediction scheme was different, it is perhaps surprising that this latter prediction was similar but slightly smaller than the previous, with both matching the observed Dimits threshold well. Although there is no guarantee that this is true in general, it does hint that the potential stabilising effect of zonal shear and conjugate mirror mode coupling may be comparable.

In §§ 4.2–4.4 and § 5.1 we saw that the choice of focusing on the Z-pinch, despite its simple geometry, was not uniformly advantageous for tertiary instability analysis. Despite the complete conservation of zonal flows, the necessary inclusion of kinetic electrons considerably complicates things. This vastly increases the zonal configuration space which complicates the optimisation scheme underlying the prediction. Although we attempted to capture this, we cannot guarantee that more stable zonal configurations can be found outside our limited selection of trial distributions. The fact that the Dimits threshold, with its emergence of unstable zonal flows, so closely coincides with the point at which we cease to discover stable profiles nevertheless hints that such configurations are very rare.

On the topic of kinetic distributions, it has already been observed by Li & Diamond (2018) that zonal flows need not have a uniformly stabilising effect even for configurations which are KH-stable. Thus, the stabilising influence of a zonal configuration can be altered significantly without altering the zonal flow profile itself. Nevertheless, the extent to which the full kinetic distribution affects tertiary (in)stability beyond its modification of the background gradients has perhaps not been adequately appreciated before in the literature. At best, it is possible, though not certain, that zonal self-organisation of the kind observed here lessens the importance of this effect by precluding such resonant instability as demonstrated in figure 11.

Regarding the tertiary instability, there have been many recent advances in its understanding. In accordance with these, we stress that the tertiary instability should primarily be thought of as a modified primary instability, differing from the former in how it extracts energy from background gradients instead of zonal flows (Zhu & Dodin 2021). Indeed, we expect the KH-instability affecting zonal flows not to play a significant role for the Z-pinch. As figure 15 indicates, it could, in principle, effectively limit zonal amplitudes. However, since it only truly arises at zonal amplitudes quite a bit greater than what is typical in simulations, it seems dubious to expect this to actually occur.

Now, although some attention has been given to features other than the tertiary instability in this article, naturally it has been the main focus. That is not to say that other nonlinear dynamics could not still be generally important for the Dimits shift. Indeed, much attention has recently been paid to non-diffusive transport in the form of avalanche bursts and solitary travelling structures (see e.g. Ivanov et al. 2020; Qi et al. 2020). Although no exhaustive search was carried out on this front, no immediately apparent signal indicating the presence of such features was observed. Allowing ourselves some liberty, we might speculate that their apparent absence is related to the fact that in the completely collisionless Z-pinch Dimits regime, drift wave amplitudes are too small for their self-interactions to enable manifestations of this kind. In conclusion, with the fuller understanding of the tertiary instability, observations generally point towards it being the dominant, sufficient contributor to the Dimits shift, at least for the collisionless Z-pinch.