On Nonlocal Energy Transfer via Zonal Flow in the Dimits Shift

TL;DR

This paper demonstrates the existence of the Dimits shift in a modified 2D Terry-Horton equation, revealing a nonlocal energy transfer mechanism via zonal flows that persists across various simplifications and is relevant for plasma physics.

Contribution

It introduces a modified 2D Terry-Horton model capturing the Dimits shift and elucidates a nonlocal energy transfer mechanism through zonal flows, applicable to more complex systems.

Findings

The Dimits shift persists in simplified models including quasilinear and four-mode truncations.

An adiabatic electron response enables efficient nonlocal energy transfer to sound-scale lengths.

The shift size matches numerical solutions and is explained by coupling between unstable and stable modes.

Abstract

The two-dimensional Terry-Horton equation is shown to exhibit the Dimits shift when suitably modified to capture both the nonlinear enhancement of zonal/drift-wave interactions and the existence of residual Rosenbluth-Hinton states. This phenomenon persists through numerous simplifications of the equation, including a quasilinear approximation as well as a four-mode truncation. It is shown that the use of an appropriate adiabatic electron response, for which the electrons are not affected by the flux-averaged potential, results in an $\boldsymbol{E}\boldsymbol{\times}\boldsymbol{B}$ nonlinearity that can efficiently transfer energy nonlocally to length scales on the order of the sound radius. The size of the shift for the nonlinear system is heuristically calculated and found to be in excellent agreement with numerical solutions. The existence of the Dimits shift for this system is then understood as an ability of the unstable primary modes to efficiently couple to stable modes at smaller scales, and the shift ends when these stable modes eventually destabilize as the density gradient is increased. This nonlocal mechanism of energy transfer is argued to be generically important even for more physically complete systems.

Equations

$${\partial }_t\zeta +\mathbf{v}\cdot \nabla \zeta -\beta {\partial }_y\phi =0,$$

$$\zeta \doteq -\delta n_{\mathrm{i}}=({\nabla }_{\perp }^2-1)\phi$$

$$\mathbf{v}\cdot \nabla f$$

$$\left\{\phi ,\zeta \right\}=\left\{\phi ,{\nabla }_{\perp }^2\phi \right\}-\left\{\phi ,\delta n_{\mathrm{e}}\right\}=\left\{\phi ,{\nabla }_{\perp }^2\phi \right\}-\left\{\phi ,\phi \right\}=\left\{\phi ,{\nabla }_{\perp }^2\phi \right\}.$$

$$\mathcal{E}\doteq \frac{1}{L_x{L}_y}{\int }_0^{L_x}dx{\int }_0^{L_y}dy\frac{1}{2}\left[({\nabla }_{\perp }\phi {)}^2+{\phi }^2\right],$$

$$\mathcal{Z}\doteq \frac{1}{L_x{L}_y}{\int }_0^{L_x}dx{\int }_0^{L_y}dy\frac{1}{2}{\zeta }^2.$$

$$-e{n}_{\mathrm{e}}{E}_{\parallel }-{\nabla }_{\parallel }{P}_{\mathrm{e}}=0.$$

$${\nabla }_{\parallel }\left(\phi -\delta n_{\mathrm{e}}\right)=0$$

$$\delta n_{\mathrm{e}}=\phi +C,$$

$$\delta n_{\mathrm{e}}=\phi -⟨\phi {⟩}_{\psi },$$

$$⟨f⟩\doteq \frac{1}{L_y}{\int }_0^{L_y}dyf.$$

$$\zeta ={\nabla }_{\perp }^2\phi -\phi +⟨\phi ⟩\doteq ({\nabla }_{\perp }^2-\hat{\alpha })\phi ,$$

$${\mathcal{E}}^{\mathrm{ZF}}$$

$${\mathcal{E}}^{\mathrm{DW}}$$

$$\delta n_{\mathrm{e}}=(1-\mathrm{i}\hat{\delta })\phi .$$

$$\zeta =({\nabla }_{\perp }^2+\mathrm{i}\hat{\delta }-1)\phi ,$$

$$\hat{\zeta }=-(k_{\perp }^2-\mathrm{i}{\hat{\delta }}_{\mathbf{k}}+1)\hat{\phi }$$

$${\mathrm{i}}\widehat{\delta}_{\boldsymbol{k}}={\mathrm{i}}(\math@atom{\pi}{\mathchoice{\hbox{$\displaystyle\pi$}}{\hbox{$\textstyle\pi$}}{\hbox{$\scriptstyle\pi$}}{\hbox{$\scriptscriptstyle\pi$}}}/2)^{1/2}\left[\omega_{\boldsymbol{k}}-\omega_{*{\mathrm{e}}}\left(1-\frac{1}{2}\eta_{\mathrm{e}}\right)\right]\Big{/}|k_{\parallel}|v_{\mathrm{the}},$$

$$\mathrm{i}{\hat{\delta }}_{\mathbf{k}}=\mathrm{i}{\delta }_0{k}_y\left(\frac{k_{\perp }^2}{1+k_{\perp }^2}-\frac{1}{2}{\eta }_{\mathrm{e}}\right),$$

$$\mathrm{i}{\hat{\delta }}_{\mathbf{k}}=\mathrm{i}{\delta }_0{k}_y,$$

$${\Gamma }_n$$

$$\delta n_{\mathrm{e}}=(1-\mathrm{i}\hat{\delta })\phi -⟨\phi {⟩}_{\psi }.$$

$$\frac{\partial \zeta }{\partial t}+\mathbf{v}\cdot \nabla \zeta =\beta \frac{\partial \phi }{\partial y}-\hat{\alpha }D\zeta ,$$

$$\hat{\zeta }=-(k_{\perp }^2-\mathrm{i}{\hat{\delta }}_{\mathbf{k}}+{\hat{\alpha }}_{\mathbf{k}})\hat{\phi }.$$

$${\gamma }_{\mathbf{k}}$$

