2.1 The (modified) Hasegawa–Mima equation

The paradigmatic model for density-gradient-driven drift waves is the Hasegawa–Mima equation (HME), which captures advection of both the background and perturbed ion gyrocentre densities caused by fluctuations of the electrostatic potential on a local segment of the outboard midplane. The system is two-dimensional in that it neglects any toroidal variation in the electrostatic potential (i.e. $k_{\Vert }\ll k_{\bot }$ ); thus the fields can be described with just radial and poloidal coordinates.

The HME is given by

where $\boldsymbol{v}=(-\unicode[STIX]{x2202}_{y}\unicode[STIX]{x1D711},\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D711})$ is the $\boldsymbol{E}\times \boldsymbol{B}$ drift velocity, $\unicode[STIX]{x1D6FD}\doteq a/L_{n}$ parameterizes the density gradient, $L_{n}^{-1}\doteq -\text{d}\,\ln n_{0}/\text{d}r$ is the background density-gradient scale length,

is the generalized vorticity derived from the gyrokinetic Poisson equation assuming an adiabatic electron response $\unicode[STIX]{x1D6FF}n_{e}=\unicode[STIX]{x1D711}$ , $\unicode[STIX]{x1D6FF}n_{e}$ is the perturbed electron density and $\unicode[STIX]{x1D711}\doteq (e\unicode[STIX]{x1D719}/T_{e})(a/\unicode[STIX]{x1D70C}_{\text{s}})$ is the dimensionless perturbation to the electrostatic potential. Here, $\unicode[STIX]{x1D6FF}n_{i}$ is the perturbed ion gyrocentre density normalized to the background density, $e$ is the unit charge, $\unicode[STIX]{x1D719}$ is the perturbed electrostatic potential and $T_{e}$ is the background electron temperature. Throughout this paper, time and space are normalized by $a/c_{\text{s}}$ and $\unicode[STIX]{x1D70C}_{\text{s}}$ , respectively, where $a$ is the minor radius, $\unicode[STIX]{x1D70C}_{\text{s}}\doteq c_{\text{s}}/\unicode[STIX]{x1D714}_{\text{ci}}$ is the sound radius and $c_{\text{s}}\doteq (ZT_{e}/m_{i})^{1/2}$ is the sound speed. Coordinates are in a local Cartesian grid where $x$ represents the radial coordinate and $y$ represents the poloidal coordinate. The operator $\unicode[STIX]{x1D6FB}_{\bot }$ denotes the gradient perpendicular to the magnetic field and acts in the $x$ and $y$ direction for this particular set-up. For any function $f$ ,

which defines the Poisson bracket $\{.\,\text{,}\,.\}$ . This nonlinearity, which arises from the $\boldsymbol{E}\times \boldsymbol{B}$ drift, consist of two components: the ion-polarization or vorticity nonlinearity $\{\unicode[STIX]{x1D711},\unicode[STIX]{x1D6FB}_{\bot }^{2}\unicode[STIX]{x1D711}\}$ and the $\boldsymbol{E}\times \boldsymbol{B}$ nonlinearity $\{\unicode[STIX]{x1D711},\unicode[STIX]{x1D6FF}n_{e}\}$ . While the ion-polarization nonlinearity has a specific form, the $\boldsymbol{E}\times \boldsymbol{B}$ nonlinearity can vary depending on which physical effects are included in the electron response. As the Poisson bracket is anticommutative and bilinear, it follows that for purely adiabatic electron response,

Thus, the $\boldsymbol{E}\times \boldsymbol{B}$ nonlinearity does not appear in the HME.

The nonlinear interaction in the HME conserves two quadratic invariants. These are the energy density

and the generalized enstrophy density

Because the HME is a two-dimensional equation with two nonlinearly conserved quantities, the system experiences a dual cascade where injected energy at some intermediate scale simultaneously causes a flow of energy to larger scales and a flow of enstrophy to smaller scales. Note that conservation of the generalized enstrophy ${\mathcal{Z}}$ is a direct result of the nonlinear interaction being in the form of a Poisson bracket and does not depend on the specific form of $\unicode[STIX]{x1D701}$ given by (2.2). This fact is important for the Terry–Horton equation, which will be discussed in § 2.2.

One shortcoming of the HME is that it does not capture the correct zonal-flow physics as seen in more complete gyrokinetic simulations (Dorland & Hammett Reference Dorland and Hammett1993; Hammett et al. Reference Hammett, Beer, Dorland, Cowley and Smith1993). This shortcoming stems from how one derives the adiabatic electron response. For the HME, this is done via parallel force balance in the electron momentum equation:

where $E_{\Vert }$ is the electric field parallel to the magnetic field, $n_{e}$ and $P_{e}$ are respectively the electron density and pressure and $\unicode[STIX]{x1D6FB}_{\Vert }$ is the gradient along the magnetic field ( $\unicode[STIX]{x1D6FB}_{\Vert }\rightarrow \unicode[STIX]{x2202}_{z}$ for HME). Assuming isothermal electrons, (2.7) leads to

in dimensionless units. This has the solution

where $C$ is constant along the magnetic field. Naïvely, one could set this constant to zero and arrive at the HME. However, to correctly determine the value of $C$ , one must use an additional constraint.

Physically, electrons that evolve adiabatically will experience no radial transport, and so the density perturbation averaged over a flux surface must vanish. Thus, taking the flux surface average of (2.9) leads to $C=-\langle \unicode[STIX]{x1D711}\rangle _{\unicode[STIX]{x1D713}}$ , or

where $\langle \ldots \rangle _{\unicode[STIX]{x1D713}}$ denotes the flux surface average which, in two dimensions, is equivalent to the zonal average $\langle \,f\rangle$ ,

Equation (2.2) is thus modified to now read

where $\widehat{\unicode[STIX]{x1D6FC}}$ is an operator that is zero when acting on zonal modes and is unity otherwise. Notice that unlike (2.4), the adiabatic electron response given by (2.10) does not disappear in the Poisson bracket (Jenko et al. Reference Jenko, Dorland, Kotschenreuther and Rogers2000), so the $\boldsymbol{E}\times \boldsymbol{B}$ nonlinearity does appear in the mHME as $\{\unicode[STIX]{x1D711},\unicode[STIX]{x1D6FF}n_{e}\}=\{\langle \unicode[STIX]{x1D711}\rangle ,\unicode[STIX]{x1D711}\}$ , the importance of which is discussed in § 5.

The Hasegawa–Mima equation that uses this form of the modified vorticity is referred to as the modified-Hasegawa–Mima equation (mHME). While the nonlinear interaction is changed by this new Poisson equation, it still conserves both the energy density ${\mathcal{E}}$ and the generalized enstrophy density ${\mathcal{Z}}$ , where $\unicode[STIX]{x1D701}$ in ${\mathcal{Z}}$ is now given by (2.12). It is helpful to define the zonal ( $k_{y}=0$ ) and drift-wave ( $k_{y}\neq 0$ ) energy densities,

where $\unicode[STIX]{x1D711}^{\prime }=\unicode[STIX]{x1D711}-\langle \unicode[STIX]{x1D711}\rangle$ . These quantities will be helpful when determining whether a system is in a zonal-flow-dominated state or a turbulence-dominated state, and will be later used in § 4.2. These definitions remain unchanged for both the Terry–Horton equation and modified-Terry–Horton equation (discussed below).

Both the HME and the mHME have been extensively studied in terms of both the modulational instability (Connaughton et al. Reference Connaughton, Nadiga, Nazarenko and Quinn2010), which concerns the instability of a background drift wave to a zonal flow, and the more general zonostrophic instability (Srinivasan & Young Reference Srinivasan and Young2012; Parker & Krommes Reference Parker and Krommes2013, Reference Parker and Krommes2014), which encompasses the instability of any statistically homogeneous steady state (including steady states of single realizations) to a zonal flow.

