2.1. Magnetic geometry

The magnetic geometry of constant magnetic curvature is chosen because it is the simplest one that enables an ITG instability by coupling the electrostatic potential and the temperature perturbations via the magnetic drift. The integration domain is positioned in the magnetic field of a line of current ($Z$-pinchFootnote ¹) at radial distance $L_B$ from the current line (see figure 1). We define the $x$ and $y$ axes as pointing radially outwards and parallel to the current, respectively. We assume $L_B \gg L_x, L_y$, where $L_x$ and $L_y$ are the ‘radial’ ($x$) and ‘poloidal’ ($y$) sizes of the domain, respectively. Here we use the terms ‘radial’ and ‘poloidal’ to reflect the intended similarity of the domain to one positioned at the outboard midplane in a tokamak geometry. In that sense, we can think of the radial $x$ coordinate as perpendicular to flux surfaces. These surfaces are parametrised by the poloidal $y$ and field-parallel $\hat {\boldsymbol {b}}$ coordinates. Here ${\boldsymbol {B}} = B \hat {\boldsymbol {b}}$ is the magnetic field and the unit vectors $\{\hat {\boldsymbol {x}}, \hat {\boldsymbol {y}}, \hat {\boldsymbol {b}}\}$ form a right-handed basis. In the two-dimensional approximation employed here, all perturbed fields depend only on $x$ and $y$. The magnetic field of the $Z$-pinch with total current $I$ is azimuthal around the current line (as shown on figure 1) and has magnitude $B(x) = 2I/cx$. The radial gradient of this field is then

(2.1)\begin{equation} \frac{1}{B} \frac{\textrm{d}B}{\textrm{d}\kern0.7pt x} = -\frac{1}{L_B}. \end{equation}

This value is constant across the domain to lowest order in $L_x/L_B \ll 1$. We define $L_B$ to be the magnetic scale length. Similarly, we can define the ITG scale length

(2.2)\begin{equation} \frac{1}{L_T} \equiv -\frac{1}{T_i} \frac{\textrm{d}T_i}{\textrm{d}\kern0.7pt x}. \end{equation}

In a tokamak, $L_B$ scales with the major radius of the device, while $L_T$ scales with the minor radius. Here we will take the limit

(2.3)\begin{equation} L_B \gg L_T, \end{equation}

equivalent to a large-aspect-ratio approximation in a tokamak geometry. We do this in order to ensure that the magnetic drift in the density equation is of the appropriate order, see (2.9). This drift is essential for the linear curvature-driven ITG instability that we aim to capture.

Figure 1.

Illustration of the $Z$-pinch magnetic geometry.

2.2. Electron response

The electron density is assumed to follow the modified adiabatic response

(2.4)\begin{equation} \frac{\delta n_e}{n_e} = \frac{e (\phi - \overline{\phi})}{T_e}, \end{equation}

taking into account the fast parallel streaming of the electrons within the flux surfaces of a tokamak, in the small-mass-ratio limit $m_e / m_i \to 0$ (Dorland & Hammett Reference Dorland and Hammett1993; Hammett et al. Reference Hammett, Beer, Dorland, Cowley and Smith1993). Here $\delta n_e$ and $n_e$ are the perturbed and equilibrium electron densities, respectively, $\phi$ is the electric potential, $T_e$ is the electron temperature and

(2.5)\begin{equation} \overline{\phi}(x) \equiv \frac{1}{L_y} \int\,\textrm{d}y \ \phi(x, y) \end{equation}

is the poloidal (zonal) spatial average of the perturbed electric potential $\phi$. We refer to zonally averaged fields as being ‘zonal’. In the two-dimensional approximation, the turbulent fields (e.g. $\phi$) are independent of $z$, and hence we do not need to integrate over the $z$ direction. The difference

(2.6)\begin{equation} \phi'(x, y) \equiv \phi(x, y) - \overline{\phi}(x) \end{equation}

is the ‘nonzonal’ part of the field.

A cautious reader has spotted that there is no way to define flux surfaces in the magnetic geometry of a $Z$-pinch, as the magnetic field lines do not describe two-dimensional surfaces, but rather close on themselves after one turn around the current axis. This problem can be fixed by demanding that the magnetic field be, in fact, sheared: ${{\boldsymbol {B}} = B_0 (\hat {\boldsymbol {z}} + x\hat {\boldsymbol {y}}/L_s)}$, where $L_s$ is the characteristic scale length of the magnetic shear. This is the magnetic field of a helimak (Gentle & He Reference Gentle and He2008). The field lines define cylindrically symmetric concentric flux surfaces and the electron parallel streaming mixes the azimuthal ($z$) and poloidal ($y$) directions. We can then take the limit $L_s \to \infty$ after performing the small-mass-ratio ($m_e/m_i \to 0$) expansion and eliminate magnetic shear from the ion equations, while retaining the flux-surface effect in the electron response.

2.3. Cold-ion limit

The cold-ion limit allows us to simplify the gyroaveraging operator that appears in GK. Its corresponding Fourier-space operator is a multiplication by the Bessel function

(2.7)\begin{equation} J_0\left(\frac{k_{\perp} v_{\perp}}{\varOmega_i}\right) = 1 - \frac{1}{4} k_{\perp}^2\rho_i^2 \frac{v_{\perp}^2}{v_{ti}^2} + {O}\left(k_{\perp}^4\rho_i^4\right). \end{equation}

The square of the ion gyroradius $\rho _i = v_{ti} m_ic/ Z eB$ and the ion temperature $T_i = m_i v_{ti}^2 / 2$ are both proportional to the square of the ion thermal speed $v_{ti}$ (here $m_i$ and $Z$ are the ion mass and charge in units of $e$, respectively). Thus, the cold-ion limit is equivalent to a long-wavelength expansion ${O}\left (k_{\perp } \rho _i\right) \ll 1$ with a finite sound radius ${O}\left (k_{\perp } \rho _s\right) \sim 1$, where $k_{\perp }$ is the perpendicular (to the magnetic field) wavenumber, $\rho _s = \rho _i / \sqrt {2\tau }$ is the sound radius and $\tau = T_i / ZT_e$ is the temperature ratio ($T_e$ is assumed finite). The sound radius $\rho _s$ is the natural normalisation for the microphysical length scales in the problem: see (2.9) and (2.10) below and their derivation in appendix A.

2.4. Model equations

We take the density and temperature moments of the electrostatic ion GK equation and adopt the high-collisionality, cold-ion, long-wavelength, large-aspect-ratio ordering

(2.8)\begin{equation} \frac{\partial_t}{\nu_i} \sim \tau \sim k_{\perp}^2\rho_i^2 \sim \frac{L_T}{L_B} \ll 1 \sim \frac{\varphi}{T}, \end{equation}

where $\varphi \equiv Z_ie\phi / T_i$ is the normalised electric potential, $T = \delta T / T_i$ is the normalised ion-temperature perturbation and $\nu _i$ is the ion–ion collision frequency. The equations that we obtain in appendix A are

(2.9)\begin{align} & \frac{\partial}{\partial t} \left( \tau \varphi' - \frac{1}{2} \rho_i^2\nabla_{\perp}^2 \varphi\right)- \frac{\rho_i v_{ti}}{L_B} \frac{\partial}{\partial y} \left( \varphi + T \right) + \frac{\rho_i v_{ti}}{2L_T} \frac{\partial}{\partial y} \left( \frac{1}{2} \rho_i^2\nabla_{\perp}^2 \varphi \right)\nonumber\\ & \qquad + \frac{1}{2} \rho_i v_{ti} \left( \left\{ \varphi, \tau\varphi' - \frac{1}{2} \rho_i^2\nabla_{\perp}^2 \varphi \right\} + \frac{1}{2}\rho_i^2 \boldsymbol{\nabla_{\perp}} \boldsymbol{\cdot} \left\{ \boldsymbol{\nabla_{\perp}}\varphi, T \right\} \right)\nonumber\\ & \quad = - \frac{1}{2} \chi \rho_i^2 \nabla_{\perp}^4 (a\varphi - bT), \end{align}

(2.10)\begin{equation} \frac{\partial T}{\partial t} + \frac{\rho_i v_{ti}}{2L_T} \frac{\partial\varphi}{\partial y} + \frac{1}{2} \rho_i v_{ti} \left\{ \varphi, T \right\} = \chi \nabla_{\perp}^2 T, \end{equation}

where the Poisson bracket is defined by

(2.11)\begin{equation} \left\{\, f, g \right\} = \hat{\boldsymbol{b}} \boldsymbol{\cdot} \left( \boldsymbol{\nabla} f \times \boldsymbol{\nabla} g \right) = \frac{\partial f}{\partial x} \frac{\partial g}{\partial y} - \frac{\partial f}{\partial y} \frac{\partial g}{\partial x}, \end{equation}

and

(2.12)\begin{equation} \chi \equiv \frac{8}{9} \sqrt{\frac{2}{{\rm \pi}}} \nu_{i}\rho_i^2 \end{equation}

is the thermal diffusivity. The numerical factor in (2.12) and the constants $a = 9/40$, $b = 67/160$ in (2.9) are specific to the Landau collision operator (see appendix A.5). These agree with more general calculations of the collisional perpendicular viscosity, gyroviscosity and collisional heat flux (Mikhailovskii & Tsypin Reference Mikhailovskii and Tsypin1971; Catto & Simakov Reference Catto and Simakov2004, Reference Catto and Simakov2005).

Let us discuss the physics content of (2.9) and (2.10). Equation (2.10) is the more obvious one – it describes the advection of the total temperature (perturbations plus equilibrium) by the $\boldsymbol {E}\times \boldsymbol {B}$ drift $\boldsymbol {V_{E}} = c\hat {\boldsymbol {b}}\times \boldsymbol {\nabla }\phi / B$, and the thermal diffusion perpendicular to the magnetic field. Indeed, (2.10) can be rewritten as

(2.13)\begin{equation} \frac{\textrm{d}}{\textrm{d}t} \left(\delta T + T_i\right) = \chi \nabla_{\perp}^2 \delta T, \end{equation}

where the advective time derivative is

(2.14)\begin{equation} \frac{\textrm{d}}{\textrm{d}t} \equiv \frac{\partial}{\partial t} + \boldsymbol{V_{E}} \boldsymbol{\cdot} \boldsymbol{\nabla}. \end{equation}

The advection of the equilibrium temperature $\boldsymbol {V_{E}} \boldsymbol {\cdot } \boldsymbol {\nabla } T_i$, which is the second term on the left-hand side of (2.10), is responsible for the injection of free energy (see § 2.7), causing the ITG instability.

Equation (2.9) describes the time evolution of the sum of perturbed ion density and the vorticity of the $\boldsymbol {E}\times \boldsymbol {B}$ drift velocity:

(2.15)\begin{equation} \frac{\delta n_i}{n_i} + \frac{1}{n_i}\int \textrm{d}^3\boldsymbol{v} \ \left(\varphi - \left\langle \left\langle \varphi \right\rangle_R \right\rangle_r\right)F_i = \tau \varphi' - \frac{1}{2} \rho_i^2\nabla_{\perp}^2 \varphi + {O}\left(k_{\perp}^4\rho_i^4 \varphi\right), \end{equation}

where $n_i$ and $\delta n_i$ are the equilibrium and perturbed ion densities, respectively. See (A 21) for more details. Thus, (2.9) can be thought of both as a perturbed-ion-density equation and as the curl of the perpendicular-momentum equation. The equality in (2.15) follows from the quasi-neutrality condition $Z\delta n_i = \delta n_e$, the electron response (2.4), the approximation (2.7) and the ordering (2.8).