$${\omega }_{\mathbf{k}}$$

$$\frac{\partial {\zeta }^{\prime }}{\partial t}$$

$$\frac{\partial U}{\partial t}$$

$$\frac{\partial {\phi }_{\mathbf{k}}}{\partial t}$$

$$\times {\delta }_{\mathbf{k},{\mathbf{k}}_1+{\mathbf{k}}_2}[{\hat{\alpha }}_{{\mathbf{k}}_2}-{\hat{\alpha }}_{{\mathbf{k}}_1}-\mathrm{i}({\hat{\delta }}_{{\mathbf{k}}_2}-{\hat{\delta }}_{{\mathbf{k}}_1})+k_{2\perp }^2-k_{1\perp }^2],$$

$${\partial }_t{\phi }_{\mathbf{p}}$$

$${\partial }_t{\phi }_{{\mathbf{r}}_{+}}$$

$${\partial }_t{\phi }_{{\mathbf{r}}_{-}}^{\ast }$$

Full text

1 Introduction

The Dimits shift is the nonlinear upshift of the critical temperature gradient for the onset of turbulent transport witnessed in simulations of collisionless tokamak plasmas (Dimits et al., 2000). This shift results from the shearing away of turbulent radial streamers by poloidal zonal flows that are generated from the so-called secondary instability (Rogers et al., 2000). The shearing of radial streamers leads to fine-scale structure, which is subsequently damped. The zonal flows are then able to persist on a longer time scale than are the turbulent eddies, since they are not Landau-damped. These residual zonal flows are called Rosenbluth-Hinton states (Rosenbluth and Hinton, 1998), named for those who were the first to show this property.

However, as the temperature gradient is further increased, the system eventually experiences a tertiary instability (Rogers et al., 2000), giving way to turbulent transport. While the qualitative aspects of the Dimits shift are understood, there is yet no complete theory that can simultaneously predict the basic features of the shift, such as its size and dependence on various physical parameters. Understanding the essential aspects of this phenomenon is critical, as it is known that zonal flows have the ability to suppress turbulence in both physical systems (Biglari et al., 1990; Carter and Maggs, 2009; Schaffner et al., 2012) and simplified ones (Ricci et al., 2006). Such a mechanism is one candidate for explaining the L-H transition seen in tokamaks (Burrell, 1997). Being able to predict the saturated level of zonal flows and their effects on turbulence is crucial for the enhancement of plasma confinement (Terry, 2000).

The Dimits shift was first witnessed in gyrokinetic simulations of electrostatic toroidal plasmas (Dimits et al., 1998), and was eventually demonstrated in nonlinear gyrofluid models with Landau-fluid closures (Beer and Hammett, 1998). This culminated in a landmark comparative study of various gyrokinetic and gyrofluid codes by Dimits et al. (2000). The shift has also been captured in some fluid models using a variety of simplifications. One such model is the minimal two-field ion-temperature-gradient (ITG) system, which retains both the ion continuity equation and an equation for the perpendicular ion temperature (Ottaviani et al., 1997; Kolesnikov and Krommes, 2005a, b). Another important model in the study of the zonal-flow/drift-wave interaction has been the two-field Hasegawa-Wakatani model (studied separately by Numata et al. 2007 and Farrell and Ioannou 2009), which is a system that includes both ion-density-gradient drift waves and non-adiabatic electron effects.

The first quantitative study on the secondary instability was by performed by Cowley et al. (1991), who considered a simplified three-field model of the slab ITG mode. It was found that the primary slab instability resulted in elongated eddies (‘streamers’) that were not well accounted for by mixing length arguments, and that these eddies were themselves subject to secondary instabilities that lead to the fragmentation of the primary streamers. Both the secondary and tertiary instabilities of the toroidal ITG mode were later studied by Rogers et al. (2000) using both numerical gyrokinetic simulations and analytical study of a simplified two-field ITG model with finite-Larmor-radius effects. It was found that the nature of the adiabatic electron response resulted in an asymmetry in the secondary and tertiary instabilities, the former being of Kelvin-Helmholtz type, whereas the latter was generally much weaker, requiring a zonal component of perpendicular temperature for instability. This was further elaborated in Jenko et al. (2000), where it was shown that the differences between the electron-temperature-gradient (ETG) and ITG modes were due to the differences between the response of the adiabatic species; the adiabatic electron response for ITG modes led to an additional ${\boldsymbol{E}}\,\boldsymbol{\times}\,{\boldsymbol{B}}$ nonlinearity which suppressed the Kelvin-Helmholtz mode in the tertiary instability, a feature that is absent in the ETG mode.

Since these pioneering works, further progress has been made into understanding the behaviour of the shift using both numerical simulation and analytical techniques, with much work focusing on the two-field Hasegawa-Wakatani and ITG systems. Numata et al. (2007) were able to show a sharp transition between steady turbulent states and relatively quiescent states dominated by zonal flows. However, their simulations had similar levels of viscous damping on both the drift-wave and zonal modes, which prevented them from making a connection to the steady Rosenbluth-Hinton states. They also made a prediction on the boundary of the shift using a stability analysis based on the Kelvin-Helmholtz instability using a simplified flow profile, though the result over-predicts the shift’s size. Farrell and Ioannou (2009) used statistical closure techniques to clearly separate the dynamics of the zonal flows and of the drift-wave modes using a closure called the second-order cumulant expansion, or CE2. Here, the equations of motion are separated into zonally averaged and fluctuating parts, with eddy-eddy self-interaction terms being dropped in the fluctuation evolution equation. (This closure is equivalent to a quasilinear system when ergodicity is assumed.) They were able to demonstrate that the essential physical effects are captured even without the eddy-eddy self-interactions. As the drift-wave and zonal equations are also separated, they were able to use a physically relevant frictional zonal damping (as would result from ion-ion collisions, see Lin et al. 1999) that was distinct from the drift-wave damping. However, their simulations were also stochastically forced to imitate the homogeneous turbulence lost upon neglecting the eddy-eddy self-interactions. This approach is invalid when studying the Dimits shift, however, as there is no turbulence in this regime, thus they were not able to observe a Dimits shift. It remains unclear as to whether the zonal flows in the Hasegawa-Wakatani model can entirely quench drift-wave turbulence as seen in the collisionless simulations of Dimits et al. (2000).