2.2 The Terry–Horton equation

Another shortcoming of the HME is the lack of irreversibility. Studies of the HME usually add dissipation by hand and are stochastically forced to drive the system. An alternative to adding forcing is to include destabilizing effects in the electron response obtained iteratively from kinetic theory. These effects materialize most simply as an additional non-Hermitian operator $\widehat{\unicode[STIX]{x1D6FF}}$ in the electron response,

This results in a modification of the Poisson equation (2.2), leading to

or

when Fourier transformed. Here, $\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{k}}$ is a real variable that depends on the wavevector $\boldsymbol{k}$ . The resulting system is named the Terry–Horton equation (THE) and has been extensively studied (Terry & Horton Reference Terry and Horton1982, Reference Terry and Horton1983). It has, however, fallen out of favour with respect to the more rigorously derived Hasegawa–Wakatani system, which has now become the standard system of equations for dealing with linearly unstable drift waves in a two-field fluid system.

The non-adiabatic part of the electron response $\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{k}}$ for the THE can be chosen to describe various types of physical mechanisms. As an example, for the untrapped collisionless electron–wave resonance (the ‘universal’ mode), one has in physical units (Tang Reference Tang1978)

where $\unicode[STIX]{x1D714}_{\ast e}\doteq ck_{y}T_{e}/eB_{0}L_{n}$ is the electron diamagnetic drift frequency, $c$ is the speed of light, $B_{0}$ is the background magnetic field magnitude and $\unicode[STIX]{x1D702}_{e}\doteq \text{d}\,\ln T_{e}/\text{d}\,\ln n_{0}$ . To proceed, $\unicode[STIX]{x1D714}_{\boldsymbol{k}}$ is taken to be the frequency of the linearized HME. Then, in dimensionless units,

where $\unicode[STIX]{x1D6FF}_{0}=(\unicode[STIX]{x03C0}/2)^{1/2}(m_{e}/m_{i})^{1/2}/|k_{\Vert }L_{n}|$ is a constant of order unity that parameterizes the parallel wavenumber $k_{\Vert }$ (as $L_{n}$ is already parameterized by $\unicode[STIX]{x1D6FD}$ ). One arrives at the original form of $\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{k}}=\unicode[STIX]{x1D6FF}_{0}k_{y}(k_{\bot }^{2}-\unicode[STIX]{x1D702}_{e}/2)$ given by Terry & Horton (Reference Terry and Horton1983) by taking the long-wavelength ( $k_{\bot }^{2}\ll 1$ ) limit. As another example, for trapped collisional electron dynamics,

where $\unicode[STIX]{x1D6FF}_{0}=\unicode[STIX]{x1D702}_{e}(\unicode[STIX]{x1D716}/2)^{1/2}(6/\unicode[STIX]{x03C0}^{1/2})/\unicode[STIX]{x1D708}_{\text{eff}}$ is a constant of order unity, $\unicode[STIX]{x1D716}$ is the inverse aspect ratio and $\unicode[STIX]{x1D708}_{\text{eff}}$ is the collisional detrapping rate. A thorough review of such mechanisms is given by Tang (Reference Tang1978).

In addition to instability, the THE equation introduces non-zero particle transport. For instance, if one takes the non-adiabatic electron response to be simply $\text{i}\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{k}}=\text{i}k_{y}$ , one obtains a non-zero particle flux

which is positive–definite. Another important feature of the THE is that it has only one nonlinearly conserved quantity, the generalized enstrophy density ${\mathcal{Z}}$ , which can alter the system’s ability to experience a dual cascade. For instance, $\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{k}}$ as given by (2.20) becomes much larger relative to $\unicode[STIX]{x1D6FB}_{\bot }^{2}$ at large scales, leading to degeneracy between ${\mathcal{E}}$ and ${\mathcal{Z}}$ . This specific form of $\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{k}}$ has been shown by Liang et al. (Reference Liang, Diamond, Wang, Newman and Terry1993) to disable the dual cascade. However, $\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{k}}$ as given by (2.19) is small relative to $\unicode[STIX]{x1D6FB}_{\bot }^{2}$ at both small and large scales. For this type of system the dual cascade is expected to remain prevalent.

3.1 Description

The model that is the focus of this article is a modified version of the THE that is designed to capture the essential zonal physics found in gyrokinetic simulations via a corrected electron response

The resulting system is hence referred to as the modified Terry–Horton equation, and is given by

where $\boldsymbol{v}=(-\unicode[STIX]{x2202}_{y}\unicode[STIX]{x1D711},\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D711})$ and

in Fourier space. Here, $D$ is a damping operator that, in Fourier space, is assumed to be even in both $k_{x}$ and $k_{y}$ . This model contains two modifications to the original THE. First, the adiabatic electron response is modified to ensure that electrons do not respond to a potential that is constant along a flux surface. This results in an enhancement of the zonal interaction between drift waves in the nonlinear term. The second modification is the appearance of the $\widehat{\unicode[STIX]{x1D6FC}}$ operator in front of the damping term, which ensures that only the drift-wave modes ( $k_{y}\neq 0$ ) are linearly damped. By doing so, the mTHE is made to model the residual Rosenbluth–Hinton zonal states (Rosenbluth & Hinton Reference Rosenbluth and Hinton1998) witnessed in the simulations performed in Dimits et al. (Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000). As a result, any state that consists purely of zonal flows is a steady-state solution to (3.2). The damping on the drift waves can then be interpreted to be related to Landau damping of the potential fluctuation. As this is a fluid model, it is agnostic to the eventual fate of the fine-scale velocity-space structure that would result in the ion distribution function.

One may also add a separate damping component to the zonal flows, as was done in Lin et al. (Reference Lin, Hahm, Lee, Tang and Diamond1999). This damping typically results in bursty behaviour involving transitions between zonally dominated states and turbulence-dominated states within what would normally be the region of the Dimits shift. While this phenomenon is interesting in its own right, it is not touched upon in this article.

A question that should be raised is whether or not adding non-adiabatic effects to the electron response will also require modification to the constant of integration when solving the parallel electron force balance. However, electrons that are close to adiabatic should still not respond to a potential that is constant along a flux surface. For the remainder of this paper, I shall assume that these non-adiabatic effects are sufficiently small so as to not affect the adiabatic component of the response. Finally, I require that the flux surface average of the non-adiabatic electron response $\langle \text{i}\widehat{\unicode[STIX]{x1D6FF}}\unicode[STIX]{x1D711}\rangle _{\unicode[STIX]{x1D713}}=0$ , which is the case for the examples noted in § 2.2. This is akin to saying that the background density gradients lie in the radial direction.

3.2 Linear properties

The eigenvalues of (3.2) for the linearization around the zero state can be readily calculated. For drift-wave modes, the growth rates and frequencies for a time dependence of $\text{e}^{\unicode[STIX]{x1D706}_{\boldsymbol{k}}t}$ where $\unicode[STIX]{x1D706}_{\boldsymbol{k}}=\unicode[STIX]{x1D6FE}_{\boldsymbol{k}}-\text{i}\unicode[STIX]{x1D714}_{\boldsymbol{k}}$ are, respectively,

For zonal modes, both are identically zero by construction, as I have neglected collisional damping. For the remainder of this article, I assume reflection symmetry in the $x$ direction for all linear terms (i.e. $\unicode[STIX]{x1D706}(k_{x},k_{y})=\unicode[STIX]{x1D706}(-k_{x},k_{y})$ ).

3.3 Quasilinear model

