Enstrophy non-conservation and the forward cascade of energy in two-dimensional electrostatic magnetized plasma turbulence

Abstract

A fluid system is derived to describe electrostatic magnetized plasma turbulence at scales somewhat larger than the Larmor radius of a given species. It is related to the Hasegawa–Mima equation, but does not conserve enstrophy, and, as a result, exhibits a forward cascade of energy, to small scales. The inertial-range energy spectrum is argued to be shallower than a $-11/3$ power law, as compared to the $-5$ law of the Hasegawa–Mima enstrophy cascade. This property, confirmed here by direct numerical simulations of the fluid system, may help explain the fluctuation spectrum observed in gyrokinetic simulations of streamer-dominated electron-temperature-gradient driven turbulence (Plunk et al., Phys. Rev. Lett., vol. 122, 2019, 035002), and also possibly some cases of ion-temperature-gradient driven turbulence where zonal flows are suppressed (Plunk et al., Phys. Rev. Lett., vol. 118, 2017, 105002).

1 Introduction

The turbulent cascade, a mechanism for the nonlinear transfer of energy across scales, is a key idea for understanding kinetic magnetized plasma turbulence. By considering simplified models, in a uniform magnetic geometry, one can obtain a theoretical prediction for the spectrum of fluctuations, valid across an ‘inertial range’ of scales, free from energy sources and sinks. Although such a theory is not able to fully describe the behaviour of realistic turbulence, which hosts instabilities, damped modes, complicated magnetic geometries, etc., it nevertheless constitutes a quantitative prediction of nonlinear behaviour of the underlying gyrokinetic equation, an equation which generally governs actual systems of practical interest. The existence of such theoretical test cases is valuable for validating the solution methods employed by gyrokinetic codes, and as a foundation for physical interpretation of the volumes of data they produce.

Here, a novel quasi-two-dimensional fluid system is derived from the electrostatic gyrokinetic system, to describe fluctuations that predominantly vary in the directions perpendicular to the mean magnetic field, i.e. in the ‘drift plane’, at scales $\ell$ larger than the Larmor radius $\rho$, corresponding to a species of interest. The notion that quasi-two-dimensional behaviour might underly magnetized plasma turbulence is intuitively justified by the fact that the magnetic guide field renders the dynamics inherently anisotropic. Furthermore, instabilities that drive electrostatic turbulence in fusion plasmas, e.g. the ion-temperature-gradient (ITG) and electron-temperature-gradient (ETG) modes, exhibit a kind of localization along the field line, accompanied by the domination of perpendicular dynamics over parallel dynamics. The fluid limit $\ell \gg \rho$ is of particular importance, because the energy of the fluctuations is predominantly found at such scales – these are the scales of importance, most directly affecting the performance of fusion devices. Furthermore, the reduction of complexity afforded by fluid limits can reveal important features of the dynamics, that do not manifest in the analysis of the general gyrokinetic equations.

Although similar systems as the one presented here have been proposed and studied in the past, most notably the Hasegawa–Mima (HM) equation (Hasegawa & Mima 1978) the present derivation takes special care in considering the consequences of the appearance of nonlinear finite-Larmor-radius (FLR) terms that appear in the dynamical equation for the electrostatic potential – i.e. the ‘vorticity’ equation. Such terms introduce a closure problem in the fluid moment hierarchy, where lower moments are coupled to ever higher ones, generally without end. This motivates the cold ion limit that underlies the HM equation, which eliminates the inconvenient terms, but is, however, not generally appropriate for application to fusion plasmas. In the present work, it is noted that the presence of these terms introduces rapid dynamics, and a multiple-scale analysis is proposed in which the fluid moment hierarchy closes at the pressure moment, without using an ad hoc closure scheme, leading to a relatively simple system involving only two fields.

What is immediately apparent is that the presence of the additional field (the pressure perturbation) breaks the nonlinear conservation of enstrophy that is famously satisfied by the HM equation, and there causes an ‘inverse cascade’ of energy to large scales. The new system, we argue, should exhibit distinct nonlinear behaviour, including a shallower energy spectrum when the effect of the nonlinear FLR terms is sufficiently strong. Direct numerical simulation of the fluid model gives some confidence in these predictions. The results of this work may help to interpret observations of turbulence in parameter regimes where the dynamics tends towards the quasi-two-dimensional limit. We discuss possible examples, including cases explored in previous gyrokinetic turbulence simulations in tokamak and stellarator geometries.

2 Equations and definitions