Analytical work has also been done using the minimal two-field toroidal ITG model.  Kolesnikov and Krommes (2005a, b) were able to calculate the Dimits shift of a four-mode truncated system (4MT) using the tools of dynamical systems. However, the size of the shift was found to be strongly dependent on the truncation number of the system. Additionally, the system failed to saturate beyond the Dimits shift. This model has also been studied under the CE2 framework by St-Onge and Krommes (2017), in which the effect of discrete-particle noise on the onset of zonostrophic instability was studied. Even for such a simple two-field system, the CE2 results in a system of five quantities (two zonal-averaged fields and three independent components of a covariance tensor with the off-diagonal component being complex) and even the simplest of linear calculations can quickly become tedious. In addition, both the Hasegawa-Wakatani and minimal two-field ITG systems are non-normal, and as a result can be made to exhibit subcritical turbulence, a phenomenon that can obfuscate the essential physics of drift waves and zonal flows.

What, then, is the simplest model that exhibits the Dimits shift with a minimal amount of physics? To be successful, the model must

contain a nonlinear boost in the interaction between zonal flows and drift waves, leading to an increase in energy transfer between these modes and a quenching of turbulence through the shearing of eddies (this ‘boost’ is made specific in section 2); 2. 2.

allow for the existence of steady states consisting purely of zonal modes. For these solutions to be proper steady states, zonal modes would necessarily be linearly undamped; 3. 3.

exhibit a finite Dimits shift. In other words, given a sufficient amount of linear drive, steady zonal states must give way to turbulence leading to a finite amount of turbulent transport; 4. 4.

have the ability for saturation beyond the end of the Dimits shift without the need for zonal damping. Otherwise, saturated states may depend on the amount of zonal damping, in contrast with the original numerical studies of Dimits et al. (2000).

The goal of this article is to show that the Terry-Horton equation can be suitably modified to exhibit all four traits. The resulting equation, henceforth referred to as the modified Terry-Horton equation (mTHE), can be used to gain insight into the fundamental aspects of the Dimits shift. In particular, it is shown that the use of an appropriate adiabatic electron response, for which the electrons are not affected by the flux-averaged potential, results in an ${\boldsymbol{E}}\,\boldsymbol{\times}\,{\boldsymbol{B}}$ nonlinearity that can transfer energy nonlocally in $k$-space via zonal shearing, depositing energy from scales at $L\gg\rho_{\mathrm{s}}$ to scales at $L\sim\rho_{\mathrm{s}}$ (here, $\rho_{\mathrm{s}}$ is the sound radius, defined in section 2.) This energy transfer is in stark contrast to the usual Kolmogorov-type cascade which is scale-by-scale. By using some simple considerations of this nonlocal transfer, the existence of the Dimits shift can be understood in terms of efficient coupling to small-scale stable modes, and the size of the shift depends on when these small-scale modes destabilize. This effect is a generic property for systems with this type of adiabatic electron response, and should be relevant even when adopting a more physically complete framework (e.g. gyrokinetics). The concepts presented here should then hopefully serve to connect the previous work of Rogers et al. (2000) with the more recent ideas of mode coupling (see, for example, Hatch et al. 2011; Makwana et al. 2014).

Before proceeding, this article focuses solely on the collisionless Dimits shift, where the zonal flows are linearly undamped. In physical systems, some amount of zonal damping must of course be present. Even so, the underlying mechanisms present in idealized models, such as the one presented herein, which capture the fundamental aspects of the Dimits shift, should remain important when additional physical effects, such as zonal damping, are included. Remarkably, it has been noted recently by Mikkelsen and Dorland (2008) that the Dimits shift can persist even in more physically complete systems. Thus, as the mechanisms that underlie the Dimits shift appear to be robust, studying them in isolation by using a minimal model can pinpoint which effects are truly necessary for the shift to exist. This is precisely the point of view adopted in this paper.

The remainder of this article is organized as follows. In Sec. 2 I describe the historically important models that lead up to the mTHE, which is described in Sec. 3. In Sec. 4 I demonstrate by direct numerical simulation that the mTHE exhibits the Dimits shift. Analytical results are presented in Sec. 5, starting with analysis of the 4MT, followed by a heuristic calculation of the size of the Dimits shift for the nonlinear system. This is done in Sec. 5.3. Finally, the work is summarized in Sec. 6.

2 Review

2.1 The (Modified) Hasegawa-Mima Equation

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

The HME is given by

[TABLE]

where ${\boldsymbol{v}}=(-\partial_{y}\varphi,\partial_{x}\varphi)$ is the ${\boldsymbol{E}}\,\boldsymbol{\times}\,{\boldsymbol{B}}$ drift velocity, $\beta\doteq a/L_{n}$ parameterizes the density gradient, $L_{n}^{-1}\doteq-\operatorname{d}\!\,\ln n_{0}/\operatorname{d}\!r$ is the background density gradient scale length,

[TABLE]