I also consider numerous approximations to (3.2). The typical first step is to decompose the modified vorticity into zonal and non-zonal components, viz. $\unicode[STIX]{x1D701}=\langle \unicode[STIX]{x1D701}\rangle +\unicode[STIX]{x1D701}^{\prime }$ . By zonally averaging (3.2) and subtracting the result from the original equation, one obtains the new system of equations

where $u^{\prime }\doteq -\unicode[STIX]{x2202}_{y}\unicode[STIX]{x1D711}^{\prime }$ , $v^{\prime }\doteq \unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D711}^{\prime }$ , $U(x)\doteq \unicode[STIX]{x2202}_{x}\langle \unicode[STIX]{x1D711}\rangle$ is the zonal velocity, $F^{\prime }\doteq \{\unicode[STIX]{x1D711}^{\prime },\unicode[STIX]{x1D701}^{\prime }\}-\langle \{\unicode[STIX]{x1D711}^{\prime },\unicode[STIX]{x1D701}^{\prime }\}\rangle$ and $\langle \langle \cdots \rangle \rangle$ denotes the total spatial average. The last term in (3.6b ) results from integrating in $x$ to arrive at the equation for zonal velocity and is chosen to ensure mathematical consistency for $k_{x}=0$ . So far, these equations are an exact description of the original equation (3.2). One arrives at the quasilinear system (QL) when the eddy–eddy interactions in the fluctuation equation (3.6a ) are neglected, which is done by simply setting $F^{\prime }=0$ .

The physical effects that the quasilinear equations neglect are the eddy–eddy self-interactions. However, they do retain the eddy–eddy interactions that act on zonal modes, which appear in the form of a Reynolds stress (first term on the right-hand side of (3.6b )), as well as a term unique to the Terry–Horton equation (second term on the right-hand side of (3.6b )). The latter term describes radial $\boldsymbol{E}\times \boldsymbol{B}$ advection of the background electron gradient for the model $\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{k}}$ given by (2.19). It has already been shown that the quasilinear system captures the essential qualitative aspects of systems that are zonally dominant (Farrell & Ioannou Reference Farrell and Ioannou2009; Parker & Krommes Reference Parker and Krommes2013), and has provided many key insights into the role of zonal-flow-catalyzed energy transfer in the stabilization of drift-wave turbulence.

3.4 Four-mode truncation

Previous work on the HME and the two-field ITG model has concentrated on low-order Galerkin truncations, focusing on the modulational instability of a single mode to calculate zonal-flow growth rates, as well as on calculating the Dimits shift using the tools of dynamical systems (Kolesnikov & Krommes Reference Kolesnikov and Krommes2005a ,Reference Kolesnikov and Krommes b ). To connect this article to this previous work, I formulate a four-mode truncation (4MT) of (3.2); in Fourier space,

where $\unicode[STIX]{x1D6FF}_{\boldsymbol{a},\boldsymbol{b}}$ denotes the Kronecker delta. In the 4MT, I retain a pure drift-wave mode $\boldsymbol{p}=(0,p_{y})$ , a pure zonal mode $\boldsymbol{q}=(q_{x},0)$ and two sidebands $\boldsymbol{r}_{\pm }=(\pm q_{x},p_{y})$ . In addition, the complex conjugate modes are retained in order to satisfy the reality condition. I also assume that $\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{r}_{-}}=\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{r}_{+}}\doteq \widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{r}}$ , which is the case for the two examples of $\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{k}}$ given in § 2.2. Finally, $\widehat{\unicode[STIX]{x1D6FC}}_{\boldsymbol{q}}$ is kept arbitrary to compare with previous results.

The resulting set of equations is

along with their complex-conjugate counterparts. Here, $\widehat{\unicode[STIX]{x1D6FF}}_{\pm }\doteq \widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{p}}\pm \widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{r}}$ and the values of $M_{\boldsymbol{k}}$ are mode-coupling coefficients of the interaction between zonal flows and drift waves and are given by It will be useful to express these coefficients in terms of their real and imaginary parts, To obtain these, I have used the fact that $\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{-}\boldsymbol{k}}=-\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{k}}$ .

One quantity that can be readily calculated for the 4MT is the point at which the system loses its ability to saturate (becomes globally unstable). For an $N$ -mode truncation, this is done by considering the phase space of modes consisting of the $2N$ real and imaginary components of the mode amplitudes and their corresponding velocities that define the equations of motion. Terry & Horton (Reference Terry and Horton1982) have shown that the necessary condition for the system to be absolutely stable is for the phase space at any arbitrary point to be volume contracting. Equivalently, the sum of the growth rates of every mode must be lesser or equal to zero ( $\sum _{\boldsymbol{k}}\unicode[STIX]{x1D6FE}_{\boldsymbol{k}}\leqslant 0$ ). This condition provides insight into a fundamental property of diagonalized partial differential equations with a Poisson-bracket-type quadratic nonlinearity and constant linear coefficients: when one linearizes around any arbitrary point in the phase space of modes, the diagonal of the resulting matrix can only be populated by linear terms. Thus, if the trace of that matrix is positive, then there exists at least one positive eigenvalue and the system is globally unstable. While this trace is usually negative once a sufficient amount of stable modes is included (e.g. modes at small scales that are viscously damped), this has interesting consequences when a severe truncation of the system can be justified. This will be seen in § 5.3.

For the 4MT considered here, the necessary condition for absolute stability is

although nonlinear mode coupling tends to approximately make this the sufficient condition as well, in the sense that only linear gradients near the necessary threshold seem to destabilize the entire system. This is seen in numerical simulations of the 4MT (see figure 3).

4.1 Set-up

The nonlinear (NL) system (3.2) and quasilinear (QL) system (3.6) are simulated pseudospectrally and dealiased on a square Cartesian grid with $L=20\unicode[STIX]{x03C0}$ and $N=256$ on each side. Time stepping is performed using third-order Adams–Bashforth with an integrating factor. The vorticity at $t=0$ is initialized with Gaussian noise of zero mean and standard deviation $5\times 10^{-3}$ , with no energy in zonal modes. The random number generator used for the initial state is initialized with the same seed for all simulations. For all simulations in this section, $D=\unicode[STIX]{x1D707}-\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}_{\bot }^{2}$ where $\unicode[STIX]{x1D707}=1$ and $\unicode[STIX]{x1D708}=10^{-2}$ . Here, friction $\unicode[STIX]{x1D707}$ is added to dissipate energy at large scales and is added as a preventative measure against a possible inverse cascade. (It is shown in § 5 that the qualitative aspects of the simulations do not depend on the choice of damping operator.) To simulate the four-mode truncation, I choose modes with $q_{x}=p_{y}=1$ (the dependence of this system on the specific values of $q_{x}$ and $p_{y}$ , as well as the damping operator $D$ , is discussed in § 5). Finally, the value of $\unicode[STIX]{x1D6FD}$ at which the Dimits shift ends is denoted as a critical density gradient $\unicode[STIX]{x1D6FD}^{\ast }$ . The size of the shift then is $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FD}\doteq \unicode[STIX]{x1D6FD}^{\ast }-\unicode[STIX]{x1D6FD}_{\text{lin}}$ .

I run separate simulations using the different $\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{k}}$ values given by (2.19) and (2.20), the former with $\unicode[STIX]{x1D6FF}_{0}=2$ and $\unicode[STIX]{x1D702}_{e}=0$ , and the latter with $\unicode[STIX]{x1D6FF}_{0}=1.5$ . To differentiate between the two systems, I denote the former by (NL1), (QL1), etc., and the latter by (NL2), (QL2), etc. For results that are insensitive to the details of $\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{k}}$ , I simply use the unnumbered abbreviations. The parameters given above correspond to a threshold of linear instability at $\unicode[STIX]{x1D6FD}_{\text{lin}}\approx 4.74$ for $\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{k}}$ given by (2.19), and $\unicode[STIX]{x1D6FD}_{\text{lin}}\approx 4.21$ for $\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{k}}$ given by (2.20).

