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. Reference Plunk, Cowley, Schekochihin and Tatsuno2010), 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 Reference Bender and Orszag1978). 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. Reference Plunk, Cowley, Schekochihin and Tatsuno2010), 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 Reference Frisch1995; Plunk et al. Reference Plunk, Cowley, Schekochihin and Tatsuno2010). 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 (Reference Plunk2019). 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 Reference Hossain, Matthaeus and Ghosh1992). 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$).