is the generalized vorticity derived from the gyrokinetic Poisson equation assuming an adiabatic electron response $\delta n_{\mathrm{e}}=\varphi$, $\delta n_{\mathrm{e}}$ is the perturbed electron density, and $\varphi\doteq(e\phi/T_{\mathrm{e}})(a/\rho_{\mathrm{s}})$ is the dimensionless perturbation to the electrostatic potential. Here, $\delta n_{\mathrm{i}}$ is the perturbed ion gyrocenter density normalized to the background density, $e$ is the unit charge, $\phi$ is the perturbed electrostatic potential, and $T_{\mathrm{e}}$ is the background electron temperature. Throughout this paper, time and space are normalized by $a/c_{\mathrm{s}}$ and $\rho_{\mathrm{s}}$, respectively, where $a$ is the minor radius, $\rho_{\mathrm{s}}\doteq c_{\mathrm{s}}/\omega_{\mathrm{ci}}$ is the sound radius, and $c_{\mathrm{s}}\doteq(ZT_{\mathrm{e}}/m_{\mathrm{i}})^{1/2}$ is the sound speed. Coordinates are in a local Cartesian grid where $x$ represents the radial coordinate and $y$ represents the poloidal coordinate. The operator $\nabla_{\perp}$ denotes the gradient perpendicular to the magnetic field and acts in the $x$ and $y$ direction for this particular setup. For any function $f$,

[TABLE]

which defines the Poisson bracket $\{\ldotp,\ldotp\}$. This nonlinearity, which arises from the ${\boldsymbol{E}}\,\boldsymbol{\times}\,{\boldsymbol{B}}$ drift, consist of two components: the ion-polarization or vorticity nonlinearity $\left\{\varphi,\nabla_{\perp}^{2}\varphi\right\}$ and the ${\boldsymbol{E}}\,\boldsymbol{\times}\,{\boldsymbol{B}}$ nonlinearity $\left\{\varphi,\delta n_{\mathrm{e}}\right\}$. While the ion-polarization nonlinearity has a specific form, the ${\boldsymbol{E}}\,\boldsymbol{\times}\,{\boldsymbol{B}}$ nonlinearity can vary depending on which physical effects are included in the electron response. As the Poisson bracket is anticommutative and bilinear, it follows that for purely adiabatic electron response,

[TABLE]

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

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

[TABLE]

and the generalised enstrophy density

[TABLE]

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

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

[TABLE]

where $E_{\parallel}$ is the electric field parallel to the magnetic field, $n_{\mathrm{e}}$ and $P_{\mathrm{e}}$ are respectively the electron density and pressure, and $\nabla_{\parallel}$ is the gradient along the magnetic field ($\nabla_{\parallel}\rightarrow\partial_{z}$ for HME). Assuming isothermal electrons, (7) leads to

[TABLE]

in dimensionless units. This has the solution

[TABLE]

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

Physically, electrons that evolve adiabatically will experience no radial transport, and so the density perturbation averaged over a flux surface must vanish. Thus, taking the flux-surface average of (9) leads to $C=-\langle{\varphi}\rangle_{\psi}$, or

[TABLE]

where $\langle{\ldots}\rangle_{\psi}$ denotes the flux-surface average which, in two dimensions, is equivalent to the zonal average $\langle{f}\rangle$,

[TABLE]

Equation (2) is thus modified to now read

[TABLE]

where $\widehat{\alpha}$ is an operator that is zero when acting on zonal modes and is unity otherwise. Notice that unlike (4), the adiabatic electron response given by (10) does not disappear in the Poisson bracket (Jenko et al., 2000), so the ${\boldsymbol{E}}\,\boldsymbol{\times}\,{\boldsymbol{B}}$ nonlinearity does appear in the mHME, the importance of which is discussed in section 5.

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

[TABLE]

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

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

2.2 The Terry-Horton Equation

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

[TABLE]

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

[TABLE]

or

[TABLE]

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

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

[TABLE]

where $\omega_{*{\mathrm{e}}}\doteq ck_{y}T_{\mathrm{e}}/eB_{0}L_{n}$ is the electron diamagnetic drift frequency, $c$ is the speed of light, $B_{0}$ is the background magnetic field magnitude, and $\eta_{\mathrm{e}}\doteq\operatorname{d}\!\,\ln T_{\mathrm{e}}/\operatorname{d}\!\,\ln n_{0}$. To proceed, $\omega_{\boldsymbol{k}}$ is taken to be the frequency of the linearized HME. Then, in dimensionless units,

[TABLE]

where \delta_{0}=(\math@atom{\pi}{\mathchoice{\hbox{\displaystyle\pi}}{\hbox{\textstyle\pi}}{\hbox{\scriptstyle\pi}}{\hbox{\scriptscriptstyle\pi}}}/2)^{1/2}(m_{\mathrm{e}}/m_{\mathrm{i}})^{1/2}/|k_{\parallel}L_{n}| is a constant of order unity that parameterizes the parallel wavenumber $k_{\parallel}$ (as $L_{n}$ is already parameterized by $\beta$). One arrives at the original form of $\widehat{\delta}_{\boldsymbol{k}}=\delta_{0}k_{y}(k_{\perp}^{2}-\eta_{\mathrm{e}}/2)$ given by Terry and Horton (1983) by taking the long-wavelength ($k_{\perp}^{2}\ll 1$) limit. As another example, for trapped collisional electron dynamics,

[TABLE]

where \delta_{0}=\eta_{\mathrm{e}}(\epsilon/2)^{1/2}(6/\math@atom{\pi}{\mathchoice{\hbox{\displaystyle\pi}}{\hbox{\textstyle\pi}}{\hbox{\scriptstyle\pi}}{\hbox{\scriptscriptstyle\pi}}}^{1/2})/\nu_{\mathrm{eff}} is a constant of order unity, $\epsilon$ is the inverse aspect ratio and $\nu_{\mathrm{eff}}$ is the collisional detrapping rate. A thorough review of such mechanisms is given by Tang (1978).

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

[TABLE]

which is positive-definite. Another important feature of the THE is that it has only one nonlinearly conserved quantity, the generalized vorticity $\zeta$, which can alter the system’s ability to experience a dual cascade. For instance, $\widehat{\delta}_{\boldsymbol{k}}$ as given by (20) becomes much larger relative to $\nabla_{\perp}^{2}$ at large scales, leading to degeneracy between $\mathcal{E}$ and $\mathcal{Z}$. This specific form of $\widehat{\delta}_{\boldsymbol{k}}$ has been shown by Liang et al. (1993) to disable the dual cascade. However, $\widehat{\delta}_{\boldsymbol{k}}$ as given by (19) is small relative to $\nabla_{\perp}^{2}$ at both small and large scales. For this type of system the dual cascade is expected to remain prevalent.

