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

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 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.


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 is difficult, requiring numerically expensive gyrokinetic simulations (Dimits et al. 2000;Lin 1998;Parker et al. 2004). Therefore, linear instability growth rates, among other things, have been used as a shortcut to expedite this analysis, forming the basis of simpler, more efficient mixing length estimates (Waltz et al. 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) that efficiently suppress radial streamers by E × B-flow shear decorrelation (Biglari et al. 1990;Lin 1998) or zonal-flow-catalysed energy transfer (Terry et al. 2018. In this range 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). Unfortunately, though the prediction of the Dimits shift would prove highly valuable, a general quantitative understanding for its exploitation has been lacking. This is because a precise, qualitative description of the upshift has proven elusive, with many competing explanations having been put forth (see e.g. Ivanov et al. 2020;Pueschel et al. 2021).
The reason why no single explanation has emerged victorious is a result of the veritable wealth of physics present within the Dimits regime, even for the simpler fluid models usually investigated. Though a full overview is beyond the present scope, a few examples are illustrative, the first of which could be categorised as nonlocal transport. Ivanov et al. (2020) observed the presence of radially propagating soliton-like coherent structures, calling them "ferdinons" due to their resemblance to those previously discovered by van Wyk et al. (2016) in fully developed gyrokinetic turbulence. In addition to these, Qi et al. (2020) also observed transient disordered radial avalanches similar to those of fully developed turbulence (McMillan et al. 2009), which are also known to coexist with such gyrokinetic coherent structures (McMillan et al. 2018). Precisely how these two are related to each other and their effect on the manifestation of a Dimits shift is unknown, but they seem to possess the same root cause while often providing a dominating transport contribution.
Next, although the time-averaged transport by definition is small in the Dimits regime, this does not exclude the possibility of transient turbulent bursts (see Berionni & Gürcan 2011), of which the aforementioned avalanches indeed is an example. Thus, a prevalent feature is predator-prey type energy oscillations between zonal flows, in the role of predator, and drift waves, as prey (Diamond et al. 1994). Now these oscillations are observed to vary qualitatively, both in amplitude, frequency, energy channels, or other ways. This holds even within a single system as the Dimits regime is traversed (Malkov et al. 2001;Berionni & Gürcan 2011;Zhu et al. 2020a). These ostensibly vanish when continuous finite transport arises around the Dimits threshold, but to what extent there exists a causal link remains unknown.
Related to the previous features, it is well-known that the strong zonal flows that characterise the Dimits regime cause drift waves and their associated turbulence to localise around those limited regions of zero E × B flow shear (Kobayashi & Rogers 2012;Kim et al. 2018Kim et al. , 2019Zhang & Krasheninnikov 2020). Naturally, the associated transport barriers with limited drift waves in between these points enhances the transport reduction within the Dimits regime. However, an associated nonlinear effect also seems to be that zonal flows eventually saturate into so called E × B-staircases (Dif-Pradalier et al. 2010;Peeters et al. 2016;Garbet et al. 2021) that may extend well above the Dimits regime, raising doubts that the breakdown of these barriers are sufficient to capture the Dimits shift. Nevertheless, though the formation and sustainability of these staircases are poorly understood, they are such a omnipresent feature that they surely must still be accounted for.
Despite the complexity the above examples 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 downward with the inclusion and emphasis of some particular piece of physics. Naturally, a clear advantage to this approach is how it connects with already existing transport modelling. Unfortunately fully nonlinear kinetic modelling is computationally expensive and it is probable that an accurate prediction of this kind will share this disadvantage. Unfortunately, in the Dimits regime of zonalflow regulated turbulence, where mean zonal shear exceeds instability growth rates, more efficient conventional quasilinear modelling fails . Thankfully, progress continues to be made on this front, with e.q. 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 the observation of , using gyrokinetic pseudospectral analysis, that nonlinear coupling is favoured between unstable modes and their conjugate side band mirror modes, even when these are not present linearly. This is because the three-wave correlation time is inversely proportional to the frequency mismatch, which is minimal for this coupling. Though this approach is promising, it may of course be that some piece of physics not captured by this scheme is in fact be of greater importance, and it remains to be seen which, of the multitude of potential extrapolation candidates, will emerge as most plausible.
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 and then proceeds from the linearised behaviour of small amplitude drift waves, in the presence of an essentially static zonal flow, to extrapolate "upward", i.e. toward the unset 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 the collision rate must be increased considerably to affect the Dimits shift significantly (Weikl et al. 2017). Additionally, however, it is not clear to what extent the more unstable Dimits range remains the same as the more stable range and if the Dimits transition can even be treated as a sharp qualitative transition. Even should this not be the case, investigations of this second kind may nevertheless still be of interest to clarify the Dimits regime, aiding those of the first kind.
Above all, the main advantage of the second approach is that the zonal flow may be treated as fixed, reducing the problem into an essentially linear one, which should be efficiently solvable. This could prove very useful for reactor optimisation, and in particular stellarators with their large space of possible configurations (Boozer 1998). Since it is well known that zonal shear quenching exhibits a complex magnetic geometry dependence (Kinsey et al. 2007;Belli et al. 2008), so should the size of the Dimits shift, opening the possibility of linearly unstable configurations that nevertheless are Dimits stable. However, a fruitful search for such configurations require a general Dimits shift prediction. Some success has been found predicting the Dimits shift for individual cases. However, owing to the complexity of a full gyrokinetic description, these cases all employ some kind of simplified fluid model, predominantly of the Hasegawa-Mima-Wakatani kind (Hasegawa & Mima 1978;Hasegawa & Wakatani 1983). Frequently, though not always (see e.g. Ivanov et al. (2020)), the analysis underlying these emphasises the so called the tertiary instability (St-Onge 2017; Zhu et al. 2020a,b). The name of this instability arises from its place in the primary-secondary-tertiary hierarchy (Kim & Diamond 2002;Diamond et al. 2005). Here the primary instability is a linear drift wave instability, then the secondary instability is that which magnifies slight drift wave inhomogeneities to give rise to zonal flows, which finally themselves are subject to the tertiary instability that causes the reemergence of drift waves. When originally studied, this instability, like the secondary instability, was thought to be essentially a plasma nonlinear Kelvin-Helmholtz-type (KH) instability . Subsequent investigations have however revealed this assumption to typically not be true under experimentally relevant conditions. Instead, the dominant tertiary instability is in fact a zonal shear modified primary instability, and so it is vital that linear terms driving the primary instability are included for its treatment (Zhu et al. 2020a;Zhu & Dodin 2021).
Though eminently fruitful in elucidating key features, a fluid model cannot encapsulate the full dynamics of full gyrokinetics. Thus a Dimits shift prediction obtained from these models must typically be modified for the kinetic case, which often is not straightforward, and then validated. 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 using the GS2-code, while investigating and highlighting various features of near-marginal turbulence, have already established that the Z-pinch indeed exhibits a Dimits shift (Kobayashi et al. 2010;. Here we will again observe similar dynamics in nonlinear simulations using the GENE code (Jenko et al. 2000), but also employ said code to perform sophisticated kinetic tertiary instability analysis and arrive at a Dimits shift prediction to match the nonlinear findings well. This paper is outlined as follows. The necessary fundamental physics of 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 follows in § 3, before the tertiary instability is presented in § 4 and its numerical implementation discussed in § 4.1, the detailed analysis of which is presented 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.

Basic Physics
The focus of the present paper will be collisionless gyrokinetics in the Z-pinch geometry. This has the considerable advantage for an exploratory investigation that, when the plasma β is low, the dominant electrostatic instability does not exhibit any parallel (along the magnetic field) component, i.e. k = 0 (Simakov et al. 2002), effectively reducing the system from being spatially 3D to just 2D. Additionally zonal flows are then linearly conserved, simplifying and emphasising the tertiary picture of zonal flows since geodesic acoustic modes (Winsor et al. 1968;Qiu et al. 2018) are not present with their accompanying zonal/drift wave predator-prey type intermittent oscillations (Miki et al. 2007;, and with the Rosenbluth-Hinton remnant fraction (Rosenbluth & Hinton 1998) being identically one. Thus this system is an excellent test bed when one wishes to extend observations obtained in simplified models (see e.g. Kolesnikov & Krommes 2005a;Numata et al. 2007;St-Onge 2017;Zhu et al. 2020a) to see how kinetic effects will modify the picture.

Z-pinch geometry
The set up for the Z-pinch can be expressed as follows in the usual cylindrical coordinates (r, θ, z) as follows: a strong current runs in the longitudinalẑ-direction which gives rise to a azimuthal magnetic field confining the plasma, assumed to have sufficiently low magnetic β so that its presence does not alter the longitudinally uniform magnetic vacuum field Like the magnetic field strength B, the large-scale plasma density n(r) and ion/electron temperatures T i (r) and T e (r), which we for simplicity henceforth will assume satisfy T i (r) = T e (r) = T (r) since a scaling argument reveals this to be the most unstable scenario (Ricci et al. 2006), only depends on the radial coordinate r, since the plasma is free to rapidly equilibrate along the magnetic field lines. In the typical gyrokinetic fashion we are then interested in the small scale fluctuations deriving their energy from the gradients of the large-scale background plasma. To encapsulate these it is advantageous to employ flux tube simulations (Candy et al. 2004), and thus we proceed to the local picture, expanding around the point r = R, setting B = B(R). Then it is advantageous to introduce the local gradient scale length where L n ∼ L r ρ i , and adopting a more suitable local flux tube coordinate system (x , y , z )x in which x ∈ (−L x /2, L x /2) is the radial coordinate, y ∈ (−L y /2, L y /2) is the bilinear "poloidal" coordinate and z is the parallel coordinate, where L x , L y R are the flux tube box widths. We will however be dropping the prime notation for future convenience and simply write (x, y, z) where it is implicitly understood that we refer to this curvilinear coordinate system and not Cartesian coordinates. 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.
In the absence of flux surfaces, since the magnetic field lines close on themselves after having encircled the central current axis once, the Z-pinch does not, strictly speaking, possess the familiar flux surfaces of toroidal geometries (Diamond et al. 2005) that constitute the "zones" in zonal flows. This is however a minor point, since these can easily be recovered via the inclusion of a small poloidal magnetic field δBŷ which twists the individually closed magnetic field line loops into a spiral outlining a "flux surface cylinder". Taking the limit δB → 0 the initial configuration is recovered while this flux surface nevertheless remains.

Gyrokinetics
The collisionless gyrokinetic equation (Catto 1978;Frieman 1982;Abel et al. 2013) of ion-scale dynamics for the present system can in Fourier space and dimensionless form be expressed as (see Plunk et al. 2014) for the component with wave number k. Here the macroscopic and microscopic spatial dimensions are normalised to R and the ion gyroradius ρ i = m i v T i / √ 2eB respectively, while the temporal dimension is normalised to the ion streaming time √ 2R/v T i . The subscript s = i/e denotes the ion/electron species with charge Z s e = ±e and macroscopic Maxwellian distribution f 0s (T ) with mean thermal velocity v T s = 2T /m s . The corresponding gyrocentre distribution is given by g sk = J 0sk δf sk R/ρ i , where δf sk is the fluctuating distribution, and the Bessel function of the first kind potential φ as ϕ = eφR/T ρ i . Naturally, the velocity and wavenumber are both split into their parallel and perpendicular components v , k , v ⊥ , k ⊥ with respect to the magnetic field.
Next in (2.4), the Fourier space Poisson bracket, representing the E × B-nonlinearity, is given by where δ k,k1+k2 is the Kronecker delta. Finally, employing a small magnetic β approximation so that ∇lnB = b · ∇b =x/R, the velocity-dependent diamagnetic and magnetic drift frequencies are given by where η = L n /L T and the local gradients L n and L T are given by (2.2).
To complement the gyrokinetic equation and provide closure, the gyrocentre distribution is related to the potential ϕ via the quasineutrality condition where Γ 0sk can be expressed in terms of the modified Bessel function I 0 through

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 (λ p k t)-time dependence, neglecting the nonlinear mode coupling through the Poisson bracket. The partial time derivative ∂ t is symbolically replaced by the primary complex frequency λ p k = γ p k − iω p k in the gyrokinetic equation (2.4), which is then multiplied by J 0sk . The resulting equation is then integrated over velocity space and inserted into the quasineutrality condition (2.8), resulting in the primary dispersion relation where the k-dependence enters through ω * s , ω ds , J 0sk , and Γ 0sk . Note the full inclusion of kinetic electrons in lieu of a simplified adiabatic response 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 presently necessary because the dominant electrostatic instability at small gradients, originally considered by Kadomtsev (1960), is the so called entropy mode. This mode, like the MHD interchange mode (Rosenbluth & Longmire 1957), leaves both the magnetic field and pressure p = nT unchanged, but modifies the entropy ln T 5/2 /p via exchange of density n and temperature T . Thus 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 electron dynamics are excluded (Ware 1962).
For the present purpose of investigating the Dimits shift, it is worth noting that the entropy mode is also occasionally referred to as the drift-temperature-gradient mode (Kesner 2000, see e.g.), despite existing also in the absence of a temperature gradient. This is because the drive term is the same as the ITG mode, and so it is unsurprising that its associated turbulence, subject to the same E × B-nonlinearity, could conceivably be ripe for the same shear-stabilisation characteristic of the Dimits shift, something already confirmed by Kobayashi & Rogers (2012). Turning our attention back to the dispersion relation (2.10), the fact that all linear terms in Equation (2.4) are proportional to k y makes it clear that the zonal potential ϕ = 1 L y Ly 0 ϕdy, (2.11) and thus the zonal E × B-flow ∂ x ϕ, are completely linearly stable, as previously mentioned, since it is composed of modes with k y = 0. Furthermore, as a consequence of us having assumed that T e /T i = 1, it is plain from (2.7) that the only two remaining independent parameters are R/L n and η, which span the configuration space.
As to the solution of (2.10), since it is an integral equation that cannot be solved analytically, it must therefore be treated numerically. Though one can, in fact, derive the equivalent entropy mode dispersion relation from a gyrofluid model, a kinetic treatment is necessary close to the marginal stability point of relevance to us, as both Ricci et al. (2006) and Kobayashi & Rogers (2012) note that the kinetic stability threshold is significantly lower than the collisionless gyrofluid value. Now, instead of explicitly solving (2.10), as will be apparent in § 4, it is convenient for us to instead 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 γ p by fitting an exponential to the system's long-term evolution.
The result of a calculation of this kind, using the code GENE (Jenko et al. 2000) with Z-pinch geometry (Bañón , can be seen in figure 1, where (a) shows the range of gyrokinetic instability thus obtained compared to 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 parameters, i.e. one that uniformly destabilises the plasma, unlike η or R/L T . In (b) and (c), on the other hand, the entropy mode growth rate is shown as a function of the wavenumber k for η = 0.25 with R/L n = 1.4 and R/L n = 1.8, which, as we will see in § 3.1, corresponds to 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 k ≈ (0, 0.6), while other modes constitute a broader, lower growth rate shoulder.

The interchange mode
So far we have just 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. Now, the first large scale MHD mode 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 et al. 1956) and possesses a very similar dispersion relation with a simple instability criterion given by (2.12) Figure 1. (a) Region of collisionless gyrofluid entropy mode instability (blue) and collisionless gyrokinetic entropy mode instability (green). The instability threshold of the latter is seen to occur at significantly lower gradients. In (b) and (c), the entropy mode (primary) growth rate γ p R/cs as a function of the radial and poloidal wavelengths kx, ky is shown for η = 0.25 and (b) R/Ln = 1.4 and (c) R/Ln = 1.8, the latter of which corresponds to the Dimits transition. With increasing R/Ln the unstable range is seen to widen significantly, but its maximum remains mostly stationary, peaked around ky = 0.6.
when electrons/ions have the same temperature (Ricci et al. 2006). For our purposes, the presence of the interchange mode at larger gradients, with its characteristic k-independent growth rate eliminates the possibility beyond this point of self-generated small-scale and small-amplitude zonal flow stabilisation. Simulations in this range are therefore pathological, 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 shearflow stabilised, for example it was predicted that the m = 1 kink mode (Newcomb & Kaufman 1961) could be stabilised (Shumlak & Hartman 1995), which was also confirmed experimentally (Shumlak et al. 2001(Shumlak et al. , 2009). These flows are however of a different type, being of significantly larger scales and externally imposed instead of self-generated. As such they affect the background equilibrium, and so we will not be considering the range given by (2.12).

Nonlinear simulations
Before proceeding further to discuss the tertiary instability it is suitable to describe how nonlinear simulations are performed, as well as some typical results. As previously mentioned, gyrokinetic simulations were performed using the GENE code (Jenko et al. 2000). Naturally, it is desirable to use as minimal a numerical implementation as possible to allow 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 space and 12 points in µ-space. These values were chosen in advance to yield γ p -values within a margin of 0.01v th /R of the converged value close above the linear instability threshold and and 0.0001v 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 was found to prove suitably minimal without compromising key results.
Though collisions were absent in the simulations, explicit numerical dissipation had to be included for stability. However, this dissipation was made not to affect zonal modes in order to eliminate the kind of predator-prey-type dynamics (see e.g.  which would complicate the tertiary picture. Naturally, with the inclusion of any zonal dissipation mechanism (collisional or not) this kind of behaviour would be recovered. Though the net effect on mean transport levels would be negligible for small values, with increasing values this would no longer be the case.
As to the explicit form of the numerical dissipation, dynamically tuned gyroLES hyperviscous dissipation of the form ν x k 4 x + ν y k 4 y (Morel et al. 2011(Morel et al. , 2012Bañón Navarro et al. 2014) was employed. This reduced the excitation of smaller scale modes and simulates the flow of energy to smaller, unresolved scales when turbulence flares up, while flexibly reducing numerical dissipation in the marginal Dimits regime and preventing it from inordinately affecting results. Indeed, this scheme has already been noted as capable of allowing, and being consistent with, a Dimits shift (Morel et al. 2011). Typical values for quiescent, zonally dominated states were ν x ∼ 10 −4 and ν y ∼ 10 −2 , both of which flared up to ∼ 0.05 in periods of significant turbulence. Some example simulations can be observed in figure 2, which depicts the box-averaged particle flux where v E =ẑ × ∇ϕ is the E × B-drift, as well as the zonal flow ∂ x ϕ over time for three different configurations: the first is below the Dimits threshold, the second far above it, and the last is only somewhat above. The dynamics are observed to vary significantly in a way consistent with previous findings (see e.g. . For small R/L n significant transport only exists transiently before an essentially stable zonal flow is established that effectively quenches drift waves and their associated transport. On the other hand, at large R/L n the zonal flow continually evolves and fails to quench turbulent transport. At intermediary values, the two types of dynamics coexist at different times, with quasi-stationary zonal flows eventually succumbing in turbulent bursts before being reestablished, a pattern which continually repeats. This is precisely the pattern, dubbed zonal flow cycling, observed for the gyrofluid limit described in Hallenbert & Plunk (2021), differing significantly only in that the cycling here is much slower.

The Dimits shift
It has been noted that tokamak turbulence around the Dimits threshold transport exhibit bursts that render time-resolved statistics frustratingly imprecise and possibly misleading (Pueschel & Jenko 2010), a feature that has continually plagued investigations of the Dimits shift, in particular when determining whether a given configuration belongs to the Dimits regime or not (for example, compare St-Onge 2017; Zhu et al. 2020b). The same holds true here, with statistics averaged over different, similar sections of a burst Figure 2. Time average particle transport rate Γ , normalised to the gyro-Bohm level ΓGB, and the zonal flow ∂xϕ over time for η = 0.25 and (a) R/Ln = 1.6, (b) R/Ln = 2.5, and (c) R/Ln = 2.0, corresponding to below, well above, and just 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 intermediary 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 state.
cycle yielding very different results. Nevertheless, keeping this in mind, employing large enough sample sizes, and focusing on features wherein this problem is lessened, a Dimits threshold can be discerned.
Thus we return to what was described in § 3, for which the time-averaged particle flux Γ as a function of R/L n for η = 0.25 can be seen in figure 3. Γ is seen to increase exponentially, but exhibits a discontinuity at R/L n ≈ 1.8 from a negligible level of ∼ 10 −2 Γ GB to the gyro-Bohm value ∼ Γ GB , which we interpret to correspond to the Dimits threshold.
This corresponds to the point where the fraction of time Θ the system exhibits turbulent bursts, instead of remaining quiescent with quasistatic zonal flows, increases from 0. To make this quantitative, we introduce Figure 3. In black: particle flux Γ as a function of the density gradient R/Ln for η = 0.25. In blue: the fraction Θ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 ± 0.3 and 4.9±0.4, discontinuously jumping between the two at R/Ln ∼ 1.8. This corresponds to that same point at which Θ0.1 becomes greater than 0 and will be identified as the Dimits threshold.
where H is the Heaviside function and t i is the time when the initial, linear growth phase seizes, i.e. Γ (t i ) is taken to be a representative level of turbulent transport. This quantity, for α = 0.1, is also plotted in figure 3. It is seen that this fraction begins to increases from 0 at the Dimits threshold. Though it only gradually increases thereafter, with fully developed turbulence (Θ α ∼ 1) only commencing well above the threshold, this clear signal was used to simply and unambiguously characterise the Dimits threshold. Now some comments on the observed exponential growth of Γ with respect to R/L n on both sides of the Dimits threshold are in order, since these appear to have differing slopes, 4 ± 0.3 and 4.9 ± 0.4 respectively. Below the Dimits threshold the zonal flow is large and almost stationary, and as we will describe in § 4, drift modes will, to lowest order, couple together quasilinearly into marginally (un)stable tertiary modes that only slowly exchange energy with the zonal flow and each other. These are very sensitive to any small perturbation of their dynamics, be it the zonal flow itself, sideband-sideband coupling, or dissipation mechanisms, even of high k x sidebands. That this is the case has been confirmed by running simulations where hyperviscosity, instead of dynamically changing, was kept constant at comparable values to those obtained through the gyroLES procedure, at which point Γ = 0 was eventually obtained for R/L n -values below 1.2−1.6, the precise value depending on ν x , ν y . One might be worried that the same holds true for the actual Dimits threshold at R/L n ≈ 1.8. This is indeed the case, but thankfully this is a much smaller effect, and Γ remains quite insensitive at and above the Dimits threshold. To achieve a comparable modification at these values, ν x , ν y must be increased much further, similar to previous collisional findings where zonal flows are also affected (Weikl et al. 2017).
Turning to the range above the Dimits threshold, the turbulent bursts appear very similar to the continuous turbulence at larger gradients. Unfortunately, this turbulence is not well-behaved, because an unbalanced inverse energy cascade inevitably leads to a pileup, as energy tries to shuffle energy to unresolved large scales, preventing a steady state from being obtained. This feature is well known to plague 2D turbulence (Kraichnan 1967;Qian 1986;Terry 2004). Thus we stress that, strictly speaking, the turbulent transport levels and their exponential increase above the Dimits threshold should not be Figure 4. Nonlinear particle flux as a function of the density gradient around the observed Dimits threshold ∼ 1.8 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 k-space resolution (red circles) (the latter two are offset for clarity), with an uncertainty calculated by splitting the full run into ten smaller pieces. Since the observed Dimits transition is seen to be consistent across configurations, the standard configuration employed is seen to be sufficient to determine the Dimits threshold.
trusted. Thankfully, the collisionless Dimits regime is not characterised by turbulence, but rather completely stable or very slowly "quasilinearly" evolving states, which are resolved. The Dimits threshold itself, which occurs when these states fail to remain uniformly quiescent, thus depends on this resolved range. This is confirmed in figure 4 which makes it clear that the discontinuous transport increase at the Dimits threshold remains persistent, and consistent, as the resolution is increased.

Zonal shear
A common quantitative measure of the stabilising effect of zonal flows for turbulence suppression is the box-averaged zonal shear magnitude |∂ 2 x ϕ| 2 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 |∂ 2 x ϕ| 2 1/2 ∼ γ p holds (Waltz et al. 1998;Kinsey et al. 2005). However, as can be seen in figure 5, the nonlinear shear significantly exceeds γ p , attaining a value of ∼ 1 in the Dimits regime. It is likely that this discrepancy can be traced to zonal dissipation, which would push the flow down to the 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 zonal shear levels above the quench rule in the Z-pinch Dimits regime that decreased with increasing collisionality.
Another heuristic estimate, which forms the basis of simple diffusivity scalings of the kind 1 − α |∂ 2 x ϕ| 2 1/2 /γ p with some α that are sometimes employed to model zonal shear quenching (Waltz et al. 1998;Kinsey et al. 2005), is that only when |∂ 2 x ϕ| 2 1/2 γ 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 |∂ 2 x ϕ| 2 1/2 increases faster than γ p above the Dimits threshold, even when accounting for the fact that the effective shearing decorrelation rate decreases when the zonal flow varies rapidly in time, like described by Hahm et al. (1999). A similar increase in ITG simulations was noted by Terry et al. (2021) as being inconsistent with the tertiary picture of the Dimits shift, seemingly under the assumption that greater zonal shear should be more x ϕ| 2 1/2 , obtained from nonlinear simulations after the transient period, as a function of the density gradient R/Ln for η = 0.25 (blue) compared to the entropy mode growth rate γ p (green), which is significantly lower. Also shown is the range where a sinusoidal zonal profile with kx = 0.4, as shown in § 4 to be approximately maximally stabilising, may be stable (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 beginning to gently increase thereafter.
stabilising. However, that this need not be the case has already been observed (see e.g. Kinsey et al. 2005;Kobayashi & Rogers 2012), and by also plotting the approximate range zonal shear amplitudes (determined as outlined in § 4 and § 5) where a sinusoidal profile can be completely stable in figure 5, it is seen to eventually be reduced to a single value, both from above and below, close to the Dimits threshold.

Localisation
A persistent feature of marginal drift waves in the presence of a zonal flow is the localisation of said modes to regions of zero zonal shear, ∂ 2 x ϕ = 0. This is in accordance with the intuitive picture that the zonal flow at these points is unable to decorrelate radial streamers, effectively only transporting them poloidally, i.e. in the y-direction, in a uniform way (Ivanov et al. 2020). However, this necessary condition is not sufficient for drift waves to congregate at such a point, since they asymmetrically respond to minima ∂ 3 x ϕ > 0 and maxima ∂ 3 x ϕ < 0: significant localisation is only observed at minima (McMillan et al. 2011;Zhu et al. 2020a). That the same holds here is very clear in figure  6, where it is additionally observed that more developed staircase states, i.e. broader rectangular shapes (Dif-Pradalier et al. 2010;Peeters et al. 2016) of nearly constant shear interspersed with sharp transitions from positive to negative values, arises as one increases R/L n to approach the Dimits threshold from below.
In a modified Hasegawa-Mima system, the formation of well-defined staircase states like these was recently conjectured by Zhang & Krasheninnikov (2020) to provide effective transport barriers of drift waves, and in particular drifting coherent structures like ferdinons that carry transport (Ivanov et al. 2020), reflecting and trapping them in the vicinity of zonal extrema. This effect was even boldly speculated to be more important than zonal shearing for transport quenching. Though possibly relevant, particularly close to marginal stability, this conclusion seems dubious for the system at hand. This is because drifting structures were not observed. Instead drift waves, though indeed Figure 6. Mean drift wave density perturbationñ as a function of the radial coordinate x in the presence of the quasistationary zonal flow ∂xϕ and zonal fluctuating temperature T for: (a) R/Ln = 1.2, (a) R/Ln = 1.5, and (c) R/Ln = 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 ∂ 2 x ϕ. Meanwhile, the drift waves increasingly localise around zonal flow minima, i.e. where the zonal shear ∂ 2 x ϕ vanishes and ∂ 3 x ϕ > 0.
localised, commonly remained at very low amplitudes with essentially no self-interactions for extended periods of time, seemingly unable to grow in the presence of strong shear. As a final note before we proceed, it is apparent from figure 6 that ion/electron temperature moments predominantly self-organise in such a way that they are out of phase with each other, something we will clarify further in § 4.3-4.4.

The tertiary instability
Having described how nonlinear simulations were performed and some typical results, we now consider a static, high amplitude zonal profile ϕ and proceed like in § 2.3 to study the tertiary instability in greater detail, since it served as the motivation for the present study. Naturally we can no longer neglect the nonlinear term, but, as is plain from (2.6), the E × B nonlinearity {Φ s , g s } k only couples modes together if they satisfy the triplet condition k = k + k . (4.1) Thus, zonal modes with 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 problem is characterised by a single poloidal wavenumber, and the discrete spectral decomposition that the vector wavenumbers given by p, p ± q, p ± 2q... , (4.2) where p = (0, p) is purely poloidal and q = (q, 0) purely radial, should be considered to constitute a combined tertiary mode e γ t t+ipy s,n a sn e inqx g sp+nq (4.3) where n runs over the integers and a n have to be determined by insertion into the gyrokinetic equation. Galerkin-truncating (4.2) after some q G = jq we thus have a finite set of modes whose radial wavenumbers can be combined into a (2j + 1)-dimensional vector (4.4) whose elements we denote by Q m , where −j m j. Then the Fourier components of other variables can also be combined into the corresponding vectors (4.5) with elementsg sm , g sm ,φ sm , and ϕ sm , for the potential and gyrocentre distributions. Letting J 0s , and Γ 0s denote the corresponding diagonal matrices, i.e. with mth entries J 0smm = J 0 ( 2(p 2 + Q 2 m )m s /m i w s ) and Γ 0smm = I 0 ((p 2 +Q 2 m )m s /m i )e −(p 2 +Q 2 m )ms/mi , we can thus write the set of coupled equations given by (2.4)  where I is the identity matrix and the matrices A s and B s are given by while the Poisson bracket is encapsulated by the matrices C 1s and C 2s with elements The procedure is now analogous to that of the primary instability, so upon multiplying (4.6) by A −1 and substituting the result into (4.7) we find that, in order for the solution to possess nonzeroφ, the tertiary dispersion relation must be satisfied. The exact solution can then recovered in the continuum limit where q G → ∞, though 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 precisely be a linear confluence of the poloidal band of primary modes coupled together by the zonal flow. Though this allows the description of a KHlike mode, like the tertiary mode is frequently 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).

Tertiary instability simulations
Of course, in general (4.10), when expanded out, becomes an intractable sum of products of integrals like (2.10) but with ever more complex integrands. Since (2.10) already must be solved numerically, the same holds true here, but as q G is increased it becomes cumbersome to do so directly. Instead the alternate, aforementioned procedure of evolving the gyrokinetic equation directly is employed, in a procedure that is essentially as described in § 2.3.
A zonal profile, i.e. the set of 2n zonal distributions g sk with potential modes ϕ 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 an n k x -mode and 2 k y -mode (k y = 0 and k y = p) Fourier decomposition 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. By the triplet condition (4.1) drift wave self-interaction are then not present, and having frozen the zonal profile the drift wave "back reaction" is also removed. Combined, this enables the convergence of the observed tertiary growth rate without drift waves nonlinearly affecting the zonal profile or each other as they grow in amplitude Before proceeding it is worth discussing the matter of numerical damping and resolution for this procedure, since this was found to be important. Performing tertiary simulations of stable profiles extracted from low R/L n -simulations, these were observed to only be marginally stable. Changing ν x by a small amount, though it might be expected to have little effect, therefore nevertheless often causes the profile to become unstable. Thus it seems that nonlinearly obtained profiles are "balanced" for tertiary sideband "resonant" coupling that is finely modulated by ν x . This might make one nervous that this numerical dissipation has a nonphysically important role when establishing stability. However, this is not the case. Modifying ν x to destabilise an initially stable profile, after restarting a nonlinear simulation it will after some transient period of increased drift waves quickly be only slightly modified to restore tertiary stability. Therefore, the specific realisation of ν x , unless massively increased, does not matter for the formation of stable profiles. The situation is entirely analogous for the numerical resolution, where the removal or inclusion of further modes can destabilise specific instances of zonal profiles, Figure 7. Tertiary growth rates of different, randomly chosen zonal profiles, obtained from nonlinear simulations with η = 0.25 and R/Ln = 1.8, as a function of the poloidal wave number 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 γ p is plotted in green for comparison. As can be seen, the quiescent profiles are only unstable to a few p-values clustered around ∼0.6, the most primary unstable point, with growth rates severely reduced to about 1% of the primary growth rates. In comparison, the turbulent profiles are much more unstable to a broader range of p-values.
but not change whether stability is expected to nonlinearly arise so long as the resolution is not taken to be patently too low.
To begin investigating the present tertiary instability using the aforementioned scheme, in figure 7 the tertiary growth rate γ t of the poloidal band p under the influence of different zonal profiles extracted from nonlinear simulations are are compared with the primary growth rate γ p . The profiles are obtained at the Dimits threshold η = 0.25 and R/L n = 1.8, where turbulent bursts first manifest. It is seen that γ t is much greater for zonal profiles of turbulent periods compared with their quasistationary counterparts that are only barely unstable with γ t ∼ 0.01γ p . However in both cases the growth rate is seen to peak around the same value p ≈ 0.6 where γ p attains its maximum.
Having observed that the 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 while decaying, markedly altering them. 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. The need for such stringent stability is presently much lesser. 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.
In spite of what was just discussed, investigating the question of "stability robustness" is still of interest. In figure 8, therefore, the maximal γ t is plotted for various zonal flow profiles, versus the rescaling factor a used to (artificially) adjust the zonal flow amplitude. Clearly, γ t increases relatively swiftly away from marginal stability as a is moved away from 1 for quasistationary profiles. On the other hand, the turbulent profiles can be made somewhat more stable for some values a = 1, but never to the same extent as the quasistationary ones. We interpret this as meaning that quasistationary flows, in Figure 8. Tertiary growth rate γ t , normalised to the primary growth rate γ p , of different zonal flow profiles obtained from nonlinear simulations 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.
being established when turbulent bursts end, are delicately pushed towards configurations of greater tertiary stability, all the way to marginal stability if possible. We also note again that increasing the zonal shear by increasing a does not necessarily imply that γ t decreases, a fact 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.

Tertiary instability of a single zonal mode
To gain some basic understanding of the tertiary instability, we now turn to the case of a single zonal mode i.e. a sinusoidal zonal profile, and in particular the 4M-reduced system obtained when q G = q in (4.4), which is the minimal system that can encapsulate zonal shear stabilisation. As is apparent from the presence of g s in C 2 the zonal flow ∂ x ϕ alone is insufficient to determine tertiary instability; the specific form of the zonal electron/ion distribution giving rise any particular zonal flow also matters by modifying the sideband coupling. Naturally, one reason for this is because g s in effect act as a modification to the background gradients. However, its influence is not restricted to this effect since the full kinetic zonal distribution g sq has to be accounted for, not just its density and temperature moments. 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 g sq , holding ϕ q fixed by employing a multiplicative factor.
Before delving into this topic, we will first discuss a way in which g s does not affect the tertiary instability. In Appendix A we show that the tertiary growth rate of a single zonal mode is unchanged under the transformation (4.11) Figure 9. Tertiary growth rate γ t at η = 0.25 and R/Ln = 1.8 for p = 0.6, normalised to the primary growth rate γ p , of a single zonal mode with wavenumber q = 0.4 for different example configurations of the form (4.12) for (a) ϕ q =3.5, α1 + α2i = (0.5, 2 + i, 0, −1), and (b) ϕ q =35, α1 + α2i = (1 + i, 0.4 + i, 0.1 + 0.25i, −1 + i). Also plotted in green is the primary growth rate γ p p+q G of the sideband mode with kx = qG. The configurations are chosen to give a wide spread of converged γ t , which occurs around qG/q ∼ 10 and qG/q ∼ 30 for ϕ q = 3.5 and ϕ q = 35 respectively, beyond which point further coupling to smaller scale sidebands becomes negligible. As is shown, the growth rate decreases as smaller scale sidebands are engaged, but not below γ p p+q G . For greater zonal amplitudes, however, both patterns are disrupted for small qG.
Since this substitution also leaves ϕ unchanged, the tertiary instability can only possess a g sq -dependence of the (species-summed) form |g sq | and g sq + g * sq ϕ 2 q /|ϕ q | 2 , i.e. Re (g sq ϕ q * ). As a consequence, it does not matter which species distribution trails the potential and which leads it. This result can be interpreted in light of the tertiary localisation highlighted in § 3.3. The part of the particle distribution g sout that is out of phase with ϕ satisfies ∂ x g sout = 0 at the points ∂ 2 x ϕ = 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 g sq in § 4.3 for the tertiary instability, it is prudent to investigate its effect. Indeed, the single zonal mode presently under consideration is suitable for this purpose. Thus, how γ t , precisely obtained with GENE as outlined in § 4.1, changes as q G is varied for some different single-mode configurations with q = 0.4, p = 0.6, and ϕ q = 3.5 (a typical nonlinear value) or ϕ q = 35, is compared with the primary sideband growth rate γ p p+q G in figure 9. For small to "normal" zonal amplitudes γ t > γ p p+q G seems to uniformly hold in accordance with the view that the tertiary instability mixes primary modes (Hallenbert & Plunk 2021). When significantly increased, γ t > γ p p+q G may however hold for small q G , a feature presumably explained by immediate subdominant stable mode coupling (see e.g Terry et al. 2006;Hatch et al. 2011), the matching of which becomes less resonant when more modes are included, i.e. when q G is increased.
Continuing, we note that the point at which γ t converges effectively constitutes a measure of how many sidebands are coupled together into a combined tertiary mode, which intuitively increases with increasing ϕ q , since in the opposite weak coupling limit ϕ q → 0 different wavenumbers decouple. For the typical zonal amplitude of figure 9 (a) this point occurs at k x ≈ 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.

Role of electron/ion response phase
Leaving aside the complicated question of how the inclusion of multiple zonal modes with their accompanying cross-phases affect the tertiary instability, we proceed to investigate how important the different zonal distribution effects uncovered in § 4.2 are, the first of which is the relative phase of the ion/electron response. To investigate this feature, both distributions are taken to possess a simple representative Maxwellian f 0sdependence, i.e.
where M is a multiplicative factor that ensures ϕ q attains the desired value ϕ qdes , and which can be calculated using quasineutrality (2.8) as Here is the cross-species phase whose influence we want to investigate. On the other hand, α 3 = α 4 = 0.2 are "effective temperatures", the choice of which will be explained in § 5. Figure 10 shows how γ t of the poloidal band p = 0.6 at the Dimits threshold η = 0.25 and R/L n = 1.8 changes as α 3 /α 4 are varied, Four different cases, corresponding to ϕ 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 γ t γ p , centred at the point (α 1 , α 2 ) = (1, 0) where the ion/electron responses are completely in phase and which spreads with increasing ϕ q . The second feature is the stabilisation γ t < γ p of out-of-phase responses with α 1 < 0 so long as α 2 is not too large.
Based on the above result and the findings of § 4, we conclude that quiescent zonal profiles should exhibit this ion/electron phase shift, which the zonal temperature profiles of figure 6 already hinted at. To firmly establish that this is the case, multiple simulations scattered closely around the point η = 0.25 and R/L n = 1.7 were conducted to extract different long-term stable zonal distributions. The result of these can be seen in figure  11, where the distribution of (a) the contribution 2q 2 |ϕ q | to the zonal shear and (b) arg(n geq /n geq ) for the first 15 zonal modes is plotted. It is striking that the q = 0.3 (whose shear contribution typically is the largest) and q = 0.4-modes indeed develops the expected phase shift arg(n geq /n geq ) ≈ π, but confusingly other modes are seen to exhibit similar "avoidance zones" centred elsewhere or appear nearly uniformly distributed.
The above picture is clarified when we selectively filter out Fourier components to define high-and low-pass-filtered zonal fields ϕ >q cutoff and ϕ <q cutoff . The tertiary modes are then simulated with GENE and the resulting growth rate γ t of three such profiles shown in 11 (c). γ t of the low-pass-filtered fields is seen to rapidly increase from ∼ 0 to greater than γ p where it then remains as q cutoff goes from 0.3 to 0.5, i.e. when the aforementioned q = 0.3 and q = 0.4-modes are excluded. For the high-pass-filtered fields, a similar trend is observed when q cutoff traverses in the opposite direction, though here Figure 10. Tertiary instability growth γ t for η = 0.25, R/Ln = 1.8, and p = 0.6 in the presence of a sinusoidal zonal profile of q = 0.4 and amplitude (a) ϕ q = 0.35, (b) ϕ q = 3.5, (c) ϕ q = 35, and (d) ϕ q = 350, as a function of the relative phase α1 + α2i = d 3 vgeq/ d 3 vgiq of the zonal electron/ion responses of (4.12). Though the stabilisation changes significantly as the zonal amplitude is varied, two noteworthy features persist: maximum instability occurs for (α1, α2) = (1, 0); the case (α1, α2) = (−1, 0), though not necessarily most stabilising, is always very stabilising. γ t also exhibits fluctuating large values also in the range 0.5γ p − 1.2γ p . Combined, it therefore seems that the phase-shifted modes contributes especially strongly to the total stabilisation compared with other modes, as one would suspect based solely on their shear. These modes can be thought of as acting as a foundation to which modes at smaller scale modes can contribute, constructively or destructively, but without which smaller scale modes fail to even approach stability; in short, they are the most important component.

Velocity space dependence
Another dependence of the tertiary we noted in § 4.2 was that the dependence upon velocity space moments and so we wish to explore how γ t depends upon d 3 vv 2n g sq , n ∈ {1, 2, 3, ...}, (4.15) for some representative g sq -functions. Naturally, we know from § 4.2 that the relevant velocity moments for the tertiary theory have the form (4.16) Figure 11. (a) The individual mode contribution 2q 2 |ϕ q | to the squared box-averaged zonal shear |∂ 2 x ϕ| 2 of different zonal quiescent zonal profiles obtained from nonlinear simulations with (η, R/Ln) closely clustered around (0.25, 1.7), three of which (red, blue, green) are highlighted. (b) the argument of the corresponding relative phase ngeq/ngiq. (c) the corresponding tertiary growth rates for p = 0.6 of the modified profiles ϕ < and ϕ > that are obtained via the removal of modes not satisfying q < q cutoff (triangles) and q > q cutoff (squares) respectively. It is apparent that the modes with q = 0.3 and q = 0.4 typically develops a phase shift of ∼ π and thus play an inordinately important role (compared to their shearing) for the profile's total stabilising ability, since γ t greatly increases around q cutoff ∼ 0.3 − 0.4 when large scale modes are excluded. When instead small scale modes are excluded, γ t instead decreases there, but continues to vary significantly and grow large again beyond this point, showing that smaller scale modes still can play an important role.
where the resonant denominator L s originates from A −1 s . Unfortunately, L s will generally depend upon λ t and v 2 , which necessitates the solution of an intractable integral equation for a rigorous determination of the velocity moments controlling the tertiary mode. We therefore focus on the general moments (4.15) to obtain a rough assessment of the importance of velocity space dependence.
Proceeding by lifting the most stable distributions of the form (4.12) in figure 10 (b) and (c), we then add an additional part g T sq = p 1 (v 2 ⊥ ) exp −v 2 , where p 1 is the third degree polynomial satisfying where p 2 is the third degree polynomial satisfying 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. Now, the resulting degree ∆γ t that the tertiary growth rate γ t is modified by the Figure 12. Tertiary growth rate modification ∆γ t when an additional part g T sq = p1(v 2 ⊥ ) exp −v 2 , with p1(v 2 ⊥ ) being the second degree v 2 ⊥ -polynomial such that g T sq has vanishing zeroth and fourth velocity moments but non-vanishing second velocity moment, is added to the initial, pre-renormalisation zonal response g n sq employed in figure 10. Specifically, for (a) and (b) (α1, α2, ϕ q ) = (−1, 0, 3.5), corresponding to the most stable point of figure 10 (b) where ϕ q = 3.5, while for (c) and (d) (α1, α2, ϕ q ) = (0, 0.5, 35), corresponding to the most stable point of figure 10 (c) where ϕ q = 35. The inclusion of the "pure temperature perturbation" g T sq is seen to predominantly destabilise the tertiary instability, and when it is stabilising its effect does not exceed ∼ 1 − 2% of the primary growth rate γ p .
presence of g T sq and g χ sq can be seen in Figures 12 and 13. It is very clear that in all cases, even when tuned correctly, the inclusion of either term could only marginally reduce γ t by some 0.01γ p . On the other hand, for other values that are not tuned the net effect of said inclusion rapidly becomes strongly destabilising. Of course we cannot guarantee that it is impossible to employ some trial function that will significantly stabilise γ t . Nevertheless a range of different exploratory investigations with different trial functions hints that the velocity space stabilisation cannot move γ t much at all below its minimum value among the configurations of figure 10.

A reduced mode tertiary instability Dimits shift prediction
At last we proceed to the motivation behind the present work: attempting 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, but introducing a 4M -subscript to explicitly show that a "reduced" 4M truncated tertiary instability is considered, the predicted Dimits Figure 13. The equivalent result to figure 12, but where instead g χ sq = p2(v 2 ⊥ ) exp −v 2 with p2(v 2 ⊥ ) being the second degree v 2 ⊥ -polynomial such that g T sq has vanishing zeroth and second velocity moments but non-vanishing fourth velocity moment. Once again, the inclusion of a higher moment perturbation is seen to predominantly destabilise the tertiary instability, with marginal stabilisation not exceeding ∼ 2%. threshold was obtained from the critical solution of the equation system involving the four system parameters p, q, ϕ q , and R/L n that is given by Here the last equation is the "most important", as its solution specifies a L n -value that is taken to correspond to the Dimits threshold. The reason why the 4M-reduction was employed essentially reduced to the fact that it greatly increases computational efficiency, a very valuable feature in increasingly complex systems. Additionally, a 4M prediction was there found to be lower than a non-truncated prediction. This was advantageous since it was observed that zonal profiles characterising that Dimits range typically were "robustly stable" to small-amplitude drift waves rather than merely tertiary stable. The lower 4M prediction could thus be thought to approximate this effect, but may not prove generally suitable.
We will now go through Equations (5.1)-(5.4) to make their meaning clear. Starting from the bottom, Equation (5.4) succinctly expresses that the Dimits threshold corresponds to a tertiary stability threshold, specifically one which is captured by the 4Mreduced system. 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 single zonal mode in the system, and whose instability threshold is taken to approximate a typical full zonal profile, is then set through Equations (5.2) and (5.3), where the former fixes its amplitude and the latter its wavelength, both so that the tertiary mode is maximally stable.
Equations (5.1)-(5.4) taken together can thus be summarised as the conditions for marginal stability, (5.4), of the most unstable tertiary mode, (5.1), assuming a maximally stabilising sinusoidal zonal profile, (5.2) and (5.3). 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, if it exists, is able to stabilise the system. Particularly, since energy will be injected at a faster rate for unstable profiles, the system will then be more free, and thus prone, to evolve into a different configuration.
Rather than accepting the above description wholesale, it is prudent to again review its applicability and validity. To this end we recall that, as indicated by Figures 2, 7, and 8, γ 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 7 makes it patently clear that 5.1 is a suitable and robust condition. That leaves just the 4M-reduction to a single zonal mode of (5.2) and (5.3), but once again, as a result of the observation of figure 11 that a few modes gives rise to most stabilisation and the evolution towards stable flows entailed by figure 8, it is hinted that a reduced mode description should prove apt.

Adapting the prediction
Now, the astute reader will clearly have identified the problem with Condition (5.2) as written: formulated for a fluid model it does not prescribe how one should choose g sq . Of course, this is in principle easily fixed by extending the optimisation to find the maximally stabilising distributions, but straightforwardly doing so entirely eliminates the simplicity and computational efficiency of the prediction with the vast possibility space. We therefore reduce our attention to some trial functions with n e + n i variable parameters α si with i ∈ 1, 2, ..., n s , the determination of whom is given by ∂γ t ∂α si = 0, (5.6) which are added to Conditions (5.1)-(5.4) to complete our system. Here we must obviously attempt to find as minimal a set as possible that we can still trust to lead us reasonably close to the optimal configuration. Now with a drive term proportional to f 0s (T ), it is natural to assume a similar Maxwellian velocity-dependence, but an initial, exploratory investigation revealed that f 0s (bT ) with b ≈ 0.2 generally yielded much lower tertiary growth rates. This proved consistent enough that, though included as an optimisation parameter, it could safely Figure 14. Growth rate γ t at the Dimits threshold η = 0.25 and R/Ln = 1.8 of the full tertiary mode, i.e. qG = 32 (see the convergence of figure 9), with poloidal wavenumber p = 0.6, normalised to the primary growth rate γ p , for a sinusoidal zonal profile ϕ k with radial wavenumber q and where the ion/electron non-adiabatic responses are correctly out of phase with each other for maximum stability, like seen in figure 10. Strong instability commences only when q > p and, unless q is large, only for large zonal flows. Below this threshold the zonal flow is uniformly stabilising, to varying degrees. Note the peculiar feature that γ t exhibits two separate local minima with respect to the zonal amplitude ϕ k for q 1.3. be left at this value with minimal impact for typical values q p and ϕ q ∼ 1 − 10. Beyond this, the analysis of § 4.3, and figure 10 in particular, makes it clear that the cross-species phase n geq /n giq is a vital parameter to include. On the other hand, § 4.4 strongly hints that, when searching for stable distributions, 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γ 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 a massively favourable result for the present prediction, and means that we need not search for further variable parameters. Thus we will employ (4.12) in (5.5) with α 3 = 0.2, thereby restricting (5.6) to just α 1 and α 2 , which should prove sufficient.
Turning to figure 14, we can see the full (not 4M-reduced) tertiary growth rate γ t , for the parameters chosen according to (5.6), as a function of the zonal amplitude ϕ q and the radial wavenumber q. The 4M-case looks qualitatively similar, except that the stabilisation at larger zonal amplitudes is much greater, with the minimum, i.e. the solution to (5.1)-(5.3) and (5.6), occurring for ϕ q around an order of magnitude larger than the full tertiary minimum at q ≈ 0.55 and ϕ q ≈ 5 in figure 14. This difference 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. The fact that typical nonlinear values, not to mention the most stable 4M-value, correspond to points far outside the valid range of a 4M-truncation, like figure 9 (b) demonstrates, furthermore indicates that caution is warranted when making the reduction in only considering direct sideband coupling. Now a stark feature is observed in figure 14 for large ϕ q as q crosses the boundary q = p. Though it is not shown, below this value there always exists stabilising zonal distributions g sq with γ t < γ p , whereas above it none cannot be found. This is the behaviour of KHlike modes (Zhu et al. 2018), with the generalised KH instability criterion that the radial wavenumber must exceed the poloidal of gyrokinetics, k 2 x > k 2 y , already well known (Kim & Diamond 2002;Diamond et al. 2005). Thus we can easily identify the unstable upper right region of 14 as KH-dominated, with an associated k x k y |ϕ|-dependent growth rate. Note however that the stabilising profiles elsewhere nevertheless only have a finite stabilising effect and fail to actually be stable, since the value of R/L n used in figure 14 corresponds to the Dimits threshold.

Prediction comparison
Finally proceeding to put the prediction into action, the γ t -minimum of a given R/L nvalue, with poloidal wavenumber p given by (5.1), is obtained through a steepest-descent walk of the parameters ϕ q , q, α 1 , and α 2 discussed in § 5.1 that simultaneously solves Equations (5.2), (5.3), and (5.6). Because of the presence of multiple undriven modes, the convergence of γ t is exceptionally slow around the Dimits threshold, and so the algorithm is commenced at R/L n well into the expected turbulent range. Having thus selected a single, minimal γ t (R/L n )-value, R/L n is then repeatedly lowered and a steepest descent search re-initiated to determine the slope which enables Newton's algorithm to be deployed for γ t (R/L n ) to find the solution to (5.4) while ameliorating the problem of slow convergence. The result of this entire procedure compared with the result of nonlinear simulations, ran for up to t = 8000 and categorised through the measure Θ 0.1 of (3.2), can be seen in figure 15, both for the 4M-prediction and when effectively no restrictive truncation is employed with q G = 32. The 4M-prediction is found to match the observed Dimits threshold remarkably well. However, as mentioned in § 5.1, the 4M solution produces high zonal amplitudes that would couple to many more sidebands if they were included. The 32M-solution does not suffer from this problem, and also does not differ much from the observed Dimits threshold. Curiously, the critical gradient of the latter case is actually lower than the former, a feature that will be discussed below. In either case, the key result is that a reduced mode tertiary instability prediction is able to accurately capture the Dimits transition.

Discussion
In this article we have performed nonlinear flux tube simulations and employed extensive associated 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 Figure 15. Characterisation of the qualitative long-term behaviour of nonlinear simulations, with unstable or quiescent zonal flows, as a function of the system configuration R/Ln and η = Ln/LT . The estimated Dimits threshold, obtained as the solution to (5.1)-(5.4) and (5.6) using both a 4M and 32M calculation, is also plotted, exhibiting remarkable agreement.
for toroidal geometries is a suitable means of accounting for the fact that only Rosenbluth-Hinton residual flows can be stable. Hopefully these can be captured by a similar trial function 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), the Dimits regime dynamics here observed differ in some key ways, which should be kept in mind while evaluating the prediction as we go from one to the other. There, the Dimits range was characterised by large but transient turbulent bursts before the complete removal of drift waves after the emergence of a robustly stable zonal profile, i.e. one whose similar profiles are also stable so that instability is not reestablished before drift waves have completely decayed. Here, only at the very lowest Dimits range do drift waves completely decay. Instead, for a majority of the Dimits range localised, low amplitude drift waves continue to be present, but 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 burst do not occur within the Dimits regime of Z-pinch gyrokinetics like those in the fluid limit, they appear as R/L n is further increased and then dominate transport before continuous turbulence eventually develops. Indeed, these bursts seem to be of the same nature as those in Hallenbert & Plunk (2021), exhibiting little change in zonal energy, rather than conventional predator-prey type zonal-drift wave oscillations (Malkov et al. 2001;. Such oscillations arise as a result of zonal collisional damping. Higher collisionalities, though 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), instead the present bursts coincide with periods of zonal profile cycling between profiles that are tertiary unstable at different points.
For computational efficiency, the present prediction uses a 4M reduction. Additionally, 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 .
It is helpful to contrast the view above with the conventional shearing picture, in which energy is continually shuffled, through shear, to smaller and smaller radial scale sidebands before being dissipated. Despite the difference in energy transfer, both are consistent with our tertiary picture, with poloidal bands exhibiting collective, coherent behaviour. Nevertheless, Terry et al. (2021), based on findings of Whelan et al. (2019) that finite β-effects initially diminishes transport despite zonal flow shearing then becoming weakened, considers it likely that direct mirror mode coupling rather than zonal shearing is the dominant stabilising mechanism within the Dimits shift. For the collisionless Zpinch 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 9 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 32Mprediction 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 transition well. Though 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. Due to the necessary inclusion of kinetic electrons and the complete conservation of zonal flows, rather than solely Rosenbluth-Hinton remnant states (Rosenbluth & Hinton 1998), the linear tertiary instability picture is considerably more complicated. Both effects vastly increase the zonal configuration space which complicates the prediction method. Though we attempted to capture this, we can of course not guarantee that more stable zonal configurations can be found outside our limited selection of trial functions. 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 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 can matter beyond its mere density/temperature moments 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 10.
Regarding the tertiary instability, there have been many recent advances in its understanding. From its introduction, it was long assumed to be a mere KH-like instability . Here we instead stress that, in accordance with recent understanding (Zhu & Dodin 2021), 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. Indeed, we expect the KH-instability affecting zonal flows not to play a significant role for the Z-pinch. As figure 14 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, though 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 threshold. 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). Though 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.
where A s = (λ + iω dsp )I + C 1s , B s = iZ s f 0s (ω * sp − ω dsp )I + C 2s , (A 2) and the elements of C 1s , C 2s , and J 0s are given by and Γ 0smn = I 0 ((p 2 + m 2 q 2 )m s /m i )e −(p 2 +m 2 q 2 )ms/mi δ mn . (A 6) Restricting the zonal profile to the single mode q (and its conjugate −q) we now set ϕ l = g sl = 0 for l = ±1. Then we can write where the real condition g −q = g * q was used. Next we introduce the transformation matrix where m is the row number. It is clearly diagonal, T T = T , but also orthonormal with determinant −1, so that T 2 = I. As is easily checked, the T -transformation replaces the elements of a matrix with their mirror like E mn ↔ (−1) m+n E −m−n (remember that the indexing runs from −j to j). By (A 3)-(A 6), we therefore have TC 1s T = C * 1s , TC 2s T = C * 2s , TJ 0s T = J 0s , and T Γ 0s T = Γ 0s . Thus (A 1) can equivalently be expressed as i.e. the tertiary dispersion relation of a single mode zonal profile is unchanged under the substitution g sq → g * sq . Of course this substitution also modifies ϕ q , but if we also include a translation of the periodic domain (which patently does not modify the tertiary instability) in our substitution like g sq → g * sq ϕ 2 q |ϕ q | 2 , (A 10) it is clear upon insertion into the quasineutrality condition (2.8) that ϕ q now also remains unchanged.