4.2 Direct numerical simulation

Figure 1 contains an animation that shows the evolution of the NL2 system (3.2) with $\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{k}}=1.5\text{i}k_{y}$ within the Dimits shift regime ( $\unicode[STIX]{x1D6FD}=5.5$ ). Additional movies with varying values of $\unicode[STIX]{x1D6FD}$ for the nonlinear system ( $\unicode[STIX]{x1D6FD}=4.5$ , $6.5$ ) as well as a movie of the quasilinear system (3.6) with $\unicode[STIX]{x1D6FD}=5.5$ are included as supplementary movies and are available at https://doi.org/10.1017/S0022377817000708, which can be accessed online. Every simulation begins with a short period of linear damping of stable drift-wave modes and linear growth of unstable ones. If the system begins below the threshold of linear instability ( $\unicode[STIX]{x1D6FD}<\unicode[STIX]{x1D6FD}_{\text{lin}}$ ) then all modes damp and the final state is the zero state; otherwise, growth of the drift-wave modes are clustered around the most unstable mode with growth rate $\unicode[STIX]{x1D6FE}=\unicode[STIX]{x1D6FE}_{\text{max}}$ . Because the perturbation level at this stage is quite small ( $\unicode[STIX]{x1D711}_{\boldsymbol{k}}\ll 1$ ), these modes are allowed to grow without any influence from the nonlinear interaction. As the drift-wave energy density ${\mathcal{E}}^{\text{DW}}$ grows and becomes of order unity, the nonlinear interaction becomes effective in transferring energy to zonal modes, resulting in fast growth in zonal energy.

Figure 2. Direct numerical simulation of the NL2 system (3.2) with $D=1-0.01\unicode[STIX]{x1D6FB}_{\bot }^{2}$ and $\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{k}}=1.5\text{i}k_{y}$ for three values of density-gradient parameter $\unicode[STIX]{x1D6FD}$ . The solid line denotes the zonal component of the energy density, while the dashed line denotes the energy density in all non-zonal modes.

Once ${\mathcal{E}}^{\text{ZF}}\sim {\mathcal{E}}^{\text{DW}}$ , the nonlinear interaction becomes important and two different scenarios can take place. One possibility is for the system to find a stable zonal state that occurs after a burst of turbulent interactions (Kolesnikov & Krommes Reference Kolesnikov and Krommes2005a ,Reference Kolesnikov and Krommes b ). Once a stable zonal spectrum is found, the system then relaxes to a pure zonal state with the non-zonal modes damped away to zero. For simulations near the end of the Dimits shift, the system can cycle through a number of zonal spectra until finding one that is ultimately stable. The time for this to happen typically increases with $\unicode[STIX]{x1D6FD}$ , though for a given realization this may not always strictly be true. The other possibility is that no such stable zonal spectrum can be reached, if it even exists. In this case the system remains in a turbulent state with finite particle flux. As an example, the nonlinear system with $\unicode[STIX]{x1D6FD}=6.5$ and $\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{k}}=\text{i}\unicode[STIX]{x1D6FF}_{0}k_{y}$ was run for a time of $t=25\,000$ without ever reaching a stable zonal state. The system with $\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{k}}$ given by (2.19), on the other hand, tends to always find a stable zonal state, thus it does not exhibit a finite Dimits shift. This situation is further discussed in the last paragraph of this section. The qualitative behaviour of both the QL system and 4MT are also similar, though the quantitative aspects of the 4MT, which will be discussed in § 5.2, are quite different.

Figure 2 shows the evolution of both zonal and drift-wave energy densities of the NL2 system for three values of $\unicode[STIX]{x1D6FD}$ versus time scaled by the respective growth rate of the fastest growing mode. The first case ( $\unicode[STIX]{x1D6FD}=4.5$ ) is slightly above the threshold for linear instability. Here, the stable zonal state is found in a short time, as can be seen in the supplemental movie DS_b4.5_NL.mp4. The second case is near the end of the Dimits shift ( $\unicode[STIX]{x1D6FD}=5.5$ ). Now the system spends significantly more time cycling through several zonal spectra until a final one is found. This is shown in figure 1. Finally, the third case is past the end of the Dimits shift ( $\unicode[STIX]{x1D6FD}=6.5$ ). For this case, turbulence persists and no stable zonal state is found (see supplemental movie DS_b6.5_NL.mp4). These are quantitatively similar for the quasilinear system, and a representative evolution of this system is shown in the supplemental movie DS_b5.5_QL.mp4.

Figure 3. Particle flux for the NL2, QL2 and 4MT2 systems with $D=1-0.01\unicode[STIX]{x1D6FB}_{\bot }^{2}$ and $\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{k}}=1.5\text{i}k_{y}$ as a function of ion-density-gradient parameter $\unicode[STIX]{x1D6FD}$ . The medium-dashed line denotes the linear threshold for instability $\unicode[STIX]{x1D6FD}_{\text{lin}}$ . The dotted-line line denotes predicted end of the Dimits shift $\unicode[STIX]{x1D6FD}^{\ast }$ given from the solution of (5.23).

Figure 3 shows the long-time turbulent flux as a function of $\unicode[STIX]{x1D6FD}$ for the NL2, QL2 and 4MT2 systems with $\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{k}}=1.5\text{i}k_{y}$ . (This plot serves the same purpose as figure 3 in Dimits et al. Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000.) As expected, the Dimits shift is observed as a distinct lack of any flux immediately beyond the point of linear instability. Once the end of the Dimits shift is reached ( $\unicode[STIX]{x1D6FD}^{\ast }\approx 5.75$ for both the NL2 and QL2 systems), turbulence is again allowed to persist and turbulent transport ensues. It is important to note that the Dimits shift is quantitatively similar between the NL and the QL systems. This is a rather profound result, as the NL system contains additional avenues of energy dissipation through the direct cascade, where turbulent eddies continuously self-interact, forming smaller turbulent eddies which eventually leads to damping. On the other hand, eddies in the QL system can only be sheared via zonal flows. However, this shearing is the dominant mechanism when the enhancement of the nonlinear zonal interaction is introduced via a physical adiabatic electron response. (I shall show what this means in § 5.) As a result, energy transfer is principally in the direction of $k_{x}$ . This has interesting consequences, the important one being that a QL-like closure, such as CE2, suffices to capture all the relevant physics needed for the Dimits shift.

It is also important to note the discrepancy of the size of the Dimits shift, as well as the saturation levels, of the 4MT. This discrepancy is simply a result of the fact that there is only a single stable mode that can accept energy. This also explains the sharp discrepancy between the saturated levels of flux between the NL/QL and 4MT; because there are far fewer stable modes, the effective damping rate is greatly reduced in the 4MT, resulting in larger levels of saturation when the system strikes a balance between energy production (which has not changed significantly for the 4MT as the most unstable mode is retained) and energy dissipation (which has). Finally, the necessary condition for absolute stability is given by (3.11). This gives a gradient threshold for saturation of $\unicode[STIX]{x1D6FD}^{\text{sat}}\approx 6.0$ for the parameters given in figure 3. This agrees well with what is observed numerically.