We assume uniform magnetic geometry, where the magnetic guide field is constant and points in the $\hat {z}$-direction. One species is assumed to be kinetic, with the other species satisfying a simple Boltzmann response model. We begin with a non-dimensional form of the gyrokinetic system (Plunk et al. 2010), normalized relative to the kinetic species: $v_{\perp }/v_{\textrm {th}}\rightarrow v$ (with $v_{\textrm {th}}= \sqrt {T/m}$, and $T$ and $m$ are the temperature and mass of the kinetic species) is the normalized perpendicular velocity and the normalized wavenumber is $k_{\perp }\rho \rightarrow k$ where the thermal Larmor radius of the kinetic species is $\rho = v_{\textrm {th}}/\varOmega _c$ and $\varOmega _c = qB/m$, B is the magnitude of the magnetic field and m is the particle mass. The two-dimensional gyrokinetic equation is written as follows in terms of the perturbed gyrocentre distribution function $g({\boldsymbol R}, v, t)$, where ${\boldsymbol R} = \hat {\boldsymbol x} X + \hat {\boldsymbol y} Y$ is the gyrocentre position:

(2.1)\begin{equation} \frac{\partial g}{\partial t} + \left\{\left\langle\varphi\right\rangle_{\boldsymbol R},\;g\right\} = \left\langle C[h]\right\rangle_{\boldsymbol R},\end{equation}

where $\left \langle C [h]\right \rangle _{\boldsymbol R}$ is the collision operator (not treated here explicitly); the Poisson bracket is $\left \{A,\;B\right \} = \hat {\boldsymbol z}\times \boldsymbol {\nabla } A \boldsymbol {\cdot } \boldsymbol {\nabla } B = \partial _x A \partial _y B - \partial _y A \partial _x B$ and the gyro-average is defined $\left \langle A ({\boldsymbol r})\right \rangle _{\boldsymbol R} = ({1}/{2{\rm \pi} })\int _0^{2{\rm \pi} } \textrm {d}\vartheta A({\boldsymbol R} + \boldsymbol {\rho }(\vartheta ))$, where the Larmor-radius vector is $\boldsymbol {\rho }(\vartheta ) = \hat {{\boldsymbol z}}\times {\boldsymbol v} = v_{\perp }(\hat {{\boldsymbol y}}\cos {\vartheta } - \hat {{\boldsymbol x}}\sin {\vartheta })$ and $\vartheta$ is the gyro-angle. (Note that the quantity inside of the collision operator is $h = g + \left \langle \varphi \right \rangle _{\boldsymbol R}F_0$. Note also that the spatial coordinate is ${\boldsymbol R}$ in the gyrokinetic equation and, formally, the spatial derivatives are to be interpreted in this variable, but for simplicity we avoid making the distinction explicit.) We mostly ignore the collision operator but note that some mechanism of coarse graining will be necessary to get sensible solutions out of the equation. Quasi-neutrality yields the electrostatic potential $\varphi ({\boldsymbol r}, t)$, where ${\boldsymbol r} = \hat {\boldsymbol x} x + \hat {\boldsymbol y} y$ is the position-space coordinate

(2.2)\begin{equation} 2{\rm \pi}\int_0^{\infty} v \,\textrm{d} v \left\langle g\right\rangle_{\boldsymbol r} = (1 + \tau)\varphi - \varGamma_0\varphi,\end{equation}

where the $g$ is implicitly assumed to be integrated over parallel velocity so that $2{\rm \pi} \int _0^{\infty } v dv$ completes the integration over three-dimensional velocity space. The angle average is defined as $\left \langle A ({\boldsymbol R})\right \rangle _{\boldsymbol r} = ({1}/{2{\rm \pi} })\int _0^{2{\rm \pi} } \textrm {d}\vartheta A({\boldsymbol r} - \boldsymbol {\rho }(\vartheta ))$, and the term $\tau \varphi$ is the adiabatic density response, with $\tau = T_i/(Z T_e)$ for the case of ion scales and $\tau = Z T_e/ T_i$ for the case of electron scales. For the ion case, this Boltzmann response might be considered reasonable if zonal flows are strongly suppressed. The operator $\varGamma _0\phi = 2{\rm \pi} \int _0^{\infty } v \,\textrm {d} v F_0(v)\left \langle \left \langle \phi \right \rangle _{\boldsymbol R}\right \rangle _{\boldsymbol r}$ is more naturally expressed in Fourier space, assuming a Maxwellian background $F_0 = \exp [-v^2/2]/(2{\rm \pi} )$, i.e. $\varGamma _0\varphi = \sum _{\boldsymbol k} \exp (\textrm {i} {\boldsymbol k}\boldsymbol {\cdot }{\boldsymbol r}) \hat {\varGamma }_0 \hat {\varphi }$, with

(2.3)\begin{equation} \hat{\varGamma}_0(k) = \int_0^{\infty} v \,\textrm{d} v \,\textrm{e}^{-v^2/2}J_0^2(kv) = I_0(k^2)\,\textrm{e}^{-k^2},\end{equation}

