Dimits transition in three-dimensional ion-temperature-gradient turbulence

Abstract

We extend our previous work on the two-dimensional (2-D) Dimits transition in ion-scale turbulence (Ivanov et al., J. Plasma Phys., vol. 86, 2020, 855860502) to include variations along the magnetic field. We consider a three-field fluid model for the perturbations of electrostatic potential, ion temperature, and ion parallel flow in a constant-magnetic-curvature geometry without magnetic shear. It is derived in the cold-ion, long-wavelength asymptotic limit of the gyrokinetic theory. Just as in the 2-D model, a low-transport (Dimits) regime exists and is found to be dominated by a quasistatic staircase-like arrangement of strong zonal flows and zonal temperature. This zonal staircase is formed and maintained by a negative turbulent viscosity for the zonal flows. Unlike the 2-D model, the three-dimensional (3-D) one does not suffer from an unphysical blow up beyond the Dimits threshold where the staircase becomes nonlinearly unstable. Instead, a well-defined finite-amplitude saturated state is established. This qualitative difference between the 2-D and 3-D models is due to the appearance of small-scale ‘parasitic’ modes that exist only if we allow perturbations to vary along the magnetic field lines. These modes extract energy from the large-scale perturbations and provide an effective enhancement of large-scale thermal diffusion, thus aiding the energy transfer from large injection scales to small dissipative ones. We show that in our model, the parasitic modes always favour a zonal-flow-dominated state. In fact, a Dimits state with a zonal staircase is achieved regardless of the strength of the linear drive, provided the system is sufficiently extended along the magnetic field and sufficient parallel resolution is provided.

1 Introduction

In our previous work (Ivanov et al. 2020), we discussed the two-dimensional (2-D) dynamics of ion-scale turbulence driven by the ion-temperature-gradient (ITG) instability in the plane perpendicular to the magnetic field. We identified the fundamental mechanism of the Dimits transition that demarcates saturation dominated by strong coherent zonal flows (ZFs) – the ‘Dimits state’ – and the strongly turbulent regime where no coherent ZFs exist. The turbulent momentum flux of turbulence sheared by ZFs – viz., whether the zonal ‘turbulent viscosity’ was positive or negative – was found to be the key to the demise of the Dimits state.

However, those findings were based on a simplified model (to which we shall here refer as the ‘2-D model’), obtained as an asymptotic, highly collisional limit of ion gyrokinetics (GK), with the additional assumption that the dynamics were 2-D. This assumption cannot be justified asymptotically. In fact, GK studies of tokamak turbulence have revealed that parallel dynamics are linked to turbulence in the perpendicular plane via the ‘critical balance’ between the nonlinear mixing time and the parallel propagation time (Barnes, Parra & Schekochihin 2011).

In this paper, we carry our work over to a more general model that is a true asymptotic limit of the GK equations by relaxing the two-dimensionality assumption to determine whether the three-dimensional (3-D) Dimits transition is governed by the same mechanism as the 2-D one. In the highly collisional limit discussed in Ivanov et al. (2020), we obtain virtually the same equations for the perturbations of ion temperature and electric potential, with the addition of parallel dynamics and of a new equation for the perturbed parallel ion flow. These three equations (to which we refer as the ‘3-D model’) describe both of the classic ITG instabilities: one mediated by compression along the magnetic field, which we shall call the slab-ITG (sITG) instability (Rudakov & Sagdeev 1961; Coppi, Rosenbluth & Sagdeev 1967; Cowley, Kulsrud & Sudan 1991); the other by magnetic curvature, which we shall call the curvature-driven ITG (cITG) instability (Pogutse 1968; Guzdar et al. 1983). Note that we shall consider only the case of zero magnetic shear.

Our numerical results indicate that the Dimits-regime dynamics of the 3-D model are essentially the same as those of the 2-D model. Namely, we find that the Dimits regime is dominated by a quasistatic staircase-like arrangement of strong ZFs that rip and suppress turbulence. This zonal staircase, reminiscent of the so-called $\boldsymbol {E}\times \boldsymbol {B}$ staircase seen in global GK simulations (Dif-Pradalier et al. 2010; Villard et al. 2013; Rath et al. 2016; Dif-Pradalier et al. 2017), slowly decays due to collisional viscosity. This viscous decay results in recurrent turbulent bursts that are triggered by localised travelling structures emerging from the ZF maxima, where they are created by a local (‘tertiary’) instability of the ZF profile. The turbulence that develops during a burst is sheared by the ZFs. Locally, the shear breaks the fundamental parity symmetry of GK turbulence (Parra, Barnes & Peeters 2011; Fox et al. 2017). This gives rise to a radial flux of poloidal momentum whose sign is controlled by the sign of the zonal shear. This momentum flux consists of two parts: the usual Reynolds stress of the $\boldsymbol {E}\times \boldsymbol {B}$ flow, which is known to generate strong ZFs (Diamond et al. 2005); and a diamagnetic contribution, which is found to oppose the Reynolds stress. The distinguishing feature of the Dimits regime is that the Reynolds stress overcomes the diamagnetic one. The zonal staircase is stable to turbulent bursts because ZF-sheared turbulence provides an effective negative viscosity for the ZFs. All of these effects are found to be qualitatively identical between the 2-D and 3-D models.

The Dimits transition to higher turbulent transport occurs when the diamagnetic stress overcomes the Reynolds one, so the effective turbulent viscosity flips its sign and the coherent ZFs that support the Dimits state become nonlinearly unstable. The 2-D model fails to reach finite-amplitude saturation in this state; instead, box-sized exponentially growing streamers emerge (Ivanov et al. 2020). While such a blow up has not been observed in prior gyrokinetic studies of turbulence in a $Z$-pinch (Ricci, Rogers & Dorland 2006; Kobayashi & Rogers 2012), it is not entirely unexpected in a 2-D fluid system. The 3-D fluid system does not suffer from such an unphysical blow up. Instead, a finite-amplitude saturated state without strong ZFs is established. This qualitative difference between the 3-D and 2-D models is due to the appearance of small-scale sITG modes, which exist only in the 3-D model and are primarily driven by the temperature perturbations associated with the large-scale 2-D perturbations (rather than by the equilibrium temperature gradient). These ‘parasitic’ modes extract energy from those large-scale perturbations and transfer it to smaller perpendicular scales where it is dissipated, thus enabling the system to achieve saturation at finite amplitudes. The idea of such parasitic modes is hardly original (see, e.g., Drake, Guzdar & Hassam 1988; Cowley et al. 1991; Rath & Sridhar 1992). We back their existence both by analytical arguments and by numerical results (§ 4.2) and show that their influence on the large-scale perturbations is to provide an effective enhancement to thermal diffusion (§ 4.2.4).

The rest of the paper is organised as follows. In § 2, we discuss the 3-D extension of the 2-D model of Ivanov et al. (2020). Detailed derivations can be found in Appendix A. Section 3 deals with the linear instability of the 3-D model. Then, in § 4, we describe the nonlinear saturated state: § 4.1 is devoted to the 3-D Dimits regime; § 4.2 to the small-scale sITG instability and to its role in both the Dimits and the strongly turbulent state. We summarise and discuss our results in § 5.

2 Collisional, cold-ion $Z$-pinch in three dimensions

2.1 Model equations

The 3-D model can be derived by following Appendix A of Ivanov et al. (2020), with the addition of the 3-D terms worked out in Appendix A of the present paper. We consider a cold-ion plasma in $Z$-pinch magnetic geometry (shown in figure 1) with magnetic scale length $L_B \equiv -\partial _x \ln B$, where the magnetic field points in the $z$ direction, $\boldsymbol {B} = B \hat {\boldsymbol {z}}$, and $x$ and $y$ are the radial and poloidal coordinates, respectively. Here, $z$ is the coordinate around the current line of the $Z$-pinch ($L_B$ times the azimuthal angle). The ITG scale length is defined as $L_T \equiv - \partial _x \ln T_i$, where $T_i$ is the equilibrium ion temperature. We also assume a large-aspect-ratio system, viz., $L_B \gg L_T$.¹

Figure 1. Illustration of the 3-D $Z$-pinch magnetic geometry.

The perturbed electron density $\delta n_e$ is assumed to obey a modified adiabatic response (Dorland & Hammett 1993; Hammett et al. 1993)

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

where $n_e$ is the equilibrium electron density, $\phi$ is the electric potential, $T_e$ is the electron temperature and

(2.2)\begin{equation} \overline{\phi}(x) \equiv \frac{1}{L_yL_z} \int \,{{\rm d}y}\,{\rm d}z \ \phi(x, y, z) \end{equation}

is the zonal (flux-surface) spatial average of the perturbed electric potential $\phi$. We refer to zonally averaged fields as ‘zonal fields’. We also define the nonzonal field $\phi ' \equiv \phi - \overline {\phi }$. Even though, strictly speaking, there are no well-defined flux surfaces in a $Z$-pinch geometry, our aim is to model a tokamak-like system, thus our definition of a flux-surface average (2.2) is an average over both $y$ and $z$. This can be rationalised by the presence either of asymptotically small, but nonzero, magnetic shear (Ivanov et al. 2020), or of asymptotically small irrational rotational transform. Note that neither of these is present in the final form of our equations.

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

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

where $\varphi \equiv Ze\phi / T_i$ is the normalised electric potential, $Ze$ is the ion charge, $T = \delta T / T_i$ is the normalised ion-temperature perturbation, $\tau = T_i / ZT_e$ is the temperature ratio, $\rho _i \equiv v_{\text {th}i} / \varOmega _i$ is the ion gyroradius given in terms of the ion thermal speed $v_{\text {th}i} \equiv \sqrt {2T_i/m_i}$ and the ion gyrofrequency $\varOmega _i \equiv ZeB/m_ic$, $m_i$ is the ion mass, and $\nu _i$ is the ion–ion collision frequency (for an exact definition of $\nu _i$, see Appendix A.1 of Ivanov et al. 2020). The resulting equations are

(2.4)\begin{align} & \frac{\partial}{\partial t} \left( \tau \varphi' - \frac{1}{2} \rho_i^2\nabla_\perp^2 \varphi\right) +\frac{\partial u_\parallel}{\partial z} -\frac{\rho_i v_{\text{th}i}}{L_B} \frac{\partial}{\partial y} \left( \varphi + T \right) + \frac{\rho_i v_{\text{th}i}}{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_{\text{th}i} \bigg( \bigg\{ \varphi, \tau\varphi' - \frac{1}{2} \rho_i^2\nabla_\perp^2 \varphi \bigg\} + \frac{1}{2}\rho_i^2 \boldsymbol{\nabla_\perp} \boldsymbol{\cdot} \left\{ \boldsymbol{\nabla_\perp}\varphi, T \right\} \bigg) \nonumber\\ & \quad ={-} \frac{1}{2} \chi \rho_i^2 \nabla_\perp^4 (a\varphi - bT), \end{align}

(2.5)\begin{align} & \frac{\partial T}{\partial t} + \frac{5}{2} \frac{\partial u_\parallel}{\partial z} + \frac{\rho_i v_{\text{th}i}}{2L_T} \frac{\partial\varphi}{\partial y} + \frac{1}{2} \rho_i v_{\text{th}i} \left\{ \varphi, T \right\} = \chi \nabla_\perp^2 T, \end{align}

(2.6)\begin{align} & \frac{\partial u_\parallel}{\partial t} + \frac{v_{\text{th}i}^2}{2}\frac{\partial(\varphi + T)}{\partial z} + \frac{1}{2} \rho_i v_{\text{th}i} \left\{ \varphi, u_\parallel \right\} = s \chi \nabla_\perp^2 u_\parallel, \end{align}

where the Poisson bracket is defined by

(2.7)\begin{equation} \left\{ f, g \right\} = \hat{\boldsymbol{b}} \boldsymbol{\cdot} \left( \boldsymbol{\nabla}_\perp f \times \boldsymbol{\nabla}_\perp 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 $\boldsymbol {\nabla }_\perp \equiv \partial _x \hat {\boldsymbol {x}} + \partial _y \hat {\boldsymbol {y}}$ denotes the gradient operator in the perpendicular plane. The values of the thermal diffusivity $\chi$ and the numerical constants $a = 9/40$, $b = 67/160$, $s = 9/10$ are determined by the collisional operator, for which we have used the linearised Landau collision integral. We have omitted the magnetic-drift terms in (2.5) and (2.6) because those are an order $L_T / L_B \sim O(k_{\perp }^2\rho _i^2) \ll 1$ smaller than the rest of the terms in their respective equations. The derivations of (2.4) and (2.5) can be found in Ivanov et al. (2020); the equation (2.6) for the evolution of the parallel flow velocity is derived in Appendix A. Note that we are yet to order the (inverse) parallel scale $k_\parallel \sim \partial _z$ and flow velocity $u_\parallel$, so we have kept parallel streaming in all three equations.

Let us discuss briefly the physics of the ‘new’ (compared with the 2-D model) terms in (2.4)–(2.6). The terms $\propto \partial _z u_\parallel$ in (2.4) and (2.5) describe the compressions and rarefactions due to the parallel ion flow. Equation (2.6) has a straightforward interpretation: the parallel flow is driven by the parallel gradient of the pressure $p = \varphi + T$, advected by the $\boldsymbol {E}\times \boldsymbol {B}$ flow $\boldsymbol {V_{E}} = c\hat {\boldsymbol {b}}\times \boldsymbol {\nabla }_\perp \phi / B$, and damped by the collisional viscosity $s\chi$.

We would like to find an ordering for $k_\parallel$ and $u_\parallel$ that allows for both sITG and cITG. The former depends on the presence of the parallel-streaming terms in (2.4) and (2.6). Thus, we require

(2.8)\begin{equation} \omega \tau \varphi \sim k_\parallel u_\parallel, \quad \omega u_\parallel \sim v_{\text{th}i}^2 k_\parallel \varphi \implies \omega^2 \tau \sim v_{\text{th}i}^2 k_\parallel^2, \end{equation}

where $\omega \sim \partial _t$ is the inverse time scale. We want to retain the curvature-driven instability, so we order $\omega \sim \rho _s \varOmega _i / L_B$, where $\varOmega _i$ is the ion gyrofrequency. Then (2.8) implies

(2.9)\begin{equation} k_\parallel \sim L_B^{{-}1}, \quad u_\parallel \sim \tau c_s \varphi, \end{equation}

where $\rho _s \equiv \rho _i / \sqrt {2\tau }$ is the sound radius and $c_s \equiv \rho _s \varOmega _i$ is the sound speed.

We now introduce the following normalisations (consistent with those that we used for our 2-D model):

(2.10)\begin{equation} \left. \begin{gathered} \hat{t} \equiv \frac{2\rho_s\varOmega_i}{L_B}t, \quad \hat{x} \equiv \frac{x}{\rho_s}, \quad \hat{y} \equiv \frac{y}{\rho_s}, \quad \hat{z} \equiv \frac{2z}{L_B}, \\ \widehat{\varphi} \equiv \frac{\tau L_B\varphi}{2\rho_s} = \frac{\tau L_B}{2\rho_s} \frac{Ze\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}, \quad \hat{u} \equiv \frac{u_\parallel}{\rho_s \varOmega_i \tau}, \\ \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{gathered} \right\} \end{equation}

All hatted quantities are ordered as $O(1)$. Dropping hats, we obtain from (2.4)–(2.6) the following equations in normalised units:

(2.11)\begin{gather} \partial_t \left( \varphi' - \nabla_{{\perp}}^2 \varphi\right) + \partial_\parallel {u} - \partial_y \left( \varphi + T \right) + \kappa_T \partial_y \nabla_{{\perp}}^2 \varphi \nonumber\\ +\left\{ \varphi, \varphi' - \nabla_{{\perp}}^2 \varphi \right\} + \boldsymbol{\nabla}_\perp \boldsymbol{\cdot} \left\{ \boldsymbol{\nabla}_\perp \varphi, T \right\} ={-}\chi \nabla_\perp^4 (a\varphi - bT), \end{gather}

(2.12)\begin{gather}\partial_tT +\kappa_T\partial_y \varphi +\left\{ \varphi, T \right\} = \chi \nabla_{{\perp}}^2 T, \end{gather}

(2.13)\begin{gather}\partial_t {u} + \partial_\parallel (\varphi + T) + \left\{ \varphi, {u} \right\} = s \chi \nabla_{{\perp}}^2 {u}, \end{gather}