The first term of (2.9) is the time derivative of (2.15). The second term, ${(\rho _i v_{ti} / L_B) \partial _y ( \varphi + T)}$, is the magnetic drift (both curvature and $\boldsymbol {\nabla } B$) of pressure perturbations. This, or rather the $\partial _y T$ part of it, is essential for the curvature-driven ITG instability. It appears in the density equation because the magnetic drift creates charge separation, and hence electrostatic potential, which is then coupled to the perturbed density via quasi-neutrality. The third term, $(-\rho _i v_{ti} / 4L_T) \partial _y ( \rho _i^2\nabla _{\perp }^2 \varphi)$, is a finite-Larmor-radius (FLR) term originating from the gyroaveraged $\boldsymbol {E}\times \boldsymbol {B}$ drift. It is the diamagnetic drift due to the equilibrium temperature gradient. The first of the nonlinear terms represents the advection of the quantity (2.15) by the $\boldsymbol {E}\times \boldsymbol {B}$ drift. The second nonlinear term $\boldsymbol {\nabla _{\perp }} \boldsymbol {\cdot } \left \{ \boldsymbol {\nabla _{\perp }}\varphi , T \right \}$ is another FLR effect, which describes the $\boldsymbol {E}\times \boldsymbol {B}$ advection of diamagnetic momentum. This term provides a crucial source of turbulent poloidal momentum flux that destabilises the ZF profiles, destroying the ZF-dominated Dimits regime. The nature of this term and its role in the Dimits transition are discussed in detail in § 4.3. Note that the nonlinear terms in (2.9) and (2.10) are equivalent to the nonlinearities appearing in the model analysed by Rogers et al. (Reference Rogers, Dorland and Kotschenreuther2000) in the limit (2.8). Finally, the collisional terms in (2.9) represent the viscous damping of the $\boldsymbol {E}\times \boldsymbol {B}$ flow and also couple the density and temperature perturbations. The latter coupling does not appear to be important for the results of this paper, but has been kept for consistency.

To prepare (2.9) and (2.10) for numerical analysis and distil important parameters, we introduce normalised variables and fields:

(2.16)\begin{equation} \left. \begin{aligned} \hat{t} & \equiv \frac{2\rho_s\varOmega_i}{L_B}t, \quad \hat{x} \equiv \frac{x}{\rho_s}, \qquad \hat{y} \equiv \frac{y}{\rho_s}, \\ \hat{\varphi} & \equiv \frac{\tau L_B\varphi}{2\rho_s} = \frac{\tau L_B}{2\rho_s} \frac{Z_ie\phi}{T_i}, \quad \hat{T} \equiv \frac{\tau L_BT}{2\rho_s} = \frac{\tau L_B}{2\rho_s} \frac{\delta T}{T_i}, \\ \kappa_T & \equiv \frac{\tau L_B}{2 L_T}, \quad \hat{\chi} \equiv \frac{L_B}{2\rho_s} \frac{\chi}{\varOmega_i \rho_s^2}, \end{aligned} \right\} \end{equation}

where $\varOmega _i = v_{ti} / \rho _i$ is the ion gyrofrequency. Dropping hats and subscripts ($\boldsymbol {\nabla }_{\perp } \mapsto \boldsymbol {\nabla }$), we obtain from (2.9) and (2.10) the following equations in normalised units:

(2.17)\begin{align} &\partial_t ( \varphi' - \nabla^2 \varphi) - \partial_y \left( \varphi + T \right) + \kappa_T \partial_y \nabla^2 \varphi +\{ \varphi, \varphi' - \nabla^2 \varphi \} + \boldsymbol{\nabla} \boldsymbol{\cdot} \{ \boldsymbol{\nabla} \varphi, T \} \nonumber\\ &\quad=-\chi \nabla^4 (a\varphi - bT), \\[-3pc]\nonumber\end{align}

(2.18)\begin{equation} \partial_tT +\kappa_T\partial_y \varphi +\left\{ \varphi, T \right\} = \chi \nabla^2 T. \end{equation}

These equations have two independent parameters: the normalised equilibrium temperature gradient $\kappa _T$ and the normalised collisionality $\chi$.Footnote ² There are two other parameters – $L_x$ and $L_y$, the domain lengths in $x$ and $y$ – but any physically relevant results must be independent of these if our equations are indeed a valid local model of the plasma. This turns out to be true for the saturated Dimits state.

We solve (2.17) and (2.18) numerically in a doubly periodic box of size $L_x$ and $L_y$ using a pseudo-spectral algorithm. We integrate the linear terms implicitly in time, while the nonlinear terms are integrated explicitly using the Adams–Bashforth three-step method. This integration scheme is similar to the one implemented in the popular GK code GS2 (Kotschenreuther, Rewoldt & Tang Reference Kotschenreuther, Rewoldt and Tang1995b; Dorland et al. Reference Dorland, Jenko, Kotschenreuther and Rogers2000).

2.5. Relationship to Hasegawa–Mima equation and related models

It is easy to see that setting $\kappa _T = 0$ effectively decouples (2.18) from (2.17) – taking an initial condition $T(t = 0) = 0$ then leads to a trivial solution $T(t) = 0$. In that case, (2.17) reduces to

(2.19)\begin{equation} \partial_t ( \varphi' - \nabla^2 \varphi) - \partial_y \varphi +\{ \varphi, \varphi' - \nabla^2 \varphi \} = -a \chi \nabla^4 \varphi, \end{equation}

which is the well-known (modified) Charney–Hasegawa–Mima (mCHM) equation (Hasegawa & Mima Reference Hasegawa and Mima1978) that includes the appropriate modified adiabatic electron response, with viscous damping. Even though we have considered a situation with no equilibrium density gradient, the magnetic drift provides a term identical to the one that would have arisen from the $\boldsymbol {E}\times \boldsymbol {B}$ advection of an inhomogeneous equilibrium density profile. This puts the model considered here in the same class of systems as those proposed by Hasegawa & Wakatani (Reference Hasegawa and Wakatani1983), Terry & Horton (Reference Terry and Horton1983) and others – all essentially extensions of the Hasegawa–Mima equation with additional physics to account for microinstabilities in the plasma.

As (2.19) is contained within the model considered in this paper, (2.17) and (2.18) should, in principle, capture the behaviour of the mCHM equation as well as additional temperature and ITG effects. There has recently been a significant effort to advance the understanding of (2.19) and its relatives (Parker & Krommes Reference Parker and Krommes2013, Reference Parker and Krommes2014; Parker Reference Parker2016; Ruiz et al. Reference Ruiz, Parker, Shi and Dodin2016; Plunk & Bañón Navarro Reference Plunk and Bañón Navarro2017; St-Onge Reference St-Onge2017; Majda, Qi & Cerfon Reference Majda, Qi and Cerfon2018; Zhu, Zhou & Dodin Reference Zhu, Zhou and Dodin2018b, Reference Zhu, Zhou and Dodin2019, Reference Zhu, Zhou and Dodin2020a,Reference Zhu, Zhou and Dodinb; Qi, Majda & Cerfon Reference Qi, Majda and Cerfon2019; Ruiz, Glinsky & Dodin Reference Ruiz, Glinsky and Dodin2019; Zhou, Zhu & Dodin Reference Zhou, Zhu and Dodin2019). The mCHM equation does capture certain important phenomena, such as the generation of ZFs through a secondary instability (see § 2.8); however, its predictive capabilities for ITG turbulence are unclear. In particular, we shall find that the break up of the Dimits state of (2.17) and (2.18) is, in an essential way, governed by the behaviour of the temperature perturbations, which are absent in (2.19) (see § 4).

2.6. Linear physics of ITG instability

Let us analyse the linear stability of (2.17) and (2.18). Dropping the nonlinear terms, we look for Fourier modes $\varphi , T \propto \exp [ (\gamma _{\boldsymbol{k}} - \textrm {i}\omega _{\boldsymbol{k}}) t + \textrm {i}{\boldsymbol {k}}\boldsymbol {\cdot }{\boldsymbol {r}}]$, where $\gamma _{\boldsymbol {k}}$ and $\omega _{\boldsymbol {k}}$ are the real growth rate and frequency, respectively. Figure 2 shows $\gamma _{{\boldsymbol {k}}}$ as a function of the wavenumber ${\boldsymbol {k}}$. Qualitatively it resembles the growth rate of toroidal ITG modes in tokamaks (Horton, Choi & Tang Reference Horton, Choi and Tang1981). This is expected because the mechanism of the toroidal ITG instability is similar to that of the two-dimensional curvature-driven ITG instability. The terms that give rise to the instability are the magnetic drift term $-\partial _y T$ in (2.17) and the $\boldsymbol {E}\times \boldsymbol {B}$ advection of the equilibrium temperature $\kappa _T\partial _y \varphi$ in (2.18). We find that the fastest growing linear modes are radially extended across the entire box, i.e. they have $k_x = 0$. Such modes are sometimes called ‘streamers’.

Figure 2.

(a) Dependence of the growth rate $\gamma _{\boldsymbol {k}}$ on $k_y$ for the streamer modes ($k_x = 0$). (b) Dependence of $\gamma _{\boldsymbol {k}}$ on $k_x$ and $k_y$ for $\kappa _T = 0.36, \chi = 0.1$. The dashed line is the boundary between stable and unstable modes ($\gamma _{\boldsymbol {k}} = 0$). In § 3, we consider nonlinear simulations with these same parameters.

The exact dispersion relation is

(2.20)\begin{align} & (\gamma_{\boldsymbol{k}} - \textrm{i}\omega_{\boldsymbol{k}})^2(1+k^2) +(\gamma_{\boldsymbol{k}} - \textrm{i}\omega_{\boldsymbol{k}})\{ -\textrm{i}k_y (1+\kappa_T k^2) + \chi k^2 [1 + (1+a)k^2]\} \nonumber\\ &\qquad +a\chi^2k^6 - \kappa_T k_y^2 - \textrm{i}k_y\chi k^2 [1+\kappa_T (1-b)k^2] = 0. \end{align}

We can get a good qualitative idea of the properties of the instability by setting $\chi =0$. Then the solution of (2.20) is

(2.21)\begin{equation} \gamma_{\boldsymbol{k}} - \textrm{i}\omega_{\boldsymbol{k}} = \frac{\textrm{i}k_y(1+\kappa_T k^2) \pm k_y\sqrt{4\kappa_T(1+k^2) - (1+\kappa_T k^2)^2}}{2(1+k^2)}, \end{equation}

so the growth rate of the unstable mode is

(2.22)\begin{equation} \gamma_{{\boldsymbol{k}}} = \frac{k_y\sqrt{(2\sqrt{\kappa_T} - 1 + \kappa_T k^2)(2\sqrt{\kappa_T} + 1 - \kappa_T k^2)}}{2(1+k^2)}. \end{equation}

To simplify further, consider $\kappa _T \gg 1 \gg \kappa _T^{-1/4} \gg k$. Then

(2.23)\begin{equation} \gamma_{{\boldsymbol{k}}} \approx k_y \sqrt{\kappa_T}. \end{equation}

This expression, with the normalisations (2.16) undone, is the well-known ‘bad-curvature-instability’ growth rate (Beer Reference Beer1995):

(2.24)\begin{equation} \gamma_{{\boldsymbol{k}}} = \varOmega_i \frac{\rho_i^2 k_y}{\sqrt{2\tau L_BL_T}}. \end{equation}

Note that there is another, physically distinct, ITG instability usually referred to as the ‘slab ITG mode’. This instability relies on coupling density and temperature through parallel-velocity perturbations, and so is naturally three-dimensional (Cowley et al. Reference Cowley, Kulsrud and Sudan1991). This mode is entirely absent from our two-dimensional model.

Now let us return to the general dispersion. An important feature of the modes described by (2.20) is the boundedness of the region of unstable wavenumbers in the ${\boldsymbol {k}}$ plane (figure 2b). This allows us to integrate (2.17) and (2.18) without the need for artificial dissipation. There are both collisionless and collisional mechanisms that lead to the suppression of the ITG instability. Let us consider these mechanisms.

2.6.1. Collisionless bounds on unstable wavenumbers

It is easy to see that, in order to be positive, the collisionless growth rate (2.22) requires $k < k_{\max ,\text {FLR}}$, where

(2.25)\begin{equation} k_{\max,\text{FLR}}^2 = \frac{1 + 2\sqrt{\kappa_T}}{\kappa_T}. \end{equation}

For $\kappa _T < 1/4$, (2.22) also gives a lower bound on the wavenumbers $k$ of the unstable collisionless modes, i.e. $k > k_{\min ,\text {FLR}}$, where