where $I_0$ is the zeroth-order modified Bessel function.

3 Fluid limit

We expand in the limit

(3.1)\begin{equation} \delta = k^2 \ll 1, \end{equation}

i.e. we assume that the scales of interest are larger than the Larmor radius of the species of interest. For electron scales, the limit is considered subsidiary to the adiabatic ion limit, so scales of the turbulence must remain much smaller than the ion Larmor scale, i.e. $\rho _e/\rho _i \ll k \ll 1$. Note that there may also be a minimum applicable $k$ imposed by the dynamics parallel to the magnetic field, but treating this explicitly is outside the scope of this work. We will only need the first two orders of the expansion in $\delta$. The gyrokinetic equation, henceforth omitting explicit collisional effects, is written as

(3.2)\begin{equation} \frac{\partial g}{\partial t} + \left\{\left(1+ \frac{v^2}{4} \nabla^2\right)\varphi,\;g\right\} \approx 0,\end{equation}

and (2.2) becomes

(3.3)\begin{equation} \tau\varphi - \nabla^2\varphi \approx 2{\rm \pi} \int_0^{\infty} v\,\textrm{d} v \left(1+ \frac{v^2}{4} \nabla^2\right) g.\end{equation}

We will denote $v$-moments of $g$ as

(3.4)\begin{equation} G_n = 2{\rm \pi} \int v\,\textrm{d} v \left(\frac{v}{2}\right)^n g. \end{equation}

3.1 Naive expansion

To give a taste for the issues that arise in the expansion, let us take an initial informal look at the moments of the gyrokinetic equation. We first examine the density moment. We include only the ostensibly dominant nonlinear terms. We stress that this equation is given only for illustrative purposes, and is not to be taken as a basis for the later derivations of the paper

(3.5)\begin{equation} \frac{\partial}{\partial t}\left(\tau \varphi -\nabla^2\varphi - \nabla^2 G_2\right) + \left\{\varphi,\;-\nabla^2 \varphi -\nabla^2 G_2 \right\} + \left\{G_2,\; -\nabla^2 \varphi\right\} = 0. \end{equation}

Note that the term $-\partial _t \nabla ^2 \varphi$, which appears in the HM equation, should be neglected here because it is formally smaller than $\partial _t \varphi$ by one power of the ordering parameter $\delta$. Likewise, the term $-\partial _t \nabla ^2 G_2$ must be considered negligible if the ordering $G_2 \sim \varphi$ and $\partial _t G_2 \sim \partial _t \varphi$ holds. The situation is, however, a bit more subtle. The above equation couples to the $v^2$ moment of $g$, $G_2$ and the equation for this and other such moments can be written, neglecting higher-order FLR terms, as

(3.6)\begin{equation} \frac{\partial G_n}{\partial t} + \left\{\varphi,\;G_n\right\} = 0.\end{equation}

What we now notice, examining these two equations, is that the density equation is driven by nonlinear terms that appear to be much smaller than those controlling the higher moments of the distribution function – that is, the equations for $G_n$ have dominant contributions from the ‘$E\times B$’ nonlinearity, while the potential evolves under the influence of terms like the ‘polarization drift’ nonlinearity, which is smaller by a factor of $\delta = k^2$. One possible resolution of this apparent imbalance is to consider $G_n$ itself to be large, as for instance in the non-resonant limit of the ITG/ETG mode, i.e. $G_n \sim \delta ^{-1} \varphi$. In this case, the polarization drift nonlinearity can be neglected, and we see the justification for retaining the additional time derivative term above, since $\partial _t\varphi \sim \partial _t \nabla ^2 G_2$. This term can be evaluated from the Laplacian of (3.6), yielding

(3.7)\begin{equation} \tau \frac{\partial \varphi}{\partial t} + \nabla^2\left\{\varphi,\;G_2\right\} + \left\{\varphi,\;-\nabla^2 G_2 \right\} + \left\{G_2,\; -\nabla^2 \varphi\right\} = 0.\end{equation}

Equations (3.6)–(3.7) demonstrate a consistent fluid limit, but cannot describe ITG or ETG turbulence in the resonant limit, where $\varphi \sim G_2$. This corresponds the more physically reasonable scenario of a modest turbulence drive – i.e. not very far from the linear critical gradient, or considering the weakly unstable, large-scale modes that dominate the turbulence spectrum. To treat this limit properly, we must account for the fact that $\varphi$ evolves much more slowly than $G_n$. Physically, it can be argued that, in a turbulent state, (3.6) will then describe rapid mixing of $G_n$ by $E\times B$ vortices, so that any initial variation along streamlines of constant $\varphi$ will decay on a fast time scale (with the help of some explicit dissipation), leaving $G_n$ to be constant along those streamlines (S.C. Cowley, private communication 2008). To account for such processes more carefully, we abandon conventional perturbation theory in favour of the method of multiple scales (see, i.e. Bender & Orszag 1978). We will henceforth disregard the equations presented here, in § 3.1, and proceed to derive equations that contain only terms justified by a set of explicitly stated ordering assumptions.