I have mentioned that the system with $\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{k}}$ given by (2.19) does not exhibit a finite Dimits shift; rather, a steady zonal state seems to be always found. Even so, two different regimes materialize, as shown in figure 4. Here, the time taken to reach a steady zonal state is defined as $\unicode[STIX]{x0394}t\doteq t_{f}-t_{i}$ , where $t_{i}$ is the first time where ${\mathcal{E}}^{\text{ZF}}={\mathcal{E}}^{\text{DW}}$ , and $t_{f}$ is the first time where ${\mathcal{E}}^{\text{DW}}/{\mathcal{E}}^{\text{ZF}}=10^{-6}$ . In the first regime where $\unicode[STIX]{x1D6FD}<7$ , steady zonal states are found within $\unicode[STIX]{x0394}t\sim 500$ . When $\unicode[STIX]{x1D6FD}>7$ , this time increases by at least a full order of magnitude. One possible explanation is that in the first regime, the zonal shearing associated with the Dimits shift is operational and zonal states that can sufficiently quench the drift wave turbulence are quickly found, while in the second regime some other channel of energy transfer is eventually established, leading to stable solutions. A potential candidate for energy transfer in this second regime could be the non-local cascade of enstrophy to small scales that is typically associated with the dual cascade. More work needs to be done to better understand this regime; that is left for future study.

5.1 Zonal growth rates of the secondary instability

At the beginning of the simulation, the most unstable drift-wave mode grows without any influence from the zonal or sideband modes in what I shall call the secondary regime. One can then solve (3.8b–d ) with $\unicode[STIX]{x1D711}_{\boldsymbol{p}}=\unicode[STIX]{x1D711}_{0}\text{e}^{(\unicode[STIX]{x1D6FE}_{\boldsymbol{p}}-\text{i}\unicode[STIX]{x1D714}_{\boldsymbol{p}})t}$ , which renders the remaining equations a linear system:

The goal is now to determine the growth rate of the zonal flow against this growing drift-wave mode.

At this point, one can make a further approximation by focusing on the region where the zonal mode grows much faster than the drift-wave mode, making the assumption $\unicode[STIX]{x1D6FE}_{\boldsymbol{p}}=\unicode[STIX]{x1D6FE}_{\boldsymbol{r}}=0$ . Alternatively one can also proceed, without any further assumptions, to derive an evolution equation for $\unicode[STIX]{x1D711}_{\boldsymbol{q}}$ ; this is done here.

To do so, the transformations $\unicode[STIX]{x1D711}_{\boldsymbol{r}_{+}}=\unicode[STIX]{x1D711}_{0}\unicode[STIX]{x1D711}_{\boldsymbol{r}_{+}}^{\prime }\text{e}^{-\text{i}\unicode[STIX]{x1D714}_{\boldsymbol{p}}t}$ and $\unicode[STIX]{x1D711}_{\boldsymbol{r}_{-}}^{\ast }=\unicode[STIX]{x1D711}_{0}^{\ast }{\unicode[STIX]{x1D711}^{\prime }}_{\boldsymbol{r}_{-}}^{\ast }\text{e}^{\text{i}\unicode[STIX]{x1D714}_{\boldsymbol{p}}t}$ are made to eliminate the rapid oscillatory behaviour of the drift-wave mode. New equations of motion are then formulated for the variables $\unicode[STIX]{x1D6EC}_{\pm }\doteq \unicode[STIX]{x1D711}_{\boldsymbol{r}_{+}}^{\prime }\pm {\unicode[STIX]{x1D711}^{\prime }}_{\boldsymbol{r}_{-}}^{\ast }$ . Defining $\unicode[STIX]{x1D6FE}_{\pm }\doteq \unicode[STIX]{x1D6FE}_{\boldsymbol{p}}\pm \unicode[STIX]{x1D6FE}_{\boldsymbol{r}}$ and $\unicode[STIX]{x1D714}_{\pm }\doteq \unicode[STIX]{x1D714}_{\boldsymbol{p}}\pm \unicode[STIX]{x1D714}_{\boldsymbol{r}}$ , equation (5.1a–c ) become

From here, it is a simple exercise to derive an ordinary differential equation for the zonal mode by combining (5.2a–c ):

where

The usual dispersion relation of the modulational instability for the HME is recovered by setting $\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{k}}=0$ and $\widehat{\unicode[STIX]{x1D6FC}}_{\boldsymbol{q}}=1$  (Connaughton et al. Reference Connaughton, Nadiga, Nazarenko and Quinn2010).

The goal is to analyse (5.3) in the asymptotic limits $t\rightarrow 0$ and $t\rightarrow \infty$ . The latter, however, is made difficult by the fact that $C$ , $D\propto |\unicode[STIX]{x1D711}_{0}^{2}|\lll 1$ . Furthermore, $C$ and $D$ can have values with quite different magnitudes relative to each other, exacerbating the situation. The limit $t\rightarrow \infty$ cannot be taken at face value then, as this asymptotic time may occur well beyond the range of validity of this approximation (that is to say, beyond the secondary regime). Thus, equation (5.3) must be analysed in the asymptotic limit , which is coarse-grained according to the relative size of the coefficients $A$ , $B$ , $C$ and $D$ , based on a given set of parameters $q_{x}$ , $p_{y}$ and $\unicode[STIX]{x1D6FD}$ .

For $t\rightarrow 0$ , the $C$ and $D$ coefficients become subdominant. The resulting equation is solved readily with solutions

and $y\sim c$ , where $c$ is a constant. This corresponds to the case where the sideband mode also grows independently with growth rate $\unicode[STIX]{x1D6FE}_{\boldsymbol{r}}$ . The zonal-flow growth rate is then $\unicode[STIX]{x1D6FE}_{+}$ , which is maximized at the largest scale of the box. This early type of growth is seen in both the simulations presented here (see figure 1) as well as the simulations of Rogers et al. (Reference Rogers, Dorland and Kotschenreuther2000).

To analyse the limit , the ansatz $\unicode[STIX]{x1D711}_{\boldsymbol{q}}\sim \text{e}^{S(t)}$ is used in (5.3), resulting in

Two cases are now considered. First, let both $C$ and $D$ be of the same order. One then finds the leading-order behaviours $\unicode[STIX]{x1D711}_{\boldsymbol{q}}\sim \exp (-Dt/C)$ and

where

These apply equally well to the case where $D\rightarrow 0$ . For the other case where $C$ is subdominant, one finds the leading-order behaviour

where

and $I_{3}$ is a cube root of unity. Both (5.7) and (5.9) represent instabilities of Kelvin–Helmholtz type ( $\unicode[STIX]{x1D6FE}\sim |q_{x}p_{y}\unicode[STIX]{x1D711}_{\boldsymbol{p}}|^{\unicode[STIX]{x1D700}}\sim |q_{x}v_{x}|^{\unicode[STIX]{x1D700}}$ where $\unicode[STIX]{x1D700}$ is a constant of order unity). This form of the growth rate agrees with the usual picture of the secondary instability given by Rogers et al. (Reference Rogers, Dorland and Kotschenreuther2000), and is roughly maximized at scales between $k_{y}^{-1}$ and the sound radius (i.e. $k_{x}\sim k_{y}$ to $k_{x}\unicode[STIX]{x1D70C}_{\text{s}}\sim 1$ for $k_{y}\unicode[STIX]{x1D70C}_{\text{s}}\ll 1$ ).

Again, it is emphasized that which asymptotic behaviour is relevant depends on the specific values of $q_{x}$ , $p_{y}$ and $\unicode[STIX]{x1D6FD}$ , as well as where the secondary regime ends. For instance, with $q_{x}=0.6$ , $p_{y}=1.3$ and $\unicode[STIX]{x1D6FD}=5$ , $C$ is subdominant for $t\lesssim 200$ , and so the second situation applies. Figure 5 shows the evolution of the drift-wave energy and the zonal energy for the 4MT with $\unicode[STIX]{x1D6FD}=5$ . Both the numerical solution of (5.3) and the scaling given by (5.9) are displayed, showing excellent agreement.