(2.26)\begin{equation} k_{\min,\text{FLR}}^2 = \frac{1 - 2\sqrt{\kappa_T}}{\kappa_T}. \end{equation}

Adding collisions re-establishes the instability at low $k$. We deem this to be an unimportant peculiarity of our model, thus we shall only consider $\kappa _T > 1/4$.

2.6.2. Collisional bounds on unstable wavenumbers

For nonzero ($\chi > 0$) collisionality, the term $a\chi ^2 k^6$ in (2.20) dominates over the ITG term $\kappa _T k_y^2$ when $k$ is large enough and gives strictly damped modes. To show this, let us simplify (2.20) by writing it as

(2.27)\begin{equation} (\gamma_{\boldsymbol{k}} - \textrm{i}\omega_{\boldsymbol{k}} + A) (\gamma_{\boldsymbol{k}} - \textrm{i}\omega_{\boldsymbol{k}} + B - \textrm{i}C) - fAB + \textrm{i}gAC = 0, \end{equation}

where

(2.28a–e)\begin{align} A = \chi k^2, \quad B = \frac{a\chi k^4}{1+k^2}, \quad C = k_y \frac{1+\kappa_T k^2}{1+k^2}, \quad f = \frac{\kappa_T k_y^2}{a\chi^2k^6}, \quad g = \frac{b \kappa_T k^2}{1 + \kappa_T k^2}. \end{align}

The instability threshold is given by $\gamma _{\boldsymbol {k}} = 0$. The real and imaginary parts of (2.27) for $\gamma _{\boldsymbol {k}} = 0$ are

(2.29)\begin{gather} -\omega_{\boldsymbol{k}}^2 -\omega_{\boldsymbol{k}} C + (1-f)AB = 0, \end{gather}

(2.30)\begin{gather} -\omega_{\boldsymbol{k}} (A+B) - (1-g)AC = 0. \end{gather}

Substituting into (2.29) the value of $\omega _{\boldsymbol {k}}$ derived from (2.30), and using $A \neq 0$, we find

(2.31)\begin{equation} g(1-g)AC^2 + BC^2(1-g) + (1-f)B(A+B)^2 = 0. \end{equation}

Since $g \in (0, 1)$,Footnote ³ a necessary condition for instability is

(2.32)\begin{equation} f > 1 \implies a\chi^2k^6 < \kappa_T k_y^2. \end{equation}

Thus, the region of unstable modes is bounded by $k < k_{\max ,\chi }$, where

(2.33)\begin{equation} k_{\max,\chi}^2 = \sqrt{\frac{\kappa_T}{a\chi^2}}. \end{equation}

2.7. Conservation laws

Equations (2.17) and (2.18) have several conservations laws describing the time evolution of quantities that would be conserved in the absence of equilibrium gradients and dissipation:

(2.34)\begin{align} &\partial_t \int\,\textrm{d}\kern0.7pt x\,\textrm{d}y \ \tfrac{1}{2} T^2 = -\kappa_T \int\,\textrm{d}\kern0.7pt x\,\textrm{d}y \ T \partial_y \varphi - \chi \int\,\textrm{d}\kern0.7pt x\,\textrm{d}y \ \left( \boldsymbol{\nabla} T \right)^2, \end{align}