3 The Modified Terry-Horton Equation

3.1 Description

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

[TABLE]

The resulting system is hence referred to as the modified Terry-Horton Equation (mTHE), and is given by

[TABLE]

where ${\boldsymbol{v}}=(-\partial_{y}\varphi,\partial_{x}\varphi)$ and

[TABLE]

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

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

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

3.2 Linear Properties

The eigenvalues of (23) for the linearization around the zero state can be readily calculated. For drift-wave modes, the growth rates and frequencies for a time dependence of ${\mathrm{e}}^{\lambda_{\boldsymbol{k}}t}$ where $\lambda_{\boldsymbol{k}}=\gamma_{\boldsymbol{k}}-{\mathrm{i}}\omega_{\boldsymbol{k}}$ are, respectively,

[TABLE]

For zonal modes, both are identically zero by construction, as I have neglected collisional damping.

3.3 Quasilinear Model

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

[TABLE]

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

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

3.4 Four-Mode Truncation

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

[TABLE]

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

The resulting set of equations is

[TABLE]

along with their complex-conjugate counterparts. Here, $\widehat{\delta}_{\pm}\doteq\widehat{\delta}_{\boldsymbol{p}}\pm\widehat{\delta}_{\boldsymbol{r}}$, and the $M_{\boldsymbol{k}}$’s are mode-coupling coefficients of the interaction between zonal flows and drift waves and are given by

[TABLE]

It will be useful to express these coefficients in terms of their real and imaginary parts,

[TABLE]

To obtain these, I have used the fact that $\widehat{\delta}_{\boldsymbol{-k}}=-\widehat{\delta}_{\boldsymbol{k}}$.

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

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

[TABLE]

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

4 Numerical results

4.1 Setup

The nonlinear (NL) system (23) and quasilinear (QL) system (27) are simulated pseudospectrally and dealiased on a square Cartesian grid with L=20\math@atom{\pi}{\mathchoice{\hbox{\displaystyle\pi}}{\hbox{\textstyle\pi}}{\hbox{\scriptstyle\pi}}{\hbox{\scriptscriptstyle\pi}}} and $N=256$ on each side. Time stepping is performed using third-order Adams-Bashforth with an integrating factor. The vorticity at $t=0$ is initialized with Gaussian noise of zero mean and standard deviation $5\times 10^{-3}$, with no energy in zonal modes. The random number generator used for the initial state is initialized with the same seed for all simulations. For all simulations in this section, $D=\mu-\nu\nabla_{\perp}^{2}$ where $\mu=1$ and $\nu=10^{-2}$. Here, friction $\mu$ is added to dissipate energy at large-scales and is added as a preventative measure against a possible inverse cascade. (It is shown in Sec. 5 that the qualitative aspects of the simulations do not depend on the choice of damping operator.) To simulate the four-mode truncation, I choose modes with $q_{x}=p_{y}=1$ (the dependence of this system on the specific values of $q_{x}$ and $p_{y}$, as well as the damping operator $D$, is discussed in Sec. 5). Finally, the value of $\beta$ at which the Dimits shift ends is denoted as a critical density gradient $\beta^{*}$. The size of the shift then is $\Delta\beta\doteq\beta^{*}-\beta_{\mathrm{lin}}$.

I run separate simulations using the different $\widehat{\delta}_{\boldsymbol{k}}$’s given by (19) and (20), the former with $\delta_{0}=2$ and $\eta_{\mathrm{e}}=0$, and the latter with $\delta_{0}=1.5$. To differentiate between the two systems, I denote the former by (NL1), (QL1), etc., and the latter by (NL2), (QL2), etc. For results that are insensitive to the details of $\widehat{\delta}_{\boldsymbol{k}}$, I simply use the unnumbered abbreviations. The parameters given above correspond to a threshold of linear instability at $\beta_{\mathrm{lin}}\approx 4.74$ for $\widehat{\delta}_{\boldsymbol{k}}$ given by (19), and $\beta_{\mathrm{lin}}\approx 4.21$ for $\widehat{\delta}_{\boldsymbol{k}}$ given by (20).

4.2 Direct Numerical Simulation

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

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

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

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

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

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

5 Analytical Results

I begin this section with a linear stability analysis of both the secondary and tertiary instabilities for the 4MT. With these results in hand, the size of the Dimits shift can be estimated for the fully nonlinear system, and the underlying mechanism behind the shift in this system can be well understood in terms of efficient energy transfer between unstable and stable modes.

5.1 Zonal growth rates of the secondary instability

At the beginning of the simulation, the most unstable drift-wave mode grows without any influence from the zonal or side-band modes in what I shall call the secondary regime. One can then solve (29b–d) with $\varphi_{\boldsymbol{p}}=\varphi_{0}{\mathrm{e}}^{(\gamma_{\boldsymbol{p}}+{\mathrm{i}}\omega_{\boldsymbol{p}})t}$, which renders the remaining equations a linear system:

[TABLE]

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

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

To do so, the transformations $\varphi_{{\boldsymbol{r}}_{+}}=\varphi_{0}\varphi^{\prime}_{{\boldsymbol{r}}_{+}}{\mathrm{e}}^{-{\mathrm{i}}\omega_{\boldsymbol{p}}t}$ and $\varphi^{*}_{{\boldsymbol{r}}_{-}}=\varphi_{0}^{*}\varphi^{\prime*}_{{\boldsymbol{r}}_{-}}{\mathrm{e}}^{{\mathrm{i}}\omega_{\boldsymbol{p}}t}$ are made to eliminate the rapid oscillatory behaviour of the drift-wave mode. New equations of motion are then formulated for the variables $\Lambda_{\pm}\doteq\varphi^{\prime}_{{\boldsymbol{r}}_{+}}\pm\varphi^{\prime*}_{{\boldsymbol{r}}_{-}}$. Defining $\gamma_{\pm}\doteq\gamma_{\boldsymbol{p}}\pm\gamma_{\boldsymbol{r}}$ and $\omega_{\pm}\doteq\omega_{\boldsymbol{p}}\pm\omega_{\boldsymbol{r}}$, (33a–c) become