3.2 Multiscale expansion

We introduce the fast and slow time variables $t_\mathrm {f}$, and $t_\mathrm {s}$, such that $\partial _{t_\mathrm {f}} \sim \left \{\varphi ,\;.\right \}$ and $\partial _{t_\mathrm {s}} \sim \left \{\nabla ^2\varphi ,\;.\right \} \sim \delta \partial _{t_\mathrm {f}}$ and expand the fields as

(3.8)\begin{equation} \varphi = \varphi^{(0)}(t_{\textrm {s}}, t_{\textrm {f}}, x, y) + \varphi^{(1)}(t_{\textrm {s}}, t_{\textrm {f}}, x, y) + \ldots, \end{equation}

(3.9)\begin{equation} G_n = G_n^{(0)}(t_\mathrm{s}, t_\mathrm{f}, x, y) + G_n^{(1)}(t_\mathrm{s}, t_\mathrm{f}, x, y) + \ldots, \end{equation}

where $\varphi ^{(m+1)}/\varphi ^{(m)} \sim {\mathcal {O}}(\delta )$, etc. We reiterate that the assumptions we have made are $\delta \ll 1$, the above multi-scale expansion, and the quasi-two-dimensional approximation, whereby variation in the direction along the magnetic field is neglected, and the non-kinetic species is assumed to follow a Boltzmann distribution, implying (2.2); no further approximations will be made in this section. We proceed to examine the moments of gyrokinetic equation, order by order. The density moment at dominant order in $\delta$ is

(3.10)\begin{equation} \frac{\partial \varphi^{(0)}}{\partial t_\mathrm{f}} = 0, \end{equation}

from which we formally establish that $\varphi ^{(0)}$ depends only on the slow time variable. At next order in $\delta$ we obtain

(3.11)\begin{equation} \tau \frac{\partial \varphi^{(0)}}{\partial t_\mathrm{s}} + \tau\frac{\partial \varphi^{(1)}}{\partial t_\mathrm{f}} - \frac{\partial}{\partial t_\mathrm{f}}{{\nabla}^2} G_2^{(0)} + \left\{\varphi^{(0)}, -{{\nabla}^2} \varphi^{(0)} -{{\nabla}^2} G_2^{(0)} \right\} + \left\{G_2^{(0)}, -{{\nabla}^2} \varphi^{(0)}\right\} = 0.\end{equation}

We introduce a time-average operator to extract the smooth-time behaviour from this equation

(3.12)\begin{equation} {\overparen{A}} = \frac{1}{{\rm \Delta} t} \int_{t_\mathrm{f}-{\rm \Delta} t/2}^{t_\mathrm{f}+ {\rm \Delta} t/2}\,\textrm{d} t_\mathrm{f}^\prime A(t_\mathrm{f}^\prime).\end{equation}

This time average extends over a period of time much longer than the short time scale (${\rm \Delta} t^{-1} \ll \left \{\varphi ,\;.\right \}$). Applying this average to (3.11), we obtain

(3.13)\begin{equation} \tau \frac{\partial \varphi^{(0)}}{\partial t_\mathrm{s}} + \left\{\varphi^{(0)},\;-\nabla^2 \varphi^{(0)} -\nabla^2 \overparen{G}_2^{(0)} \right\} + \left\{\overparen{G}_2^{(0)},\; -\nabla^2 \varphi^{(0)}\right\} = 0.\end{equation}

The dominant-order equation for $G_n$ is

(3.14)\begin{equation} \frac{\partial G_n^{(0)}}{\partial t_\mathrm{f}} + \left\{\varphi,\;G_n^{(0)}\right\} = 0,\end{equation}

from which, upon time averaging, we conclude that $\overparen {G}_n^{(0)}$ is constant along closed streamlines of constant $\varphi$. Informally, we will say that $\overparen {G}_n^{(0)}$ is a function of $\varphi$, although it can be multi-valued. For $n = 2$ we adopt the notation

(3.15)\begin{equation} \overparen{G}_2^{(0)} = \chi(\varphi, t_\mathrm{s}). \end{equation}

Note that, formally, we must exclude special points and lines where $\boldsymbol {\nabla } \varphi = 0$ (o-points, and the ‘separatrices’ that include x-points) but these should occupy negligible volume in the $x$–$y$ plane. At the next order, we will obtain the smooth evolution of $G_n$,