Unfortunately, it is difficult to extend this analysis to the NL or QL systems in order to determine the dominant zonal mode, as the cumulative effect of the asymptotic behaviour of every mode is somewhat unclear. Even then, such a prediction would only predict a dominant zonal mode during the initial secondary stage. There is no a priori reason why this mode would remain the dominant one once the fully nonlinear interaction stage comes to an end. More sophisticated approaches, such as the wave kinetic equation (Diamond et al. Reference Diamond, Itoh, Itoh and Hahm2005; Parker Reference Parker2016; Ruiz et al. Reference Ruiz, Parker, Shi and Dodin2016), have been used in the past to calculate zonal growth rates with some success. However, this approach assumes a homogeneous background of drift-wave turbulence (Krommes & Kim Reference Krommes and Kim2000), which is not the case in this regime, and so the usefulness of the wave kinetic equation for calculating zonal growth rates of the secondary instability within the Dimits shift regime is at this point uncertain.

Figure 5. Typical energy evolution for the 4MT system as described in § 5.1 with $q_{x}=0.6$ , $p_{y}=1.3$ , $D=1-0.01\unicode[STIX]{x1D6FB}_{\bot }^{2}$ , $\unicode[STIX]{x1D6FD}=5$ and $\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{k}}$ given by (2.19) with $\unicode[STIX]{x1D6FF}_{0}=2$ and $\unicode[STIX]{x1D702}_{e}=0$ . The solid line denotes the drift-wave mode energy while the dot-dashed line denotes the zonal mode energy. The numerical solution to (5.3) is plotted with points, with the t→‘∞’ asymptotic behaviour being shown as a dotted line.

In order to appreciate the effect of the physical adiabatic electron response given by (2.10), it is instructive to reconsider the above analysis in the limit of no growth or dispersion, retaining only the nonlinear interactions (G. Hammett, private communication). This is akin to taking the strong turbulence limit. In addition, the Terry–Horton term $\widehat{\unicode[STIX]{x1D6FF}}$ is set to zero. Then (5.1) immediately leads to the zonal-flow growth rate

where the $\widehat{\unicode[STIX]{x1D6FC}}_{\boldsymbol{k}}$ operator has been kept arbitrary. This growth rate is plotted in figure 6 using the HME and mHME with $p_{y}=0.3$ and $\unicode[STIX]{x1D711}_{0}=1$ . In addition, the same growth rate using the two-dimensional incompressible Navier–Stokes equation (NSE) is calculated by taking $\widehat{\unicode[STIX]{x1D6FC}}_{\boldsymbol{k}}=0$ identically. For $q_{x}$ , $p_{y}\ll 1$ , physical adiabatic electron response under the partial time derivative of (2.1) leads to an enhancement of zonal-flow growth by a factor of $1/(q_{x}p_{y})$ . A similar level of growth is also seen in the NSE. It is this correction to the partial time derivative that led to the idea that zonal flows would grow at the largest scales, forming the idea of what was called the ‘mesoscale’ between microturbulence and equilibrium scales (Diamond et al. Reference Diamond, Itoh, Itoh and Hahm2005). However, when the physical adiabatic electron response is taken into account in the Poisson bracket (resulting in an $\boldsymbol{E}\times \boldsymbol{B}$ nonlinearity), the range of zonal-flow scales in which energy is deposited is greatly increased, down to the scale given by

In tokamak plasmas, $p_{y}\ll 1$ , and so this results in a non-local transfer of energy from the primary drift wave to zonal flows at scales of order ${\sim}\unicode[STIX]{x1D70C}_{\text{s}}$ , which is at odds with the conventional wisdom that energy flows to zonal flows at the largest scales (Diamond et al. Reference Diamond, Itoh, Itoh and Hahm2005). This was previously seen in both Rogers et al. (Reference Rogers, Dorland and Kotschenreuther2000) and Guzdar, Kleva & Chen (Reference Guzdar, Kleva and Chen2001). When physical electron response is inconsistently taken into account under the partial time derivative but not in the Poisson bracket (as was done in Holland et al. Reference Holland, Diamond, Champeaux, Kim, Gurcan, Rosenbluth, Tynan, Crocker, Nevins and Candy2003), the zonal-flow growth rate is greatly reduced, which I have added to figure 6 using a dashed line. While the above analysis of the zonal-flow growth rate is based on a simplified four-mode truncation, the cutoff given by (5.12) is also observed in statistical closures of the mHME, including CE2 and a geometrical optics reduction of CE2 (Parker Reference Parker2016). The importance of non-local energy transfer effected by the $\boldsymbol{E}\times \boldsymbol{B}$ nonlinearity is further elaborated upon in the next section.

Figure 6. Zonal growth rate as a function of the zonal-flow radial wavenumber $q_{x}$ over a large-amplitude drift-wave background with poloidal wavenumber $p_{y}=0.3$ and $\unicode[STIX]{x1D711}_{0}=1$ . The solid line denotes the growth rate using a standard two-dimensional incompressible Navier–Stokes equation (NSE, i.e. $\widehat{\unicode[STIX]{x1D6FC}}$ is identically zero). The dotted line denotes the growth rate for the HME, while the dashed-dotted line denotes the growth rate for the mHME, both with $\unicode[STIX]{x1D6FD}=0$ . The dashed line denotes the mHME when the electron response is modified only within the partial time derivative, as studied by Holland et al. (Reference Holland, Diamond, Champeaux, Kim, Gurcan, Rosenbluth, Tynan, Crocker, Nevins and Candy2003). (From G. Hammett, private communication.)

5.2 Zonal-mode stability analysis of the tertiary instability

The stability analysis of the previous section can be repeated to study the stability of a zonal mode to a drift-wave perturbation, otherwise known as the tertiary instability. For a background zonal state $\unicode[STIX]{x1D711}_{\boldsymbol{q}}=\unicode[STIX]{x1D711}_{0}$ of the 4MT, the linearized equations of motion are

The eigenvalues are quickly found: $\unicode[STIX]{x1D706}=(\unicode[STIX]{x1D6FE}_{\boldsymbol{r}}-\text{i}\unicode[STIX]{x1D714}_{\boldsymbol{r}})$ and

The destabilizing root is the positive branch which has the real component

where

Immediately, certain terms appear with physical importance. For instance, $\unicode[STIX]{x1D714}_{-}$ represents the modulational part of the dispersion relation that is due to dispersive effects, while $\unicode[STIX]{x1D6FE}_{+}$ is the coupling of the linear growth rates by the nonlinear interaction. Equation (5.15) contains a wealth of information that must be carefully parsed.

When the outer discriminant of (5.15) vanishes, the real part of both the unstable and stable eigenvalues becomes $\text{Re}(\unicode[STIX]{x1D706}_{+})=\unicode[STIX]{x1D6FE}_{+}/2$ . This corresponds to a special situation where both the drift-wave and sideband modes become maximally coupled, and occurs when $\unicode[STIX]{x1D6E9}=0$ and $\unicode[STIX]{x1D6FA}\leqslant 0$ . Here, both modes work together in tandem to create a coupled mode with an effective growth rate that is the average of the individual growth rates. This coupled mode will be referred to as the maximally coupled mode (MCM). In principle, MCMs can exist for arbitrary drift-wave/sideband pairings where the primary drift-wave mode contains a non-zero radial component $p_{x}$ , though in this section I only consider the most unstable mode with $p_{x}=0$ .