[TABLE]

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

[TABLE]

where

[TABLE]

The usual dispersion relation of the modulational instability for the HME is recovered by setting $\widehat{\delta}_{\boldsymbol{k}}=0$ and $\widehat{\alpha}_{\boldsymbol{q}}=1$ (Connaughton et al., 2010).

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

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

[TABLE]

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

To analyze the limit $t\rightarrow``\infty"$, the ansatz $\varphi_{\boldsymbol{q}}\sim e^{S(t)}$ is used in (35), resulting in

[TABLE]

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

[TABLE]

where

[TABLE]

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

[TABLE]

where

[TABLE]

and $I_{3}$ is a cube root of unity. Both (39) and (41) represent instabilities of Kelvin-Helmholtz type ($\gamma\sim|q_{x}p_{y}\varphi_{\boldsymbol{p}}|^{\varepsilon}\sim|q_{x}v_{x}|^{\varepsilon}$ where $\varepsilon$ is a constant of order unity). This form of the growth rate agrees with the usual picture of the secondary instability given by Rogers et al. (2000), and is roughly maximized at scales between $k_{y}^{-1}$ and the sound radius (i.e. $k_{x}\sim k_{y}$ to $k_{x}\rho_{\mathrm{s}}\sim 1$).

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

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

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

[TABLE]

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

[TABLE]

In tokamak plasmas, $p_{y}\ll 1$, and so this results in a nonlocal transfer of energy from the primary drift wave to zonal flows at scales of order $\sim\rho_{\mathrm{s}}$, which is at odds with the conventional wisdom that energy flows to zonal flows at the largest scales (Diamond et al., 2005). This was previously seen in both Rogers et al. (2000) and Guzdar et al. (2001). When physical electron response is inconsistently taken into account under the partial time derivative but not in the Poisson bracket (as was done in Holland et al. 2003), the zonal flow growth rate is greatly reduced. This is denoted in figure 6 by a dashed line. While the above analysis of the zonal flow growth rate is based on a simplified four-mode truncation, the cutoff given by (44) is also observed in statistical closures of the mHME, including CE2 and a geometrical-optics reduction of CE2 (Parker, 2016). The importance of nonlocal energy transfer effected by the ${\boldsymbol{E}}\,\boldsymbol{\times}\,{\boldsymbol{B}}$ nonlinearity is further elaborated upon in the next section.

5.2 Zonal mode stability analysis of the tertiary instability

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

[TABLE]

The eigenvalues are quickly found: $\lambda=(\gamma_{\boldsymbol{r}}-{\mathrm{i}}\omega_{\boldsymbol{r}})$ and

[TABLE]

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

[TABLE]

where

[TABLE]

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

When the outer discriminant of (47) vanishes, the real part of both the unstable and stable eigenvalues becomes ${\mathrm{Re}}(\lambda_{+})=\gamma_{+}/2$. This corresponds to a special situation where both the drift-wave and sideband modes become maximally coupled, and occurs when $\Theta=0$ and $\Omega\leq 0$. Here, both modes work together in tandem to create a coupled mode with an effective damping rate that is the average of the individual growth rates. This coupled mode will be referred to as the Maximally-Coupled Mode (MCM). In principle, MCMs can exist for arbitrary drift-wave/sideband pairings where the primary drift-wave mode contains a non-zero radial component $p_{x}$, though in this section I only consider the most unstable mode with $p_{x}=0$.

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

[TABLE]

This predicts the end of the Dimits shift for the system considered in Sec. 4.2 to be $\beta^{*}=5.2$, in excellent agreement with the simulation result in figure 3.

In order to develop MCMs, the condition $\Omega\leq 0$ must be satisfied. In order for zonal flows to be stabilizing, the condition

[TABLE]

must be met. It is important to realize that $M_{\boldsymbol{k}}\sim q_{x}p_{y}$ times a factor of order unity. This term then clearly represents a Kelvin-Helmholtz-type destabilization of the zonal mode, and can become dominant for sufficiently large zonal amplitudes and negative values of the mode-coupling coefficient. However, this term only needs to be comparable to $\gamma_{-}^{2}-\omega_{-}^{2}$ in order to spoil the coupling between drift-wave modes. This happens for a sufficiently large value of $q_{x}$ (denoted by $q_{x}^{*}$) and signifies the smallest scale where the last MCM is formed.

The other condition that must be satisfied in order for MCMs to develop is $\Theta=0$. There then exists a zonal amplitude that is most stable, given by

[TABLE]

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

In the introduction, four ingredients in a successful model of the Dimits shift were highlighted, one of which was the necessity for a nonlinear boost in the interaction between zonal flows and drift waves. This boost is provided by the $\widehat{\alpha}$ operator, which modifies the ion density in both the partial derivative and the Poisson bracket of the HME. While the former effect has been appreciated in the literature (Diamond et al., 2005), the latter has largely gone unnoticed, and in some cases is entirely neglected (e.g., Holland et al. 2003; Makwana et al. 2012). To appreciate the importance of the ${\boldsymbol{E}}\,\boldsymbol{\times}\,{\boldsymbol{B}}$ nonlinearity, consider the mode-coupling coefficients of the (m)HME where $M_{\boldsymbol{p}}^{{\mathrm{Im}}}M_{\boldsymbol{r}}^{\mathrm{Im}}=0$ (equation (30) with $\widehat{\delta}_{\boldsymbol{k}}=0$) with a background drift-wave at wavenumber $p_{y}$. In order for energy to flow from an unstable drift-wave mode to a stable one through a zonal flow, the zonal flow itself must be stable to the Kelvin-Helmholtz instability. Otherwise, energy is fed back from the zonal flow to the unstable drift-wave mode. For a zonal flow to be stabilizing, the condition $M_{\boldsymbol{p}}^{{\mathrm{Re}}}M_{\boldsymbol{r}}^{\mathrm{Re}}>0$ must be met, leading to destabilizing zonal flows at scales given by