(3.16)\begin{equation} \frac{\partial G_n^{(0)}}{\partial t_\mathrm{s}} + \frac{\partial G_n^{(1)}}{\partial t_\mathrm{f}} + \left\{\varphi^{(0)},\;G_n^{(1)}\right\} + \left\{\varphi^{(1)},\;G_n^{(0)}\right\} + \left\{\nabla^2 \varphi^{(0)},\;G_{n+ 2}^{(0)}\right\} = 0.\end{equation}

To this equation we apply two averages, the time average, and also an average along streamlines of constant $\varphi$. To define this average, we introduce a coordinate $s$ which parameterizes these streamlines and satisfies $\hat {\boldsymbol z}\times \boldsymbol {\nabla } \varphi \boldsymbol {\cdot } \boldsymbol {\nabla } s = 1$ for convenience. Then we define

(3.17)\begin{equation} \left\langle A\right\rangle_{s} = \frac{\oint \textrm{d} s A(s)}{\oint \textrm{d} s}. \end{equation}

The integral over $s$ is either closed in the sense that the streamlines are closed, or effectively closed by periodic boundary conditions. The second term of (3.16) is annihilated by the time average. Noting that $\left \{F(\varphi ),\; A\right \} = \partial _s(A F')$, the third term is zero under the $s$-average, as is the last term after time average, since $\overparen {G}_{n + 2}^{(0)}$ is a function of $\varphi ^{(0)}$ by (3.14). This is a crucial cancellation since the fluid moment hierarchy is consequently shown to be closed.

It is convenient to now introduce notation for the part of a field that varies on the fast time scale, i.e. the ‘fluctuating part’, complementing the mean component defined by (3.12)

(3.18)\begin{equation} \tilde{A} = A - \overparen{A}.\end{equation}

What results from the double average of (3.16) can then be expressed

(3.19)\begin{equation} \left\langle\frac{\partial \overparen{G}_n^{(0)}}{\partial t_\mathrm{s}}\right\rangle_{s} + \left\langle\overparen{\left\{\tilde{\varphi}^{(1)},\;\tilde{G}_n^{(0)}\right\}}\right\rangle_{s} = 0.\end{equation}

To evaluate the nonlinear term of (3.19), we must obtain dynamical equations for the fluctuating fields $\tilde {\varphi }^{(1)}$ and $\tilde {G}_n^{(0)}$. These come from (3.11) and (3.14), respectively. The fluctuating part of (3.14) is

(3.20)\begin{equation} \frac{\partial \tilde{G}_n^{(0)}}{\partial t_\mathrm{f}} + \left\{\varphi^{(0)},\;\tilde{G}_n^{(0)}\right\} = 0.\end{equation}

Subtracting (3.13) from (3.11), and using the Laplacian of (3.20) to evaluate $\partial _{t_\mathrm {f}}\nabla ^2 \tilde {G}_n^{(0)}$, we find

(3.21)\begin{equation} \tau\frac{\partial \tilde{\varphi}^{(1)}}{\partial t_\mathrm{f}} + \nabla^2\left\{\varphi^{(0)},\;\tilde{G}_2^{(0)}\right\} + \left\{\varphi^{(0)},\;-\nabla^2 \tilde{G}_2^{(0)} \right\} + \left\{\tilde{G}_2^{(0)},\; -\nabla^2 \varphi^{(0)}\right\} = 0.\end{equation}

Finally, noting that $\left \langle \partial _{t_\mathrm {s}}\varphi ^{(0)}\right \rangle _{s} = 0$ (from (3.13)), we obtain, from (3.19), an expression determining the explicit time dependence of $\chi$

(3.22)\begin{equation} \left(\frac{\partial \chi}{\partial t_\mathrm{s}}\right)_{\varphi} + \left\langle\overparen{\left\{\tilde{\varphi}^{(1)},\;\tilde{G}_2^{(0)}\right\}}\right\rangle_{s} = 0,\end{equation}

where the partial time derivative is taken at constant $\varphi ^{(0)}$. To summarize, (3.20) and (3.21) are the fast-time equations that determine $\tilde {\varphi }^{(1)}$ and $\tilde {G}_n^{(0)}$, which can be substituted into the slow-time equation (3.22) for $\chi$, and coupled with the following equation (a repetition of (3.13) written in terms of $\chi$) to close the system

(3.23)\begin{equation} \tau \frac{\partial \varphi^{(0)}}{\partial t_\mathrm{s}} + \left\{\varphi^{(0)},\;-\nabla^2 \varphi^{(0)} -\nabla^2 \chi \right\} + \left\{\chi,\; -\nabla^2 \varphi^{(0)}\right\} = 0.\end{equation}

Noting that $\overparen {\varphi } \approx \varphi ^{(0)}$ and $\tilde {\varphi } \approx \tilde {\varphi }^{(1)}$, we can, without introducing ambiguity, simply drop the superscripts in what follows.

The final system, (3.20)–(3.23), has some noteworthy features. First, (3.22) has the appearance of a heat transport equation, where the flux is carried by the rapidly varying pressure perturbation $\tilde {G}_2$ and the small amplitude fluctuating potential $\tilde {\varphi }$. Note also how (3.20) and (3.21) bear a strong resemblance to the fluid system given by (3.6)–(3.7), where a similar ordering is satisfied, namely $\tilde {\varphi } \ll \tilde {G}_2$.

4 Decaying turbulence

We will avoid the complications introduced by the instabilities that physically drive turbulence, and instead now consider decaying turbulence. (We note that a linear instability could be added to this fluid system using the non-resonant limit of the toroidal branch of the ITG or ETG mode, but this would require some care to maintain consistency with the ordering assumptions, as discussed in § 3.1.) Let us consider periodic boundary conditions, and include explicit dissipation using a fourth-order hyperviscosity term. Without drive terms, (3.20) implies the rapid decay of $\tilde {G}_n$ to zero, implying $(\partial _{t_\mathrm {s}}\chi )_{\varphi } = 0$ (i.e. it only depends on the time via its dependence on $\varphi$). The variation of $\chi$ in $\varphi$ (or more formally, its variation between distinct lines of constant $\varphi$) can be considered as an initial condition of our calculation. We need only then solve a single equation, which, neglecting superscripts for order and the subscripts of the time variable $t_\mathrm {s}$, becomes

(4.1)\begin{equation} \tau \frac{\partial \varphi}{\partial t} + \left\{\varphi,\;-\nabla^2 \varphi -\nabla^2 \chi \right\} + \left\{\chi,\; -\nabla^2 \varphi\right\} = \nu_4 \nabla^4 \varphi.\end{equation}

The electrostatic energy

(4.2)\begin{equation} E = \frac{\tau}{2}\int \textrm{d} x \,\textrm{d} y \varphi^2, \end{equation}

is conserved by the nonlinearity for arbitrary $\chi$, which can be verified by multiplying the equation by $\varphi$ and integrating over the $x$–$y$ domain. Note that the resulting integral of the final nonlinear term of (4.1) can be rewritten as $- \int \textrm {d} x \,\textrm {d} y {\boldsymbol v}_E\boldsymbol {\cdot }\boldsymbol {\nabla }(\varphi \chi ^\prime \nabla ^2\varphi )$, with ${\boldsymbol v}_E = \hat {\boldsymbol z}\times \boldsymbol {\nabla }\varphi$, which is zero using $\boldsymbol {\nabla }\boldsymbol {\cdot }{\boldsymbol v}_E = 0$ and periodicity.

Another quantity of interest is the enstrophy, which we will define here as

(4.3)\begin{equation} Z = \frac{\tau}{2}\int \textrm{d} x \,\textrm{d} y |\boldsymbol{\nabla}\varphi|^2. \end{equation}

The enstrophy balance equation is found by multiplying (4.1) by $-\nabla ^2\varphi$ and integrating over $x$ and $y$. Note that the presence of $\chi$ in the equation breaks enstrophy conservation if $\chi$ is a nonlinear function of $\varphi$. The nonlinear invariance of $Z$ is associated with the inverse cascade of energy in Hasegawa–Mima turbulence. We thus expect to recover the spectra corresponding to the potential limit of the Hasegawa–Mima equation (i.e. where $\nabla ^2 \varphi \ll \varphi$; see Plunk et al. 2010), if $\chi$ is small, and qualitatively different cascade when $\chi$ is sufficiently large.

The Hasegawa–Mima spectra can be derived in the rough ‘phenomenological’ style, in terms of the fluctuation amplitude at scale $\ell$, denoted $\varphi _\ell$, by assuming constancy of the nonlinear flux of its nonlinear invariants (see e.g. Frisch 1995; Plunk et al. 2010). For the forward cascade, i.e. at scales smaller than the scale of energy injection, the enstrophy flux, denoted $\varepsilon _Z$, is assumed constant (independent of scale $\ell$), which is expressed as follows:

(4.4)\begin{equation} \varepsilon_Z = \tau_{\mathrm{NL}}^{-1} \ell^{-2} \varphi_\ell^2 \sim \varphi_\ell^3\ell^{-6}, \end{equation}

with $\tau _{\mathrm {NL}}(\ell )$ denoting the nonlinear turnover time. This leads to the scaling $\varphi _\ell \sim \ell ^2\varepsilon _Z^{1/3}$, implying a one-dimensional energy spectrum of $E(k) \sim k^{-5}$. The constancy of the scale-by-scale flux of energy, expected for the inverse cascade at scales larger than the injection scale, is expressed as

(4.5)\begin{equation} \varepsilon_E = \tau_{\mathrm{NL}}^{-1} \varphi_\ell^2 \sim \varphi_\ell^3\ell^{-4}, \end{equation}

implying $\varphi _\ell \sim \ell ^{4/3}\varepsilon _E^{1/3}$ and a spectrum $E(k) \sim k^{-11/3}$.

Because the additional nonlinear terms of (4.1), henceforth called the ‘$\chi$ nonlinearity’, formally break enstrophy conservation, we expect that if they are sufficiently strong, the inverse cascade should be eliminated and the forward cascade of $Z$ replaced with a direct cascade of $E$. If this flux is carried by the HM nonlinearity, one might expect to observe the spectrum $E(k) \sim k^{-11/3}$, as suggested by Plunk (2019). On the other hand, balancing the $\chi$ nonlinearity with the HM nonlinearity, scale by scale, implies the linear relation $\chi _\ell \sim \varphi _\ell$, i.e. $\chi \propto \varphi$, which would imply that enstrophy is actually a nonlinear invariant, preventing the forward cascade of $E$. For this reason, we may expect to observe an energy spectrum distinct from $k^{-11/3}$, whose steepness depends on the relationship between $\chi _\ell$ and $\varphi _\ell$, which itself depends on details of the turbulence.

Providing a definitive prediction of this relationship is beyond the scope of the present work, but a power law seems to be a reasonable possibility to explore, i.e. $\chi _\ell \sim \varphi _\ell ^\alpha$. Note that any super-linear scaling $\alpha > 1$ should lead to a spectrum shallower than $k^{-11/3}$, while a sub-linear scaling $\alpha < 1$ would imply $\chi$ is not analytic in $x$ and $y$. The fluctuating fields $\tilde {G}_2$ and $\tilde {\varphi }$ could be especially active in regions of low $E\times B$ shear (see (3.20)), causing local extrema in the function $\chi$, via (3.22), so that a quadratic relationship prevails in such regions, $\chi _\ell \sim \varphi _\ell ^2$. Whether or not this seems plausible, assuming a simple nonlinear relationship will allow us to make the discussion now more concrete; qualitatively similar conclusions should apply for all $\alpha > 1$. Let us consider the following form for $\chi$:

(4.6)\begin{equation} \chi(\varphi) = \frac{\lambda}{2} \varphi^2.\end{equation}

The nonlinear energy flux by the $\chi$ terms is then expressed as $\varepsilon _E \sim \lambda \varphi _\ell ^4 \ell ^{-4}$, implying $\varphi _\ell \sim \ell (\varepsilon _E/\lambda )^{1/4}$, and the corresponding energy spectrum

(4.7)\begin{equation} E(k) \sim k^{-3}. \end{equation}

This spectrum should prevail in cases where the $\chi$-nonlinearity dominates (e.g. large $\lambda$). At sufficiently low $\lambda$, one expects a return to the HM behaviour, implying $E(k) \sim k^{-5}$ for the forward cascade.

Some sort of hybrid behaviour may also be possible, although the broad scale range needed for clear observation of this may be not be present for realistic conditions encountered in fusion plasmas. One might argue that, because the amplitude of fluctuations $\varphi _\ell$ is generally expected to decrease as scales do, the cubic nonlinearity should be dominant at large scales, and subdominant at small scales. Thus, for sufficiently large $\lambda$, the energy cascade scaling $\varphi _\ell \sim \ell (\varepsilon _E/\lambda )^{1/4}$ should hold from the injection scale, down to a transition scale, which can be found by balancing the HM nonlinearity with the $\chi$-nonlinearity, i.e. $\varphi _\ell ^2\ell ^{-4} \sim \lambda \varphi _\ell ^3 \ell ^{-4}$. Defining the outer scale $\ell _\mathrm {o}$ as the scale of energy injection (or initial energy containing scale), and $\varphi _{\mathrm {o}} = \varphi _{\ell _\mathrm {o}}$, we can write the $\varphi _\ell$ scaling as $\varphi _\ell \sim (\ell /\ell _\mathrm {o}) \varphi _{\mathrm {o}}$, so that the above balance occurs at the ‘transition’ scale $\ell _t \sim \ell _\mathrm {o} / (\lambda \varphi _{\mathrm {o}})$. Thus, if $\lambda \varphi _{\mathrm {o}} \gtrsim 1$ one might expect $E \sim k^{-3}$ scaling for $\ell _{\mathrm {o}}^{-1} < k < \ell _t^{-1}$ followed by $E \sim k^{-5}$ for $k > \ell _t^{-1}$.

4.1 Direct numerical simulations

To explore the behaviour of the model, (4.1), and test the theoretical predictions, we perform direct numerical simulations, assuming the simple quadratic form of $\chi (\varphi )$ in (4.6). This introduces a nonlinearity that is cubic in $\varphi$, which can be treated pseudo-spectrally using a padding factor of $2$ for dealiasing; higher order nonlinearities require additional padding (Hossain, Matthaeus & Ghosh 1992). The boundary conditions for the simulations are periodic in $x$ and $y$, and $\tau = 1$ for all simulations.

Figure 1 compares the simulation results with the theoretical scaling laws. All simulations are initialized with randomly phased fluctuation amplitude of $\varphi \sim 1$ around $k = 1$, falling off exponentially at higher $k$. Note that although the model assumes $k \ll 1$ there is no conflict in using $k > 1$ for the simulations, as scaling symmetries of the model allow the results to be reinterpreted for $k \ll 1$. The spectrum found for the $\lambda = 0$ case is roughly consistent with the theoretical power law $k^{-5}$ expected for the potential limit of the HM equation. We note that similar results (not shown here) are encountered for $\lambda \lesssim 0.1$. At larger $\lambda$, the breaking of enstrophy conservation is indeed observed in the time trace of $Z$, as the energy fills in the spectrum at large $k$. For the case labelled $\lambda \rightarrow \infty$ in figure 1, the spectrum seems consistent with the theoretical $k^{-3}$ prediction at scales smaller than the injection scale. Note that this limit is obtained by actually setting $\lambda = 1$ and simply removing the HM nonlinearity (i.e. the first term of (4.1)) from the equation, as can be formally justified by rescaling (4.1) in the limit $\lambda \rightarrow \infty$. Similar behaviour is observed for $\lambda \gtrsim 1$. Intermediate values of $\lambda$ show intermediate behaviour.

Figure 1. Comparison of spectra exhibited by HM-type system ($\lambda = 0$) and our two-dimensional turbulence model ($\lambda \rightarrow \infty$).

One example is shown in figure 2, which seems to show evidence of a transition scale between the two theoretical power laws, giving some support to the predictions of a hybrid scenario described theoretically in the previous section. A more extensive set of simulations would be needed to test the predictions in detail, for instance the dependence of the transition scale $\ell _t$ on system parameters. We would like to generally stress that the results of all of the numerical simulations presented here come at a very modest computational expense, and larger scale computational effort, especially using a gyrokinetic code, could offer a more extensive test of the conclusions of this work.

Figure 2. Energy spectrum for a case of intermediate strength of $\chi$ nonlinearity ($\lambda = 0.5$).

5 Discussion

A novel fluid system has been derived to describe the behaviour of certain classes of quasi-two-dimensional electrostatic magnetized plasma turbulence. A possible application is to describe the energy cascade in cases of streamer-dominated ETG turbulence (note the spectrum, noted to be close to $k^{-11/3}$, in figure 5 of Plunk 2019), where the nonlinear stability of elongated turbulent eddies is believed to stem from the two-dimensional character of the dominant instabilities, e.g. the absence of sufficient variation of the mode structure in the direction along the magnetic field (Jenko & Dorland 2002). This turbulent state is, however, sensitive to magnetic geometry, and seems to vanish when, for instance, the global magnetic shear is varied in such a way as to induce stronger parallel electron flow to the ETG mode. The ensuing dynamics then depends on kinetic physics involving the parallel streaming term, absent from two-dimensional models. A second possible application of the present model might be to describe ITG turbulence in cases where the zonal flows are suppressed. One candidate is a case observed with simulations of the HSX stellarator having surprisingly steep fluctuation spectra (Plunk, Xanthopoulos & Helander 2017), found to be close to $k^{-10/3}$.

Although the presented model has limited application, it fills a significant gap in present theories describing gyrokinetic turbulence cascades, as it accounts for the essential nonlinear terms that arise when the cold ion approximation is invalid. These terms, it is found, alter the conservative properties of the nonlinearity, with significant consequences on the cascade, so that, even in the two-dimensional limit, the inverse cascade of energy can be shut down. The numerical simulations confirm that the size of the pressure perturbation ($\chi$) can control the cascade type, and HM-like behaviour can be recovered if it is sufficiently small. This may underlie the slow secular growth of large-scale zonal flows (Guttenfelder & Candy 2011) and other coherent structures (Nakata et al. 2010) in simulations of ETG turbulence, and the related appearance of a Dimits shift phenomenon in near-marginal cases (Colyer et al. 2017).