The upper bound of the Dimits shift for the 4MT can then be quickly found as the solution to the equation $\unicode[STIX]{x1D6FE}_{+}=0$ , or

This predicts the end of the Dimits shift for the system considered in § 4.2 to be $\unicode[STIX]{x1D6FD}^{\ast }=5.2$ , in excellent agreement with the simulation result in figure 3. Interestingly, the magnitude of the background zonal flow $\unicode[STIX]{x1D711}_{0}$ does not appear in this stability criterion.

For zonal flows to be stabilizing in for MCMs to develop, the condition $\unicode[STIX]{x1D6FA}\leqslant 0$ must be satisfied. This requires

It is important to realize that $M_{\boldsymbol{k}}\sim q_{x}p_{y}$ times a factor of order unity. This term then clearly represents a Kelvin–Helmholtz-type destabilization of the zonal mode, and can become dominant for sufficiently large zonal amplitudes when the left-hand side of (5.19) becomes negative. However, this term only needs to be comparable to $\unicode[STIX]{x1D6FE}_{-}^{2}-\unicode[STIX]{x1D714}_{-}^{2}$ in order to spoil the coupling between drift-wave modes. This happens for a sufficiently large value of $q_{x}$ (denoted by $q_{x}^{\ast }$ ) and signifies the smallest scale where the last MCM is formed.

The other condition that must be satisfied in order for MCMs to develop is $\unicode[STIX]{x1D6E9}=0$ . There then exists a zonal amplitude that is most stable, given by

The stability of such a mode has been verified numerically. A few special cases are noted. First, if $M_{\boldsymbol{p}}^{\text{Im}}=M_{\boldsymbol{r}}^{\text{Im}}=0$ and $\unicode[STIX]{x1D714}_{-}\neq 0$ (as would be the case if linear instability was introduced by hand to the (m)HME without modifying the Poisson equation) then maximal coupling can only be achieved with $\unicode[STIX]{x1D711}_{0}\rightarrow \infty$ , though typically the frequency mismatch $\unicode[STIX]{x1D714}_{-}$ is quite small. Second, if $M_{\boldsymbol{p}}^{\text{Im}}=M_{\boldsymbol{r}}^{\text{Im}}=\unicode[STIX]{x1D714}_{-}=0$ , then any amount of zonal amplitude is either stabilizing or destabilizing, depending on the sign of $M_{\boldsymbol{p}}^{\text{Re}}M_{\boldsymbol{r}}^{\text{Re}}-M_{\boldsymbol{p}}^{\text{Im}}M_{\boldsymbol{r}}^{\text{Im}}$ . Once the discriminant becomes purely imaginary, additional zonal amplitude ceases to make any further difference in terms of stability. This happens when $\unicode[STIX]{x1D6FA}=0$ . Finally, in the (m)HME limit, $\unicode[STIX]{x1D6E9}=0$ and $\unicode[STIX]{x1D6FA}=-\unicode[STIX]{x1D714}_{-}^{2}-8|\unicode[STIX]{x1D711}_{0}|^{2}M_{\boldsymbol{p}}^{\text{Re}}M_{\boldsymbol{r}}^{\text{Re}}$ . Clearly, the effect of frequency mismatch between modes is a stabilizing influence, and is related to the requirement of resonant triads in the theory of wave turbulence where the size of the nonlinearity is asymptotically small compared to the linear terms.

In the Introduction, four ingredients in a successful model of the Dimits shift were highlighted, one of which was the necessity for a nonlinear boost in the interaction between zonal flows and drift waves. This boost is provided by the $\widehat{\unicode[STIX]{x1D6FC}}$ operator, which modifies the ion continuity equation in both the partial derivative and the Poisson bracket of the HME. While the former effect has been appreciated in the literature (Diamond et al. Reference Diamond, Itoh, Itoh and Hahm2005), the latter has largely gone unnoticed, and in some cases is entirely neglected (e.g. Holland et al. Reference Holland, Diamond, Champeaux, Kim, Gurcan, Rosenbluth, Tynan, Crocker, Nevins and Candy2003, Makwana, Terry & Kim Reference Makwana, Terry and Kim2012). To appreciate the importance of the $\boldsymbol{E}\times \boldsymbol{B}$ nonlinearity, consider the criterion given by (5.19) using the mode-coupling coefficients of the (m)HME where $M_{\boldsymbol{p}}^{\text{Im}}M_{\boldsymbol{r}}^{\text{Im}}=0$ ((3.9) with $\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{k}}=0$ ) with a background drift wave at wavenumber $p_{y}$ . In order for energy to flow from an unstable drift-wave mode to a stable one through a zonal flow, the zonal flow itself must be stable to the Kelvin–Helmholtz instability. Otherwise, energy is fed back from the zonal flow to the unstable drift-wave mode. For a zonal flow to be stabilizing, the condition $M_{\boldsymbol{p}}^{\text{Re}}M_{\boldsymbol{r}}^{\text{Re}}>0$ must be met, leading to destabilizing zonal flows at scales given by

for the HME ( $\widehat{\unicode[STIX]{x1D6FC}}_{\boldsymbol{q}}=1$ ). This is equivalent to the stability threshold given by Diamond et al. (Reference Diamond, Itoh, Itoh and Hahm2005), where the hydrodynamic limit was considered. However, the $\widehat{\unicode[STIX]{x1D6FC}}$ operator in the mHME (viz. (3.9) with $\widehat{\unicode[STIX]{x1D6FC}}_{\boldsymbol{q}}=0$ ) leads to zonal flows that are uniformly stable against the Kelvin–Helmholtz instability down to the scale $1/q_{x}^{\ast }$ given by

Typically, $p_{y}\ll 1$ in toroidal plasmas, leading to $q_{x}^{\ast }$ and $p_{y}$ that are of disparate scales. This is a rather surprising result, as stable zonal flows at scales of the order of ${\sim}\unicode[STIX]{x1D70C}_{\text{s}}$ can be an efficient catalyst of non-local energy transfer, which is in stark contrast to the scale-by-scale transfer in a Kolmogorov-type cascade, whose transfer rate is set by the local eddy turnover time. Additionally, the magnitude of $M_{\boldsymbol{p}}^{\text{Re}}M_{\boldsymbol{r}}^{\text{Re}}$ in the mHME is increased by a factor of $p_{y}^{-2}$ for $p_{y}\ll 1$ compared to that of the HME, resulting in a much stronger stabilizing effect.Footnote ¹

5.3 Heuristic calculation of the Dimits shift

I now generalize the above results for the 4MT in order to heuristically calculate the Dimits shift for the fully nonlinear system. One plausible way to do so is to assume that the zonal interaction couples the most unstable drift-wave mode to every sideband mode within the interaction range $q_{x}^{\ast }$ given in the previous section. Then the threshold for stability is not given by $\unicode[STIX]{x1D6FE}_{+}=0$ , but rather $(2\unicode[STIX]{x03C0}/L_{x})\sum _{\text{MCM}}\unicode[STIX]{x1D6FE}_{\boldsymbol{k}}\approx \int _{0}^{q_{x}^{\ast }}\unicode[STIX]{x1D6FE}_{\boldsymbol{k}}(q^{\prime },p_{y}^{\ast })\,\text{d}q^{\prime }=0$ , where the sum is over all modes with $k_{y}=p_{y}^{\ast }$ in the range $|k_{x}|\leqslant q_{x}^{\ast }$ . Here, I have used the idea of phase-space contraction that was given in § 3.4, but now only consider modes that are well-coupled by zonal shearing to the primary (most unstable) drift-wave mode. Passing into the continuum here is equivalent to taking the large box limit ( $L_{x}\rightarrow \infty$ ).