[TABLE]

for the HME ($\widehat{\alpha}_{\boldsymbol{q}}=1$). This is equivalent to the stability threshold given by Diamond et al. (2005), where the hydrodynamic limit was considered. However, the $\widehat{\alpha}$ operator in the mHME (viz. (30) with $\widehat{\alpha}_{\boldsymbol{q}}=0$) leads to zonal flows that are stable against the Kelvin-Helmholtz instability down to the scale set by

[TABLE]

Typically, $p_{y}\ll 1$ in toroidal plasmas, leading to $q_{x}$ and $p_{y}$ that are of disparate scales. This is a rather surprising result, as stable zonal flows at scales on the order $\sim\rho_{\mathrm{s}}$ can be an efficient catalyst of nonlocal energy transfer, which is in stark contrast to the scale-by-scale transfer in a Kolmogorov-type cascade, whose transfer rate is set by the local eddy turnover time. Additionally, the magnitude of $M_{\boldsymbol{p}}^{{\mathrm{Re}}}M_{\boldsymbol{r}}^{\mathrm{Re}}$ in the mHME is increased by a factor of $p_{y}^{-2}$ for $p_{y}\ll 1$ compared to that of the HME, resulting in a much stronger stabilizing effect. 111This, however, comes at the price of an increased destabilizing effect by the Kelvin-Helmholtz instability for scales with wavenumber $q_{x}^{2}\gg 1+p_{y}^{2}$. Luckily, the amplitudes of the zonal flows at these scales tend not to be large.

5.3 Heuristic calculation of the Dimits shift

I now generalize the above results for the 4MT in order to heuristically calculate the Dimits shift for the fully nonlinear system. One plausible way to do so is to assume that the zonal interaction couples the most unstable drift-wave mode to every sideband mode within the interaction range given in the previous section. Then the threshold for stability is not given by $\gamma_{+}=0$, but rather (2\math@atom{\pi}{\mathchoice{\hbox{\displaystyle\pi}}{\hbox{\textstyle\pi}}{\hbox{\scriptstyle\pi}}{\hbox{\scriptscriptstyle\pi}}}/L_{x})\sum_{{\mathrm{MCM}}}\gamma_{\boldsymbol{k}}\approx\int_{0}^{q_{x}^{*}}\gamma_{\boldsymbol{k}}(q^{\prime},p_{y}^{*})\operatorname{d}\!q^{\prime}=0, where the sum is over all modes with $k_{y}=p_{y}^{*}$ in the range $|k_{x}|\leq q_{x}^{*}$. Here, I have used the idea of phase space contraction that was given in section 3.4, but now only consider modes that are well-coupled by zonal shearing to the primary (most unstable) drift-wave mode. Passing into the continuum here is equivalent to taking the large box limit ($L_{x}\rightarrow\infty$).

The above calculation differs fundamentally from the linear stability analyses of both Rogers et al. (2000) and Numata et al. (2007). While the latter set of calculations focus on the stability of simplified background zonal flow profiles (resulting in localized instabilities), here the details of the background profile are unimportant; rather, the background profile is assumed to be whatever it needs to be in order to maximally couple the primary drift-wave mode to stable ones, and that it is such that Kelvin-Helmholtz modes are not excited. The criteria for stability then is not determined by the flow profile, but rather by the zonal flow’s ability to transfer enough energy to quench the primary instability. Thus, the stability of only a small subset of modes needs to be considered, those that are maximally coupled to the primary drift wave.

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

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

[TABLE]

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

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

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

One reason why this heuristic calculation works well for this system is that there are only two competing mechanisms for energy transfer: the direct cascade (which is local), and the zonal shearing (which is nonlocal). While the nonlocal mechanism can be taken into account in more complete systems, added complexity comes with additional competing forms of energy transfer. One such mechanism is the dual cascade, which has the ability to transfer both energy to large scales (where large-scale sinks may reside) and enstrophy nonlocally to small scales. This process is dominated by the largest eddies of the system, and underlies the typical picture of shearing by large-scale zonal flows (Diamond et al., 2005). More complex models also have an increased number of fields, leading to multiple branches (see the recent work of Makwana et al. 2012, 2014). Typically, one branch is unstable in a region of $k$-space while the others are purely stable for all $k$. Coupling energy between branches, then, can lead to an efficient source of energy dissipation. This is exacerbated when moving to a kinetic system, where there could be an arbitrary number of stable branches. Finally, additional nonlinearities can also arise, such as finite Larmor-radius terms. Thus, one has to take into account all these effects as they arise, and determine their relative impact on stability accordingly.

6 Discussion and Conclusion

I have shown through direct numerical simulation that the Terry-Horton equation can be made to exhibit the Dimits shift with two suitable modifications. First, an appropriate adiabatic electron response is added to ensure that electrons do not respond to a potential that is constant along a flux surface. Secondly, zonal modes are made to be explicitly undamped, thus capturing the residual Rosenbluth-Hinton states seen in gyrokinetic simulations. The Dimits shift is shown to persist even after various simplifications are made. Analytical progress was made on a four-mode truncation of the system, focusing on the behaviour of the zonal mode growth during the secondary instability, and calculating an upper bound on the end of the Dimits shift.