(2.35)\begin{align} &\partial_t \int\,\textrm{d}\kern0.7pt x\,\textrm{d}y \ \tfrac{1}{2} \left[\varphi'^2 + (\boldsymbol{\nabla} \varphi)^2\right] \nonumber\\ &\qquad = -\int\,\textrm{d}\kern0.7pt x\,\textrm{d}y \ T \partial_y \varphi - \chi \int\,\textrm{d}\kern0.7pt x\,\textrm{d}y \ (\nabla^2 \varphi) \left( a\nabla^2 \varphi - b \nabla^2 T \right) , \end{align}

(2.36)\begin{align} &\partial_t \int\,\textrm{d}\kern0.7pt x\,\textrm{d}y \ \left[ \tfrac{1}{2} \varphi'^2 + T \varphi' + \tfrac{1}{2} \left(\boldsymbol{\nabla} T + \boldsymbol{\nabla} \varphi\right)^2 \right] \nonumber\\ &\qquad = -\chi \int\,\textrm{d}\kern0.7pt x\,\textrm{d}y \ \left[ \left(\boldsymbol{\nabla} \varphi'\right) \boldsymbol{\cdot} \left(\boldsymbol{\nabla} T\right) + a \left( \nabla^2 \varphi \right)^2 \right. \nonumber\\ &\left.\qquad+(a + 1 - b) \left( \nabla^2 \varphi \right)\left( \nabla^2 T \right) + (1-b)\left( \nabla^2 T \right)^2 \right]. \end{align}

These conservation laws can be deduced directly from (2.17) and (2.18): e.g. (2.34) is obtained by multiplying (2.18) by $T$ and integrating over $x$ and $y$. They are also particular cases of the conservation laws of the GK equation. The conservation of the variance of $T$, given by (2.34), is the lowest-order version of the GK free-energy budget. The other two conservation laws, (2.35) and (2.36), can be derived from the conservation of the two-dimensional GK invariant (see Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009; Plunk et al. Reference Plunk, Cowley, Schekochihin and Tatsuno2010). This invariant is a function of velocity in the GK formalism. The model presented here is based on only two velocity moments of the distribution function, namely density and temperature, and so the two-dimensional invariant yields two independent conservation laws. More specifically, (2.35) is a generalisation of the ‘electrostatic GK invariant’. The derivations of the three invariants of our system directly from the corresponding GK invariants can be found in appendix B.

Equations (2.34) and (2.35) imply that a steady saturated state, i.e. $\partial _t = 0$ for all averaged quantities, can be achieved only if appropriate balance between injection and dissipation terms is established:

(2.37)\begin{equation} \frac{\chi}{\kappa_T} \int \frac{\textrm{d}\kern0.7pt x\,\textrm{d}y}{L_xL_y} \ \left( \boldsymbol{\nabla} T \right)^2 = \chi \int \frac{\textrm{d}\kern0.7pt x\,\textrm{d}y}{L_xL_y} \ (\nabla^2 \varphi) \left( a\nabla^2 \varphi - b \nabla^2 T \right) = Q, \end{equation}

where the total radial heat flux $Q$ isFootnote ⁴

(2.38)\begin{equation} Q = -\frac{1}{L_xL_y}\int\,\textrm{d}\kern0.7pt x\,\textrm{d}y \ T \partial_y \varphi. \end{equation}

Thus, a saturated state would necessarily have a net positive ‘turbulent’ (or ‘anomalous’) heat flux $Q > 0$. Note that the first term on the right-hand side of (2.34), which represents injection of free energy, is $\kappa _T Q$. The turbulent heat flux is enabled by the turbulence excited by the ITG instability.

Note as well that the linearly unstable modes have a positive radial heat flux. Indeed, from (2.38),

(2.39)\begin{equation} Q = \sum_{{\boldsymbol{k}}} \textrm{i}k_y T_{{\boldsymbol{k}}} \varphi_{{\boldsymbol{k}}}^* = \sum_{{\boldsymbol{k}}} \textrm{i} k_y |\varphi_{{\boldsymbol{k}}}|^2 \frac{T_{{\boldsymbol{k}}}}{\varphi_{{\boldsymbol{k}}}}. \end{equation}

The relative phase of the temperature and potential perturbations can be obtained from (2.18):

(2.40)\begin{equation} \frac{T_{{\boldsymbol{k}}}}{\varphi_{{\boldsymbol{k}}}} = \frac{-\textrm{i}k_y\kappa_T}{\gamma_{\boldsymbol{k}} - \textrm{i}\omega_{\boldsymbol{k}} + \chi k^2}, \end{equation}

where $\gamma _{\boldsymbol {k}} - \textrm {i}\omega _{\boldsymbol {k}}$ is the solution of the dispersion relation (2.20). Then

(2.41)\begin{equation} Q = \sum_{{\boldsymbol{k}}} \kappa_T k_y^2 |\varphi_{\boldsymbol{k}}|^2 \frac{\gamma_{\boldsymbol{k}} + \chi k^2}{(\gamma_{\boldsymbol{k}} + \chi k^2)^2 + \omega_{\boldsymbol{k}}^2} > 0 \end{equation}

for the unstable modes, which have $\gamma _{\boldsymbol {k}} > 0$.

Finally, the third conservation law (2.36) has some peculiar properties. First, neither the conserved quantity on the left-hand side nor the dissipation rate on the right-hand side is sign-definite. Second, all of the evolution is dissipative, i.e. this invariant is not injected by any equilibrium gradients and is constant in time if $\chi = 0$.

2.8. Secondary instability

Before we delve into the study of nonlinear saturation, let us show how ZFs can be generated from the linearly unstable ITG modes. Consider the stability of a streamer mode with $k_x = 0, k_y = q$ (the ‘primary’ mode) to infinitesimal ‘secondary’ perturbations:

(2.42)\begin{gather} \varphi = (\varphi_q \textrm{e}^{\textrm{i} q y } + \text{c.c.}) + \delta \varphi (x, y), \end{gather}

(2.43)\begin{gather} T = (T_q \textrm{e}^{\textrm{i} q y } + \text{c.c.}) + \delta T (x, y). \end{gather}

A common way of analysing the secondary instability is to take a Galerkin truncation by considering only four Fourier modes $(k_x, k_y) = \{ (0, q), (p, 0), (p, \pm q) \}$ and their complex conjugates: the $(0, q)$ mode is the primary streamer in (2.42) and (2.43) and the others are

(2.44)\begin{gather} \delta \varphi = (\delta\varphi_+ \textrm{e}^{\textrm{i} q y } + \delta\varphi_- \textrm{e}^{- \textrm{i} q y } + \delta\varphi_0) \textrm{e}^{\textrm{i}px} \textrm{e}^{\gamma_2 t} + \text{c.c.}, \end{gather}

(2.45)\begin{gather} \delta T = (\delta T_+ \textrm{e}^{\textrm{i} q y } +\delta T_- \textrm{e}^{- \textrm{i} q y } + \delta T_0) \textrm{e}^{\textrm{i}px} \textrm{e}^{\gamma_2 t} + \text{c.c.}, \end{gather}

where $p$ is the radial wavenumber of the secondary perturbations, $\delta \varphi _0$ and $\delta T_0$ are the ZF and zonal temperature and $\delta \varphi _{\pm }$ and $\delta T_{\pm }$ are known as ‘sidebands’. Substituting all this into (2.17) and (2.18) and linearising the nonlinear terms for $\delta \varphi \ll \varphi _q$ and $\delta T \ll T_q$, we obtain a closed set of equations. In order to keep things simple, we drop the linear terms in (2.17) and (2.18) – this is valid when the amplitude of the primary mode is large enough, so that interactions with it are more important for the evolution of $\delta \varphi$ and $\delta T$ than the effects of the equilibrium gradients and collisions. Observe that, due to the structure of the Poisson bracket (2.11), all nonlinear terms are proportional to $pq$. Defining for convenience $\gamma _2 \equiv \sqrt {2}pq |\varphi _q| \hat {\gamma }_2$, we obtain the following equation for $\hat {\gamma }_2$:

(2.46)\begin{equation} (\hat{\gamma}_2^2 + U)(\hat{\gamma}_2^2 + V) = W, \end{equation}

where

(2.47)\begin{equation} \left. \begin{aligned} U & = 1 + \frac{q^2{\text{Re}{(T_q / \varphi_q)}}}{1+p^2+q^2},\\ V & = \frac{p^2 {\text{Im}{(T_q / \varphi_q)}}^2 + p^2 \left[1 + {\text{Re}{(T_q / \varphi_q)}}\right]^2 - \left(1+q^2\right)\left[1 + {\text{Re}{(T_q / \varphi_q)}}\right]}{1 + p^2 + q^2},\\ W & = \frac{p^2q^2}{(1+p^2+q^2)^2} \left[|T_q|^2 / |\varphi_q|^2 + 2 {\text{Re}{(T_q / \varphi_q)}}\right] \left[1 + {\text{Re}{(T_q / \varphi_q)}}\right]. \end{aligned} \right\} \end{equation}

We see that the growth rate $\gamma _2$ of the secondary instability depends both on the amplitudes of the primary fields $\varphi _q$ and $T_q$ and on their relative phase.

2.8.1. No temperature perturbation

If we set $T_q = 0$, i.e. ignore the temperature perturbation of the primary streamer, (2.46) gives the well-known dispersion relation for the secondary instability of the modified Hasegawa–Mima model (Rogers et al. Reference Rogers, Dorland and Kotschenreuther2000; Strintzi & Jenko Reference Strintzi and Jenko2007):

(2.48)\begin{equation} \gamma_2^{\text{HM}} = pq |\varphi_q| \sqrt{\frac{2(1+q^2-p^2)}{1+p^2+q^2}}. \end{equation}

This form of the secondary instability has long been associated with the strong ZFs observed numerically in ITG turbulence (Hammett et al. Reference Hammett, Beer, Dorland, Cowley and Smith1993).Footnote ⁵ We will show that the inclusion of the temperature perturbations can introduce qualitative and quantitative changes, and even suppress the secondary instability completely.

2.8.2. Long-wavelength limit

To simplify (2.46), we can consider the long-wavelength limit $p \ll 1$. Then (2.46) gives

(2.49)\begin{equation} \left[\hat{\gamma}_2^2 + 1 + \frac{q^2{\text{Re}{(T_q / \varphi_q)}}}{1+q^2}\right][\hat{\gamma}_2^2 - 1 - {\text{Re}{(T_q / \varphi_q)}}] = {O}(p^2) \approx 0. \end{equation}

Thus, there are two independent branches of the secondary instability with instability conditions given by ${\text {Re}{(T_q / \varphi _q)}} + 1/q^2 < -1$ and ${\text {Re}{(T_q / \varphi _q)}} > -1$, respectively. The second branch is a modified form of the long-wavelength Hasegawa–Mima secondary instability (2.48)Footnote ⁶:

(2.50)\begin{equation} \gamma_2 = pq |\varphi_q| \sqrt{2\left(1+{\text{Re}{\frac{T_q}{\varphi_q}}}\right)}. \end{equation}

We observe that (2.50) relies only on a handful of the nonlinear terms in (2.17) and (2.18). Substituting (2.42) and (2.44) into (2.17) and taking the limit $p \ll 1$ gives us the following equations for $\delta \varphi _0, \delta \varphi _+$ and $\delta \varphi _-$:

(2.51)\begin{gather} \hat{\gamma}_2 \delta\varphi_0 = \frac{1}{|\varphi_q|\sqrt{2}} \left[ \delta\varphi_+(\varphi_q^* + T_q^*) - \delta\varphi_-(\varphi_q + T_q) \right], \end{gather}

(2.52)\begin{gather} \hat{\gamma}_2 \delta\varphi_+ = \frac{1}{|\varphi_q|\sqrt{2}} \varphi_q \delta\varphi_0, \end{gather}

(2.53)\begin{gather} \hat{\gamma}_2 \delta\varphi_- = -\frac{1}{|\varphi_q|\sqrt{2}} \varphi_q^* \delta\varphi_0. \end{gather}

Substituting (2.52) and (2.53) into (2.51) yields precisely (2.50). We do not consider the equations for the temperature perturbations because $\delta T = 0$ is a consistent solution and it corresponds to (2.50). The terms on the right-hand side of (2.52) and (2.53) arise from the zonal advection term $\left \{ \overline {\varphi }, \varphi ' - \nabla ^2\varphi ' \right \}$ in (2.17) and represent the tilting of the primary streamer by the ZF. The terms on the right-hand side of (2.51) are the poloidal $\boldsymbol {E}\times \boldsymbol {B}$ and diamagnetic flows caused by the interaction of the primary mode and the two sidebands $(k_x, k_y) = (p, \pm q)$. The quantity ${\text {Re}{(T_q / \varphi _q)}}$ controls the response of the primary mode to the zonal perturbation: ${\text {Re}{(T_q / \varphi _q)}} > -1$ yields an unstable ZF, while ${\text {Re}{(T_q / \varphi _q)}} < -1$ results in a stable, oscillatory perturbation.

Let us consider the collisionless case ($\chi = 0$), where analytical progress is possible, and ask for what values of $\kappa _T$ are the two modes described by (2.49) unstable. Let us take $(k_x, k_y) = (0, q)$ to be the linear mode with the largest growth rate. We then define the critical gradient $\kappa _T^{\text {sec}}$ for the long-wavelength secondary instability of the fastest-growing streamer as the value of $\kappa _T$ at which ${\text {Re}{(T_q / \varphi _q)}} = -1$. Using (2.21) and the relationship (2.40) between $T_q$ and $\varphi _q$, we obtain

(2.54)\begin{equation} {\text{Re}{\frac{T_q}{\varphi_q}}} = -\frac{1 + \kappa_T q^2}{2}. \end{equation}

To determine $q$, we seek the maximum of $\gamma _{{\boldsymbol {k}}}$, as given by (2.22) for $k_x = 0$ and $k_y = q$. We find

(2.55)\begin{equation} \frac{\partial \gamma_{{\boldsymbol{k}}}}{\partial q} \propto \kappa_T^2 q^6 + 3\kappa_T^2q^4 - q^2 - 4\kappa_T + 1 = 0, \end{equation}

where the equality holds for the most unstable mode. As an equation for $q^2$, (2.55) is a cubic with only one positive solution for $\kappa _T > 1/4$. Substituting that solution into (2.54), we find ${\text {Re}{(T_q/\varphi _q)}}$ as a function of $\kappa _T$. This relationship is given in figure 3. In particular, we obtain that ${\text {Re}{(T_q/\varphi _q)}} = -1$ at $\kappa _T = 1$, as can indeed be verified analytically from (2.55) and (2.54), and ${\text {Re}{(T_q/\varphi _q)}} < -1$ for $\kappa _T > 1$. We also find that ${\text {Re}{(T_q / \varphi _q)}} + 1/q^2 > -1$ always. Thus, for $\kappa _T > 1$, the most unstable collisionless ITG mode is stable to the secondary perturbations. Note that (2.54) depends crucially on the diamagnetic drift $\kappa _T\partial _y \nabla ^2\varphi$ in (2.17). If we do not include the diamagnetic drift, we find that ${\text {Re}{(T_q / \varphi _q)}} = -1/2$ regardless of $\kappa _T$ and $q$, and thus the collisionless secondary instability is never quenched.

Figure 3.

Temperature-gradient dependence of ${\text {Re}{(T_{q}/\varphi _{q})}}$ (solid) and ${\text {Re}{(T_{q}/\varphi _{q})}} + 1/q^2$ (dashed) for the most unstable collisionless ($\chi = 0$) mode. We find that ${\text {Re}{(T_{q}/\varphi _{q})}} < -1$ for $\kappa _T > 1$ and ${\text {Re}{(T_{q}/\varphi _{q})}} + 1/q^2 > -1$ for all $\kappa _T$. The secondary instability is present only for $\kappa _T < 1$.

2.8.3. General case

Let us go back to the general secondary dispersion relation (2.46). Its solution is

(2.56)\begin{equation} \hat{\gamma}_2^2 = \frac{-(U+V) \pm \sqrt{(U-V)^2 + 4W}}{2}, \end{equation}

where

(2.57)\begin{equation} U + V = \frac{p^2\left[1 + {\text{Im}{(T_q / \varphi_q)}}^2\right] + p^2 \left[1 + {\text{Re}{(T_q / \varphi_q)}}\right]^2 - {\text{Re}{(T_q / \varphi_q)}}}{1+p^2+q^2}. \end{equation}

We can use the primary dispersion relation (2.20) to show that ${\text {Re}{(T_{\boldsymbol {k}} / \varphi _{\boldsymbol {k}})}} < 0$, and hence $U+V>0$, for any unstable primary mode with wavenumber $\boldsymbol {k}$. Indeed, the real part of (2.40) is

(2.58)\begin{equation} {\text{Re}{\frac{T_{\boldsymbol{k}}}{\varphi_{\boldsymbol{k}}}}} = \frac{k_y\kappa_T\omega_{\boldsymbol{k}}}{|\gamma_{\boldsymbol{k}} -\textrm{i}\omega_{\boldsymbol{k}} + \chi k^2|^2} < 0 \end{equation}

if $k_y \omega _{\boldsymbol {k}} < 0$. Let us show that this is true. For $k \ll 1$ and $\kappa _T > 1/4$ (the reasons for the latter are discussed at the end of § 2.6.1), the dispersion relation (2.20) gives simply $k_y\omega _{\boldsymbol {k}} = -k_y^2 / 2 < 0$. Since the solutions to (2.20) are continuous functions of $\boldsymbol {k}$, if $k_y\omega _{\boldsymbol{k}}$ changes sign and becomes positive, then $\omega _{\boldsymbol {k}} = 0$ somewhere. However, if we set $\omega _{\boldsymbol {k}} = 0$, the imaginary part of (2.27) gives $\gamma _{\boldsymbol {k}} = (g - 1)A < 0$. Therefore, $k_y\omega _{\boldsymbol {k}}$ cannot change sign within the region of linear instability and so $k_y\omega _{\boldsymbol {k}} < 0$ for all linearly unstable modes.

We now consider the solution (2.56) assuming that the relationship between $\varphi _q$ and $T_q$ is given by (2.40) with $k_y = q$ and $\gamma _{{\boldsymbol {k}}}$ and $\omega _{{\boldsymbol {k}}}$ corresponding to the most unstable mode. This gives us $\gamma _2$ as a function of $\kappa _T$, $\chi$ and $p$. Figure 4 shows the real part of $\gamma _2$ maximised over $p$ for each pair of equilibrium parameters $\kappa _T$ and $\chi$, and the wavenumber $p_{\max }$ at which that maximum is attained. Let us discuss this figure. There are three distinct regions:

1. (i) $\kappa _T < \kappa _T^{\text {sec}}$, where $\kappa _T^{\text {sec}}$ is defined as the value of $\kappa _T$ where ${\text {Re}{(T_q / \varphi _q)}} = -1$; in this region, ${\text {Re}{(T_q / \varphi _q)}} > -1$. Additionally, $UV - W < 0$ for $p = p_{\max }$, so $\hat {\gamma }_2^2$ given by (2.56) is real and positive. The instability exists for arbitrarily small values of $p$ (i.e. for an arbitrarily long wavelength of the ZF). Increasing the temperature gradient $\kappa _T$ towards $\kappa _T^{\text {sec}}$ has a dramatic effect on the secondary instability of the most unstable mode: it diminishes both the growth rate and the region of zonal wavenumbers that go unstable. On the line $\kappa _T = \kappa _T^{\text {sec}}$, $\hat {\gamma }_2$ is purely imaginary and there are no growing secondary modes, just like in the long-wavelength analysis of § 2.8.2. Indeed, substituting ${\text {Re}{(T_q / \varphi _q)}} = -1$ in (2.47), we obtain $W = 0$ and $U, V > 0$. Then, by (2.56), $\hat {\gamma }_2^2 = -U$ or $-V$. Figure 5 ($\kappa _T = 0.7, 1.1, 1.5$) shows $\gamma _2$ versus $p$ in region (i).

2. (ii) $\kappa _T > \kappa _T^{\text {sec}}$. Increasing $\kappa _T$ past $\kappa _T^{\text {sec}}$ changes the fastest-growing secondary mode discontinuously. The fastest-growing secondary mode now has $UV - W > 0$ and $(U-V)^2 + 4W < 0$. Hence the $\hat {\gamma }_2^2$ given by (2.56) is complex. In this region, there is always $\hat {\gamma }_2$ with a positive real part. The peak-growth wavenumber $p_{\max }$ changes discontinuously across the ${\text {Re}{(T_q / \varphi _q)}} = -1$ line. For $\kappa _T > \kappa _T^{\text {sec}}$, the secondary instability does not extend to arbitrarily small $p$ (figure 5, $\kappa _T = 1.9, 2.3$), consistent with the discussion of the long-wavelength secondary instability in § 2.8.2.

3. (iii) $\kappa _T > \kappa _T^{\text {sec}}$, but now $(U-V)^2 + 4W > 0$ for all values of $p$, so $\hat {\gamma }_2^2$ given by (2.56) is real and negative. The location of this region of stability depends on the value of ${\text {Im}{(T_q / \varphi _q)}}$, as well as ${\text {Re}{(T_q / \varphi _q)}}$, and does not have a simple analytic form like the boundary between regions (i) and (ii).

Figure 4.

(a) Secondary-instability growth rate (2.56) of the most unstable streamer mode, maximised over all values of $p$. The growth rate vanishes on the ${\text {Re}{(T_q / \varphi _q)}} = -1$ curve (shown in black). (b) Radial wavenumber $p_{\max }$ at which the maximum growth rate shown in panel (a) is attained. A discontinuity in the most unstable wavenumber across the ${\text {Re}{(T_q / \varphi _q)}} = -1$ curve is evident. The absolute-stability region (iii), visible in the bottom right of both panels, where $\gamma _2 = 0$, is the region where $(U-V)^2 + 4W$ is always positive and $\gamma _2$ is purely imaginary.

Figure 5.

(a) Secondary-instability growth rate for $\chi = 0.1$ and a number of values of $\kappa _T$ versus zonal (radial) wavenumber, as given by (2.46). The primary mode $\varphi _q, T_q$ is taken to be the most unstable one in every case. (b) ${\text {Re}{(T_{\boldsymbol{k}}/\varphi_{\boldsymbol{k}})}}$ for the fastest-growing linear mode ${\boldsymbol {k}} = (0, q)$ for $\chi = 0.1$ as a function of $\kappa _T$. The dashed lines correspond to the same values of $\kappa _T$ as in panel (a).

This analysis of the secondary instability suggests that the system will fail to generate ZFs at a high enough $\kappa _T$. In what follows, we will indeed find that the zonally dominated Dimits regime ceases to exist when the temperature gradient exceeds a certain threshold, $\kappa _T > \kappa _T^c$. However, the naïve guess $\kappa _T^c \approx \kappa _T^{\text {sec}}$, as given by the secondary-instability threshold of the most unstable streamer, does not yield satisfactory agreement with the observed threshold for the Dimits regime (see § 4.4). The secondary-instability picture is incomplete because we must take into account not only whether ZFs can be generated by the ITG modes, but also whether the strong ZFs that support the Dimits regime are resilient to nonzonal perturbations. We shall pick up this topic in § 4.

2.9. Tertiary instability

To study the stability of a zonal state, we consider infinitesimal ITG perturbations over a background of strong ZF and zonal temperature:

(2.59)\begin{equation} \left. \begin{aligned} \varphi & = \overline{\varphi} + \varphi', \ \varphi' \ll \overline{\varphi},\\ T & = \overline{T} + T', \ T' \ll \overline{T}, \end{aligned} \right\} \end{equation}

and linearise (2.17) and (2.18) to obtain evolution equations for $\varphi '$ and $T'$. We refer to the ITG modes governed by these linearised equations as ‘tertiary modes’ and to their linear instability as the ‘tertiary instability’ (in truth, this is just the primary ITG instability but for an equilibrium state modified by the zonal fields). We will discover that this instability can seed turbulent perturbations in the Dimits regime, but is not solely responsible for the transition to strong turbulence (see §§ 3 and 4). Further discussion of the tertiary instability has been exiled to appendix C.

3.1. Shear zones

3.1.1. Suppression of turbulence

The zonal staircase is arranged in such a way that it efficiently suppresses turbulence in the shear zones via strong ZF shear. We find that this shear, $S \equiv \partial _x^2 \overline {\varphi }$, satisfies $S \gtrsim \gamma _{\max }$, where $\gamma _{\max }$ is the largest linear ITG growth rate determined from the dispersion relation (2.20). The notion that ITG turbulence requires comparable $S$ and $\gamma _{\max }$ to be suppressed by shear is known as the ‘quench rule’ (Waltz, Kerbel & Milovich Reference Waltz, Kerbel and Milovich1994; Waltz, Dewar & Garbet Reference Waltz, Dewar and Garbet1998; Kinsey, Waltz & Candy Reference Kinsey, Waltz and Candy2005; Kobayashi & Rogers Reference Kobayashi and Rogers2012). Quantitatively, this is supported by figure 12, which shows that the time- and space-averaged (over the shear zones) shear satisfies $S \approx 2 \gamma _{\max }$ over a range of simulation parameters.Footnote ⁷ Note that the particular snapshots of zonal profiles seen in figure 8 suggest $S \approx \gamma _{\max }$. However, the time-averaged $S$ is larger due to the variation of shear over time (see also figure 13).

Figure 12.

The ZF shear $S$ (time- and space-averaged over the shear zones) versus maximum linear ITG growth rate $\gamma _{\max }$ (as given by 2.20). The data are taken from a number of simulations over a range of parameters: $\kappa _T \in [0.16, 7.29]$ and $\chi \in [0.1, 10]$. The error bars represent the smallest and the largest values of the spatially averaged ZF shear for each simulation. The best-fit line (dashed) is $S \approx 2 \gamma _{\max }$.

Figure 13.

Heat flux $Q$ (orange) and zonal shear $\partial _x^2 \overline {\varphi }$ (black) versus time for SimL. The highlighted (in green) sections of the zonal shear correspond to the quiescent periods. They are identified as those in which $Q$ is smaller than a threshold value (dashed black line), defined as $60\,\%$ of the time-averaged $Q$. The average decay rate of the zonal shear in the thus-identified quiescent periods is $\gamma \approx -3.5\times 10^{-4}$ and the decay rate given by (3.6) is $\gamma _{s} \approx - 3.6\times 10^{-4}$.

3.1.2. Decay of ZFs

Let us study the viscous decay of the ZFs. The equation for the evolution of the ZFs is given by the zonal part of (2.17):

(3.1)\begin{equation} \partial_t \partial_x^2 \overline{\varphi} = \partial_x^2 \overline{\partial_x \varphi \partial_y \left( \varphi + T \right)} + \chi \partial_x^4 \left(a \overline{\varphi} - b \overline{T}\right). \end{equation}

Integrating (3.1) once yields

(3.2)\begin{equation} \partial_t \overline{u}_y = \partial_x \left[\,\overline{\left(\partial_x \varphi\right) \partial_y \left( \varphi + T \right)} + \chi \partial_x^2 \left(a \overline{\varphi} - b \overline{T}\right)\right] = -\partial_x \left( \varPi_t + \varPi_d \right), \end{equation}

where the ZF velocity is $\overline {u}_y \equiv \partial _x \overline {\varphi }$ and we have identified the turbulent, $\varPi _t \equiv - \overline {\left (\partial _x \varphi \right) \partial _y \left ( \varphi + T \right)}$, and diffusive, $\varPi _d \equiv - \chi \partial _x^2 (a \overline {\varphi } - b \overline {T})$, radial fluxes of poloidal momentum. The integration constant in (3.2) is zero because both sides of the equation are exact derivatives with respect to $x$ and our domain is periodic. Integrating (3.2) once more yields a term that is not necessarily an exact derivative – the turbulent momentum flux $\varPi _t$:

(3.3)\begin{equation} \partial_t \overline{\varphi} + \varPi_t + \varPi_d = \varPi, \end{equation}

where the integration constant $\varPi = (1/L_x) \int _0^{L_x}\,\textrm {d}x \ \varPi _t$ is the total box-averaged poloidal momentum flux. However, (2.17) and (2.18) are invariant under the symmetry

(3.4)\begin{equation} x \mapsto -x, \ y \mapsto y, \ \varphi \mapsto -\varphi, \ T \mapsto - T. \end{equation}

Under this symmetry, $\varPi \mapsto -\varPi$, a property of our model inherited from GK (Parra, Barnes & Peeters Reference Parra, Barnes and Peeters2011). Therefore, assuming that the volume-averaged solutions to (2.17) and (2.18) respect (3.4), we conclude that $\varPi = 0$. This is confirmed by our numerical solutions. Thus, the right-hand side of (3.3) vanishes.

During the quiescent periods of the Dimits-state evolution (i.e. between turbulent bursts), the turbulent momentum flux in the shear zones is negligible compared to the diffusive momentum flux, $\varPi _t \ll \varPi _d$. This is a consequence of the suppression of the $k_y \neq 0$ ITG modes by the zonal shear.Footnote ⁸ We also find that the zonal temperature gradient $\partial _x \overline {T}$ is approximately constant in the quiescent shear zones (see figure 8), so $\partial _x^2 \overline {T} = 0$. Therefore, (3.3) becomes

(3.5)\begin{equation} \partial_t \overline{\varphi} = a\chi\partial_x^2 \overline{\varphi}. \end{equation}

This is a diffusion equation governing the viscous decay of the ZFs with a collisional viscosity $a\chi$. As figure 13 shows, quiescent periods of low heat flux and, thus, low levels of nonzonal perturbations are correlated with the periods of decay of the zonal shear. We find that, despite the ever-present turbulence in the convection zones, where (3.5) does not hold, the decay rate of the zonal shear is closely approximated by the viscous decay rate of the longest-wavelength ZF that comprises the zonal staircase, i.e.

(3.6)\begin{equation} \gamma_{s} = -a \chi \left(\frac{2 {\rm \pi} n}{L_x}\right)^2, \end{equation}

where $n$ is the number of periods of the zonal staircase in the domain of radial size $L_x$.

Let us now discuss what the ZF periodicity is.

3.1.3. Scale of ZFs

In general, increasing/decreasing the radial extent of the integration domain by a factor increases/decreases the number $n$ of shear zones by the same factor. This suggests that the characteristic length scale of the staircase, i.e. the time-averaged radial separation of ZF extrema, is determined by a box-size-independent mechanism (see further discussion in § 5). Ascertaining definitively whether this is the case is made difficult by the numerically observed time scales of convergence, which are at least an order of magnitude larger than the longest linear time scales, i.e. than the box-scale diffusion time $t_{\text {box}} = (L_x / 2{\rm \pi})^2 / a\chi$ (see figure 14). Note that the long-time evolution of the zonal profile and its length scale is not accompanied by a significant change in the average heat flux. In fact, the latter appears to reach saturation on time scales comparable to the box-scale diffusion time. Therefore, it is reasonable to trust the numerical values of the box-averaged heat flux (e.g. those shown in figure 6), even though we could not be certain that the zonal profiles have reached ultimate saturation.

Figure 14.

(a) Heat flux, raw (black) and rolling-averaged with a window of $6000$ time units (orange) for $\kappa _T = 0.36$, $\chi = 0.1$, $L_x = 100$, $L_y = 100$ and $125\times 125$ Fourier modes. This is the same as SimL, but we have used small-amplitude white noise as an initial condition. The diffusion time for the box-scale ZF is $t_{\text {box}} = (L_x / 2{\rm \pi})^2 / a\chi \approx 11\,200$. Convergence to a box-sized ZF occurs on a very long time scale ($>50t_{\text {box}}$). (b) Same as panel (a), but for $\kappa _T = 4$, $\chi = 1$, $L_x = 200$ and $L_y = 200$. Here $t_{\text {box}} = (L_x / 2{\rm \pi})^2 / a\chi \approx 4500$. (c) Same as panel (b), but with a single-peak triangular ZF as an initial condition. It is evident that the ZF profile does not converge to a single peak.

As we increase collisionality and thus move away from the near-marginal regime and into the collisionality-independent regime (the plateau seen in figure 6b), the ZFs become more dynamic – they can merge, split and drift, as shown in figure 15(b). Here we focus on the near-marginal regime and the transition to strong turbulence, so this higher-collisionality regime will not be studied.

Figure 15.

Same as figure 11, but for $\chi = 1.5$. The diffusion time for the box-scale ZF is $t_{\text {box}} = (L_x / 2{\rm \pi})^2 / a\chi \approx 750$. The locations of the ZF extrema of the staircase drift significantly over times comparable to $t_{\text {box}}$.

Even though the zonal staircase arises naturally from white-noise initial conditions for both the zonal and the nonzonal fields, its shape suggests initialising the ZFs with a ‘triangular’ pattern. We find that this helps achieve more quickly a ‘less noisy’ and more symmetric final state, which is easier to handle both numerically and analytically. Of course, we do not know in advance how many ‘steps’ the staircase will ‘choose’ to have in the saturated state, so their number for the ‘triangular’ initial condition is just an informed guess. Most results in this paper are from simulations that used such a triangular ZF initial condition, including SimL and SimH. Notable exceptions are figures 6, 12 and 23, where we used data from many simulations, some with white-noise initial conditions and others with ‘triangular’ ones.

3.2. Convection zones

The convection zones located at the extrema of the ZFs contain localised patches of ITG turbulence and have a high radial turbulent heat conductivity (see figures 8 and 9). The imposed equilibrium temperature gradient is flattened in the convection zones and slightly steepened in the shear zones. This results in a staircase-like radial temperature profile, shown in figure 16. The turbulence in the convection zones is driven by a tertiary instability, localised by the zonal shear. In the low-collisionality, near-marginal regime, which we consider to be the most important (see the footnote 2 on p. 9), there is a qualitative difference between the way in which the tertiary instability operates at the ZF maxima and minima. A similar difference exists in both the Hasegawa–Mima equation (Zhu et al. Reference Zhu, Zhou and Dodin2018b) and GK (McMillan et al. Reference McMillan, Hill, Bottino, Jolliet, Vernay and Villard2011).

Figure 16.

Time-averaged normalised total temperature (SimL) relative to the absolute temperature $T_R$ at the right edge of the domain: $(\tau L_B / 2 \rho _s)(T_i - T_R + \overline {\delta T}) = (L_x-x)\kappa _T + \overline {T}$. A strong flattening of the gradient is visible around the ZF minima at $x \approx 25$ and $75$, a weaker one around the ZF maximum at $x \approx 50$.

3.2.3. Scale of convection zones

The width of the convection zones can be characterised by the quantity

(3.7)\begin{equation} \delta \equiv \sqrt{\frac{\partial_x \overline{\varphi}}{\partial_x^3 \overline{\varphi}}}. \end{equation}

Figure 17 shows that $\delta$ does not depend very strongly on either $\kappa _T$ or $\chi$, except far from the marginal state, where collisionality appears to smooth out the gradients in the convection zones and thus increase $\delta$. This suggests that near the Dimits threshold, $\delta$ is an ${O}\left (1\right)$ quantity in the normalised units of (2.17) and (2.18), i.e. it is equal to a few times the sound radius $\rho _s$.

Figure 17.

Numerically determined time-averaged values of $\delta$ (see 3.7) at the global maximum (a) and global minimum (b) of the ZFs in the simulations whose heat flux is shown in figure 6. We find no significant variation of $\delta$ with $\chi$ or $\kappa _T$, except at large $\chi$.

3.3. Ferdinons

After the ZF has decayed sufficiently to weaken its ability to suppress perturbations, vortex-like propagating structures are spawned from the ZF maxima and drift radially through the shear zones. Strikingly similar structures – ferdinons – have been observed in GK simulations with imposed background flow shear (van Wyk et al. Reference van Wyk, Highcock, Schekochihin, Roach, Field and Dorland2016, Reference van Wyk, Highcock, Field, Roach, Schekochihin, Parra and Dorland2017). Figures 8 and 9 ($t = 136$) show a particular instance of the launching of ferdinons (see also figure 10).

As the ferdinons smash into the turbulent modes in the convection zones at the ZF minima, more structures are produced and a burst of turbulence ensues (see figures 8c and 9c). These travelling structures, as well as the resulting turbulence, cause a significant spike in the box-averaged heat flux (see figure 10). As figure 8 shows, they do not carry a significant ZF perturbation. They are created and propagate even if the ZF is held artificially constant in the numerical simulations, but the zonal temperature is left to evolve according to (2.18). In other words, a localised ZF perturbation is not an essential part of these structures.

Figure 18 shows that ferdinons consist of a vortex dipole and a strong temperature perturbation trapped in one of the vortices of this dipole. There are ferdinons carrying both positive (‘hot’) and negative (‘cold’) temperature perturbations. Hot ferdinons drift towards the cooler (right) side of the domain, while the cold ones drift in the opposite direction, towards the hotter (left) side (see also figure 10). Figures 18(a) and 18(b) demonstrate that the direction of the drift does not depend on the sign of the zonal shear. Net flow circulation around the ferdinons is also independent of the sign of the shear – it is always anticlockwise for hot and clockwise for cold ones (see figure 18c).

Figure 18.

(a) A close-up view (from figure 9b) of the temperature and electrostatic-potential perturbations of a hot ferdinon in a region of negative zonal shear, $S = \partial _x^2 \overline {\varphi } < 0$. The arrows represent the local nonzonal $\boldsymbol {E}\times \boldsymbol {B}$ velocity ${\boldsymbol {u}}' = \hat {\boldsymbol {z}}\times \boldsymbol {\nabla }\varphi '$ (in arbitrary units). (b) A hot ferdinon in a region of positive zonal shear, $S = \partial _x^2 \overline {\varphi } > 0$. (c) A cold ferdinon in a region of positive zonal shear, $S = \partial _x^2 \overline {\varphi } > 0$.

Note that the ferdinons that emerge in our simple ITG model bear a striking qualitative resemblance to the avalanches reported by Villard et al. (Reference Villard, Angelino, Bottino, Brunner, Jolliet, McMillan, Tran and Vernay2013) in global GK simulations, namely they propagate both inwards and outwards, but always with a positive heat flux, and originate from the local maxima of the ZF. Simple soliton solutions have already been proposed as a model for GK avalanches (McMillan et al. Reference McMillan, Jolliet, Tran, Villard, Bottino and Angelino2009; McMillan, Pringle & Teaca Reference McMillan, Pringle and Teaca2018). Vortex-dipole solitons called ‘modons’ have been investigated in Hasegawa–Mima-like models of turbulence (Horton & Hasegawa Reference Horton and Hasegawa1994). We do not yet know how and whether any of these are related to the ferdinons that we observe.

Let us discuss what we expect the ferdinon solution to be. Numerically, we find that the existence and propagation of these structures depend crucially on the two ITG-drive terms in (2.17) and (2.18), as well as on the nonlinear terms. In particular, the poloidal localisation of these structures is due to the nonzonal–nonzonal interactions. Indeed, we have found that (2.17) and (2.18) with the nonzonal–nonzonal nonlinear terms taken out (in what is sometimes referred to as the ‘quasi-linear approximation’; see Srinivasan & Young (Reference Srinivasan and Young2012)) do not have ferdinon solutions. However, the quasi-linear system does have soliton solutions that are not localised poloidally, but rather appear to have a definite poloidal wavenumber $k_y$. These solutions might be related to those described by McMillan et al. (Reference McMillan, Jolliet, Tran, Villard, Bottino and Angelino2009) and Zhou, Zhu & Dodin (Reference Zhou, Zhu and Dodin2020). Models have been proposed for structure formation in a sheared flow that rely on the tilting of turbulence by shear and a nonzero group velocity to produce moving structures (McMillan et al. Reference McMillan, Pringle and Teaca2018; Zhou et al. Reference Zhou, Zhu and Dodin2020). The radial group velocity is (at least in the Hasegawa–Mima-related models) proportional to the product $k_x k_y$,Footnote ⁹ which acquires a definite sign in the presence of flow shear (§ 4.3). However, we observe ferdinons moving in both radial directions in regions of definite zonal shear and, thus, definite radial group velocity. Therefore, at the moment, we consider it unlikely that the propagation of ferdinons can be explained using such group-velocity arguments. We leave the detailed investigation of ferdinon generation and propagation for future work.

Understanding ferdinons and their properties can also put an upper bound on the radial scale of the ZF. Indeed, our numerical simulations show that the ZFs can have a well-defined radial scale smaller than the box size (§ 3.1.3). This scale could perhaps be estimated via a causality argument – assuming that ferdinons, and, thus, turbulence, can only propagate a finite radial distance in a region of self-consistently evolving zonal shear, then an infinitely wide shear zone cannot be sustained for long. Note that finite-lifetime ferdinons over a dynamic ZF background, with which they can interact and gain or lose energy, are not in contradiction with the infinite-lifetime ferdinons seen by van Wyk et al. (Reference van Wyk, Highcock, Schekochihin, Roach, Field and Dorland2016), where a constant flow shear was imposed, and thus the shear profile was unable to react to the presence of ferdinons.

Once ferdinons are generated and turbulence develops in the shear zones, our analysis of the viscous decay of the zonal staircase in § 3.1.2, which ignored the turbulent momentum flux, is no longer valid. Instead, we must focus on the effect of the turbulence on the ZFs. We find that the turbulence in the shear zones has a restoring effect on the zonal staircase in the Dimits regime, whereas beyond the Dimits threshold, it inhibits staircase formation.

4.1. Turbulent momentum flux

In order to investigate the way in which ZF profiles are formed and maintained during turbulent periods, let us ask the following question: does ITG turbulence in the shear zones (produced in bursts) have a definite effect on the ZFs, and is it to oppose or to feed them? We shall find that this sheared turbulence enhances the ZFs in the Dimits regime and destroys them beyond the Dimits transition.

The ZF evolution (3.3) is

(4.1)\begin{equation} \partial_t \overline{\varphi} + \varPi_t + \varPi_d = 0. \end{equation}

In § 3.1.2, we discussed the effect of the diffusive momentum flux $\varPi _d$, namely, the viscous decay of the ZFs. For the rest of this section, we focus on the effects of turbulence by examining the turbulent momentum flux $\varPi _t$.

It is evident from (4.1) that ZF saturation requires

(4.2)\begin{equation} \left\langle \varPi_t(t, x) + \varPi_d(t, x) \right\rangle_{\Delta t} \approx 0 \end{equation}

to be satisfied at every radial location $x$, where $\left \langle\, f(t) \right \rangle _{\Delta t} \equiv (1 / \Delta t) \int _{\Delta t} dt\,f(t)$ is a time average in the saturated state over a time $\Delta t$ longer than the typical evolution time of the ZF (e.g. longer than the duration of turbulent bursts if the saturated state is bursty). Recall that $\varPi _d \approx -a\chi \partial _x^2\overline {\varphi } = -a\chi S$ in the shear zones (see § 3.1.2). Therefore, within a shear zone with nearly constant (in time and in space) zonal shear, (4.2) tells us that the time-averaged turbulent momentum flux in that shear zone must have a definite value,Footnote ¹⁰ determined by the local zonal shear. Thus, in the saturated state, $\varPi _t$ will be correlated with $S$.

To quantify this correlation, let us multiply both sides of (4.1) by $S = \partial _x^2 \overline {\varphi }$ and integrate across the radial extent of the domain. We find

(4.3)\begin{equation} \int_0^{L_x}\,\textrm{d}\kern0.7pt x \ \left(\partial_x^2 \overline{\varphi}\right) \partial_t \overline{\varphi} + \int_0^{L_x}\,\textrm{d}\kern0.7pt x \ \varPi_t S + \int_0^{L_x}\,\textrm{d}\kern0.7pt x \ \varPi_d S = 0. \end{equation}

Since $\varPi _d = -\chi \partial _x^2 (a\overline {\varphi } - b\overline {T}) = -a\chi S + b\chi \partial _x^2\overline {T}$, we have

(4.4)\begin{equation} \int_0^{L_x}\,\textrm{d}\kern0.7pt x \ \varPi_d S = -a\chi \int_0^{L_x}\,\textrm{d}\kern0.7pt x \ S^2 + b\chi \int_0^{L_x}\,\textrm{d}\kern0.7pt x \ S \partial_x^2\overline{T} \approx -a\chi \int_0^{L_x}\,\textrm{d}\kern0.7pt x \ S^2, \end{equation}

where we have assumed that the second term is negligible because the main contribution to $S$ comes from the shear zones, where $\partial _x^2 \overline {T} \approx 0$ (see also the discussion leading to 3.5). Therefore, after integrating by parts the first term in (4.3) and time averaging the resulting equation, we find

(4.5)\begin{equation} \left\langle \int_0^{L_x}\,\textrm{d}\kern0.7pt x \ \left( \varPi_t S - a\chi S^2 \right) \right\rangle_{\Delta t} = - \left\langle \tfrac{1}{2} \partial_t \int_0^{L_x}\,\textrm{d}\kern0.7pt x \ \left(\partial_x \overline{\varphi}\right)^2 \right\rangle_{\Delta t} \approx 0. \end{equation}

This gives a prediction for the effective ‘turbulent viscosity’ in the shear zones:

(4.6)\begin{equation} \nu_t \equiv -\frac{ \left\langle \int_0^{L_x}\,\textrm{d}\kern0.7pt x \ \varPi_t S \right\rangle_{\Delta t}}{ \left\langle \int_0^{L_x}\,\textrm{d}\kern0.7pt x \ S^2 \right\rangle_{\Delta t}} \approx -a\chi. \end{equation}

Relation (4.6) is, of course, corroborated by numerical simulations (see figure 19).

Figure 19.

Comparison of estimated turbulent viscosity $\nu _t$ and the collisional viscosity $a\chi$ for the simulations from figure 6. We see that $\nu _t \approx -a\chi$.

4.2. Sign reversal of the turbulent momentum flux at the Dimits threshold

An important consequence of (4.6) is that, in a shear zone, the sign of the turbulent momentum flux must coincide with the sign of the zonal shear. Therefore, if, for certain parameters, sheared turbulence has a momentum flux with a sign opposing that of the local shear, saturation cannot be achieved. We shall see that this is exactly what happens beyond the Dimits threshold.

Let us investigate how turbulence responds to an imposed static zonal profile. We solve (2.17) and (2.18) numerically with an imposed static triangular ZF pattern (i.e. we do not evolve the ZFs at all), in a box of size $L_x = L_y = 100$ and $169\times 169$ Fourier modes, for a range of parameters around the Dimits transition. The chosen ZF pattern is shown in figure 20(a) and is adjusted for every simulation so that the value of the zonal shear in the shear zones matches the largest ITG growth rate for that simulation. The chosen radial scale of the ZF ($= 100$) is inspired by the typical ZF scale that we observe in the low-collisionality regime, and is held fixed as we vary $\chi$ and $\kappa _T$. Then we calculate the effective turbulent viscosity $\nu _t$ associated with the turbulent momentum flux.

Figure 20.

(a) The artificial triangular zonal pattern used to generate panel (b). The zonal shear in the shear zones is chosen to be equal in absolute value to the largest ITG growth rate (represented by the dashed green line). (b) The effective turbulent viscosity $\nu _t$, as defined by (4.6), for the static triangular ZF profile given in panel (a) (coloured data points). The black line represents the numerically established Dimits threshold.

As figure 20(b) shows, we find a negative turbulent viscosity $\nu _t$ (and, thus, a positive correlation between local zonal shear and turbulent momentum flux) in the Dimits regime and a positive $\nu _t$ beyond it (thus, a negative correlation). Let us denote by $\kappa _T^{\text {static}}$ the temperature gradient at which $\nu _t$ reverses its sign. The designation ‘static’ reflects the fact that this is a numerical result for ITG turbulence with an artificially imposed static ZF profile. We find that the value of $\kappa _T^{\text {static}}$ is insensitive to the exact shape of the ZF profile, and, most importantly, it nearly perfectly coincides with the Dimits threshold, i.e. $\kappa _T^{\text {static}} \approx \kappa _T^c$.

Thus, in the Dimits regime, shear zones are resilient because, when the zonal shear there decays due to viscosity and turbulence is thus unleashed, this turbulence acts to reinforce the ZFs and the zonal shear in the shear zones is restored to its turbulence-suppressing level. Beyond the Dimits regime, the zonal staircase cannot be sustained because both turbulence and collisional viscosity act to flatten out the ZFs.

4.3. Reynolds stress and diamagnetic stress

Let us analyse what causes the turbulent momentum flux $\varPi _t$ to reverse its sign at the Dimits transition. We split $\varPi _t = \varPi _{\varphi } + \varPi _{T}$ and define

(4.7a,b)\begin{equation} \varPi_{\varphi} \equiv - \overline{(\partial_x \varphi) (\partial_y \varphi)} = \overline{u_y u_x}, \quad \varPi_{T} \equiv - \overline{(\partial_x \varphi) (\partial_y T)} = \overline{u_y w_x}, \end{equation}

where ${\boldsymbol {u}} = (-\partial _y \varphi , \partial _x \varphi)$ is the $\boldsymbol {E}\times \boldsymbol {B}$ flow and ${\boldsymbol {w}} = (-\partial _y T, \partial _x T)$ is the diamagnetic flow. Here $\varPi _{\varphi }$ is the radial flux of the poloidal momentum due to the Reynolds stress of the $\boldsymbol {E}\times \boldsymbol {B}$ flow and the ‘diamagnetic stress’ $\varPi _{T}$ is a contribution to the momentum flux that physically arises due to the advection of the poloidal diamagnetic flow $w_y = \partial _x T$ by the radial $\boldsymbol {E}\times \boldsymbol {B}$ flow $u_x = -\partial _y \varphi$.Footnote ¹¹ To see this, let us take the zonal average of (2.18) and differentiate once with respect to $x$, to obtain an equation for the zonal diamagnetic flow:

(4.8)\begin{equation} \partial_t \overline{w}_y + \partial_x \overline{u_x w_y} - \partial_x \varPi_T = \chi \partial_x^2 \overline{w}_y. \end{equation}

The evolution of the zonal poloidal $\boldsymbol {E}\times \boldsymbol {B}$ flow is described by (3.2), which can be recast as

(4.9)\begin{equation} \partial_t \overline{u}_y + \partial_x \overline{u_x u_y} + \partial_x \varPi_{T} = -\partial_x \varPi_d. \end{equation}

Added together, (4.8) and (4.9) describe the advection of the total poloidal flow ($\boldsymbol {E}\times \boldsymbol {B}$ + diamagnetic), $v_y \equiv u_y + w_y$, by the radial $\boldsymbol {E}\times \boldsymbol {B}$ flow:

(4.10)\begin{equation} \partial_t \overline{v}_y + \partial_x \overline{u_x v_y} = \text{dissipative terms}. \end{equation}

This makes physical sense because the diamagnetic flow is not a real flow and thus cannot advect anything.

The numerical solutions of (2.17) and (2.18) reveal that $\varPi _{\varphi }$ and $\varPi _{T}$ are in competition: on average, $\varPi _{\varphi }$ has the same sign as the zonal shear $S = \partial _x^2\overline {\varphi }$, while $\varPi _{T}$ has the opposite sign. This is evident in figure 21. Equation (4.1) then tells us that $\varPi _{\varphi }$ feeds the ZFs by increasing the zonal potential $\overline {\varphi }$ in the shear zones of negative zonal shear (where $\overline {\varphi }$ is concave) and decreasing it in the shear zones of positive zonal shear (where $\overline {\varphi }$ is convex), whereas $\varPi _{T}$ relaxes the ZFs by opposing $\varPi _{\varphi }$. Their combined effect either steepens or relaxes the ZF velocity $\overline {u}_y = \partial _x\overline {\varphi }$ at the turning points $\partial _x^2 \overline {\varphi } = 0$, depending on which stress is larger. Figure 22 is an illustration of this. This competition is crucial for the ability of the ZF to reconstitute itself after a turbulent burst and thus sets the threshold for the Dimits regime.

Figure 21.

(a) Time-averaged momentum fluxes $\varPi _{\varphi }$, $\varPi _{T}$ and $\varPi _t = \varPi _{\varphi } + \varPi _{T}$ for saturated ITG turbulence over a fixed zonal background. Note the correlation between the signs of the various fluxes and the zonal shear $S$: the sign of $\varPi _t$ coincides with that of $\varPi _{\varphi }$ and $S$ and opposes the sign of $\varPi _{T}$. This reflects that the temperature gradient is lower than the Dimits threshold, $\kappa _T = 0.36 < \kappa _T^c \approx 1$. The ZF profile used here was extracted from SimH at $t = 10$, but reduced by a factor of $0.8$ in order to allow ITG turbulence to develop in the shear zones. (b) Same as panel (a), but with $\kappa _T = 1.21$. The sign of $\varPi _t$ now opposes the sign of $\varPi _{\varphi }$ and $S$ and coincides with the sign of $\varPi _{T}$. This reflects that $\kappa _T = 1.21 > \kappa _T^c \approx 1$. The extracted ZF is augmented by a factor of $2$ to account for the increased ITG growth rate (due to the larger $\kappa _T$). This is necessary for the turbulence to saturate at numerically feasible amplitudes. Note that saturation is possible only because we have fixed the ZF profile. If the ZF is left to evolve according to (3.2), the poloidal momentum generated by the nonzonal perturbations flattens it and the system fails to reach a finite-amplitude saturated state (see § 4.5).

Figure 22.

An illustration of the enhancing and suppressing effects of $\varPi _{\varphi }$ and $\varPi _{T}$ on the ZF. The black curve shows the $\overline {\varphi }(x)$ profile taken from SimH at $t = 10$. The ZF extrema are the locations where $S = \partial _x^2 \overline {\varphi } = 0$ (marked by dashed lines).

In order to assess what decides the outcome of this competition (i.e. the relative size of $\varPi _{\varphi }$ and $\varPi _{T}$), let us consider how zonal shear affects ITG turbulence. For this purpose, consider a shear zone of radial extent $d$ with a constant zonal shear $S = \partial _x^2 \overline {\varphi }$ throughout it. We can then perform the usual shearing-box change of variables $(t, x, y) \mapsto (\tilde {t}, \tilde {x}, \tilde {y})$, where

(4.11a–c)\begin{equation} \tilde{t} = t, \quad \tilde{x} = x, \quad \tilde{y} = y - Stx. \end{equation}

This coordinate transformation eliminates the spatially inhomogeneous zonal-advection terms $(\partial _x \overline {\varphi })\partial _y = Sx\partial _y$ in (2.17) and (2.18). Consider a Fourier mode ${\varphi , T \propto \exp (\textrm {i} \tilde {k}_x \tilde {x} + \textrm {i} \tilde {k}_y \tilde {y})}$ in this shearing frame. In the laboratory frame $(t, x, y)$, this mode has the form $\varphi , T \propto \exp (\textrm {i} k_x x + \textrm {i} k_y y)$, where

(4.12a,b)\begin{equation} k_x = \tilde{k}_x - St \tilde{k}_y, \quad k_y = \tilde{k}_y. \end{equation}

Thus, the ZF shear introduces an effective drift of the laboratory-frame radial wavenumber.Footnote ¹² The direction of this drift is given by the sign of $S$: $S > 0$ gives rise to an anticorrelation of $k_x$ and $k_y$, i.e. $k_xk_y < 0$, whereas for $S < 0$, $k_xk_y > 0$. Integrating the effect of $\varPi _{\varphi }$ over the sheared region, we obtain

(4.13)\begin{equation} \frac{1}{\textrm{d}}\int\,\textrm{d}\kern0.7pt x \ \varPi_{\varphi} = -\frac{1}{\textrm{d}L_y}\int\,\textrm{d}\kern0.7pt x\,\textrm{d}y\ \left(\partial_x \varphi\right) \left(\partial_y \varphi\right) = -\sum_{{\boldsymbol{k}}} k_xk_y |\varphi_{\boldsymbol{k}}|^2. \end{equation}

Therefore, on average, $\varPi _{\varphi }$ has the same sign as $S$, and, thus, feeds the ZFs that generate the shear zones.Footnote ¹³

We can write a similar expression for the diamagnetic stress:

(4.14)\begin{equation} \frac{1}{\textrm{d}}\int\,\textrm{d}\kern0.7pt x \ \varPi_{T} = -\sum_{{\boldsymbol{k}}} k_xk_y |\varphi_{\boldsymbol{k}}|^2 {\text{Re}{\frac{T_{\boldsymbol{k}}}{\varphi_{\boldsymbol{k}}}}}. \end{equation}

Then the total turbulent momentum flux integrated over a shear region is

(4.15)\begin{equation} \frac{1}{\textrm{d}}\int\,\textrm{d}\kern0.7pt x \ \varPi_t = -\sum_{{\boldsymbol{k}}} k_xk_y |\varphi_{\boldsymbol{k}}|^2 \left(1 + {\text{Re}{\frac{T_{\boldsymbol{k}}}{\varphi_{\boldsymbol{k}}}}} \right). \end{equation}

Recall that we already encountered the quantity ${\text {Re}{(T_{\boldsymbol {k}}/\varphi _{\boldsymbol {k}})}}$ when dealing with the secondary instability in § 2.8. There we found that ${\text {Re}{(T_{\boldsymbol {k}}/\varphi _{\boldsymbol {k}})}} < 0$ for all linearly unstable modes. Thus, linear theory predicts that $\varPi _{T}$ and $\varPi _{\varphi }$ are anticorrelated due to the negative sign of ${\text {Re}{(T_{\boldsymbol {k}}/\varphi _{\boldsymbol {k}})}}$. Now let us perform a more detailed analysis of the linear modes and attempt to construct a model for the Dimits threshold based on it.

4.4. Dimits threshold from linear physics

Using our knowledge of ITG perturbations in a region of uniform ZF shear, and of the turbulent momentum flux produced by them, we can make a heuristic linear-physics-based estimate for the Dimits threshold $\kappa _T^c$. In view of (4.15), it is given by the temperature gradient at which the ‘relevant’ ITG modes have ${\text {Re}{(T_{\boldsymbol {k}}/\varphi _{\boldsymbol {k}})}} = -1$. By ‘relevant’ we mean those ITG modes that dominate the turbulence in the shear zones. It is tempting to assume that these modes would be the most unstable modes in the system. This, however, cannot be the case because the most unstable modes are the radial streamers with $k_x = 0$, but the zonal shear that we find is comparable in magnitude to the largest ITG growth rate ($S \sim \gamma _{\max }$), and, therefore, is bound to break these streamers. Following the discussion in § 4.3, we may assume that the typical ITG modes in sheared turbulence satisfy $k_x \sim \alpha k_y$, where $\alpha \sim S\tau _{\text {nl}}$ characterises how tilted the mode is, $\tau _{\text {nl}}$ being the nonlinear correlation time of the turbulence in the shear zones. If $\tau _{\text {nl}}^{-1} \sim \gamma _{\max } \sim S$, then $\alpha \sim 1$.

Thus, we assume that the relevant modes are tilted with $k_x =\alpha k_y$, where $\alpha \sim 1$ is an unknown tilt parameter that depends on the structure of the turbulence. We then look for the temperature gradient $\kappa _T$ at which the fastest-growing ITG mode with $k_x = \alpha k_y$ satisfies ${\text {Re}{(T_{\boldsymbol {k}}/\varphi _{\boldsymbol {k}})}} = -1$. This yields a prediction for the Dimits threshold in the $(\kappa _T, \chi)$ plane that we refer to as the ‘fastest-mode approximation’. Note that there is no a priori reason to assume that $\alpha$ is itself not a function of $\kappa _T$.

4.4.1. High-collisionality limit

We can take the $\chi \to \infty$ limit of the fastest-mode approximation analytically. We formally order $\chi \sim \kappa _T$, as suggested by figure 23. We then use the dispersion relation (2.20) to find the growth rate $\gamma _{\boldsymbol {k}}$ and real frequency $\omega _{\boldsymbol {k}}$ of the fastest mode with $k_x =\alpha k_y$. In § 2.6.2, we showed that a mode of wavenumber $\boldsymbol {k}$ is unstable if and only if

(4.16)\begin{equation} \kappa_T k_y^2 > a\chi^2k^6 = a\chi^2(1 + \alpha^2)^3k_y^6, \end{equation}

so all unstable modes with $k_x = \alpha k_y$, where $\alpha \sim {O}\left (1\right)$, satisfy

(4.17)\begin{equation} k_y^4 < \frac{\kappa_T}{a\chi^2(1 + \alpha^2)^3} \sim {O}\left(\kappa_T^{-1}\right). \end{equation}

Similarly, using the results in § 2.6.1 for the FLR bounds on the region of unstable wavenumbers, we find $k_y^4 \sim {O}\left (\kappa _T^{-1}\right)$ in the limit $\kappa _T \to \infty$. Thus, both mechanisms that bound the region of instability (and hence restrict the largest ITG growth rate) lead to the same scaling for the unstable wavenumbers. Therefore, the wavenumber of the most unstable mode must also satisfy $k_y \sim \kappa _T^{-1/4}$. Applying the ordering $\chi \sim \kappa _T \gg 1$ and $k_x = \alpha k_y \sim \kappa _T^{-1/4}$ to (2.28a–e), we find

(4.18a–e)\begin{equation} A \sim {O}\left(\kappa_T^{1/2}\right), \quad B \sim {O}\left(1\right), \quad C \sim {O}\left(\kappa_T^{1/4}\right), \quad f \sim {O}\left(1\right), \quad g \sim {O}\left(1\right). \end{equation}

The unstable solution of (2.27) is

(4.19)\begin{equation} \gamma_{\boldsymbol{k}} - \textrm{i} \omega_{\boldsymbol{k}} = \frac{-A-B+\textrm{i}C + \sqrt{(A-B+\textrm{i}C)^2 + 4fAB - 4\textrm{i}gAC}}{2}. \end{equation}

After expanding it using (4.18a–e), we find

(4.20)\begin{equation} \gamma_{\boldsymbol{k}} \sim {O}\left(1\right) \ll \omega_{\boldsymbol{k}} \sim {O}\left(\kappa_T^{1/4}\right) \ll \chi k^2 \sim {O}\left(\sqrt{\kappa_T}\right), \end{equation}

(4.21)\begin{equation} \omega_{\boldsymbol{k}} = -(1-g)C + {O}\left(\kappa_T^{-1/4}\right) = -\kappa_T (1-b)(1+\alpha^2)k_y^3 + {O}\left(\kappa_T^{-1/4}\right). \end{equation}

Therefore, (2.58) gives

(4.22)\begin{equation} {\text{Re}{\frac{T_{\boldsymbol{k}}}{\varphi_{\boldsymbol{k}}}}} = \frac{k_y\kappa_T\omega_{\boldsymbol{k}}}{|\gamma_{\boldsymbol{k}} -\textrm{i}\omega_{\boldsymbol{k}} + \chi k^2|^2} \approx \frac{k_y\kappa_T\omega_{\boldsymbol{k}}}{\chi^2k^4} \approx -\frac{1-b}{1+\alpha^2} \frac{\kappa_T^2}{\chi^2}. \end{equation}

Thus, the large-temperature-gradient fastest-mode approximation of the Dimits threshold is a straight line in the $(\kappa _T, \chi)$ plane, given by

(4.23)\begin{equation} \kappa_T^c \approx \chi \sqrt{\frac{1 + \alpha^2}{1 - b}}, \end{equation}

a posteriori confirming the ordering $\chi \sim \kappa _T$. The numerically determined Dimits threshold is indeed close to a straight line. Fitting the slope of that line to (4.23) yields $\alpha \approx 3$. Comparison of the prediction for the Dimits threshold for this value of $\alpha$, as well as $\alpha = 0$, which corresponds to the threshold for the secondary instability of a primary streamer (as discussed in § 2.8), can be found in figure 23. The convergence is slow ($\propto \kappa _T^{-1/4}$), hence the sizeable discrepancy for the values of $\kappa _T$ shown there, but we consider the asymptotic result to be sound.

Figure 23.

Comparison of numerical data with the analytical estimate for the threshold of the Dimits regime (§ 4.4). The black points represent the numerically observed $\kappa _T^c$. The other two curves correspond to the parameters for which the fastest-growing mode with $k_x = \alpha k_y$ satisfies ${\text {Re}{(T_{\boldsymbol {k}}/\varphi _{\boldsymbol {k}})}} = -1$. The case $\alpha = 0$ corresponds to the threshold for the suppression of the secondary instability of the fastest streamer ($k_x = 0$) mode (see § 2.8), and the curve with $\alpha = 3$ asymptotes to the numerically determined slope of the Dimits threshold as $\kappa _T \to \infty$. Its asymptote is represented here by the dotted line and is given by (4.23).

4.5. Beyond the Dimits regime

Beyond the Dimits threshold ($\kappa _T > \kappa _T^c$), our two-dimensional system fails to reach saturation on a scale smaller than the domain size – perturbations grow exponentially and the box-sized streamer ($k_x = 0, k_y = 2{\rm \pi} /L_y$) eventually dominates the spectrum. Figures 24 and 25 show that the large-scale, coherent ZFs that comprise the zonal staircase are quickly destroyed and never reappear. This is consistent with the illustration in figure 22 and the discussion in § 4.3. For $\kappa _T > \kappa _T^c$, if a shear zone of coherent zonal shear (like the ones we observe in the Dimits regime) were formed, the turbulent stress $\varPi _t$ would flatten out the ZF profile. Any coherent zonal shear is thus the harbinger of its own demise due to the momentum flux of the tilted turbulent eddies. The nonzonal perturbations grow exponentially, and so do the ZFs, but the latter are now dominated by small-scale time-incoherent zonal modes that are unable to quench the instability.

Figure 24.

Time evolution of the heat flux (a) and ZF velocity $\overline {u}_y$ (b) beyond the Dimits regime. The ZF amplitude in the panel (b) normalised to a maximum of $1$ at each time. The initial conditions are the same as those for SimH, but with an augmented temperature gradient $\kappa _T = 1.21 > \kappa _T^c \approx 1$ and a lower resolution of $337\times 167$ Fourier modes due to the numerical cost of simulating the blow-up regime. The ZFs are quickly destroyed and large-scale ZFs never reappear, while the box-averaged heat flux $Q$ grows exponentially. This state is eventually dominated by a streamer with a poloidal scale equal to that of the integration domain (see figure 25).

Figure 25.

Snapshots of temperature perturbations in the blow-up state beyond the Dimits threshold (supplementary movies, entitled as movie\_fig25.mp4). The amplitudes are normalised to a maximum of $1$ at each time. The initial conditions are the same as those for SimH, but with an augmented temperature gradient $\kappa _T = 1.21 > \kappa _T^c \approx 1$ and a lower resolution of $337\times 167$ Fourier modes due to the numerical cost of simulating the blow-up regime.

The lack of saturation beyond the Dimits regime in two dimensions is not surprising. Gyrokinetic simulations have shown that, beyond the Dimits regime in saturated three-dimensional ITG turbulence, the ITG frequency at the injection (‘outer’) scale of the perpendicular plane is balanced by the parallel propagation time – the turbulence is in ‘critical balance’ (Barnes et al. Reference Barnes, Parra and Schekochihin2011):

(4.24)\begin{equation} \omega^{o} \sim v_{\parallel} k_{\parallel}, \end{equation}

where $v_{\parallel }$ is some appropriate speed of parallel propagation (e.g. the ion thermal speed $v_{ti}$) and $\omega ^{o} \propto k_y$ is the ITG frequency at the outer scale, proportional to $k_y$ by (2.20). In a tokamak, the smallest allowed value of $k_{\parallel }$ is $k_{\parallel } \sim L_{\parallel }^{-1}$, where $L_{\parallel }$ is the parallel connection length of the device. Thus, a parallel length scale is enforced by the magnetic geometry. The poloidal outer scale ($k_y$) then follows by (4.24) and the radial outer scale is enforced by zonal shearing ($k_x \sim k_y$). The two-dimensional approximation can be obtained as the $k_{\parallel } \to 0$ limit of the three-dimensional system. In this case, (4.24) implies that $k_y \to 0$, in agreement with the blow up dominated by the box-sized streamer that we observe beyond the Dimits threshold. Thus, the two-dimensional approximation is fundamentally inadequate as a description of fully developed ITG turbulence.Footnote ¹⁴ However, we have shown that ZF-mediated saturation and the Dimits transition are captured by a two-dimensional model. Of course, it is an outstanding task (left for future work) to confirm that the physics of the two-dimensional Dimits transition remains (qualitatively) valid in three dimensions.