where (2.12) has lost its parallel-streaming term because it is $O(\tau )$ smaller than the other terms, and we use $\partial _\parallel \equiv \partial _z$. These equations have two independent parameters: the normalised equilibrium temperature gradient, $\kappa _T$; and the normalised collisionality, $\chi$. There are three other parameters, $L_x$, $L_y$, and $L_{\parallel }$, that are the domain sizes in $x$, $y$ (in units of $\rho _s$), and $z$ (in units of $L_B/2$), respectively. We have already seen that the physics of the 2-D model is independent of $L_x$ and $L_y$ (Ivanov et al. 2020), and that will be true for the 3-D model as well, so the interesting one is $L_{\parallel }$. As we shall later see, the saturated state is independent of $L_{\parallel }$ if $L_{\parallel }$ is large enough, but if it is not, it will play a nontrivial role. Even though the $Z$-pinch geometry imposes a natural $L_{\parallel }$, viz., $L_{\parallel } = 4{\rm \pi}$ (dimensionally this is $2{\rm \pi} L_B$), we will not limit ourselves to that. By considering $L_{\parallel }$ as an independent parameter, we are able to model a shearless flux tube with constant magnetic drifts, periodic boundary conditions, and connection length $L_{\parallel }$. Varying $L_{\parallel }$ in our model is akin to varying the connection length $2{\rm \pi} qR$ in toroidal geometry, where $q$ is the safety factor and $R$ is the major radius.

2.2 Conservation laws

The 2-D cold-ion $Z$-pinch system has three nonlinear invariants (Ivanov et al. 2020). One is the gyrokinetic free energy, while the other two result from the so-called ‘general 2-D invariants’ of GK (Schekochihin et al. 2009). The conservation law of free energy for the 3-D equations (2.11)–(2.13) is equivalent (modulo the integration domain) to that of the 2-D equations. It reads

(2.14)\begin{equation} L_xL_yL_{{\parallel}}\partial_t W \equiv \partial_t \int \,{\rm d}^3\boldsymbol{r} \ \frac{1}{2} T^2 ={-}\kappa_T \int \,{\rm d}^3\boldsymbol{r} \ T \partial_y \varphi - \chi \int \,{\rm d}^3\boldsymbol{r} \ \left( \boldsymbol{\nabla}_\perp T \right)^2. \end{equation}

The first term on the right-hand side of (2.14) is proportional to the nondimensionalised radial heat flux

(2.15)\begin{equation} Q ={-}\frac{1}{L_xL_yL_z}\int \,{\rm d}^3\boldsymbol{r} \ T \partial_y \varphi, \end{equation}

whereas the second one is the collisional thermalisation.

Surprisingly, upgrading from two dimensions to three dimensions does not eliminate both of the other two 2-D invariants. One of them survives, and the following conservation law holds even in three dimensions:

(2.16)\begin{align} L_xL_yL_{{\parallel}}\partial_t I & \equiv \partial_t \int \,{\rm d}^3\boldsymbol{r} \ \left[ \frac{1}{2} \left(\varphi' + T'\right)^2 + \frac{1}{2} \overline{T}^2 + \frac{1}{2} \left(\boldsymbol{\nabla}_\perp T + \boldsymbol{\nabla}_\perp \varphi\right)^2 + \frac{1}{2} u^2 \right] \nonumber\\ & ={-}\kappa_T \int \,{\rm d}^3\boldsymbol{r} \ T \partial_y \varphi -\chi \int \,{\rm d}^3\boldsymbol{r} \ \left[ \left(\boldsymbol{\nabla}_\perp \varphi'\right) \boldsymbol{\cdot} \left(\boldsymbol{\nabla}_\perp T\right) + \left( \boldsymbol{\nabla}_\perp T \right)^2 \right. \nonumber\\ & \left.\quad +\,a \left( \nabla_{{\perp}}^2 \varphi \right)^2 +(a + 1 - b) \left( \nabla_{{\perp}}^2 \varphi \right)\left( \nabla_{{\perp}}^2 T \right) + (1-b)\left( \nabla_{{\perp}}^2 T \right)^2 + s\left(\boldsymbol{\nabla}_\perp {u}\right)^2 \right]. \end{align}

As expected, one recovers a corresponding 2-D conservation law by setting $u = 0$ and excluding $z$ from the integration (see § 2.7 of Ivanov et al. 2020).

Later, it will prove useful to discuss the spectra of $W$ and $I$. For this, we write $W = \sum _{\boldsymbol {k}} W_{\boldsymbol {k}}$ and $I = \sum _{\boldsymbol {k}} I_{\boldsymbol {k}}$, where we have defined

(2.17)\begin{gather} W_{\boldsymbol{k}} \equiv \frac{1}{2} |T_{\boldsymbol{k}}|^2, \end{gather}

(2.18)\begin{gather}I_{\boldsymbol{k}} \equiv \frac{1}{2} \left( |\varphi'_{\boldsymbol{k}} + T'_{\boldsymbol{k}}|^2 + |\overline{T}_{k_x}|^2 + k_{{\perp}}^2 |\varphi_{\boldsymbol{k}} + T_{\boldsymbol{k}}|^2 + |{u}_{\boldsymbol{k}}|^2 \right). \end{gather}

Here the $\boldsymbol {k}$ subscript denotes Fourier components, defined for any field $\varphi (\boldsymbol {r})$ as

(2.19)\begin{equation} \varphi(\boldsymbol{r}) = \sum_{\boldsymbol{k}} \varphi_{\boldsymbol{k}} {e}^{{i}\boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{r}}, \end{equation}

and (2.17) and (2.18) follow from Parseval's theorem.

3 Linear ITG instabilities

Equations (2.11)–(2.13) support two distinct types of linear instability, viz., cITG and sITG. The former was studied by Ivanov et al. (2020) and describes the linearly unstable 2-D modes. In order to investigate the stability of the 3-D modes, we drop the nonlinear terms in (2.11)–(2.13) and look for Fourier modes $\varphi, T, {u} \propto \exp \left ( -{i}\omega _{\boldsymbol {k}} t + {i} \boldsymbol {k} \boldsymbol {\cdot } \boldsymbol {r} \right )$, where ${\textrm {Re}{\left (\omega _{\boldsymbol {k}}\right )}}$, ${\textrm {Im}{\left (\omega _{\boldsymbol {k}}\right )}}$, and $\boldsymbol {k} = (k_x, k_y, k_\parallel )$ are the real frequency, growth rate, and wavenumber of the mode, respectively. The dispersion relation can be written as

(3.1)\begin{equation} ({-{i}\omega_{\boldsymbol{k}}} + s {k_\perp^2}) D_{2D} + \frac{k_\parallel^2}{1 + {k_\perp^2}} \left( {-{i}\omega_{\boldsymbol{k}}} + \chi {k_\perp^2} - {i}\kappa_T k_y\right) = 0, \end{equation}

where the 2-D dispersion relation is given by

(3.2)\begin{gather} D_{2D} \equiv ({-{i}\omega_{\boldsymbol{k}}} + A) ({-{i}\omega_{\boldsymbol{k}}} + B - iC) - fAB + igAC = 0, \end{gather}

(3.3)\begin{gather}A = \chi k_{{\perp}}^2, \quad B = \frac{a\chi k_{{\perp}}^4}{1+k_{{\perp}}^2}, \quad C = k_y \frac{1+\kappa_T k_{{\perp}}^2}{1+k_{{\perp}}^2}, \quad f = \frac{\kappa_T k_y^2}{a\chi^2k_{{\perp}}^6}, \quad g = \frac{b \kappa_T k_{{\perp}}^2}{1 + \kappa_T k_{{\perp}}^2}. \end{gather}

An example of the solutions of (3.1) is given in figure 2. It is evident that (3.1) is too complicated for a general analytical solution. Thus, we will limit our discussion here to several important limits.

Figure 2. A visualisation of the linear growth rate, ${\textrm {Im}{\left (\omega _{\boldsymbol {k}}\right )}}$, given by (3.1), for $\kappa _T = 1$ and $\chi = 0.1$. (a) The linear growth rate in the $k_\parallel = 0$ plane. This is the 2-D cITG instability that we dealt with in Ivanov et al. (2020). (b) The linear growth rate in the $k_x = 0$ plane (where it is largest). The solid black lines denote the marginal modes with ${\textrm {Im}{\left (\omega _{\boldsymbol {k}}\right )}} = 0$. The dotted lines outline the region of unstable collisionless ($\chi = 0$), pure-slab ($L_B^{-1} = 0$) modes, given by (3.14).

3.1 Stable waves

Setting $\kappa _T = 0$ and $\chi = 0$ eliminates the linear instability and damping. The dispersion relation (3.1) reduces to

(3.4)\begin{equation} (1+{k_\perp^2})\omega_{\boldsymbol{k}}^2 + k_y \omega_{\boldsymbol{k}} - k_\parallel^2 = 0, \end{equation}

with solutions

(3.5)\begin{equation} \omega_{\boldsymbol{k}} = \frac{-k_y \pm \sqrt{k_y^2 + 4 k_\parallel^2(1+{k_\perp^2})}}{2(1 + {k_\perp^2})}. \end{equation}

In the limit $k_\parallel \ll k_y$ (which, in terms of dimensional wavenumbers, corresponds to $k_\parallel L_B \ll k_y \rho _s$), we find the familiar 2-D drift waves $\omega _{\boldsymbol {k}} = -k_y / (1 + {k_\perp ^2})$ that result from the magnetic drift. The opposite limit, $k_\parallel \gg k_y$, corresponds to equally familiar ion sound waves, modified by ion finite-Larmor-radius (FLR) effects: ${\omega _{\boldsymbol {k}} = k_\parallel / \sqrt {1 + {k_\perp ^2}}}$. We can undo the normalisations (2.10) to verify that this is the usual dispersion for the ion sound waves in terms of the dimensional wavenumbers $k_\parallel$ and $k_{\perp }$:

(3.6)\begin{equation} \omega_{\boldsymbol{k}} ={\pm} \frac{c_s k_\parallel}{\sqrt{1 + {k_\perp^2} \rho_s^2}}. \end{equation}

We now briefly recap the 2-D cITG instability before turning to the $k_\parallel \neq 0$ sITG.

3.2 Curvature-driven ITG modes

3.2.1 Instability in two dimensions

The dispersion relation for the unstable 2-D ($k_\parallel = 0$) modes is (3.2). These modes were studied carefully in Ivanov et al. (2020); let us recap some important points.

The 2-D modes exist at large perpendicular scales, viz., $k_{\perp } < \min \{k_{\perp, \max,\textrm {FLR}}, k_{\perp, \max, \chi }\}$, where the collisionless and collisional cutoffs are given by

(3.7a,b)\begin{equation} k_{{\perp}, \max,\text{FLR}}^2 = \frac{1+2\sqrt{\kappa_T}}{\kappa_T}, \quad k_{{\perp}, \max, \chi}^2 = \sqrt{\frac{\kappa_T}{a \chi^2}}, \end{equation}

respectively. As shown in Ivanov et al. (2020), the Dimits threshold in two dimensions satisfies $\kappa _T \sim \chi$. Here, however, we shall be interested in the strongly driven limit of $\kappa _T \gg \chi$ and $\kappa _T \gg ~1$, for which a saturated state exists only in three dimensions. In this limit, (3.7a) tells us that the cITG modes exist at (and below) wavenumbers

(3.8)\begin{equation} k_{{\perp}} \sim k_{{\perp}, \max,\text{FLR}} \sim \kappa_T^{{-}1/4} \ll 1. \end{equation}

Solving (3.2) shows that these modes also satisfy

(3.9)\begin{equation} {\text{Re}{\left(\omega_{\boldsymbol{k}}\right)}} \sim {\text{Im}{\left(\omega_{\boldsymbol{k}}\right)}} \sim \kappa_T^{1/4}. \end{equation}

3.2.2 $k_\parallel \neq 0$ corrections

Let us now see how $k_\parallel \neq 0$ affects the strongly driven modes at the curvature-driven scales $k_{\perp } \sim \kappa _T^{-1/4} \ll 1$. At these large perpendicular scales, the effects of collisions are negligible, so we may set $\chi = 0$. Note that the scaling $k_{\perp } \sim \kappa _T^{-1/4} \ll 1$ implies $\kappa _T k_{\perp }^2 \sim \sqrt {\kappa _T} \gg 1$, in which case the dispersion (3.1) becomes

(3.10)\begin{equation} \omega_{\boldsymbol{k}}\left[\omega_{\boldsymbol{k}} \left(\omega_{\boldsymbol{k}} + \kappa_Tk_{{\perp}}^2k_y\right) + \kappa_T k_y^2\right] = k_\parallel^2 \left(\omega_{\boldsymbol{k}} + \kappa_T k_y\right), \end{equation}

where the dispersion relation for the curvature-driven 2-D modes is the expression in the square brackets on the left-hand side. Using the results in § 3.2.1, we can estimate that for these modes, the left-hand and right-hand sides of (3.10) satisfy

(3.11a,b)\begin{equation} \omega_{\boldsymbol{k}}\left[\omega_{\boldsymbol{k}} \left(\omega_{\boldsymbol{k}} + \kappa_Tk_{{\perp}}^2k_y\right) + \kappa_T k_y^2\right] \sim \kappa_T^{3/4}, \quad k_\parallel^2 \left(\omega_{\boldsymbol{k}} + \kappa_T k_y\right) \sim k_\parallel^2 \kappa_T^{3/4}. \end{equation}

We can then conclude that for $k_\parallel \ll 1$, the solutions are essentially 2-D, i.e., the dispersion relation (3.10) is well-approximated by (3.2), whereas $k_\parallel \gg 1$ is expected to introduce qualitative changes to the modes. Let us now investigate the $k_\parallel \gg 1$ sITG instability.

3.3 Collisionless sITG modes

Let us investigate the linear instability of (2.11)–(2.13) in the absence of the magnetic- gradient term $-\partial _y (\varphi + T)$ in (2.11). We shall see shortly when this is appropriate. For now, we limit ourselves to the collisionless ($\chi = 0$) regime (see also § 3.4 and Appendix C). Then, (3.1) becomes

(3.12)\begin{equation} \left(\hat{\omega}_{\boldsymbol{k}}^2 - \frac{\hat{k}_{{\parallel}}^2}{1 + k_{{\perp}}^2} \right) \left(\hat{\omega}_{\boldsymbol{k}} + 1\right) = \frac{2k_{{\perp}}^2\hat{\gamma}_{\boldsymbol{k}}^2\hat{\omega}_{\boldsymbol{k}}^2}{1 + k_{{\perp}}^2}, \end{equation}

where we have defined $\omega _{\boldsymbol {k}} \equiv \kappa _T k_y \hat {\omega }_{\boldsymbol {k}}$, $k_\parallel \equiv \kappa _T k_y \hat {k}_{\parallel }$, and $\hat {\gamma }_{\boldsymbol {k}}^2 \equiv 1/2k_{\perp }^2$.² The last of these may seem like an inconvenience now, but will make the following analysis more easily generalisable for our needs in § 4. Since (3.12) is a real cubic in $\hat {\omega }_{\boldsymbol {k}}$, it either has three real solutions, so all linear modes are stable waves, or one real and two complex solutions, in which case one of the complex solutions has a positive imaginary part and thus corresponds to a linearly unstable mode. It can be shown (see Appendix B) that (3.12) has complex solutions if and only if $\hat {\gamma }_{\boldsymbol {k}}^2 > 0$ and $\hat {k}_{\parallel }^2 \in (\hat {k}_{\parallel,-}^2, \hat {k}_{\parallel,+}^2)$, where

(3.13)\begin{equation} \hat{k}_{{\parallel},\pm}^2 = \frac{k_{{\perp}}^4 + 10k_{{\perp}}^2\hat{\gamma}_{\boldsymbol{k}}^2(1+k_{{\perp}}^2)+k_{{\perp}}^2(4-k_{{\perp}}^2\hat{\gamma}_{\boldsymbol{k}}^2) + 2 \pm k_{{\perp}}\hat{\gamma}_{\boldsymbol{k}}\left[4(1+k_{{\perp}}^2) + k_{{\perp}}^2\hat{\gamma}_{\boldsymbol{k}}^2\right]^{3/2}}{2(1+k_{{\perp}}^2)}. \end{equation}

Substituting $\hat {\gamma }_{\boldsymbol {k}}^2 = 1/2k_{\perp }^2$ yields

(3.14)\begin{equation} \hat{k}_{{\parallel},\pm}^2 = \frac{9}{8(1+k_{{\perp}}^2)} \left[ \frac{8k_{{\perp}}^4}{9} + 4k_{{\perp}}^2 + 3 \pm 3 \left(1 + \frac{8k_{{\perp}}^2}{9}\right)^{3/2} \right]. \end{equation}

The marginal modes, i.e., those on the boundary between unstable and oscillatory modes, are given by $\hat {k}_{\parallel }^2 =\hat {k}_{\parallel,\pm }^2$; these are shown in figure 3. We now consider two distinct asymptotic limits of (3.14): $k_{\perp } \ll 1$ and $k_{\perp } \gg 1$.

Figure 3. Linear growth rates for the sITG instability (without the magnetic drift) as a function of parallel ($k_\parallel$) and poloidal ($k_y$) wavenumbers for $k_x = 0$, $\kappa _T = 1$, $\chi = 0$. The growth rate along the $k_\parallel = \kappa _T k_{\perp } k_y$ line converges to $\kappa _T / \sqrt {2} \approx 0.7$ for large $k_\parallel$. The solid black lines are the instability boundary given by (3.14).

3.3.1 Large-scale sITG instability: $k_{\perp } \ll 1$ modes

To lowest order in $k_{\perp } \ll 1$, (3.12) simplifies to

(3.15)\begin{equation} \hat{\omega}_{\boldsymbol{k}}^3 - \hat{k}_{{\parallel}}^2 \hat{\omega}_{\boldsymbol{k}} - \hat{k}_{{\parallel}}^2 = 0. \end{equation}

This is the well-known sITG dispersion relation without FLR effects and in the absence of a density gradient (Cowley et al. 1991). In this limit, the instability boundaries (3.14) become

(3.16a,b)\begin{equation} \hat{k}_{{\parallel},-} = \frac{2}{3\sqrt{3}}k_{{\perp}}^3 + O\left(k_{{\perp}}^5\right), \quad \hat{k}_{{\parallel},+} = \frac{3\sqrt{3}}{2} + \frac{\sqrt{3}}{4}k_{{\perp}}^2 + O\left(k_{{\perp}}^4\right). \end{equation}

For small $\hat {k}_{\parallel }$, the linearly unstable solution of (3.15) is $\hat {\omega }_{\boldsymbol {k}} \approx |\hat {k}_{\parallel }^{2/3}| (-1 + {i}\sqrt {3})/2$. Thus, the linear growth rate for small $\hat {k}_{\parallel }$, or $k_\parallel \ll \kappa _T k_y$, is

(3.17)\begin{equation} {\text{Im}{\left(\omega_{\boldsymbol{k}}\right)}} \approx \frac{\sqrt{3}}{2} \left(\kappa_T k_yk_\parallel^2 \right)^{1/3}. \end{equation}

This is the most widely recognised expression for the sITG growth rate at long wavelengths, however, it is not the fastest-growing mode at $k_y \ll 1$. From (3.16a,b), we know that (3.15) has ${\textrm {Im}{(\hat {\omega }_{\boldsymbol {k}})}} > 0$ solutions up to $\hat {k}_{\parallel } = O(1)$, or $k_\parallel \sim \kappa _T k_y$. The growth rate of these modes evidently satisfies ${\textrm {Im}{(\hat {\omega }_{\boldsymbol {k}})}} = O(1)$, or ${\textrm {Im}{(\omega _{\boldsymbol {k}})}} \sim \kappa _T k_y$.

We are now able to confirm that neglecting the magnetic drift in deriving (3.12), and hence (3.15), was appropriate. As we saw in § 3.2.1, the strongly driven ($\kappa _T \gg 1$) 2-D cITG modes satisfy $k_y \sim \kappa _T^{-1/4}$. At these wavenumbers, the sITG modes exist at scales $k_\parallel \sim \kappa _T^{3/4} \gg 1$ (as expected and assumed) and have a growth rate ${\textrm {Im}{(\omega _{\boldsymbol {k}})}} \sim \kappa _T k_y \sim \kappa _T^{3/4}$, which is asymptotically larger than the growth rate ${\textrm {Im}{(\omega _{\boldsymbol {k}})}}\sim \kappa _T^{1/4}$ of the cITG modes.

3.3.2 Small-scale sITG instability: $k_{\perp } \gg 1$ modes

Expanding (3.13) for $k_{\perp } \gg 1$ and using $\hat {\gamma }_{\boldsymbol {k}}=O\left (k_{\perp }^{-1}\right )$, we find

(3.18)\begin{equation} \hat{k}_{{\parallel},\pm} = k_{{\perp}} (1 \pm 2\hat{\gamma}_{\boldsymbol{k}}) + O\left(k_{{\perp}}^{{-}1}\right). \end{equation}

Therefore, at small perpendicular scales, the sITG is localised at $\hat {k}_{\parallel } = \pm k_{\perp }$, or, equivalently, at

(3.19)\begin{equation} k_\parallel \approx \pm \kappa_T k_y k_{{\perp}}. \end{equation}

In terms of the dimensional wavenumbers, (3.19) tells us that this instability is localised at $k_\parallel L_B/2 \approx \pm \kappa _T k_yk_{\perp } \rho _s^2$, or, equivalently, $k_\parallel L_T \approx \pm k_yk_{\perp }\rho _i^2$. For $\hat {\gamma }_{\boldsymbol {k}}^2 = 1/2k_{\perp }^2$, (3.18) is

(3.20)\begin{equation} \hat{k}_{{\parallel},\pm} = k_{{\perp}} \pm \sqrt{2} + O\left(k_{{\perp}}^{{-}1}\right), \end{equation}

which implies, for the unstable modes,

(3.21)\begin{equation} \left\lvert \hat{k}_{{\parallel}} - k_{{\perp}} \right\rvert = \left\lvert \frac{k_\parallel}{\kappa_T k_y} - k_{{\perp}} \right\rvert < \sqrt{2}. \end{equation}

We can also find the $k_y$ width of the region of instability at fixed $k_\parallel$. Substituting $k_y = k_\parallel / \kappa _T k_{\perp } + \delta k_y$ into (3.21) and expanding for $\delta k_y \ll k_y$, we find

(3.22)\begin{equation} |\delta k_y| < \sqrt{2}\frac{k_{{\perp}} k_y}{|k_{{\perp}}^2 + k_y^2|} \leq \frac{\sqrt{2}}{2}. \end{equation}

To find the growth rate, consider the $k_{\perp } \gg 1$ limit of (3.12) and set $\hat {k}_{\parallel } = \pm k _{\perp } + \delta \hat {k}_{\parallel }$ and $\hat {\omega }_{\boldsymbol {k}} = -1 + \delta \hat {\omega }_{\boldsymbol {k}}$, where $\delta \hat {\omega }_{\boldsymbol {k}} \sim \delta \hat {k}_{\parallel } / k_{\perp } \sim O(k_{\perp }^{-1}) \ll 1$. Keeping terms of order up to $O(k_{\perp }^{-2})$, (3.12) becomes

(3.23)\begin{equation} \left(\delta \hat{\omega}_{\boldsymbol{k}} \pm \frac{\delta \hat{k}_{{\parallel}}}{k_{{\perp}}} \right) \delta \hat{\omega}_{\boldsymbol{k}} + \hat{\gamma}_{\boldsymbol{k}}^2 \approx 0 \implies \delta \hat{\omega}_{\boldsymbol{k}} \approx{-}\frac{\delta \hat{k}_{{\parallel}}}{2k_{{\perp}}} \pm \sqrt{\frac{\delta \hat{k}_{{\parallel}}^2}{4k_{{\perp}}^2} - \hat{\gamma}_{\boldsymbol{k}}^2}. \end{equation}

Thus, in agreement with (3.21), the instability exists only for $|\delta \hat {k}_{\parallel }| < 2k_{\perp }\hat {\gamma }_{\boldsymbol {k}}$, and its growth rate is

(3.24)\begin{equation} {\text{Im}{\left(\hat{\omega}_{\boldsymbol{k}}\right)}} \approx \sqrt{\hat{\gamma}_{\boldsymbol{k}}^2 - \frac{\delta \hat{k}_{{\parallel}}^2}{4k_{{\perp}}^2}}. \end{equation}

The maximum growth rate is then achieved for $\delta \hat {k}_{\parallel } = 0$, i.e., at $k_\parallel =\pm \kappa _T k_yk_{\perp }$, and is given by ${\textrm {Im}{\left (\hat {\omega }_{\boldsymbol {k}}\right )}} \approx \hat {\gamma }_{\boldsymbol {k}}$. Since $\hat {\gamma }_{\boldsymbol {k}}^2 = 1/2k_{\perp }^2$, this is

(3.25)\begin{equation} {\text{Im}{\left(\omega_{\boldsymbol{k}}\right)}} \approx \frac{\kappa_T k_y}{\sqrt{2} k_{{\perp}}} \leq \frac{\kappa_T}{\sqrt{2}}. \end{equation}

The characteristics of the small-scale sITG instability are summarised in figure 3.

Note that for $\kappa _T \gg 1$, the linear growth rate (3.25) of the small-scale sITG modes scales as $O(\kappa _T)$ and, therefore, dominates both the curvature-driven modes (§ 3.2.1) and the large-scale slab modes (§ 3.3.1). This time-scale separation will prove critical for the saturation of strongly driven turbulence (see § 4.2.2).

Finally, an important feature of the $k_{\perp } \gg 1$ sITG modes is the approximate relation $T \approx - \varphi$, or equivalently, $p / \varphi \ll 1$, where $p = \varphi + T$ is the perturbed pressure. Indeed, using (2.12) and (3.23), we find

(3.26)\begin{equation} \frac{T_{\boldsymbol{k}}}{\varphi_{\boldsymbol{k}}} = \frac{\kappa_T k_y}{\omega_{\boldsymbol{k}}} = \frac{1}{\hat{\omega}_{\boldsymbol{k}}} ={-}1 + \frac{\delta \hat{k}_{{\parallel}}}{2 k_{{\perp}}} - {i}\sqrt{\hat{\gamma}_{\boldsymbol{k}}^2 - \frac{\delta \hat{k}_{{\parallel}}^2}{4k_{{\perp}}^2}} + O\left(k_{{\perp}}^{{-}2}\right) \end{equation}

for the modes with ${\textrm {Im}{\left (\hat {\omega }_{\boldsymbol {k}}\right )}} > 0$. Thus, these modes generally have $\,p_{\boldsymbol {k}} / \varphi _{\boldsymbol {k}} \sim O(k_{\perp }^{-1}) \ll 1$, while the most unstable of them ($\delta \hat {k}_{\parallel }=0$) satisfy $\,p_{\boldsymbol {k}} / \varphi _{\boldsymbol {k}} \sim O(k_{\perp }^{-2})$ and ${\textrm {Re}{\left (T_{\boldsymbol {k}} / \varphi _{\boldsymbol {k}}\right )}} = -1 + O(k_{\perp }^{-2})$. This relationship between $T_{\boldsymbol {k}}$ and $\varphi _{\boldsymbol {k}}$ will allow us to identify the sITG modes in the saturated state, and will prove useful in understanding their role in maintaining the Dimits state (see §§ 4.2.1 and 4.2.5).

3.4 Mechanism of the small-scale sITG instability

The analysis in § 3.3.2 is somewhat physically opaque. To get a better grasp of the small-scale sITG modes, we can consider the problem from a slightly different angle. Let us subtract the Laplacian $\nabla _{\perp }^2$ of (2.12) from (2.11) and rewrite the linear part of the system (2.11)–(2.13) as

(3.27)\begin{gather} \partial_t \left( \varphi' - \nabla_{{\perp}}^2 \varphi\right) + \partial_\parallel {u} - \partial_y p + \kappa_T \partial_y \nabla_{{\perp}}^2 \varphi + \chi \nabla_\perp^4 \left[(a+b)\varphi - bp\right] = 0, \end{gather}

(3.28)\begin{gather}-\partial_t \nabla_{{\perp}}^2 p + \partial_\parallel {u} - \partial_y p + \chi \boldsymbol{\nabla}_\perp^4(1-b)p ={-}\partial_t \varphi' - \chi \boldsymbol{\nabla}_\perp^4 (a + b - 1)\varphi, \end{gather}

(3.29)\begin{gather}\partial_t {u} + \partial_\parallel p - s \chi \nabla_{{\perp}}^2 {u} = 0, \end{gather}

where $p = \varphi + T$ is the pressure perturbation. Let us first concentrate on the $\chi = 0$ case, viz.,

(3.30)\begin{gather} \partial_t \left( \varphi' - \nabla_{{\perp}}^2 \varphi\right) + \partial_\parallel {u} - \partial_y p + \kappa_T \partial_y \nabla_{{\perp}}^2 \varphi = 0, \end{gather}

(3.31)\begin{gather}-\partial_t \nabla_{{\perp}}^2 p + \partial_\parallel {u} - \partial _{y} p ={-}\partial_t \varphi', \end{gather}

(3.32)\begin{gather}\partial_t {u} + \partial_\parallel p = 0. \end{gather}

Observe that the term $\partial _t \varphi '$ on the right-hand side of (3.31) is asymptotically small in the $k_{\perp } \gg 1$ limit. Indeed, had we approximated $\partial _t (1+k_{\perp }^2)\varphi \approx \partial _t k_{\perp }^2\varphi$ in (2.11), as we should have done for $k_{\perp } \gg 1$, the right-hand side of (3.31) would have been zero. In this approximation, (3.31) and (3.32) decouple from (3.30). Their dispersion relation coincides with the $k_{\perp } \gg 1$ limit of (3.5), so (3.31) and (3.32) describe two propagating waves, independent of $\kappa _T$. Let us call these two modes ‘pressure waves’.³ The third mode is a $p = {u} = 0$ wave, described by (3.30); its frequency in the $k_{\perp } \gg 1$ limit is $\omega _{\boldsymbol {k}} = -\kappa _T k_y$. We shall call this a ‘diamagnetic wave’ because the restoring force comes from the diamagnetic-drift term $\kappa _T \partial _y\nabla _{\perp }^2\varphi$ in (3.30).

Since the diamagnetic and pressure waves have, in general, disparate frequencies, the small coupling term $-\partial _t \varphi '$ in (3.31) can indeed be neglected. However, if the frequencies of these modes happen to coincide, i.e., if they are in resonance, the small coupling term can no longer be neglected. Using (3.5) for the frequency of the pressure waves and $\omega _{\boldsymbol {k}} = -\kappa _T k_y$ for the diamagnetic wave, we find that such a resonance occurs when $k_\parallel = \kappa _T k_{\perp } k_y$, assuming $k_{\perp } \gg 1$. Thus, the instability condition (3.19) for collisionless small-scale sITG modes is the resonance condition for the two types of linear modes in the system, viz., pressure waves and diamagnetic waves.

Let us now restore $\chi \neq 0$. Then, (3.28) shows that for $a + b \neq 1$ (as is generally the case), there is a second coupling mechanism, via the term $\chi \boldsymbol {\nabla }_\perp ^4 (a + b - 1)\varphi$. For $\omega _{\boldsymbol {k}} \sim \kappa _T k_y$, this term is comparable to the collisionless-coupling term $\partial _t \varphi '$ when $\omega _{\boldsymbol {k}} \sim \kappa _T k_y \sim \chi k_{\perp }^4$, i.e., when

(3.33)\begin{equation} k_{{\perp}} \sim \left(\frac{\kappa_T}{\chi}\right)^{1/3} \equiv k_{\chi}, \end{equation}

assuming $k_{\perp } \sim k_y$. We find that $k_{\chi }$ is the perpendicular scale at which the collisionless results of § 3.3.2 are no longer valid as the effects of finite $\chi$ can no longer be neglected. Naïvely, one might expect that for $k_{\perp } > k_{\chi }$, collisions will act to damp the sITG instability. However, this turns out not to be the case, and, in fact, the coupling term $\chi \boldsymbol {\nabla }_\perp ^4 (a + b - 1)\varphi$ can mediate a new collisional ITG instability ($\chi$ITG) for $k_{\perp } \gtrsim k_{\chi }$ in the absence of the collisionless coupling term $\partial _t \varphi '$. However, it turns out that in order for $\chi$ITG to be non-negligible compared with sITG, very large temperature gradients are required, viz., $\kappa _T/\chi \gtrsim 830$. Numerically, we shall not investigate such large gradients, so the $\chi$ITG instability will not be relevant for us. The detailed treatment of the $\chi$ITG instability has been relegated to Appendix C.

4 Nonlinear states of low and high transport

We now proceed to study the nonlinear saturated state of (2.11)–(2.13). We solve these equations using an enhanced version of the code used in Ivanov et al. (2020), whereby (2.11)–(2.13) are solved using a pseudospectral algorithm in a triply periodic box of dimensions $L_x$, $L_y$, and $L_{\parallel }$. The linear terms are integrated 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 1995; Dorland et al. 2000). As the 3-D model has no dissipation terms that depend on $k_\parallel$, we usually include small (compared with the collisional dissipation) parallel hyperviscosity of the form $\nu k_\parallel ^4$. It is incorporated in the equations by replacing $\partial _t \mapsto \partial _t + \nu k_\parallel ^4$ for all three fields in the model. The value of $\nu$ is typically chosen to give a maximum parallel hyperviscosity of $10\,\%$ of $\chi k_{\perp, \textrm {largest}}^2$, where $k_{\perp, \textrm {largest}}$ is the largest $k_{\perp }$ included the simulation. This form of hyperviscosity effectively subtracts $\nu k_\parallel ^4$ from the growth rate of every mode, but does not alter the linear mode structure, i.e., it does not influence the ratio of Reynolds to diamagnetic stresses given by ${\textrm {Re}{\left (T_{\boldsymbol {k}} / \varphi _{\boldsymbol {k}}\right )}}$ (see §§ 4.1 and 4.2.5). Thus, it dissipates energy without perturbing the saturated state either towards or away from the Dimits regime.

Recall that the 2-D model has two distinct nonlinear states: a Dimits regime, where saturation is achieved with the aid of strong ZFs that quench the cITG instability by shearing the perturbations it produces; and a blow-up regime, where no finite-amplitude saturation is achieved, but amplitudes continue to grow exponentially indefinitely (or at least until numerical efforts become futile). This unphysical blow up is arguably the main limitation of the 2-D model, and there are good reasons to believe that it is a consequence of the $k_\parallel = 0$ restriction (see § 4.5 of Ivanov et al. 2020). This will indeed be corroborated below as we find that the 3-D model is able to saturate for all values of $\kappa _T$ and $\chi$ that we have investigated numerically.

At low collisionality (which can be argued to be the most relevant case, at least for core turbulence, see Ivanov et al. 2020), the Dimits regime of the 3-D model is strikingly similar to its 2-D counterpart. The saturated state is dominated by quasistatic triangular ZFs that break up the radial domain into regions (shear zones) of constant zonal shear, where turbulence is sheared and thus suppressed (see figure 4a–c). Localised patches of turbulence remain present at the turning points of the ZFs, where the zonal shear vanishes.

Figure 4. Example of instantaneous radial profiles of perturbations in the 3-D Dimits state for $L_x = L_y = 80$: (a) ZF; (b) zonal shear; (c) zonal temperature gradient. Example of instantaneous radial profiles in strong turbulence: (d) ZF; (e) zonal shear; (f) zonal temperature gradient. The dotted green lines in (b,e) are the largest linear growth rates for the respective simulations. The dotted orange lines in (c,f) show the value of $\kappa _T$, which is equal to minus the normalised equilibrium temperature gradient. Just as in the 2-D case, the zonal shear in the Dimits state is determined by the largest linear growth rate. Strongly turbulent ZFs do not have regions of coherent shear.

Periodically, when viscosity has eroded the ZFs and their ability to suppress turbulence has diminished, turbulent bursts are triggered. Just as in two dimensions, these bursts are foreshadowed by an instability located at the ZF maxima and by the appearance of localised travelling structures produced by this instability (‘ferdinons’, discovered by van Wyk et al. 2016, 2017 in GK simulations with external flow shear). An example of a turbulent burst in the 3-D model is shown in figure 5. It is visually indistinguishable from a burst in two dimensions when viewed as a cross-section in the $(x, y)$ plane. We shall discuss the 3-D structure of the Dimits regime in detail in § 4.1.

Figure 5. Snapshots of the perturbed nonzonal (a) temperature $T'$, (b) potential $\varphi '$, (c) pressure $p' = \varphi ' + T'$, and (d) parallel velocity ${u}'$ in the 3-D Dimits state. The colour scale is relative to the maximum absolute amplitude in each panel (given in the panels’ titles). We see that ferdinons carry a ${u}$ perturbation, as well as $T$ and $\varphi$ perturbations. A more detailed view of one of the ferdinons is shown in figure 7. These snapshots are from the same simulation as figure 4(a–c).

The crucial qualitative change in physics that allowing 3-D perturbations brings about is the sITG instability. Recall that the collisionless small-scale sITG modes live at wavenumbers up to $k_{\chi } \sim (\kappa _T/\chi )^{1/3}$ (see § 3.4). This is in stark contrast with the behaviour of the 2-D cITG modes whose cutoff wavenumber (3.7a) scales as $k_{\perp, \textrm {2D cutoff}} \sim \kappa _T^{-1/4}$ (see also § 2.6.1. of Ivanov et al. 2020). Moreover, the maximal growth rate of the sITG modes (3.25) scales as ${\textrm {Im}{\left (\omega _{\boldsymbol {k}}\right )}} \sim \kappa _T$, while that of the cITG modes (3.9) satisfies ${\textrm {Im}{\left (\omega _{\boldsymbol {k}, \textrm {2D}}\right )}} \sim \kappa _T^{1/4}$. This implies that there is a natural scale separation between slow, large-scale curvature-driven modes and fast, small-scale sITG modes. Crucially, this scale separation allows small-scale turbulence to be driven both by the equilibrium gradients and by the gradients associated with the large-scale 2-D modes (which are themselves generated by the cITG instability). In fact, as we shall see in § 4.2, the latter type of driving dominates in the saturated state to such an extent that the equilibrium temperature gradient can be turned off for the $k_\parallel \neq 0$ modes and the saturated state remains largely unchanged. In other words, the sITG modes are ‘parasitic’ modes, a type of 3-D ‘secondary’ instability of the 2-D cITG modes.

Most importantly, in the Dimits state, the small-scale instability can be shown always to favour strong, coherent ZFs. It does so in two ways: by providing an effective positive thermal diffusion for the large-scale modes that would otherwise destabilise the ZFs in 2-D (see § 4.2.4); and by generating momentum transport that is beneficial for the ZFs (i.e., a negative turbulent viscosity for the ZF, see § 4.2.5). This makes the 3-D Dimits state much more resilient than the 2-D one. In fact, we find that the 3-D system stays in a Dimits state regardless of the values of the parameters $\kappa _T$ and $\chi$, provided the domain is ‘sufficiently 3-D’, i.e., provided $L_{\parallel }$ is large enough and that our numerical simulations have sufficient parallel resolution to resolve the sITG modes (see § 4.3).

We now recap the physical mechanism that gives rise to the Dimits regime and also discuss any qualitative and quantitative changes that the 3-D physics brings about. Then, in § 4.2, we turn to the small-scale sITG instability and its consequences for the saturated state. Finally, in § 4.3, we examine the circumstances that can prevent the system from establishing a Dimits state and force it into the strongly turbulent regime.

4.1 Dimits regime

4.1.1 The 2-D picture

Recall that the 2-D Dimits transition is a sharp transition from a finite-amplitude saturated state with strong ZFs to a ‘blow-up’ state dominated by ever-growing streamers (Ivanov et al. 2020). The key to understanding this is the equation for the zonal electrostatic potential

(4.1)\begin{equation} \partial_t \overline{\varphi} + \Pi_\varphi + \Pi_T + \Pi_\chi = 0, \end{equation}

where

(4.2)\begin{equation} \Pi_\varphi \equiv{-} \overline{(\partial_x \varphi) (\partial_y \varphi)}, \quad \Pi_T \equiv{-} \overline{(\partial_x \varphi) (\partial_y T)}, \quad \Pi_\chi \equiv{-} \chi \partial_x^2 \left(a \overline{\varphi} - b \overline{T}\right) \end{equation}

are the Reynolds, diamagnetic, and diffusive stresses, respectively. Equation (4.1) describes how the ZFs are generated or eroded by turbulence (via the Reynolds and diamagnetic stresses, depending on their sign) and damped by collisional viscosity. We then consider a region of nearly constant zonal shear (a ‘shear zone’) of radial width $d$ and find that the integral of the total turbulent stress $\Pi _t = \Pi _\varphi + \Pi _T$ over such a region can be written as

(4.3)\begin{equation} \frac{1}{d}\int \,{{\rm d} x} \Pi_t ={-}\sum_{\boldsymbol{k}} k_xk_y |\varphi_{\boldsymbol{k}}|^2 \left[1 + {\text{Re}{\left(\frac{T_{\boldsymbol{k}}}{\varphi_{\boldsymbol{k}}}\right)}} \right]. \end{equation}

Thus, the effect of the mode with wavenumber $\boldsymbol {k}$ on the ZFs depends on the ratio ${\textrm {Re}{\left (T_{\boldsymbol {k}}/\varphi _{\boldsymbol {k}}\right )}}$. Namely, ${\textrm {Re}{\left (T_{\boldsymbol {k}}/\varphi _{\boldsymbol {k}}\right )}} < -1$ implies that the mode will destabilise the ZFs, while ${\textrm {Re}{\left (T_{\boldsymbol {k}}/\varphi _{\boldsymbol {k}}\right )}} > -1$ means that the mode will reinforce the ZFs. This observation is based on the fact that sheared (by the ZFs) turbulence is ‘tilted’ and the sign of $k_x k_y$ is correlated with the sign of the zonal shear in each shear zone. In Ivanov et al. (2020), we derived a simple estimate for the Dimits threshold at large $\kappa _T$ that was based on applying these ideas to the linear modes of the 2-D system. More generally, we argued that the Dimits transition occurred at the threshold of a nonlinear version of the secondary instability – when sheared by ZFs, turbulence either reinforced these flows and thus a Dimits state was maintained (the Reynolds stress won), or it failed to do so (the diamagnetic stress won) and saturation had to be reached via a different route that did not rely on zonal shear. In the 2-D case, no such alternative route for finite-amplitude saturation existed. This description of how a ZF-dominated state was maintained was demonstrated to be accurate by calculating the turbulent viscosity

(4.4)\begin{equation} \nu_t \equiv{-}\frac{ \langle \int_0^{L_x}\,{{\rm d} x} \ \Pi_t S \rangle_{\varDelta t}}{ \langle \int_0^{L_x}\, {{\rm d} x} \ S^2 \rangle_{\varDelta t}} \end{equation}

in numerical simulations with an imposed static ZF profile; here $\langle \dots \rangle _{\varDelta t}$ is a saturated-state time average and $S \equiv \partial _x^2 \overline {\varphi }$ is the zonal shear. Essentially, $\nu _t$ is a measure of the correlation between the turbulent stress $\Pi _t$ and the zonal shear $S$. We found that $\nu _t < 0$ on the Dimits side of the threshold, indicating that sheared turbulence was feeding the ZFs, which were shearing it. Accordingly, we also found $\nu _t > 0$ beyond the threshold, implying that the turbulent stress was actively suppressing the ZFs.

4.1.2 The influence of $L_{\parallel }$ on the Dimits state

Taking the limit $L_{\parallel } \to 0$ effectively restricts our model equations (2.11)–(2.13) to two dimensions, and thus their saturated Dimits state converges to that of the 2-D model. In figure 6(a,b), we show what happens to the turbulent heat flux $Q$ with increasing $L_{\parallel }$ for two cases: far below the 2-D Dimits threshold ($\kappa _T = 0.36$, $\chi = 0.1$), where turbulent bursts dominate the 2-D state; and closer to it ($\kappa _T = 0.8$, $\chi = 0.1$), where the bursts start to overlap in time. As expected, if $L_{\parallel }$ is small enough, we recover the 2-D results. As $L_{\parallel }$ increases, $Q$ converges in a monotonic way to a definite 3-D value that is smaller than the 2-D heat flux. Figure 6(c–f) show that, for larger values of $L_{\parallel }$, the turbulent bursts become more frequent, but shorter in duration and lower in amplitude. There are two effects responsible for this: parallel localisation of turbulence; and the development of the ‘parasitic’ small-scale sITG modes.

Figure 6. (a) Dependence of the time-averaged heat flux $Q$ on the parallel size of the box $L_{\parallel }$ for $\chi = 0.1$, $L_x = L_y = 80$, and $\kappa _T=0.36$. The orange dotted line shows the time-averaged heat flux for the 2-D state ($L_{\parallel } = 0$). (b) Same as (a), but with $\kappa _T = 0.8$. (c–f) Time evolution of the heat flux $Q$ for $\kappa _T = 0.8$, $\chi = 0.1$, $L_x = 80$, $L_y = 80$, and four different values of $L_{\parallel }$ (notated on each panel). As $L_{\parallel }$ increases, the turbulent bursts become more frequent and less violent, and the time-averaged $Q$ drops.

Parallel localisation is inevitable because the turbulent nonzonal modes cannot propagate information infinitely quickly along the field lines. As we increase $L_{\parallel }$ away from 0, we see elongated nonzonal modes that eventually lose the ability to stay coherent along the field lines if $L_{\parallel }$ is large enough. Figure 7 shows that the typical Dimits-state ferdinons are not true 2-D structures and develop a finite parallel extent if the parallel size of the box allows it. This is in contrast with the ZFs, which do stay perfectly coherent along the entire domain regardless of $L_{\parallel }$. This puts the ZFs at an advantage because the turbulent stresses in (4.1) are parallel averages and so a turbulent burst that is localised to a fraction $\varDelta L _{\parallel }/L_{\parallel }$ of the parallel extent of the box has its turbulent stress diminished by a factor of $\varDelta L_{\parallel }/L_{\parallel }$. As we increase $L_{\parallel }$, every such localised burst provides a smaller restoring ‘kick’ to the ZFs and so it takes less time for the ZFs to decay to a level that permits the development of a new burst. The turbulent heat flux $Q$ is also a spatial average of the turbulent fields and it too is diminished for a localised burst. Thus, we expect smaller, more frequent bursts, and this is precisely what is observed. Note that the ability of ZFs to communicate infinitely fast along the field lines is a consequence of the asymptotic limit of small mass ratio and the modified adiabatic electron response (2.1), which is itself due to the assumed infinitely fast parallel electron streaming. Therefore, the inclusion of kinetic electron effects in the equations would lead to qualitative changes for a large enough $L_{\parallel }$. Naturally, this is outside the scope of this work, but is certainly an important consideration for real devices.

Figure 7. Snapshots of the 3-D temperature perturbations associated with a ferdinon. The plots in each row are cross-sections in different planes at the same $t$ taken from simulations that have the same $\kappa _T = 0.36$, $\chi =0.1$, $L_x = L_y = 80$, but (a) $L_{\parallel } = 32$, (b) $L_{\parallel } = 64$, and (c) $L_{\parallel } = 256$. The black dashed lines visualise the intersections of the cross-sectional planes. As we increase $L_{\parallel }$, turbulence loses the ability to stay coherent along the parallel extent of the box and the bursts become localised in $z$.

Secondly, we find small-scale sITG modes that feed off the perpendicular temperature gradients associated with the ferdinons. The presence of this 3-D small-scale ‘parasitic’ instability can be detected via the parallel velocity ${u}$ because the latter is only involved in the 3-D sITG modes and not in the 2-D cITG modes. Figure 8 shows an example of a ferdinon that is ‘infected’ with such small-scale sITG instabilities. As we shall discuss in § 4.2, the small-scale instability leads to an effective increase in thermal diffusion, and thus an increase in the effective damping at large scales that reduces the large-scale temperature perturbations. This additional damping likely contributes to the reduced $Q$ of the 3-D Dimits state. It also enables saturation at finite amplitude when the Dimits state is broken (§ 4.3).

Figure 8. Snapshots of perturbed nonzonal (a) temperature, (b) potential, (c) pressure, and (d) parallel velocity at fixed $z$ in the Dimits state with parameters $\kappa _T = 0.8$, $\chi = 0.1$, $L_x = L_y = 80$, $L_{\parallel } = 1$, parallel hyperviscosity $\nu = 2.4\times 10^{-8}$, and Fourier-space resolution $(n_x, n_y, n_z) = (171, 171, 21)$. The colour scale is relative to the maximum absolute amplitude in each panel (given in the panel's title). Small-scale sITG modes driven by the gradients of the ferdinon are evident in panel (d).

4.2 The parasitic sITG instability and its role in the saturated state

4.2.1 Numerical evidence

Let us now address the small-scale sITG instability. This instability exists only in the 3-D model and is the most important distinction between it and its 2-D counterpart. It is the presence of this instability that enables, in three dimensions and with finite $L_{\parallel }$, the existence of a strongly turbulent saturated state, i.e., one in which there are no strong, coherent ZFs (the zonal profiles of such a state are shown in figure 4d–f). We find that the most distinctive feature of this state is the concentration of pressure perturbations at perpendicular scales that are much larger than the typical (small) scales for the perturbations in $\varphi$ and $T$ (or, to be more precise, the absence of pressure perturbations in the small-scale structure present in $\varphi$ and $T$); this is manifest in figure 9.

Figure 9. Snapshots of perturbed nonzonal (a) temperature, (b) potential, (c) pressure and (d) parallel velocity at fixed $z$ in the strongly turbulent state with parameters $\kappa _T = 3$, $\chi = 0.05$, $L_x = L_y = 80$, $L_{\parallel } = 1$, parallel hyperviscosity $\nu = 1.5\times 10^{-10}$, and Fourier-space resolution $(n_x, n_y, n_z) = (285, 285, 83)$. The colour scale is relative to the maximum absolute amplitude in each panel (given in the panel's title). Time-averaged spectra from the same simulation are shown in figure 10.

Figure 10. Time-averaged spectra (a) $W_{\boldsymbol {k}}$ and (b) $I_{\boldsymbol {k}}$, defined by (2.17) and (2.18), respectively, in the strongly turbulent state with parameters $\kappa _T = 3$, $\chi = 0.05$, $L_x = 80$, $L_y = 80$, and $L_{\parallel } = 1$. The solid black lines demarcate the region of linear instability for $k_x = 0$, and the red dashed line is $k_\parallel = \kappa _T k_{\perp }^2$, where the collisionless modes with largest growth rate reside (see § 3.3.2). We can see that the largest contributions to the two conserved quantities are offset from the region of linear instability. The dotted black line denotes the peak $k_{\perp, I}(k_\parallel )$ of $I_{\boldsymbol {k}}$ at fixed $k_\parallel$. Zonal profiles and cross-sectional snapshots from the same simulation are shown in figures 4 and 9, respectively. The spectra of the saturated state with the same parameters, but with $\kappa _T$ set to $0$ for all $k_\parallel \neq 0$ modes, are given in (c,d). Turning off the equilibrium gradient for the 3-D modes does not alter the spectra noticeably. Snapshots from this modified simulation are shown in figure 11.

In § 3.3.2, we showed that the smallness of the pressure perturbations (compared with the perturbations of the electrostatic potential and temperature) was characteristic of the small-scale ($k_{\perp } \gg 1$) sITG instability, see (3.26). However, the small-scale structure that we see in figure 9 is not produced by the equilibrium-driven instability. In fact, the equilibrium-driven sITG instability is inconsequential in the saturated state. To show this, we ran artificially modified simulations where $\kappa _T$ was set to $0$ for all modes with $k_\parallel \neq 0$ (this is straightforward to do in our spectral code). This removed the equilibrium-driven linear instability from all 3-D ($k_\parallel \neq 0$) modes. Examining the spectra of the two conserved quantities $W_{\boldsymbol {k}}$ and $I_{\boldsymbol {k}}$ (see § 2.2), we see that turning off the equilibrium temperature gradient for the 3-D modes has no noticeable effect on the structure of turbulence (see figure 10). As figure 11 shows, the modified simulations are also visually indistinguishable from the unmodified ones shown in figure 9. The total heat flux $Q$ changes by approximately 20%–30%, likely due to the loss of radial-symmetry breaking for the 3-D modes, which are now free to transport heat in either direction equally, so, on average, they have zero radial heat flux. The nonlinear interactions between the 2-D ($k_\parallel = 0$) modes cannot produce the 3-D modes that we see in the modified simulations. Therefore, these 3-D modes must be produced by a ‘parasitic’ sITG instability of the 2-D fields (into which energy is injected by the equilibrium gradient).

Figure 11. Same as figure 9, but for a modified simulation, i.e., with $\kappa _T$ set to zero for the $k_\parallel \neq 0$ modes. Visually, the saturated state is identical to that shown in figure 9.

Furthermore, the spectra of $W_{\boldsymbol {k}}$ and $I_{\boldsymbol {k}}$ measured in regular simulations are inconsistent with the region of linear instability of the dispersion relation (3.1). Namely, figure 10(a,b) show that $W_{\boldsymbol {k}}$ and $I_{\boldsymbol {k}}$ of the linearly unstable modes of (3.1) are orders of magnitude smaller than the corresponding spectral peaks of the two conserved quantities. We can quantify this by using the turbulent spectra to determine the ‘dominant’ perpendicular scale as a function of the parallel scale $k_\parallel$. As figure 10(a,b) show, this is the scale at which $I_{\boldsymbol {k}}$ peaks and the dependence of $W_{\boldsymbol {k}}$ on $k_{\perp }$ changes from flat to steeply declining. To extract this scale, we define $k_{\perp, I} (k_\parallel )$ as the $k_{\perp }$ that maximises $I_{\boldsymbol {k}}$ at a fixed $k_\parallel$. Figure 12(a) shows that $k_{\perp, I} (k_\parallel )$ lies outside of the region of linear instability of (3.1). Thus, the 3-D structure of the saturated state is not produced by the linear sITG instability driven by the equilibrium gradient.

Figure 12. (a) Comparison of the location of the spectral peak $k_{\perp, I} (k_\parallel )$ of $I_{\boldsymbol {k}}$ (blue line) and the location of the peak of the growth rate of the collisionless linear instability driven by the equilibrium gradient (orange dashed line), given by $k_\parallel = \kappa _T k_{\perp }^2$. The black curve circumscribes the region of linear instability, i.e., all ${\textrm {Im}{\left (\omega _{\boldsymbol {k}}\right )}} \geq 0$ solutions to (3.1) are inside it and outside of it, all solutions satisfy ${\textrm {Im}{\left (\omega _{\boldsymbol {k}}\right )}} \leq 0$. (b) Comparison of the equilibrium temperature gradient $\kappa _T$ and the ‘effective’ temperature gradient $\kappa _T^\textrm {eff}(k_\parallel ) \equiv k_\parallel / k_{\perp, I}^2$. The data is from the same simulation as shown in figure 9. The spectra of this simulation are given in figure 10. This rough estimate of $\kappa _T^\textrm {eff}$ being approximately 5–10 times larger than $\kappa _T$ is consistent with the calculated growth rate of the parasitic small-scale instability (see figure 13).

In § 3.3.2, we showed that the equilibrium-driven sITG instability is localised at $k_\parallel \approx \kappa _T k_{\perp }^2$. A similar relationship holds for $k_{\perp, I} (k_\parallel )$, viz., $k_\parallel \approx \kappa _T^\textrm {eff} k_{\perp, I}^2$, where $\kappa _T^\textrm {eff}$ can be thought of as an effective temperature gradient. Figure 12(b) shows that this $\kappa _T^\textrm {eff}$ is several times larger than the equilibrium temperature gradient. As we shall see shortly, $\kappa _T^\textrm {eff}$ is actually the gradient of the large-scale 2-D temperature perturbations.

4.2.2 Scale-separated equations for cITG and sITG modes

The numerical analysis above leads us to believe that the 3-D structure of the saturated state is a consequence of an instability driven not by the equilibrium gradient $\kappa _T$, but rather by the gradients of the 2-D perturbations. Let us attack on the analytical front. As we discussed at the start of § 4, the 3-D sITG modes are naturally scale-separated from the 2-D cITG modes both in wavenumber and in frequency. We introduce the parallel average

(4.5)\begin{equation} \left\langle f \right\rangle_\parallel{\equiv} \int \frac{dz'}{L_z} f(z'). \end{equation}

This average allows us to split (2.11)–(2.13) into separate equations for the slow 2-D modes governed by the cITG instability at large perpendicular scales ($k_{\perp } \ll 1$), and for the fast 3-D sITG modes, which live at small perpendicular scales ($k_{\perp } \gg 1$).⁴ We define the small-scale perturbations as $\widetilde {f} \equiv f - \left \langle f \right \rangle _\parallel$. The large-scale equations are then

(4.6)\begin{gather} \partial_t \left( \left\langle \varphi' \right\rangle_\parallel{-} \nabla_{{\perp}}^2 \left\langle \varphi \right\rangle_\parallel\right) - \partial_y \left( \left\langle \varphi \right\rangle_\parallel{+} \left\langle T \right\rangle_\parallel \right) + \kappa_T \partial_y \nabla_{{\perp}}^2 \left\langle \varphi \right\rangle_\parallel \nonumber\\ +\left\{ \left\langle \varphi \right\rangle_\parallel, \left\langle \varphi' \right\rangle_\parallel{-} \nabla_{{\perp}}^2 \left\langle \varphi \right\rangle_\parallel \right\} + \boldsymbol{\nabla}_\perp \boldsymbol{\cdot} \left\{ \boldsymbol{\nabla}_\perp \left\langle \varphi \right\rangle_\parallel, \left\langle T \right\rangle_\parallel \right\} + \chi \nabla_\perp^4 \left(a\left\langle \varphi \right\rangle_\parallel{-} b\left\langle T \right\rangle_\parallel\right) \nonumber\\ ={-}\left\langle \left\{ \widetilde{\varphi}, \widetilde{\varphi}' - \nabla_{{\perp}}^2 \widetilde{\varphi} \right\} - \boldsymbol{\nabla}_\perp \boldsymbol{\cdot} \left\{ \boldsymbol{\nabla}_\perp \widetilde{\varphi}, \widetilde{T} \right\} \right\rangle_\parallel, \end{gather}

(4.7)\begin{gather}\partial_t\left\langle T \right\rangle_\parallel{+}\kappa_T\partial_y \left\langle \varphi \right\rangle_\parallel{+}\left\{ \left\langle \varphi \right\rangle_\parallel, \left\langle T \right\rangle_\parallel \right\} - \chi \nabla_{{\perp}}^2 \left\langle T \right\rangle_\parallel = {-}\!\left\langle \left\{ \widetilde{\varphi}, \widetilde{T} \right\} \right\rangle_\parallel, \end{gather}

(4.8)\begin{gather}\partial_t \left\langle {u} \right\rangle_\parallel{+} \left\{ \left\langle \varphi \right\rangle_\parallel, \left\langle {u} \right\rangle_\parallel \right\} - s \chi \nabla_{{\perp}}^2 \left\langle {u} \right\rangle_\parallel = {-}\!\left\langle \left\{ \widetilde{\varphi}, \widetilde{{u}} \right\} \right\rangle_\parallel. \end{gather}

The right-hand sides of (4.6)–(4.8) represent the influence of the 3-D sITG modes on the large-scale fields. Subtracting (4.6)–(4.8) from (2.11)–(2.13), we find the small-scale equations:

(4.9)\begin{gather} \partial_t \left( \widetilde{\varphi}' - \nabla_{{\perp}}^2 \widetilde{\varphi}\right) + \partial_\parallel \widetilde{{u}} - \partial_y (\widetilde{\varphi} + \widetilde{T}) + \kappa_T \partial_y \nabla_{{\perp}}^2 \widetilde{\varphi} + \widetilde{\left\{ \widetilde{\varphi}, \widetilde{\varphi}' - \nabla_{{\perp}}^2 \widetilde{\varphi} \right\} } + \boldsymbol{\nabla}_\perp \boldsymbol{\cdot} \widetilde{\Big\{ \boldsymbol{\nabla}_\perp \widetilde{\varphi}, \widetilde{T} \Big\}} \nonumber\\ + \chi \nabla_\perp^4 \left(a\widetilde{\varphi} - b\widetilde{T}\right) ={-}\left[\Big\{ \left\langle \varphi \right\rangle_\parallel, \widetilde{\varphi}' - \nabla_{{\perp}}^2 \widetilde{\varphi} \Big\} + \Big\{ \widetilde{\varphi}, \left\langle \varphi' \right\rangle_\parallel{-} \nabla_{{\perp}}^2 \left\langle \varphi \right\rangle_\parallel \Big\} \right.\nonumber\\ \left.+ \boldsymbol{\nabla}_\perp \boldsymbol{\cdot} \Big\{ \boldsymbol{\nabla}_\perp \left\langle \varphi \right\rangle_\parallel, \widetilde{T} \Big\} + \boldsymbol{\nabla}_\perp \boldsymbol{\cdot} \Big\{ \boldsymbol{\nabla}_\perp \widetilde{\varphi}, \left\langle T \right\rangle_\parallel \Big\}\right], \end{gather}

(4.10)\begin{gather}\partial_t\widetilde{T} +\kappa_T\partial_y \widetilde{\varphi} +\widetilde{\Big\{ \widetilde{\varphi}, \widetilde{T} \Big\}} - \chi \nabla_{{\perp}}^2 \widetilde{T} ={-}\left[\Big\{ \left\langle \varphi \right\rangle_\parallel, \widetilde{T} \Big\} + \Big\{ \widetilde{\varphi}, \left\langle T \right\rangle_\parallel \Big\}\right], \end{gather}

(4.11)\begin{gather}\partial_t \widetilde{{u}} + \partial_\parallel \left(\widetilde{\varphi} + \widetilde{T}\right) + \widetilde{\left\{ \widetilde{\varphi}, \hat{u} \right\}} - s \chi \nabla_{{\perp}}^2 \hat{u} ={-}\left[\Big\{ \left\langle \varphi \right\rangle_\parallel, \widetilde{{u}} \Big\} + \Big\{ \widetilde{\varphi}, \left\langle {u} \right\rangle_\parallel \Big\}\right]. \end{gather}

In order to simplify the following analysis, we shall assume both temporal and spatial scale separation between (4.6)–(4.8) and (4.9)–(4.11), i.e., that the large-scale fields are constant in time in (4.9)–(4.11) and that the spatial and temporal scales of (4.9)–(4.11) are short compared with the respective scales of (4.6)–(4.8). In particular, we shall assume that the perpendicular scales of the 2-D modes are sufficiently large for the derivatives of their gradients to be ignored. This assumption turns out to be equivalent to $k_{\perp } q_\perp \ll 1$, where $k_{\perp }$ and $q_\perp$ are the typical perpendicular wavenumbers of the 2-D and parasitic modes, respectively. According to (3.8) and (3.33), the linearly unstable modes satisfy $k_{\perp } \sim \kappa _T^{-1/4}$ and $q_\perp \sim (\kappa _T/\chi )^{1/3}$ in the limit $\kappa _T \gg \chi$. Therefore, the condition $k_{\perp } q_\perp \ll 1$ is equivalent to $\chi \gg \kappa _T^{1/4}$. Additionally, recall that the Dimits threshold in two dimensions is found at $\kappa _T \sim \chi$ (Ivanov et al. 2020) and that, as we showed in § 4.1, the 2-D Dimits regime is qualitatively unchanged when we include 3-D effects. Thus, for the remainder of this section, we shall consider the limit

(4.12)\begin{equation} \kappa_T \gg \chi \gg \kappa_T^{1/4} \gg 1, \end{equation}

which puts us beyond the 2-D Dimits transition (i.e., in two dimensions, such a state blows up). Importantly, we limit our analysis to $\kappa _T/\chi \ll 830$, in which case the $\chi$ITG instability can safely be neglected (see Appendix C). The limit (4.12) then allows us to simplify (4.6)–(4.8) and (4.9)–(4.11) significantly, and thus to describe the interplay between 2-D and parasitic modes analytically. These analytical results agree with our simulations, even though the latter do not strictly conform to (4.12).

4.2.3 Parasitic sITG instability

First, we investigate the small-scale sITG instability in the presence of large-scale 2-D modes. Linearising (4.9)–(4.11) in the limit (4.12), we obtain

(4.13)\begin{gather} \left(\partial_t + \left\langle \boldsymbol{V_{E}} \right\rangle_\parallel \boldsymbol{\cdot} \boldsymbol{\nabla}_\perp \right) \left( \widetilde{\varphi}' - \nabla_{{\perp}}^2 \widetilde{\varphi}\right) + \partial_\parallel \widetilde{{u}} - \partial_y (\widetilde{\varphi} + \widetilde{T}) \nonumber\\ \quad +\, \boldsymbol{\kappa}_T \boldsymbol{\cdot} \boldsymbol{\nabla}_\perp \nabla_{{\perp}}^2 \widetilde{\varphi} + \boldsymbol{\kappa}_n \boldsymbol{\cdot} \boldsymbol{\nabla}_\perp \widetilde{\varphi} ={-} \chi \nabla_\perp^4 (a\widetilde{\varphi} - b\widetilde{T}), \end{gather}

(4.14)\begin{gather}\left(\partial_t + \left\langle \boldsymbol{V_{E}} \right\rangle_\parallel \boldsymbol{\cdot} \boldsymbol{\nabla}_\perp \right)\widetilde{T} +\boldsymbol{\kappa}_T \boldsymbol{\cdot} \boldsymbol{\nabla}_\perp \widetilde{\varphi} = \chi \nabla_{{\perp}}^2 \widetilde{T}, \end{gather}

(4.15)\begin{gather}\left(\partial_t + \left\langle \boldsymbol{V_{E}} \right\rangle_\parallel \boldsymbol{\cdot} \boldsymbol{\nabla}_\perp \right) \widetilde{{u}} + \partial_\parallel \left(\widetilde{\varphi} + \widetilde{T}\right) = s \nabla_{{\perp}}^2 \widetilde{{u}}, \end{gather}

where the ‘local-equilibrium’ quantities

(4.16)\begin{equation} \left\langle \boldsymbol{V_{E}} \right\rangle_\parallel \equiv \hat{\boldsymbol{z}} \times \boldsymbol{\nabla}_\perp \left\langle \varphi \right\rangle_\parallel, \quad \boldsymbol{\kappa}_n \equiv{-} \hat{\boldsymbol{z}} \times \boldsymbol{\nabla}_\perp \left\langle \varphi' \right\rangle_\parallel, \quad \boldsymbol{\kappa}_T \equiv \kappa_T \hat{\boldsymbol{y}} - \hat{\boldsymbol{z}} \times \boldsymbol{\nabla}_\perp \left\langle T \right\rangle_\parallel \end{equation}

are the $\boldsymbol {E}\times \boldsymbol {B}$ advecting flow, the local density gradient, and the total local temperature gradient (large-scale perturbation plus equilibrium), respectively. We assume that $|\boldsymbol {\kappa }_T|\sim |\boldsymbol {\kappa }_n|$ (see § 4.2.4). Note that only the nonzonal electrostatic potential $\varphi '$ gives rise to a density perturbation – this is a consequence of the modified adiabatic electron response (2.1). Note also that we have ignored the large-scale 2-D parallel flow $\left \langle {u} \right \rangle _\parallel$. Since $\left \langle {u} \right \rangle _\parallel$ is not involved in any linear instability, the only way it could be driven is via the small-scale response, viz., the right-hand side of (4.8). In Appendix E, we show that a small initial $\left \langle {u} \right \rangle _\parallel$ will decay under the influence of growing small-scale modes. Accordingly, in our numerical simulations, we find that $\left \langle u \right \rangle _\parallel$ is many orders of magnitude smaller than the other two 2-D fields and is irrelevant for the saturated state.

Ignoring collisions (i.e., setting $\chi = 0$) and taking the gradients of the large-scale fields to be constant over the small scales at which (4.13)–(4.15) hold, we can investigate the small-scale linear instability in a way analogous to what we did in § 3. In particular, we shall focus on the $k_\parallel \sim k_{\perp }^2 \gg 1$ regime analysed in § 3.3.2. We look for Doppler-shifted solutions to (4.13)–(4.15) of the form $\widetilde {\varphi }_{\boldsymbol {k}}, \widetilde {T}_{\boldsymbol {k}}, \widetilde {{u}}_{\boldsymbol {k}} \propto \exp {[-{i}(\omega _{\boldsymbol {k}} + \left \langle \boldsymbol {V_{E}} \right \rangle _\parallel\ \boldsymbol{\cdot }\ \boldsymbol {k} ) t + {i}\boldsymbol {k} \boldsymbol{\cdot} \boldsymbol {r}]}$. Note that we ignore the shear in the $\boldsymbol {E}\times \boldsymbol {B}$ flow $\left \langle \boldsymbol {V_{E}} \right \rangle _\parallel$. We also ignore the magnetic-drift term $-\partial _y(\widetilde {\varphi } + \widetilde {T})$ in (4.13) because it is subdominant for the sITG modes with $k_{\perp } \gg 1$. The resulting dispersion relation for these modes is

(4.17)\begin{equation} \left(\omega_{\boldsymbol{k}}^2 - \frac{k_\parallel^2}{1 + k_{{\perp}}^2}\right) \left( \omega_{\boldsymbol{k}} + \boldsymbol{\kappa}_T \boldsymbol{\cdot} \boldsymbol{k} \right) = \frac{\omega_{\boldsymbol{k}}^2}{1 + k_{{\perp}}^2} \left[ \left(\boldsymbol{\kappa}_n + \boldsymbol{\kappa}_T\right) \boldsymbol{\cdot} \boldsymbol{k} \right]. \end{equation}

Since (4.13)–(4.15) describe real fields, (4.17) must be invariant under $\boldsymbol {k} \mapsto -\boldsymbol {k}$ and $\omega _{\boldsymbol {k}} \mapsto -\omega _{\boldsymbol {k}}^*$. We may, therefore, assume that $\boldsymbol {\kappa }_T \boldsymbol{\cdot} \boldsymbol {k} > 0$ without loss of generality. Repeating the arguments of § 3.3.2, we define $\omega _{\boldsymbol {k}} \equiv \hat {\omega }_{\boldsymbol {k}}\boldsymbol {\kappa }_T \boldsymbol{\cdot} \boldsymbol {k}$ and $k_\parallel \equiv \hat {k}_{\parallel } \boldsymbol {\kappa }_T \boldsymbol{\cdot} \boldsymbol {k}$. Then (4.17) turns out to be formally the same as our old dispersion relation (3.12), but now with

(4.18)\begin{equation} \hat{\gamma}_{\boldsymbol{k}}^2 = \frac{(\boldsymbol{\kappa}_n+\boldsymbol{\kappa}_T) \boldsymbol{\cdot} \boldsymbol{k}}{2k_{{\perp}}^2\boldsymbol{\kappa}_T \boldsymbol{\cdot} \boldsymbol{k}}. \end{equation}

Thus, the results of § 3.3.2 carry over to the parasitic instability described by (4.17). In particular, the sITG instability exists if $\hat {\gamma }_{\boldsymbol {k}}^2 > 0$, i.e., if $(\boldsymbol {\kappa }_n+\boldsymbol {\kappa }_T) \boldsymbol{\cdot} \boldsymbol {k}$ and $\boldsymbol {\kappa }_T \boldsymbol{\cdot} \boldsymbol {k}$ have the same sign, and is localised to $k_\parallel \approx \pm \boldsymbol {\kappa }_T \boldsymbol{\cdot} \boldsymbol {k} k_{\perp }$. Its growth rate is given by

(4.19)\begin{equation} {\text{Im}{\left(\omega_{\boldsymbol{k}}\right)}} \approx \text{Re} \sqrt{\frac{\boldsymbol{\kappa}_T \boldsymbol{\cdot} \hat{\boldsymbol{k}} \left(\boldsymbol{\kappa}_n + \boldsymbol{\kappa}_T\right) \boldsymbol{\cdot} \hat{\boldsymbol{k}}}{2}}, \end{equation}

where $\hat {\boldsymbol {k}} = \boldsymbol {k} / k_{\perp }$. As expected, this is the same as (3.25) if $\boldsymbol {\kappa }_n = 0$ and $\boldsymbol {\kappa }_T = \kappa _T \hat {\boldsymbol {y}}$. In figure 13(b), we show the maximum growth rate obtained from the numerical solution of the full (with collisionality and magnetic curvature turned back on) dispersion relation (3.1) with the addition of the local temperature and density gradients of the large-scale fields. As expected from the numerical analysis in § 4.2.1, the small-scale instability driven by the large-scale gradients is significantly (${\approx }5$ times in this case) stronger than the equilibrium-driven instability. This is consistent with the estimate of the effective temperature gradient $\kappa _T^\textrm {eff}$ for the sITG instability that we showed in figure 12(b).

Figure 13. (a) Snapshot of the 2-D temperature perturbation $\left \langle T \right \rangle _\parallel$ in the $(x, y)$ plane. The data is taken from the same $\kappa _T = 3$, $\chi = 0.05$ simulation that we showed in figure 9. The 2-D temperature perturbations lack the small-scale structure that was seen in figure 9(a), confirming that the parallel average (4.5) removes small-scale perpendicular structure. (b) Small-scale growth rate in the $(x, y)$ plane. This plot is obtained by finding the maximum growth rate of the full (including collisionality and magnetic curvature) dispersion relation (3.1) with the addition of the local temperature and density gradients of the large-scale fields at every point. For this simulation, $\kappa _T = 3$, and so the largest collisionless growth rate, given by (3.25), is $\kappa _T / \sqrt {2} \approx 2.1$. It is thus evident that the influence of the gradients of the large-scale fields dominates over that of the equilibrium gradient $\kappa _T$ by a factor of $5$. The ‘effective’ $\kappa _T^\textrm {eff}$ that we estimated for the same simulation in figure 12(b) is, indeed, a factor of 5–10 larger that the equilibrium gradient $\kappa _T$.

Note that if $\boldsymbol {\kappa }_n \neq 0$, (4.19) implies that modes with $\boldsymbol {\kappa }_T \boldsymbol{\cdot} \boldsymbol {k} \left (\boldsymbol {\kappa }_n + \boldsymbol {\kappa }_T\right ) \boldsymbol{\cdot} \boldsymbol {k} < 0$ are linearly stable. Is there a $\boldsymbol {\kappa }_n$ that quenches the sITG instability for all $\boldsymbol {k}$? Suppose $\boldsymbol {\kappa }_n\nparallel \boldsymbol {\kappa }_T$; then we can choose $\hat {\boldsymbol {k}} \boldsymbol{\cdot} \boldsymbol {\kappa }_n=0$, but $\hat {\boldsymbol {k}} \boldsymbol{\cdot} \boldsymbol {\kappa }_T\neq 0$. By (4.19), any such $\hat {\boldsymbol {k}}$ is an unstable mode. Therefore, to stabilise all modes, we require $\boldsymbol {\kappa }_n\parallel \boldsymbol {\kappa }_T$. In this case, it is evident that $\boldsymbol {\kappa }_n \boldsymbol{\cdot} \hat {\boldsymbol {k}} = (\boldsymbol {\kappa }_n \boldsymbol{\cdot} \boldsymbol {\kappa }_T)(\boldsymbol {\kappa }_T \boldsymbol{\cdot} \hat {\boldsymbol {k}})/|\boldsymbol {\kappa }_T|^2$. Therefore, in order to quench the sITG instability for all $\boldsymbol {k}$, we need

(4.20)\begin{equation} \boldsymbol{\kappa}_n \parallel \boldsymbol{\kappa}_T, \quad \frac{\boldsymbol{\kappa}_n \boldsymbol{\cdot} \boldsymbol{\kappa}_T}{|\boldsymbol{\kappa}_T|^2} \leq{-}1. \end{equation}

We shall now show that the effect of the growing small-scale modes on the large-scale 2-D fields can be expressed as an enhanced thermal diffusivity for the latter.

4.2.4 Anomalous heat flux due to parasitic sITG modes

We expect that the growth of small-scale sITG modes, which are driven by the gradients associated with the large-scale fluctuations, will check the growth of the amplitudes of the driving large-scale fields. This is an intuitive consequence of the conservation laws (2.14) and (2.16). As the parasitic instability is driven by the nonlinear terms that conserve $W=\sum _{\boldsymbol {k}} W_{\boldsymbol {k}}$ and $I =\sum _{\boldsymbol {k}} I_{\boldsymbol {k}}$, an excitation of parasitic small-scale modes should show up as a sink in the large-scale equations. Let us now calculate explicitly the influence of small-scale sITG modes on the large-scale modes and show that this is indeed true. This influence is represented by the terms of the form $\left \langle \left \{ ., . \right \} \right \rangle _\parallel$ on the right-hand sides of (4.6)–(4.8).

First, consider the temperature equation (4.7). The relevant term is

(4.21)\begin{equation} -\Big\langle \Big\{ \widetilde{\varphi}, \widetilde{T} \Big\} \Big\rangle_\parallel = {-}\boldsymbol{\nabla}_\perp \boldsymbol{\cdot} \Big\langle \Big[ (\hat{\boldsymbol{z}}\times\boldsymbol{\nabla}_\perp\widetilde{\varphi}) \widetilde{T} \Big] \Big\rangle_\parallel{\equiv} -\boldsymbol{\nabla}_\perp \boldsymbol{\cdot} \langle\widetilde{\boldsymbol{Q}}\rangle_\parallel, \end{equation}

where $\langle \widetilde {\boldsymbol {Q}}\rangle _\parallel$ is the turbulent heat flux associated with the small-scale modes. Let us compute it quasilinearly (i.e., assuming that $\langle \widetilde {\boldsymbol {Q}}\rangle _\parallel$ is determined by the most unstable small-scale modes), assuming scale separation. As stated in § 4.2.2, we imagine that the small-scale equations are solved in an infinitesimal (compared with the large scales) box, and thus the parallel average is equivalent to an average over such a small-scale box,

(4.22)\begin{equation} \langle\widetilde{\boldsymbol{Q}}\rangle_\parallel = {-}\sum_{\boldsymbol{q}} {i} \hat{\boldsymbol{z}}\times\boldsymbol{q}\widetilde{\varphi}_{-\boldsymbol{q}} \widetilde{T}_{\boldsymbol{q}} \approx{-}\sum_{\boldsymbol{q}} \hat{\boldsymbol{z}}\times\hat{\boldsymbol{q}} \sqrt{\frac{(\boldsymbol{\kappa}_T + \boldsymbol{\kappa}_n) \boldsymbol{\cdot} \hat{\boldsymbol{q}}}{2\boldsymbol{\kappa}_T \boldsymbol{\cdot} \hat{\boldsymbol{q}}}} |\widetilde{\varphi}_{\boldsymbol{q}}|^2, \end{equation}

where $\hat {\boldsymbol {q}} = \boldsymbol {q} / q_\perp$ and we have assumed that the sum is dominated by the wavenumbers $\boldsymbol {q}$ corresponding to the largest linear growth rate of the parasitic sITG instability, and so have replaced $\widetilde {T}_{\boldsymbol {q}}/\widetilde {\varphi }_{\boldsymbol {q}}$ with the collisionless expression (3.26) for the modes with $k_\parallel = \boldsymbol {\kappa }_T \boldsymbol{\cdot} \boldsymbol {k}k_{\perp }$ that maximise this growth rate. Note that the small-scale fields $\widetilde {T}_{\boldsymbol {q}}$ and $\widetilde {\varphi }_{\boldsymbol {q}}$, and thus $\langle \widetilde {\boldsymbol {Q}}\rangle _\parallel$ itself, depend implicitly on the position variable of the large-scale equations (4.6)–(4.8).

In order to verify that $\langle \widetilde {\boldsymbol {Q}}\rangle _\parallel$ does indeed damp the large-scale temperature perturbations $\left \langle T \right \rangle _\parallel$, we multiply (4.7) by $\left \langle T \right \rangle _\parallel$ and integrate over space to find

(4.23)\begin{align} & \partial_t \int \,{\rm d}^3\boldsymbol{r} \frac{1}{2} \left\langle T \right\rangle_\parallel^2 + \text{linear terms} = \int \,{\rm d}^3\boldsymbol{r} \langle\widetilde{\boldsymbol{Q}}\rangle_\parallel \boldsymbol{\cdot} \boldsymbol{\nabla}_\perp \left\langle T \right\rangle_\parallel \nonumber\\ & \quad\approx{-} \int \,{\rm d}^3\boldsymbol{r} \sum_{\boldsymbol{q}} (\hat{\boldsymbol{z}}\times\hat{\boldsymbol{q}}) \boldsymbol{\cdot} \boldsymbol{\nabla}_\perp \left\langle T \right\rangle_\parallel \sqrt{\frac{(\boldsymbol{\kappa}_T + \boldsymbol{\kappa}_n) \boldsymbol{\cdot} \hat{\boldsymbol{q}}}{2\boldsymbol{\kappa}_T \boldsymbol{\cdot} \hat{\boldsymbol{q}}}} |\widetilde{\varphi}_{\boldsymbol{q}}|^2 \nonumber\\ & \quad={-} \int \,{\rm d}^3\boldsymbol{r} \sum_{\boldsymbol{q}} \sqrt{\frac{\boldsymbol{\kappa}_T \boldsymbol{\cdot} \hat{\boldsymbol{q}}(\boldsymbol{\kappa}_T + \boldsymbol{\kappa}_n) \boldsymbol{\cdot} \hat{\boldsymbol{q}}}{2}} |\widetilde{\varphi}_{\boldsymbol{q}}|^2 ={-} \int \,{\rm d}^3\boldsymbol{r} \sum_{\boldsymbol{q}} {\text{Im}{\left(\omega_{\boldsymbol{k}}\right)}} |\widetilde{\varphi}_{\boldsymbol{q}}|^2, \end{align}

where ${\textrm {Im}{\left (\omega _{\boldsymbol {k}}\right )}}$ is the sITG growth rate (4.19). Thus, the linearly unstable small-scale modes have a sign-definite effect on $\left \langle T \right \rangle _\parallel$: they provide additional dissipation.

The heat flux (4.22) depends on $\left \langle T \right \rangle _\parallel$ in a nontrivial way. Let us quantify its influence on $\left \langle T \right \rangle _\parallel$ by working out its direction as a function of $\boldsymbol {\kappa }_T$. Let us assume that $\langle \widetilde {\boldsymbol {Q}}\rangle _\parallel$ is dominated by the fastest-growing sITG modes, and let their wavevector direction be $\hat {\boldsymbol {q}}_\textrm {max}$, so $\langle \widetilde {\boldsymbol {Q}}\rangle _\parallel$ is parallel to $\hat {\boldsymbol {z}}\times \hat {\boldsymbol {q}}_\textrm {max}$. In figure 14, we illustrate the influence on $\boldsymbol {\kappa }_T$ of the contribution to $\langle \widetilde {\boldsymbol {Q}}\rangle _\parallel$ from the most unstable small-scale modes. As expected, we find that the turbulent heat flux due to the small-scale modes pushes the large-scale gradient $\boldsymbol {\kappa }_T$ towards the linearly stable configuration (4.20).

Figure 14. This plot shows the direction in which the heat flux (4.22) of the most unstable small-scale mode ($\hat {\boldsymbol {q}} = \hat {\boldsymbol {q}}_{\textrm {max}}$) pushes the temperature gradient $\boldsymbol {\kappa }_T = -\hat {\boldsymbol {z}}\times \boldsymbol {\nabla }_\perp \left \langle T \right \rangle _\parallel$. We have chosen a coordinate system in which the large-scale density gradient is $\boldsymbol {\kappa }_n = (0, 1)$, denoted by the green arrow. The red line shows the values of $\boldsymbol {\kappa }_T$ for which the sITG instability has zero growth rate, according to (4.20). The black arrows represent the direction of $-\hat {\boldsymbol {z}}\times \langle \widetilde {\boldsymbol {Q}}\rangle _\parallel$. We see that $\langle \widetilde {\boldsymbol {Q}}\rangle _\parallel$ pushes the large-scale temperature gradient $\boldsymbol {\kappa }_T$ towards the linearly stable region.

Now consider (4.6), the evolution equation for $\left \langle \varphi \right \rangle _\parallel$. The relevant nonlinear terms are

(4.24)\begin{align} & \left\langle \left\{ \widetilde{\varphi}, \widetilde{\varphi}' - \nabla_{{\perp}}^2 \widetilde{\varphi} \right\} + \boldsymbol{\nabla}_\perp \boldsymbol{\cdot} \left\{ \boldsymbol{\nabla}_\perp \widetilde{\varphi}, \widetilde{T} \right\} \right\rangle_\parallel = \boldsymbol{\nabla}_\perp \boldsymbol{\cdot} \left\langle \left\{ \boldsymbol{\nabla}_\perp \widetilde{\varphi}, \widetilde{p} \right\} \right\rangle_\parallel \nonumber\\ & \quad= \boldsymbol{\nabla}_\perp \boldsymbol{\nabla}_\perp \boldsymbol{:} \sum_{\boldsymbol{q}} (\hat{\boldsymbol{z}}\times\boldsymbol{q})\boldsymbol{q} \left(1 + {\text{Re}{\frac{\widetilde{T}_{\boldsymbol{q}}}{\widetilde{\varphi}_{\boldsymbol{q}}}}}\right) |\widetilde{\varphi}_{\boldsymbol{q}}|^2 \equiv \boldsymbol{\nabla}_\perp \boldsymbol{\nabla}_\perp \boldsymbol{:} \langle\widetilde{\Pi}\rangle_\parallel. \end{align}

The collisionless calculations of § 3.3.2 are straightforward to generalise for the collisionless parasitic small-scale instability. They yield the same relation for ${\textrm {Re}{(\widetilde {T}_{\boldsymbol {q}} / \widetilde {\varphi }_{\boldsymbol {q}})}}$, viz., $1 + {\textrm {Re}{(\widetilde {T}_{\boldsymbol {q}} / \widetilde {\varphi }_{\boldsymbol {q}})}} = O(q_\perp ^{-2})$. However, as we shall discuss in § 4.2.5, the presence of nonzero $\chi$ alters this to $1 + {\textrm {Re}{(\widetilde {T}_{\boldsymbol {q}} / \widetilde {\varphi }_{\boldsymbol {q}})}} = O(q_\perp ^{-1})$. Assuming therefore that the dominant parasitic modes satisfy $1 + {\textrm {Re}{(\widetilde {T}_{\boldsymbol {q}}/\widetilde {\varphi }_{\boldsymbol {q}})}} \lesssim O(q_\perp ^{-1})$, we find $\boldsymbol {\nabla }_\perp \boldsymbol {\nabla }_\perp \boldsymbol {:} \langle \widetilde {\Pi }\rangle _\parallel \lesssim k_{\perp }^2 q_\perp |\hat {\varphi }|^2$. However, (4.22) implies $\boldsymbol {\nabla }_\perp \boldsymbol{\cdot} \langle \widetilde {\boldsymbol {Q}}\rangle _\parallel \sim k_{\perp } |\hat {\varphi }|^2$, and so

(4.25)\begin{equation} \frac{\boldsymbol{\nabla}_\perp \boldsymbol{\nabla}_\perp \boldsymbol{:} \langle\widetilde{\Pi}\rangle_\parallel}{\boldsymbol{\nabla}_\perp \boldsymbol{\cdot} \langle\widetilde{\boldsymbol{Q}}\rangle_\parallel} \lesssim k_\perp q_\perp\sim O\left[\left(\frac{\kappa_T^{1/4}}{\chi}\right)^{1/3}\right] \ll 1, \end{equation}

in line with the assumption on scales formulated at the end of § 4.2.2. Therefore, assuming that $\left \langle \varphi ' \right \rangle _\parallel \sim \left \langle T \right \rangle _\parallel$⁵ and that they evolve on the same time scale, we conclude that the main effect of the small-scale modes is to provide a feedback to the large-scale temperature in the form of the additional heat flux $\langle \widetilde {\boldsymbol {Q}}\rangle _\parallel$.

We can thus summarise the equations that govern the evolution of $\left \langle \varphi ' \right \rangle _\parallel$ and $\left \langle T \right \rangle _\parallel$ as

(4.26)\begin{gather} \partial_t \left\langle \varphi' \right\rangle_\parallel{-} \partial_y \left( \left\langle \varphi' \right\rangle_\parallel{+} \left\langle T' \right\rangle_\parallel \right) + \kappa_T \partial_y \nabla_{{\perp}}^2 \left\langle \varphi' \right\rangle_\parallel \nonumber\\ +\left\{ \left\langle \varphi \right\rangle_\parallel, \left\langle \varphi' \right\rangle_\parallel{-} \nabla_{{\perp}}^2 \left\langle \varphi \right\rangle_\parallel \right\} ' + \boldsymbol{\nabla}_\perp \boldsymbol{\cdot} \left\{ \boldsymbol{\nabla}_\perp \left\langle \varphi \right\rangle_\parallel, \left\langle T \right\rangle_\parallel \right\} ' + \chi \nabla_\perp^4 \left(a\left\langle \varphi' \right\rangle_\parallel{-} b\left\langle T' \right\rangle_\parallel\right) = 0, \end{gather}

(4.27)\begin{gather}\partial_t\left\langle T \right\rangle_\parallel{+}\kappa_T\partial_y \left\langle \varphi \right\rangle_\parallel{+}\left\{ \left\langle \varphi \right\rangle_\parallel, \left\langle T \right\rangle_\parallel \right\} - \chi \nabla_{{\perp}}^2 \left\langle T \right\rangle_\parallel = {-}\left\langle \left\{ \widetilde{\varphi}, \widetilde{T} \right\} \right\rangle_\parallel, \end{gather}

where the influence of small-scale fields appears only in the temperature equation (4.27). The system of (4.13)–(4.15) and (4.26)–(4.27) respects the conservation of the two conserved quantities described in § 2.2; this is shown in Appendix F.

The above reasoning does not apply to the ZFs. Indeed, the equation for $\overline {\varphi }$ is

(4.28)\begin{equation} \partial_t \overline{\varphi} - \overline{\partial_x \left\langle \varphi \right\rangle_\parallel \partial_y \left\langle \varphi + T \right\rangle_\parallel} - \partial_x^2(a\overline{\varphi} - b\overline{T}) = \overline{\partial_x \widetilde{\varphi} \partial_y (\widetilde{\varphi} + \widetilde{T})} = \overline{\widetilde{\Pi}_{xx}}. \end{equation}

This shows that the small-scale stress $\widetilde {\Pi }$ influences the zonal electrostatic potential $\overline {\varphi }$ more strongly (by a factor of $k_{\perp }^{-2}$) than it does the nonzonal $\left \langle \varphi ' \right \rangle _\parallel$. This is a consequence of the electron response (2.1) and the asymptotically smaller ‘inertia’ (i.e., the factor in front of the time derivative) $\propto k_{\perp }^2$ of the ZFs compared with the ‘inertia’ $\propto (1+k_{\perp }^2)$ of the nonzonal $\varphi '$. Thus, the right-hand side of (4.28) cannot be ignored. In fact, as we showed in § 4.1, the addition of 3-D effects, and hence of parasitic modes, has a profound impact on the stability of the Dimits-state ZFs, viz., the momentum flux $\widetilde {\Pi }_{xx}$ extends the Dimits state to higher temperature gradients than the 2-D system allows. Let us show why this is the case.

4.2.5 Turbulent stress due to parasitic sITG modes

In Ivanov et al. (2020), we obtained a prediction for the critical gradient $\kappa _T^{c, \textrm {2D}} (\chi )$ above which a Dimits state with strong ZFs could not be sustained. This prediction was based on considerations of the ratio ${\textrm {Re}{\left (T_{\boldsymbol {k}}/\varphi _{\boldsymbol {k}}\right )}}$ for the linear modes with largest growth rate. As explained in § 4.1.1, this ratio determines the balance of Reynolds and diamagnetic stresses for an individual Fourier mode: if ${\textrm {Re}{\left (T_{\boldsymbol {k}}/\varphi _{\boldsymbol {k}}\right )}} > -1$, then the Reynolds stress is larger and the mode favours a Dimits state, otherwise its diamagnetic stress is larger and the mode helps suppress the coherent ZFs needed for the Dimits state. In two dimensions, this ratio is sensitive to both the temperature gradient $\kappa _T$ and the collisionality $\chi$, and thus an appropriate balance between these two parameters is required in order to have ${\textrm {Re}{\left (T_{\boldsymbol {k}}/\varphi _{\boldsymbol {k}}\right )}} > -1$ for the dominant modes and thus to keep the system in the Dimits state. In particular, for $\kappa _T \gg 1$, the Dimits threshold is given by $\kappa _T / \chi = \textrm {const}$.

Let us adopt a similar approach for the fastest-growing small-scale sITG modes located at $k_\parallel \approx \boldsymbol {\kappa }_T \boldsymbol{\cdot} \boldsymbol {k} k_{\perp }$. Equation (3.26) tells us that these modes satisfy ${\textrm {Re}{\left (T_{\boldsymbol {k}}/\varphi _{\boldsymbol {k}}\right )}} = -1 + O\left (k_{\perp }^{-2}\right )$ for $k_{\perp } \gg 1$. Therefore, to lowest order, the sITG modes are Dimits-marginal, i.e., their Reynolds and diamagnetic stresses balance out. This means that the lowest-order collisionless calculations of § 3.3.2 are insufficient for our needs. While we can extend these calculations to $O\left (k_{\perp }^{-2}\right )$, ignoring $\chi$ is, in fact, an unacceptable oversimplification. As we are about to see, nonzero $\chi$ provides a $O\left (k_{\perp }^{-1}\right )$ correction to ${\textrm {Re}{\left (T_{\boldsymbol {k}}/\varphi _{\boldsymbol {k}}\right )}}$ and hence renders any collisionless higher-order corrections irrelevant.

As discussed in the beginning of § 4, the dominant collisionless sITG modes are found at $k_{\perp }^3 \sim \kappa _T / \chi$. It turns out that at those scales, it is collisional effects that determine ${\textrm {Re}{\left (T_{\boldsymbol {k}}/\varphi _{\boldsymbol {k}}\right )}}$. The details of the relevant calculation can be found in Appendix D. We find that for $\chi$ ordered as $\kappa _T \sim \chi k_{\perp }^3$, the most unstable small-scale sITG modes are still located at $k_\parallel = \boldsymbol {\kappa }_T \boldsymbol{\cdot} \boldsymbol {k}k_{\perp }$ and satisfy

(4.29)\begin{equation} \frac{T_{\boldsymbol{k}}}{\varphi_{\boldsymbol{k}}} ={-}1 - \sqrt{-\hat{\gamma}_{\boldsymbol{k}}^2 + \frac{{i}(a + b - 1)\chi k_{{\perp}}^2}{2\boldsymbol{\kappa}_T \boldsymbol{\cdot} \boldsymbol{k}}} + O\left(k_{{\perp}}^{{-}2}\right), \end{equation}

where we take the branch of the square root with a positive imaginary part. Since $\hat {\gamma }_{\boldsymbol {k}}^2 > 0$ and $\boldsymbol {\kappa }_T \boldsymbol{\cdot} \boldsymbol {k} > 0$ (as stipulated after (4.17)), we find that the sign of the real part of the square root in (4.29) is set by the sign of $a + b - 1$. Plugging in the numerical values $a = 9/40$ and $b = 67/160$, we find $a + b - 1 < 0$, hence the square root has a negative real part and ${\textrm {Re}{\left (T_{\boldsymbol {k}}/\varphi _{\boldsymbol {k}}\right )}} > -1$. Thus, collisionality always pushes the otherwise Dimits-marginal small-scale sITG modes to side with the Reynolds stress and reinforce the ZFs.⁶ This is evident in figure 15.

Figure 15. (a) Linear growth rate and (b) the ratio ${\textrm {Re}{\left (T_{\boldsymbol {k}}/\varphi _{\boldsymbol {k}}\right )}}$ of the most unstable ($k_x=0$) modes versus $k_\parallel$ and $k_y$ for $\kappa _T = 1$ and $\chi = 0.1$. The green dashed line is ${\textrm {Re}{\left (T_{\boldsymbol {k}}/\varphi _{\boldsymbol {k}}\right )}} = -1$. The black dashed line is the location of the largest collisionless growth rate $k_\parallel = \kappa _T k_y^2$. While the green and black lines would coincide to $O(k_{\perp }^{-2})$ for the collisionless modes, we see that the addition of collisions shifts the linearly unstable modes towards the Dimits-favourable ${\textrm {Re}{\left (T_{\boldsymbol {k}}/\varphi _{\boldsymbol {k}}\right )}} > -1$ ratio.

The sensitivity of ${\textrm {Re}{\left (T_{\boldsymbol {k}}/\varphi _{\boldsymbol {k}}\right )}}$ to the numerical factors $a$ and $b$ allows us to carry out a simple test of the above theory. We pick a simulation that is in the Dimits state in three dimensions, but above the 2-D Dimits threshold, i.e., has $\kappa _T > \kappa _T^{c, \textrm {2D}}$. We restart this simulation, but set $a = 1$ for all nonzonal modes. Linearly, this increases nonzonal viscosity and reduces growth rates, without affecting zonal physics. Naïvely, one might expect that with an increased damping of the turbulence, the Dimits state should become ‘stronger’. However, such reasoning does not take into account the structure of the 3-D modes and the change in the balance of Reynolds and diamagnetic stresses stemming from the change of the sign of $a + b - 1$. Indeed, in this numerical experiment, we discover that the Dimits regime is destroyed and strong turbulence sets in, just as the analysis above predicts. This is clear evidence that the most consequential role of collisionality for the Dimits regime of (2.11)–(2.13) is not to dissipate turbulent energy, but rather to regulate the turbulent stress via the ratio ${\textrm {Re}{\left (T_{\boldsymbol {k}}/\varphi _{\boldsymbol {k}}\right )}}$. This also suggests, for future analysis of the Dimits transition in different models of ITG turbulence, that the Dimits threshold may prove to be sensitive to the details of dissipation effects on the unstable modes, especially if, in the absence of collisions, these modes are Dimits-marginal, i.e., if they satisfy ${\textrm {Re}{(T_{\boldsymbol {k}}/\varphi _{\boldsymbol {k}})}} \approx -1$.

Let us also note that in the simple case of sITG modes in slab geometry, a more general calculation that includes kinetic effects is possible. In Appendices G and H, we derive the kinetic sITG dispersion relation and the kinetic equivalent of (4.1). Then, applying the ideas developed in Ivanov et al. (2020) and in this work, we show that ZF-driving small-scale sITG modes are not limited to the cold-ion limit and could play a role outside of the realm of simple fluid approximations. This, of course, can be conclusively confirmed only by appropriate GK simulations.

4.3 Breaking the Dimits state

Recall that the 2-D critical gradient $\kappa _T^{c, \textrm {2D}}$ was found to be an increasing function of $\chi$. Naïvely, this makes sense on the basis of ‘more dissipation means less turbulence’: one expects that one should be able to compensate for an increase in the drive $\kappa _T$ by an appropriate increase in $\chi$ and thus keep the system in the Dimits state. However, this simple picture is false. Collisionality and drive are important for maintaining the Dimits state not because they provide dissipation and injection of energy, but rather because they determine the ratio ${\textrm {Re}{\left (T_{\boldsymbol {k}}/\varphi _{\boldsymbol {k}}\right )}}$ for the linearly unstable modes. In two dimensions, this ratio is sensitive to both $\kappa _T$ and $\chi$; however, this is not the case in three dimensions as the small-scale sITG modes always favour the Dimits state. First, their turbulent momentum flux was shown to satisfy ${\textrm {Re}{\left (T_{\boldsymbol {k}}/\varphi _{\boldsymbol {k}}\right )}} \approx -1$, with collisions pushing this ever so slightly in the Dimits-stable direction of ${\textrm {Re}{\left (T_{\boldsymbol {k}}/\varphi _{\boldsymbol {k}}\right )}} > -1$ (see § 4.2.5). Secondly, they provide an effective thermal diffusion for the large-scale $\left \langle T \right \rangle _\parallel$, which in turn reduces the absolute value of $\left \langle T_{\boldsymbol {k}} \right \rangle _\parallel / \left \langle \varphi _{\boldsymbol {k}} \right \rangle _\parallel$ and partially suppresses the tendency of large-scale modes to destroy the Dimits state (§ 4.2.4). Our numerical simulations show that the combination of the mode structure of the small-scale instability and its influence on large-scale modes proves to be enough to keep the system in the Dimits regime regardless of $\kappa _T$ and $\chi$. As $\kappa _T$ increases beyond the 2-D Dimits threshold $\kappa _T^{c, \textrm {2D}}$, the 2-D modes flip the sign of their turbulent momentum flux and start eroding the ZFs, but the small-scale sITG modes are able to provide enough ZF drive to maintain the Dimits state (see figure 16). Figure 17(a) illustrates the Dimits saturation mechanism.

Figure 16. Dependence of the turbulent viscosity (4.4) on the temperature gradient for $\chi = 0.1$ and $L_{\parallel } = 1$. The 2-D Dimits regime ends at $\kappa _T^{c, \textrm {2D}} \approx 1$. In 3-D simulations, the 2-D modes eventually reverse their turbulent viscosity (red), but the 3-D sITG modes continue to feed the ZFs through a negative turbulent viscosity (blue). The data is taken from simulations with fixed ZF profiles.

Figure 17. Schematic of the flow of energy in (a) the Dimits regime, characterised by strong, turbulence-shearing, staircase ZFs and (b) strong turbulence, where no such ZFs can be generated or sustained. In the Dimits regime, the equilibrium gradients (EG) inject energy into large-scale modes via the 2-D cITG instability. These can then drive ZFs via the secondary instability (see § 2.8 of Ivanov et al. 2020) and small-scale perturbations via the parasitic sITG instability (see § 4.2). In the 2-D Dimits regime ($\kappa _T < \kappa _T^{c, \textrm {2D}}$), the curvature-driven large-scale modes generate a negative turbulent viscosity on the ZFs and hence reinforce the Dimits state. For $\kappa _T > \kappa _T^{c, \textrm {2D}}$, the 2-D modes erode the ZFs, but the ZF drive of the parasitic modes sustains the ZFs (see § 4.2.5). On the other hand, if a Dimits state cannot be achieved, the energy injected into the large-scale modes is transferred to small scales via the parasitic sITG instability, whence it cascades to even smaller, linearly stable scales where it is taken out of the system.

However, the small-scale sITG modes are able to maintain the Dimits state only if the 3-D system is ‘3-D enough’. Namely, if we restrict the system in $z$ by either squeezing it to a small $L_{\parallel }$ (see § 4.3.1) or by cutting off large-$k_\parallel$ modes (see § 4.3.2), we can break the ZF-dominated Dimits regime and push the system into a strongly turbulent state. The former of these methods can be deemed ‘physical’ in the sense that real system can be geometrically limited along the magnetic field, e.g., by magnetic shear. The latter is but a numerical artefact in our cold-ion system; however, parallel transport processes, which were ordered out in (2.11)–(2.13), do provide a large-$k_\parallel$ cutoff for the sITG instability (see § 4.3.2).

Regardless of how the Dimits state is broken, amplitudes remain finite. The small-scale instability is able to extract energy efficiently from the large-scale ($k_{\perp } \ll 1$) fields, into which the cITG instability inputs energy, and to dump it into the small scales $k_{\perp } \gg 1$ of the sITG instability, whence it cascades to smaller scales, where dissipation can take it out of the system. Figure 17(b) shows the flow of energy in the strongly turbulent state.

4.3.1 Effect of parallel system size on the Dimits state

Figure 18(a) shows a typical example of the dependence of the saturated turbulent heat flux $Q$ on the parallel size of the box $L_{\parallel }$ for parameters $\kappa _T$ and $\chi$ that lie beyond the 2-D Dimits regime. For such parameters, the $L_{\parallel } = 0$ system does not reach finite-amplitude saturation. For $L_{\parallel }$ large enough, $Q$ is independent of $L_{\parallel }$, just as it was for parameters that were within the 2-D Dimits threshold (see § 4.1). As $L_{\parallel }$ is decreased, the ZFs break up and the system enters a strongly turbulent state. In figure 18(a), this happens for $L_{\parallel } < 1$. As $L_{\parallel }$ approaches zero, $Q$ starts to increase rapidly, signifying the approach to the 2-D state, where a blow up occurs.

Figure 18. Dependence of the saturated turbulent heat flux $Q$ on (a) the parallel size of the box $L_{\parallel }$ and (b) the largest parallel Fourier mode $k_{\parallel, \textrm {max}}$ that is included in the simulation.

Therefore, for each pair of values of $\kappa _T$ and $\chi$, there exists a critical $L_\parallel ^{c}$ such that the system is in the Dimits state for $L_{\parallel } > L_\parallel ^{c}$ and in the strongly turbulent regime for $L_{\parallel } < L_\parallel ^{c}$. It is clear that $L_\parallel ^{c} = 0$ if $\kappa _T <\kappa _T^{c, \textrm {2D}}$, i.e., if $\kappa _T$ and $\chi$ are such that the 2-D system is able to reach saturation. The dependence of $L_\parallel ^{c}$ on $\kappa _T$ and $\chi$ for $\kappa _T > \kappa _T^{c, \textrm {2D}}$ is not known at this point, due to the numerical cost of resolving simultaneously both the large $k_\parallel$ of the small-scale modes (see § 4.3.2) and the box-sized $k_\parallel \sim L_{\parallel }^{-1}$.

4.3.2 Effect of parallel resolution on the Dimits state

The scale separation between the large-scale cITG modes and the small-scale sITG modes increases the numerical cost of solving (2.11)–(2.13). When the parallel resolution, i.e., the largest $k_\parallel$ in the simulation, is too small, the Dimits state is destroyed numerically and the system is pushed into a strong-turbulence regime for parameters for which a Dimits state would have existed if given sufficient parallel resolution. This is shown in figure 18(b). Empirically, we have found that a good rule of thumb is ‘not to chop the leaves’ of the instability, i.e., to make sure that the wavenumbers that lie within the unstable ‘leaves’ at $k_\parallel \sim \kappa _T k_{\perp }^2$ (see figure 2) are fully included in the simulation.⁷ This, however, rapidly increases the numerical cost of the simulations. Recall that according to (3.33), the collisionless sITG instability satisfies $k_{\perp } \sim (\kappa _T/\chi )^{1/3}$. Therefore, for a fixed $\chi$, the dimensional $k_\parallel$ of the unstable modes is given by

(4.30)\begin{equation} k_\parallel L_B \sim \kappa_T k_{{\perp}}^2 \sim \frac{\kappa_T^{5/3}}{\chi^{2/3}}. \end{equation}

The number of Fourier modes required to resolve a simulation properly then scales as $\kappa _T^{5/3}$, in addition to scaling linearly with $L_{\parallel }$. This quickly renders numerical efforts futile, even for a fluid code.

Of course, the infinitely extending ‘leaves’ of the instability in our 3-D model will, in reality, be ‘chopped off’ by phenomena that have been ordered out of our equations by (2.3). For example, (2.3) orders out the parallel thermal diffusion (Braginskii 1965), but we can nonetheless estimate the dimensional $k_\parallel$ at which this effect will become important. This is the parallel scale at which the collisional heat conduction rate $v_{\textrm {th}i}^2k_\parallel ^2 / \nu _i$ becomes comparable to $\partial _t \sim c_s / L_B$ in (2.11)–(2.13), which happens at

(4.31)\begin{equation} k_\parallel L_B \sim \sqrt{\frac{L_B}{\sqrt{\tau}\lambda_\text{mfp}}}, \end{equation}

where $\lambda _\textrm {mfp} = v_{\textrm {th}i} / \nu _i$ is the mean-free path. In our ordering (2.3), $\lambda _\textrm {mfp} / L_B \sim \tau ^{3/2}$, so we find that the collisional heat conduction comes into play at $k_\parallel L_B \sim 1/\tau$. Formally, this is outside of the regime $k_\parallel L_B \sim 1$ assumed in (2.11)–(2.13), but physically, we conclude that the Dimits regime could be broken if the collisional cutoff (4.30) is superseded by the Braginskii scale (4.31), i.e., if

(4.32)\begin{equation} \frac{L_B}{L_T} \gtrsim \left(\frac{L_B}{\lambda_\text{mfp}}\right)^{7/10}\tau^{{-}11/20}, \end{equation}

where we used $\kappa _T \sim \tau L_B/L_T$ and $\chi \sim L_B\tau ^{3/2}/\lambda _\textrm {mfp}$. In a real fusion device, this condition will not be very difficult to reach, but, in fact, the more relevant mechanism for limiting the parallel wavenumber of the sITG instability is parallel streaming rather than collisional heat conduction. In Appendix G, we show that this too imposes a limit on the parallel wavenumber that is $O(\tau ^{-1})$ too large to be included in our ordering of $k_\parallel$. Namely, the sITG cutoff is given by

(4.33)\begin{equation} k_\parallel^{({c})} L_B = \frac{L_B}{2\sqrt{\rm \pi}(1+\tau)L_T}, \end{equation}

which supersedes the collisional cutoff (4.30) if

(4.34)\begin{equation} \frac{L_T}{\lambda_\text{mfp}} \lesssim \tau (1+\tau)^{3/2}. \end{equation}

Again, such a regime is entirely plausible for a real fusion device.

We conclude that in a more realistic physical regime than the one assumed in the derivation of our model equations (2.11)–(2.13), the behaviour (or even existence) of parasitic sITG modes may be influenced by parallel thermal diffusion or parallel streaming in a way that breaks the Dimits regime at large enough temperature gradients.

5 Discussion

Following our analysis of the Dimits regime and its threshold in the 2-D model of Ivanov et al. (2020), we have been able to extend both our model and our understanding of ITG turbulence to three dimensions. The important qualitative features of the 2-D Dimits state, viz., strong coherent ZFs with patch-wise constant shear, turbulent bursts and localised travelling structures survive the inclusion of 3-D physics largely unchanged (see § 4.1). The ZFs are generated and destroyed by the Reynolds and diamagnetic stresses of sheared ITG turbulence, respectively. If the Reynolds stress is larger, coherent ZFs can be maintained and the system settles into a low-transport Dimits state. Otherwise, a strongly turbulent, high-transport state arises in which saturation occurs unaided by ZFs. In the 2-D model, the ratio of Reynolds to diamagnetic stress is sensitive to the equilibrium parameters – the temperature gradient $\kappa _T$ and the ion collisionality $\chi$ – and thus an appropriate balance of the two is required in order to keep the system within the Dimits regime. With the inclusion of parallel physics, however, the stresses are modified by the 3-D-exclusive sITG instability, which is found always to favour the ZFs (see § 4.2.5). Unless 3-D physics is restricted either by a small parallel box size (§ 4.3.1) or by insufficient numerical resolution (§ 4.3.2), the sITG instability is able to tip the stress balance in the Reynolds direction and a Dimits state is established regardless of the values of the equilibrium parameters.

This 3-D sITG instability is found to be scale-separated from the 2-D cITG instability (see § 4.2.2). In the absence of collisions, the former exists at arbitrarily small perpendicular and parallel scales, while the latter is confined to large scales. This scale separation allows for sITG modes that are predominantly driven not by the equilibrium gradients but rather by the local gradients of large-scale fields, which are themselves driven by the equilibrium gradients (i.e., the sITG instability is parasitic). The nonlinear energy transfer from large-scale to small-scale modes that results from the sITG instability is found to have the form of an effective large-scale thermal diffusion (see § 4.2.4). The combination of this thermal diffusion and the favourable turbulent stress of the small-scale modes are what makes the 3-D Dimits state much more resilient than its 2-D counterpart.

The fact that the Dimits state is governed by essentially the same physical mechanisms in both the 2-D and 3-D cold-ion $Z$-pinch systems gives us not only hope that one day we could understand the Dimits regime of full-blown GK, but also a solid foundation of numerical and analytical work upon which to build such an undertaking. Although there is some numerical evidence of important similarities between these simple systems and GK, e.g., the ferdinon structures seen both by us and by van Wyk et al. (2016, 2017) in their GK simulations of an experimentally realistic configuration, there is still much unknown. The details of the Dimits state in our 3-D model depend on certain peculiar features of cold-ion physics. It is the cold-ion approximation that permits the parasitic small-scale sITG instability that underlies the main differences between the 2-D and 3-D models. As this is only one asymptotic limit of GK, it is difficult to extrapolate any quantitative predictions. However, it is important to note that the kinetic, $\tau \sim 1$, dispersion relation also predicts a collisionless sITG instability at arbitrarily large $k_{\perp }$ (see Appendix G), as was already established by Smolyakov, Yagi & Kishimoto (2002). Just as in the cold-ion fluid model, these sITG modes appear to favour a ZF-dominated state (Appendix H). Thus, the appearance of parasitic modes is not necessarily limited to our cold-ion model and, in certain regimes, could also be a feature of low-collisionality GK. This may also require a careful investigation of GK collisions along the lines of Frei, Hoffmann & Ricci (2022). All of this, combined with the fact that the nature of the Dimits state in the 2-D and 3-D models is essentially the same, encourages us to carry our ideas over into the vastly more complex world of GK. At this point, it is unknown whether the Reynolds–diamagnetic stress competition is also behind the Dimits transition in GK. One of the prominent alternative ideas is the primary–secondary–tertiary scenario, first proposed by Rogers, Dorland & Kotschenreuther (2000). Recently, there have been a number of publications discussing the applicability of this paradigm to both fluid and kinetic models (St-Onge 2017; Zhu, Zhou & Dodin 2018, 2020a, b; Hallenbert & Plunk 2021). Note that, as we showed in Ivanov et al. (2020), the Dimits transition that we observe cannot be explained by the tertiary instability of ZFs. It is possible that the nature of the transition to high transport in realistic GK simulations is, in fact, not as clear-cut as it is in the simple models, but is rather a combination of both mechanisms, viz., the competition between the stresses and a tertiary instability.

Another important feature that our model lacks is magnetic shear. It is well known that this can have a significant effect on both the linear instabilities and turbulence levels in realistic-geometry GK simulations (Kinsey, Waltz & Candy 2006). Notably, much effort today is being devoted to spherical tokamak designs, which can have large values of field-line-averaged magnetic shear combined with nontrivial variations in the local shear. Therefore, we consider the addition of magnetic shear to our analytical and numerical models to be a key direction for future work.