Using this model, important insight was gained on the underlying mechanism by which the Dimits shift operates. Specifically, it was shown that the inclusion of the appropriate adiabatic electron response results in a ${\boldsymbol{E}}\,\boldsymbol{\times}\,{\boldsymbol{B}}$ nonlinearity that can transfer energy nonlocally from large scales to scales of order $\sim\rho_{\mathrm{s}}$. The nonlocal transfer as a result of nonadiabatic electron response in the ${\boldsymbol{E}}\,\boldsymbol{\times}\,{\boldsymbol{B}}$ nonlinearity has been known for some time (Liang et al., 1993). However, while the nonlocal contribution from adiabatic electrons in the ${\boldsymbol{E}}\,\boldsymbol{\times}\,{\boldsymbol{B}}$ nonlinearity has been appreciated in some areas of the literature (see, for example, Rogers et al. 2000; Guzdar et al. 2001), it has gone unnoticed in other areas. As an example, two studies (Holland et al., 2003; Makwana et al., 2012) that used a two-field gyrofluid model of the ITG instability neglected the ${\boldsymbol{E}}\,\boldsymbol{\times}\,{\boldsymbol{B}}$ nonlinearity altogether, and thus this effect was entirely absent. Similarly, the nonlinear contribution from proper adiabatic electron response is not mentioned in the work of Diamond et al. (2005); Itoh et al. (2006), while the linear contribution under the partial time derivative is considered. An important implication of this correction is that it results in a broader range of scales where zonal flows are stable to the Kelvin-Helmholtz instability, thus leading to a more favourable stability criterion than that given in the review article by Diamond et al. (2005), where the hydrodynamic limit was considered. It is important to realize that, as modes arising from the primary and secondary instabilities are usually at scales larger than the sound radius, the adiabatic ${\boldsymbol{E}}\,\boldsymbol{\times}\,{\boldsymbol{B}}$ nonlinearity will at least be as large as, if not larger than, the ion-polarization nonlinearity, and so it is never appropriate to neglect this term in the study of ITG modes.

It is worthwhile to clarify a result given in the work of Rogers et al. (2000). There, it was stated the the Kelvin-Helmholtz instability is absent in the tertiary instability, and a finite component of $T_{\perp}$ was needed for instability. As an example, they considered the stability of a zonal flow with $k_{x}\rho_{\mathrm{s}}=0.25$, which was shown to be unconditionally stable to the Kelvin-Helmholtz instability (consistent with the analysis of section 5.2.) However, it must be understood that modes with $k_{x}\rho_{\mathrm{s}}>1$ can go unstable to the Kelvin-Helmholtz instability, even without a component of $T_{\perp}$. While these small scale modes are not relevant when discussing the usual ion-polarization-driven forward cascade which is dominated by zonal flows at the largest scales, it is rather crucial for the behaviour of the adiabatic ${\boldsymbol{E}}\,\boldsymbol{\times}\,{\boldsymbol{B}}$ nonlinearity, which depends on zonal flows at a broad range of scales. At which scales these flows ultimately go unstable, then, determines the effective interaction range in $k$-space of this nonlinearity.

It is also stated in the work of Rogers et al. (2000) that the addition of background gradients generally overpowers the instability given by $T_{\perp}$. This is interesting in the context of equation 47, which states that (at least in the four-mode truncation of the mTHE) the stability criterion regarding the end of the Dimits shift is ultimately dictated by the linear terms, regardless of the size of the zonal flow. Thus, neglecting the linear terms in a stability analysis of the Dimits shift regime may never be justified, even when other avenues of destabilization are present.

Future work should focus on extending this analysis to more physically complete (e.g., gyrokinetic) systems to see if one can derive the size of shift originally seen for the Cyclone Base Case in Dimits et al. (2000). Other quantities of interest are the saturation levels and spectra of the zonal flows resulting from the secondary instability, and the saturation levels of turbulent transport beyond the shift. While this work was built upon a very idealized model, the ${\boldsymbol{E}}\,\boldsymbol{\times}\,{\boldsymbol{B}}$ nonlinearity that arises from the non-adiabatic electron response to zonal modes is present in most relevant models of tokamak plasmas, including gyrokinetics and higher order gyrofluid closures. The nonlocal transfer of energy catalyzed by zonal flows, then, is an effect that should remain important in more complete systems. Future research should aim to confirm this.

Finally, a better estimate of the size of the Dimits shift for the mTHE could be derived using more rigorous methods (at the level of Krommes and Kim 2000). While it is known that the energy transfer at large scales is dominantly to zonal models, if this can be shown to be the case down to scales with wavenumber $q_{x}$, defined using $q_{x}^{2}\rho_{\mathrm{s}}^{2}\sim 1+p_{y}^{2}\rho_{\mathrm{s}}^{2}$ where $p_{y}$ is the wavenumber of the most unstable mode, then the heuristic calculations presented here could be better justified. This can be done, in principle, by using statistical closures to study the nonlinear mechanisms. However, as Parker and Krommes (2014) have pointed out, it is necessary to begin with an inhomogeneous closure in order that one can consider inhomogeneous symmetry-breaking perturbations (zonal flows) to a state of homogeneous turbulence. (As was discussed by St-Onge and Krommes 2017, here ‘turbulence’ below the point of zonostrophic instability refers to homogeneous noise due to discrete particles.) Not only does an inhomogeneous closure allow for symmetry breaking, it contains all of the physical effects involved in destabilizing those flows and allowing for a transition from the Dimits-shift regime to states of fully developed turbulence. Because the general structure of an inhomogeneous closure is necessarily complicated, carrying out such a program to completion represents a significant challenge for the future.

Acknowledgements.

I would like to thank J. A. Krommes, G. W. Hammett, and T. Stoltzfus-Dueck for many instructive discussions. I would also like to thank J. B. Parker and J. Squire for help with the numerics, and M. W. Kunz for his steadfast and unwavering support of this research. I would like to thank the anonymous referees for constructive comments that led to a much improved presentation. Finally, I would like to thank the editor B. Dorland for offering additional background and historical context on the subject of this paper. This work was supported by U.S. DoE contract DE-AC02-09CH11466.