The physical picture is as follows: as zonal flows are generated from the secondary instability, drift waves become sheared, resulting in a direct transfer of energy to smaller scales. As energy flow is principally horizontal in $k$  space (i.e. energy is not transferred between bands with different $k_{y}$ ), the relevant value of $p_{y}$ that determines stability is that given by the most unstable drift-wave mode, denoted by $p_{y}^{\ast }$ . This transfer can be done efficiently to a cluster of modes down to the scale given by $q_{x}^{\ast }$ , where the effect of zonal flows ceases to be stabilizing. Thus the critical gradient that signifies the end of Dimits shift will be roughly given when this cluster of modes goes unstable (becomes volume expanding, see Terry & Horton Reference Terry and Horton1982).

Quantitatively, this situation can be described as a system of four equations,

and four unknowns $q_{x}^{\ast }$ , $p_{y}^{\ast }$ , $|\unicode[STIX]{x1D711}_{0}|^{2}$ and $\unicode[STIX]{x1D6FD}^{\ast }$ . The first equation relates $p_{y}^{\ast }$ of the most unstable drift-wave mode to the instability parameter $\unicode[STIX]{x1D6FD}^{\ast }$ . The second equation determines the critical gradient where the cluster of coupled modes becomes unstable (volume expanding) and relates $q_{x}^{\ast }$ , $p_{y}^{\ast }$ and $\unicode[STIX]{x1D6FD}^{\ast }$ . Finally, the last two equations leverage the calculations of the 4MT to approximately determine the range at which MCMs can be formed and relates $|\unicode[STIX]{x1D711}_{0}|^{2}$ , $q_{x}^{\ast }$ , $p_{y}^{\ast }$ and $\unicode[STIX]{x1D6FD}^{\ast }$ . As a somewhat rougher approximation, one can alternatively use $q_{x}^{\ast 2}=1+p_{y}^{\ast 2}$ in place of (5.23c,d ), eliminating the need of $|\unicode[STIX]{x1D711}_{0}|^{2}$ . Finally, for systems where the growth rates are rather difficult to come by, one could also use the approximate condition $\unicode[STIX]{x1D6FE}(0,p_{y}^{\ast })+\unicode[STIX]{x1D6FE}(q_{x}^{\ast },p_{y}^{\ast })=0$ , where the coupling between only the most unstable mode and the sideband at scale $k_{x}=q_{x}^{\ast }$ is considered. While this forgoes the phase-space argument of Terry & Horton (Reference Terry and Horton1982) (see § 3.4), it does provide a simple and reasonable estimate of the size of the shift.

Numerical solution of (5.23a–d ) with $\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{k}}=1.5k_{y}$ and $D=1-0.01\unicode[STIX]{x1D6FB}_{\bot }^{2}$ yields $\unicode[STIX]{x1D6FD}^{\ast }\approx 5.79$ and results in a Dimits shift size of $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FD}\approx 1.55$ . This is within a per cent of the value obtained by the direct numerical simulation in § 4.2 ( $\unicode[STIX]{x1D6FD}^{\ast }\gtrsim 5.75$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FD}\approx 1.5$ ). Using the same methodology, the calculated shift for the system parameters given in the caption of figure 4 is roughly $\unicode[STIX]{x1D6FD}^{\ast }\approx 6.77$ . While the Dimits shift for this system is infinite, this agrees well with the end of the first regime at $\unicode[STIX]{x1D6FD}\approx 7$ . These are fairly good estimates, considering they are the result of a straightforward calculation with only a few nonlinear considerations. Notice that this analysis is insensitive to the details of the underlying model; it applies equally well to any Hasegawa–Mima-like equation, and should be generalizable to more complicated ones, provided they retain the $\boldsymbol{E}\times \boldsymbol{B}$ nonlinearity. Additionally, it also works with any generic damping operator. I stress that this is only an estimate; it is neither an upper bound, as additional avenues of efficient energy transfer may exist, nor is it a lower bound, as the existence of a stable state does not guarantee it is physically realizable.

Figure 7. Parameter scan of the nonlinear system with $\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{k}}$ given by (2.20). (a) Uses $D_{\boldsymbol{k}}=1+0.01k_{\bot }^{2}$ , while (b) uses $D_{\boldsymbol{k}}=0.3|k_{y}|$ . Bold dots denote system that end in steady zonal states while crosses denote systems that end in turbulent states. The solid line marks the linear stability threshold, while the dashed line denotes the predicted end of the Dimits shift $\unicode[STIX]{x1D6FD}^{\ast }$ given from the solution of (5.23).

Figure 7 shows a parameter scan of the system with $\widehat{\unicode[STIX]{x1D6FF}}_{\boldsymbol{k}}=\text{i}\unicode[STIX]{x1D6FF}_{0}k_{y}$ over both $\unicode[STIX]{x1D6FF}_{0}$ and $\unicode[STIX]{x1D6FD}$ . Figure 7(a) uses frictional and viscous damping with $D_{\boldsymbol{k}}=1+0.01k_{\bot }^{2}$ , while figure 7(b) uses a damping operator appropriate for a Landau fluid closure with $D_{\boldsymbol{k}}=0.3|k_{y}|$ . These simulations are run for a maximum time of $t=10\,000$ . Bold dots denote simulations that end in steady zonal states while crosses denote systems that end with turbulence. The solid line shows the linear stability threshold while the dashed line shows the threshold calculated from (5.23a–d ). This plot shows excellent agreement between the estimation and the observed numerical boundary of the Dimits shift. Surprisingly, in (a) the estimate also reasonably predicts the boundary at $\unicode[STIX]{x1D6FF}_{0}\approx 2$ .

The above calculation differs fundamentally from the linear stability analyses of both Rogers et al. (Reference Rogers, Dorland and Kotschenreuther2000) and Numata et al. (Reference Numata, Ball and Dewar2007). While the latter set of calculations focus on the stability of simplified background zonal-flow profiles (resulting in localized instabilities), here the details of the background profile are unimportant; rather, the background profile is assumed to be whatever it needs to be in order to maximally couple the primary drift-wave mode to stable ones, and that it is such that Kelvin–Helmholtz modes are not excited. The criteria for stability then is not determined by the flow profile, but rather by the zonal flow’s ability to transfer enough energy to quench the primary instability. Thus, the stability of only a small subset of modes needs to be considered, those that are maximally coupled to the primary drift wave. This approach views the energy transfer in the Dimits shift mainly as a one-step process, and neglects any type of cascade which extends beyond the interaction range of the most unstable mode.

While the non-local mechanism can be taken into account in more complete systems, added complexity comes with additional competing forms of energy transfer. One such mechanism is the dual cascade, which has the ability to transfer both energy to large scales (where large-scale sinks may reside) and enstrophy to small scales. This process is dominated by the largest eddies of the system, and underlies the typical picture of shearing by large-scale zonal flows (Diamond et al. Reference Diamond, Itoh, Itoh and Hahm2005). More complex models also have an increased number of fields, leading to multiple branches (see the recent work of Makwana et al. Reference Makwana, Terry and Kim2012, Reference Makwana, Terry, Pueschel and Hatch2014). Typically, one branch is unstable in a region of $k$ -space while the others are purely stable for all $k$ . Coupling energy between branches, then, can lead to an efficient source of energy dissipation. This is exacerbated when moving to a kinetic system, where there could be an arbitrary number of stable branches. Additional nonlinearities can also arise, such as finite Larmor-radius terms, that can affect the energy transfer process. Finally, other modes of destabilization may exist, such as the $T_{\bot }$ -driven instability found by Rogers et al. (Reference Rogers, Dorland and Kotschenreuther2000). Thus, one has to take into account all these effects as they arise, and determine their relative impact on stability accordingly.