Hostname: page-component-586b7cd67f-t7czq Total loading time: 0 Render date: 2024-11-20T18:25:31.804Z Has data issue: false hasContentIssue false

Scale invariance and critical balance in electrostatic drift-kinetic turbulence

Published online by Cambridge University Press:  25 July 2023

T. Adkins*
Affiliation:
Rudolf Peierls Centre for Theoretical Physics, University of Oxford, Oxford OX1 3PU, UK Merton College, Oxford OX1 4JD, UK United Kingdom Atomic Energy Authority, Culham Centre for Fusion Energy, Abingdon OX14 3DB, UK
P.G. Ivanov
Affiliation:
Rudolf Peierls Centre for Theoretical Physics, University of Oxford, Oxford OX1 3PU, UK
A.A. Schekochihin
Affiliation:
Rudolf Peierls Centre for Theoretical Physics, University of Oxford, Oxford OX1 3PU, UK Merton College, Oxford OX1 4JD, UK
*
Email address for correspondence: [email protected]
Rights & Permissions [Opens in a new window]

Abstract

The equations of electrostatic drift kinetics are observed to possess a symmetry associated with their intrinsic scale invariance. Under the assumptions of spatial periodicity, stationarity, and locality, this symmetry implies a particular scaling of the turbulent heat flux with the system's parallel size, from which its scaling with the equilibrium temperature gradient can be deduced under some additional assumptions. This macroscopic transport prediction is then confirmed numerically for a reduced model of electron-temperature-gradient-driven turbulence in slab geometry. The system realises this scaling through a turbulent cascade from large to small perpendicular spatial scales. The route of this cascade through wavenumber space (i.e. the relationship between parallel and perpendicular scales in the inertial range) is shown to be determined by a balance between nonlinear-decorrelation and parallel-dissipation timescales. This type of ‘critically balanced’ cascade, which maintains a constant energy flux despite the presence of parallel dissipation throughout the inertial range (as well as order-unity dissipative losses at the outer scale) is expected to be a generic feature of plasma turbulence. The outer scale of the turbulence, on which the turbulent heat flux depends, is determined by the breaking of drift-kinetic scale invariance due to the existence of large-scale parallel inhomogeneity (the parallel system size).

Type
Research Article
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright
Copyright © The Author(s), 2023. Published by Cambridge University Press

1 Introduction

In many plasmas, energy is injected into the system on some large, system-specific macroscale (the ‘outer scale’). In order for such a system to reach a steady state, this energy must be dissipated. The usual route to this dissipation in kinetic plasmas is via a turbulent cascade of this energy to fine scales in both position and velocity space, where it is eventually thermalised by collisions (the ‘inner scale’). Given that there is often a large separation between these inner and outer scales (such as in, e.g. astrophysical systems, where energy is often injected by magnetohydrodynamic (MHD) instabilities), many studies of plasma turbulence are able to consider the dynamics of this turbulent cascade separately from the specific mechanisms of injection, simply assuming that there is some energy arriving from large scales that needs to be processed (see, e.g. Schekochihin et al. (Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009) and references therein).

There are, however, a variety of plasma systems for which such a scale separation is not a priori obvious. This is often due to the existence of gradients associated with an equilibrium (whether gravitational or magnetic) that, thermodynamically speaking, provide sources of free energy for unstable, microscale perturbations that can engender a turbulent cascade well below the usual macroscopic outer scale. In fact, the most (linearly) unstable perturbations in such systems often occur at the smallest scales. This is the case in tokamaks, in which the turbulent heat and particle transport is dominated by the (microscale) instabilities driven by the gradients of the plasma pressure between the inner core of the tokamak and its edge. The most important of these instabilities are the ion-temperature gradient (ITG) (see, e.g. Waltz Reference Waltz1988; Cowley, Kulsrud & Sudan Reference Cowley, Kulsrud and Sudan1991; Kotschenreuther et al. Reference Kotschenreuther, Dorland, Beer and Hammett1995a) and electron-temperature gradient (ETG) ones (see, e.g. Liu Reference Liu1971; Lee et al. Reference Lee, Dong, Guzdar and Liu1987; Dorland et al. Reference Dorland, Jenko, Kotschenreuther and Rogers2000; Jenko et al. Reference Jenko, Dorland, Kotschenreuther and Rogers2000). The relationship between the macroscopic scales associated with the plasma equilibrium and the microscopic scales on which turbulent fluctuations grow – and how the interaction between these two scales determines the heat and particle transport properties of the confined plasma – remains a topic of active research and great consequence.

In this paper, we consider electrostatic, drift-kinetic plasma turbulence – applicable to many regimes of tokamak operation – with a particular focus on the connection between its macroscopic transport properties and microscale dynamics. In the presence of constant perpendicular equilibrium gradients, it is observed that the equations of electrostatic drift kinetics possess a symmetry associated with their intrinsic scale invariance, in both the collisionless and collisional limits. We then show that this symmetry implies a particular scaling of the turbulent heat flux with equilibrium-scale quantities, in particular the parallel system size, provided one can assume spatial periodicity, stationarity (that the system has reached a statistical steady state), and locality (that the heat flux is independent of the system's perpendicular size, as it should be for any valid local model of plasma turbulence, provided its perpendicular size is large enough). This macroscopic transport prediction is then confirmed numerically in the context of an electron-scale, collisional model of electrostatic turbulence driven by the ETG instability in slab geometry. The choice to focus on ETG-driven turbulence was motivated, in part, by the fact that, despite significant recent progress (see, e.g. Ren et al. Reference Ren, Belova, Gorelenkov, Guttenfelder, Kaye, Mazzucato, Peterson, Smith, Stutman and Tritz2017; Hatch et al. Reference Hatch, Kotschenreuther, Mahajan, Merlo, Field, Giroud, Hillesheim, Maggi, Perez von Thun and Roach2019; Parisi et al. Reference Parisi, Parra, Roach, Giroud, Dorland, Hatch, Barnes, Hillesheim, Aiba and Ball2020; Guttenfelder et al. Reference Guttenfelder, Groebner, Canik, Grierson, Belli and Candy2021, Reference Guttenfelder, Battaglia, Belova, Bertelli, Boyer, Chang, Diallo, Duarte, Ebrahimi and Emdee2022; Chapman-Oplopoiou et al. Reference Chapman-Oplopoiou, Hatch, Field, Frassinetti, Hillesheim, Horvath, Maggi, Parisi, Roach and Saarelma2022; Parisi et al. Reference Parisi, Parra, Roach, Hardman, Schekochihin, Abel, Aiba, Ball, Barnes and Chapman-Oplopoiou2022; Field et al. Reference Field, Chapman-Oplopoiou, Connor, Frassinetti, Hatch, Roach and Saarelma2023), the saturation of such turbulence remains significantly less well-understood than its ITG cousin.

Further consideration of the microscale dynamics of our system of equations reveals that this heat flux scaling is enabled by a critically-balanced, Kolmogorov (Reference Kolmogorov1941) style cascade of energy from large to small spatial scales. The (approximately) constant flux of energy is that which survives the parallel dissipation present at the largest scales due, in our model, to thermal conduction. The existence of this parallel dissipation is also shown to play a key role in determining the saturated state of the system, limiting the cascade of free energy in wavenumber space. The outer scale of the turbulence is found to be determined by the breaking of the drift-kinetic scale invariance due to the existence of some large-scale parallel inhomogeneity, viz., the parallel system size, rather than by the smallest scales on which the ETG instability's growth rate peaks. It is thus the largest scales that are the most important in determining the saturated amplitudes to which the fluctuations grow, and the resultant turbulent transport. This is the first detailed demonstration of a critically balanced cascade in a temperature-gradient-driven plasma system since Barnes, Parra & Schekochihin (Reference Barnes, Parra and Schekochihin2011) proposed such a cascade for ITG turbulence.

The rest of this paper is organised as follows. In § 2, the scaling of the turbulent heat flux with parallel system size is derived from considerations of the scale invariance of the electrostatic drift-kinetic system of equations. Our model system of fluid equations is introduced in § 3, and the aforementioned heat flux scaling is verified in § 4. The dynamics of the inertial range are considered extensively in § 5, including the free-energy budget (§ 5.1), the existence of a constant-flux cascade and dynamical critical balance (§ 5.2), the nature of the outer scale (§ 5.3), the two-dimensional (2D) $(k_\perp, k_\parallel )$ spectra (§ 5.4), and the perpendicular isotropy in wavenumber space (§ 5.5). Lastly, we summarise our results and generic conclusions in § 6, and discuss the limits of their applicability to plasma systems in which finite-Larmor-radius (FLR) or electromagnetic effects are thought to be important.

2 Electrostatic drift-kinetic scale invariance

For systems adequately described by electrostatic drift kinetics, the heat flux through some volume $V$ is given by

(2.1a,b)\begin{equation} Q = \sum_s Q_s, \quad Q_s = n_{0s} T_{0s}\int \frac{\mathrm{d}^3 {\boldsymbol{r}}}{V} \left( {\boldsymbol{v}}_E {\boldsymbol{\cdot}} {\boldsymbol{\nabla}} x \right) \frac{\delta T_s}{T_{0s}}, \end{equation}

where ${\boldsymbol {\nabla }} x$ is the (radial) direction of the equilibrium gradients, $n_{0s}$ and $T_{0s}$ are the equilibrium density and temperature, respectively, of species $s$, $\delta T_s$ is the corresponding temperature perturbation [see (A17)], and

(2.2)\begin{equation} {\boldsymbol{v}}_E = \frac{c}{B_0} {\boldsymbol{b}}_0 \times {\boldsymbol{\nabla}} \phi, \end{equation}

is the ${\boldsymbol {E}}\times {\boldsymbol {B}}$ drift velocity due to the perturbed electrostatic potential $\phi$, $B_0$ and ${\boldsymbol {b}}_0$ being the magnitude and direction of the equilibrium magnetic field, respectively. In what follows, $L_{n_s}$ and $L_{T_s}$ denote the characteristic scale lengths associated with the gradients of the equilibrium density and temperature, respectively, while the equilibrium energy scale of species $s$ is set by its thermal speed $v_{{\rm th}s} = \sqrt {2T_{0s}/m_s}$, with $m_s$ being the particle mass.

In appendix A, we show that, for constant perpendicular equilibrium gradients, the electrostatic, drift-kinetic system of equations is invariant under a particular one-parameter transformation. Under this transformation, the perturbed temperature and electrostatic potential transform as, for any $\lambda$,

(2.3)\begin{equation} \delta \tilde{T}_s = \lambda^2 \delta T_s(x/\lambda^2, y/\lambda^2, z/\lambda^{2/\alpha}, t/\lambda^2), \quad \tilde{\phi} = \lambda^2 \phi(x/\lambda^2, y/\lambda^2, z/\lambda^{2/\alpha}, t/\lambda^2). \end{equation}

Here, $x$, $y$ and $z$ are the radial, binormal and parallel coordinates, respectively, the tildes indicate the transformed fields, and $\alpha = 1,2$ in the collisionless and collisional limits, respectively. We have assumed that the collisional limit corresponds to the case where the frequency of the perturbations $\omega$ is comparable to rate of thermal conduction, but much smaller than $\nu _{s s'}$, the characteristic collision frequency between species $s$ and $s'$, viz., $\omega \sim (k_\parallel v_{{\rm th}s})^2/\nu _{s s'} \ll \nu _{ss'}$, as in Braginskii (Reference Braginskii1965) (where $k_\parallel$ is the characteristic wavenumber of the perturbations along the direction of the equilibrium magnetic field). Mathematically, the existence of the symmetry (2.3) is a consequence of the scale invariance of electrostatic drift kinetics: in the absence of finite-Larmor-radius effects associated with $\rho _s$ – the Larmor radius of species $s$, manifest in the gyroaverages and the resultant Bessel functions appearing in gyrokinetics (see, e.g. Abel et al. Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013) – there is no intrinsic perpendicular scale in the system, with nothing to distinguish any perpendicular scale from any other.

Under (2.3), and noting the presence of the perpendicular derivative in (2.2), the heat flux (2.1) transforms as

(2.4)\begin{equation} \tilde{Q}_s = \lambda^2 Q_s. \end{equation}

Now suppose that our original solutions for $\delta T_s$ and $\phi$ were periodic in $x$, $y$ and $z$ with domain sizes $L_x$, $L_y$, and $L_\parallel$, respectively. Then, the transformed solutions $\delta \tilde {T}_s$ and $\tilde {\phi }$ are still periodic in $x$, $y$ and $z$, except with domain sizes $\lambda ^2 L_x$, $\lambda ^2 L_y$, and $\lambda ^{2/\alpha } L_\parallel$, implying that

(2.5)\begin{equation} \tilde{Q}_s(\lambda^2 L_x, \lambda^2 L_y, \lambda^{2/\alpha} L_\parallel, t/\lambda^2) = \lambda^{2} Q_s( L_x, L_y, L_\parallel, t). \end{equation}

The heat flux will, of course, depend on other parameters of the system, e.g. equilibrium gradients and collisionality. These, however, remain unchanged under the transformation (by construction), and so we did not write them explicitly in (2.5). In a strongly magnetised (gyrokinetic) plasma, structures generated by the turbulent fluctuations are ordered comparable to the equilibrium scales in the parallel direction ($k_\parallel ^{-1} \sim L_\parallel \sim L_{T_s}$), but remain microscopic in the perpendicular direction ($k_\perp ^{-1} \sim \rho _s$). This means that, as the perpendicular domain size $L_\perp$ (ordered as $L_\perp \sim L_x \sim L_y \sim \rho _s$) is increased, there must come a point at which the turbulence, and the resultant heat flux, become independent of the perpendicular domain size; if this were not the case, then the heat flux would diverge as $L_\perp /\rho _s \rightarrow \infty$, implying that drift kinetics is not a valid local model of the plasma. We thus assume that the heat flux is independent of the perpendicular domain size, viz., independent of $L_x$ and $L_y$. We also assume that the heat flux is independent of time, in the sense that it has been able to reach a statistical steady state. Then, given that $\lambda$ can be chosen arbitrarily, (2.5) directly implies that

(2.6)\begin{equation} Q_s \propto L_\parallel^{\alpha}, \end{equation}

where once again $\alpha = 1,2$ in the collisionless and collisional limits, respectively. Physically, $L_\parallel$ can be thought of either as a measure of a quantity analogous to the connection length $2{\rm \pi} q R$ in tokamak geometry (where $q$ is the safety factor and $R$ the major radius) or, in the absence of any other gradients, as a proxy for the temperature-gradient scale length. The latter follows from dimensional analysis: without loss of generality,

(2.7)\begin{equation} \frac{Q_s}{Q_{\text{gB}s}} = \left( \frac{L_\parallel}{L_{T_s}} \right)^\alpha G \left( \nu_{*s}, \frac{L_{n_s}}{L_{T_s}}, \frac{R}{L_{T_s}}, \ldots \right), \end{equation}

where $Q_{\text {gB}s} = n_{0s} T_{0s} v_{{\rm th}s} (\rho _s/L_{T_s})^2$ is the ‘gyro-Bohm’ heat flux, $G$ is an unknown function, $\nu _{*s} = (L_T/v_{{\rm th}s}) \sum _{s '} \nu _{s s'}$ is the normalised collisionality, and ‘$\dots$’ stands for other equilibrium parameters on which the heat flux can depend, normalised, wherever a scale is required, using the temperature-gradient scale length $L_{T_s}$. If the dependence of $G$ on these other parameters can be ignored in (2.7) – due either to the absence of other gradients in, e.g. slab geometry, or the system being driven far above marginality where such dependences are typically weak – then the scaling of the heat flux with $L_{T_s}$ follows directly from its dependence on $L_\parallel$.

This is perhaps a surprising result. Under the assumptions that the system is spatially periodic, that it is able to reach a statistical steady state (stationarity), and that the heat flux is independent of the system's perpendicular size, as it should be for any valid local model of a plasma (spatial locality), the scale invariance of electrostatic drift kinetics enforces the scaling (2.6), which is a non-trivial prediction about the scaling of the heat flux with equilibrium parameters. The key physics question, then, is how the system organises itself in order to obey this scaling. Namely, it has to find a way to process the free energy injected by the equilibrium gradients at a steady rate (stationarity) and to choose a spatial scale independent of $L_\perp$ (locality). In the remainder of this paper, we investigate a particular example of a system that should exhibit this scaling, being derived in an asymptotic limit of electrostatic drift kinetics, and find that a critically balanced, Kolmogorov (Reference Kolmogorov1941)-style cascade of energy from large to small spatial scales is the dynamical means by which the formal constraint imposed by (2.3) is realised.

3 Collisional fluid model

Fluid models are capable of providing remarkable insight about the dynamics of more general physical systems, while retaining the advantage of being (comparatively) simple to handle both numerically and analytically (e.g. Cowley et al. Reference Cowley, Kulsrud and Sudan1991; Newton, Cowley & Loureiro Reference Newton, Cowley and Loureiro2010; Ivanov et al. Reference Ivanov, Schekochihin, Dorland, Field and Parra2020; Ivanov, Schekochihin & Dorland Reference Ivanov, Schekochihin and Dorland2022). The ITG and ETG instabilities in tokamaks rely on destabilisation mechanisms that are fundamentally fluid (i.e. they are not resonant instabilities) and, even in kinetic regimes, they tend to be described adequately by fluid closures (Hammett & Perkins Reference Hammett and Perkins1990; Hammett, Dorland & Perkins Reference Hammett, Dorland and Perkins1992; Dorland & Hammett Reference Dorland and Hammett1993; Hammett et al. Reference Hammett, Beer, Dorland, Cowley and Smith1993; Beer & Hammett Reference Beer and Hammett1996; Snyder, Hammett & Dorland Reference Snyder, Hammett and Dorland1997). Here we consider an electron-scale, collisional (${\nu _{*e} \gg 1}$) fluid model of electrostatic turbulence driven by the electron-temperature gradient. Despite its simplicity, we expect that many of the results reported below are qualitatively applicable to more general turbulent plasma systems.

We take the local plasma equilibrium to be that of conventional slab gyrokinetics (see, e.g. Howes et al. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2006). The (homogeneous) equilibrium magnetic field is in the ${\boldsymbol {b}}_0 = \hat {{\boldsymbol {z}}}$ direction; perturbations to both its direction and magnitude are assumed negligible, consistent with the electrostatic limit. The electric field is then related to the electrostatic potential $\phi$ by ${\boldsymbol {E}} = - {\boldsymbol {\nabla }} \phi$, and has no mean part. The equilibrium profile of the electron temperature $T_{0e}$ varies radially, with the scale length

(3.1)\begin{equation} L_T^{{-}1} \equiv L_{T_e}^{{-}1} ={-} \frac{1}{T_{0e}} \frac{\mathrm{d} T_{0e}}{\mathrm{d} x}, \end{equation}

which is assumed to be constant over the domain of the system. The equilibrium gradients of density and ion temperature are assumed to be negligibly small. The omission of an equilibrium density gradient means that our system will be unstable for any finite $L_T^{-1}$; the resultant dynamics can thus be considered to apply to a plasma driven strongly above marginality (unlike the cases considered in, e.g. Guttenfelder et al. Reference Guttenfelder, Groebner, Canik, Grierson, Belli and Candy2021, Hatch et al. Reference Hatch, Michoski, Kuang, Chapman-Oplopoiou, Curie, Halfmoon, Hassan, Kotschenreuther, Mahajan and Merlo2022, Chapman-Oplopoiou et al. Reference Chapman-Oplopoiou, Hatch, Field, Frassinetti, Hillesheim, Horvath, Maggi, Parisi, Roach and Saarelma2022 and Field et al. Reference Field, Chapman-Oplopoiou, Connor, Frassinetti, Hatch, Roach and Saarelma2023, where a strong dependence of the heat flux on the equilibrium density gradient was identified).

3.1 Moment equations

With this local equilibrium, we derive, in appendix B, evolution equations for the density ($\delta n_{e}$), parallel velocity ($u_{\parallel e}$), and temperature ($\delta T_{e}$) perturbations of the electrons:

(3.2)\begin{align} &\frac{\mathrm{d} }{\mathrm{d} t} \frac{\delta n_{e}}{n_{0e}} + \frac{\partial u_{{\parallel} e}}{\partial z} = 0, \end{align}
(3.3)\begin{align}&\frac{\nu_{ei}}{c_1} \frac{u_{{\parallel} e}}{v_{{\rm th}e}} ={-} \frac{v_{{\rm th}e}}{2} \frac{\partial}{\partial z} \left[\frac{\delta n_{e}}{n_{0e}} - \varphi + \left( 1 + \frac{c_2}{c_1} \right) \frac{\delta T_{e}}{T_{0e}} \right], \end{align}
(3.4)\begin{align}&\frac{\mathrm{d}}{\mathrm{d} t} \frac{\delta T_{e}}{T_{0e}} + \frac{2}{3} \frac{\partial}{\partial z} \frac{\delta q_e}{n_{0e} T_{0e}} + \frac{2}{3} \left(1 + \frac{c_2}{c_1} \right) \frac{\partial u_{{\parallel} e}}{\partial z} ={-} \frac{\rho_e v_{{\rm th}e}}{2 L_T} \frac{\partial \varphi}{\partial y}. \end{align}

Let us discuss what these equations represent.

Equation (3.2) is the familiar continuity equation. It describes the advection of the density perturbation by the ${\boldsymbol {E}}\times {\boldsymbol {B}}$ motion (2.2) of the electrons,

(3.5)\begin{equation} \frac{\mathrm{d}}{\mathrm{d} t} = \frac{\partial}{\partial t} + {\boldsymbol{v}}_E{\boldsymbol{\cdot}} {\boldsymbol{\nabla}}_\perp{=} \frac{\partial}{\partial t} + \frac{\rho_e v_{{\rm th}e}}{2} \left\{ \varphi, \ldots\right\},\quad \varphi = \frac{e \phi}{T_{0e}}, \end{equation}

and their compression or rarefaction due to the perturbed electron flow $u_{\parallel e} {\boldsymbol {b}}_0$ parallel to the equilibrium magnetic-field direction.

This flow velocity is determined (instantaneously) from a balance between the electron–ion frictional force – proportional to the electron–ion collision frequency $\nu _{ei}$ [see (B4)] and appearing on the left-hand side of (3.3) – and the forces on the right-hand side of (3.3): the parallel pressure gradient, the electrostatic part of the parallel electric field, and the collisional ‘thermal forces’ (Braginskii Reference Braginskii1965; Helander & Sigmar Reference Helander and Sigmar2005), proportional to $c_2/c_1$, that arise due to the velocity dependence of the collision frequency associated with the Landau collision operator [see (A4)]. The order-unity constants $c_1$, $c_2$ and $c_3$ [the latter appearing in (3.6)] arise from the inversion of said operator, and depend on the magnitude of the ion charge $Z$ [see (B38)]: – e.g. for $Z=1$, $c_1 \approx 1.94$, $c_2 \approx 1.39$ and $c_3 \approx 3.16$, in agreement with Braginskii (Reference Braginskii1965).

The temperature perturbation in (3.4) is advected by the local ${\boldsymbol {E}}\times {\boldsymbol {B}}$ flow (2.2), again according to (3.5), and is locally increased (or decreased) by compressional heating (or rarefaction cooling) due to $u_{\parallel e}$, as well as by the perturbed parallel collisional heat flux,

(3.6)\begin{equation} \frac{\delta q_e}{n_{0e} T_{0e}} ={-} c_3\frac{ v_{{\rm th}e}^2}{2 \nu_{ei}} \frac{\partial}{\partial z} \frac{\delta T_{e}}{T_{0e}}, \end{equation}

caused by the gradient of the temperature perturbation along the equilibrium magnetic field direction. The term on the right-hand side of (3.4) is the familiar linear drive (advection of the equilibrium temperature profile by the perturbed ${\boldsymbol {E}} \times {\boldsymbol {B}}$ flow) responsible for extracting free energy from the equilibrium temperature gradient $L_T^{-1}$, defined in (3.1).

Finally, the electron-density perturbation is related to the non-dimensionalised potential $\varphi$ via quasineutrality:

(3.7)\begin{equation} \frac{\delta n_e}{n_{0e}} ={-} \bar{\tau}^{{-}1} \varphi, \quad \bar{\tau} = \frac{\tau}{Z}, \end{equation}

where $\tau = T_{0i}/T_{0e}$ is the ratio of the ion to electron equilibrium temperatures. This describes an adiabatic ion response at electron scales: at scales much smaller than their Larmor radius $\rho _i$, ions can be viewed as motionless rings of charge, and their density response is Boltzmann.

Given (3.3), (3.6) and (3.7), we can contract our system to two evolution equations written entirely in terms of the electrostatic potential and temperature perturbations:

(3.8)\begin{align} &\frac{\partial}{\partial t} \bar{\tau}^{{-}1} \varphi - \frac{c_1 v_{{\rm th}e}^2}{2 \nu_{ei}} \frac{\partial^2 }{\partial z^2}\left[\left(1 + \frac{1}{\bar{\tau}} \right)\varphi - \left( 1 + \frac{c_2}{c_1} \right) \frac{\delta T_{e}}{T_{0e}} \right] = 0, \end{align}
(3.9)\begin{align} & \frac{\mathrm{d}}{\mathrm{d} t} \frac{\delta T_{e}}{T_{0e}} + \frac{2}{3 } \frac{c_1 v_{{\rm th}e}^2}{2 \nu_{ei}} \frac{\partial^2}{\partial z^2} \left\{\left( 1 + \frac{1}{\bar{\tau}} \right) \left( 1 + \frac{c_2}{c_1} \right) \varphi - \left[\frac{c_3}{c_1} + \left( 1 + \frac{c_2}{c_1} \right)^2 \right] \frac{\delta T_{e}}{T_{0e}} \right\} \nonumber\\ &\quad ={-} \frac{\rho_e v_{{\rm th}e}}{2 L_T}\frac{\partial \varphi}{\partial y}. \end{align}

Note that, due to the Boltzmann density response (3.7), the advection term in (3.8) has vanished, leaving a purely linear relationship between $\varphi$ and $\delta T_{e}$. The only nonlinearity left in the system is thus the advection of $\delta T_{e}$ by the ${\boldsymbol {E}}\times {\boldsymbol {B}}$ flow in (3.9). Note that we find finite-amplitude nonlinear saturation in our numerical simulations despite this absence of any nonlinearity in (3.8); see § 5.5 for further discussion. If finite magnetic drifts (associated with an inhomogeneous equilibrium magnetic field) are included in our model, however, we find that simulations fail to saturate; this is discussed in § 6 and appendix C.

3.2 Scale invariance

Given that (3.8) and (3.9) were derived in an asymptotic subsidiary limit of drift kinetics (see appendix B.3), they are necessarily invariant under the transformation (2.3) with $\alpha = 2$, and must therefore exhibit the scaling (2.6) of the heat flux with parallel system size – this is confirmed numerically in § 4.2.

Physically, this scale invariance is a consequence of the fact that (3.8) and (3.9) are valid within the wavenumber range (see appendix B.1),

(3.10)\begin{equation} \sqrt{\beta_e} \ll k_\parallel L_T \ll 1,\quad \beta_e \frac{\lambda_{ei}}{L_T} \ll k_\perp \rho_e \ll \frac{\lambda_{ei}}{L_T}, \end{equation}

i.e. at perpendicular scales much smaller that those at which electromagnetic effects become important (the ‘flux-freezing scale’; see Adkins et al. Reference Adkins, Schekochihin, Ivanov and Roach2022), but much larger than those on which one encounters the effects of electron thermal diffusion due to the finite Larmor motion of the electrons (Hardman et al. Reference Hardman, Parra, Chong, Adkins, Anastopoulos-Tzanis, Barnes, Dickinson, Parisi and Wilson2022; Adkins Reference Adkins2023) – both of these bring in a special perpendicular scale that would break the drift-kinetic scale invarianceFootnote 1 . In other words, (3.8) and (3.9) describe physics on scales

(3.11)\begin{equation} k_\parallel L_T \sim \sqrt \sigma ,\quad k_\perp \rho_\perp \sim 1, \quad \rho_\perp = \frac{\rho_e}{\sigma} \frac{L_T}{\lambda_{ei}}, \end{equation}

where $\sigma$ is, formally, some arbitrary constant satisfying

(3.12)\begin{equation} \beta_e \ll \sigma \ll 1. \end{equation}

The fact that it should be arbitrary follows from the fact that there is no special scale within the wavenumber ranges (3.11). Our normalisation of perpendicular and parallel wavenumbers in (3.8)–(3.9) will thus also be arbitrary, up to the definition of $\sigma$.

3.3 Collisional slab ETG instability

Let us now summarise briefly the linear stability properties of the system of equations (3.8)–(3.9). Linearising and Fourier-transforming these equations, one obtains the dispersion relation:

(3.13)\begin{align} & \omega^2 + \left[ 1+ \bar{\tau} + \frac{2}{3} \left(1 + \frac{c_2}{c_1}\right)^2 + \frac{2}{3} \frac{c_3}{c_1} \right] {\rm i} {\omega_\parallel} \omega \nonumber\\ &\quad - \left[\frac{2}{3}(1+\bar{\tau}) \frac{c_3}{c_1} {\omega_\parallel} + \left(1 + \frac{c_2}{c_1} \right) {\rm i}\omega_{*e}\bar{\tau} \right] {\omega_\parallel} = 0, \end{align}

where we have introduced the characteristic parallel and perpendicular frequencies

(3.14)\begin{equation} {\omega_\parallel} = c_1 \frac{\left(k_\parallel v_{{\rm th}e} \right)^2}{2 \nu_{ei}},\quad \omega_{*e} = \frac{k_y \rho_e v_{{\rm th}e}}{2 L_T}. \end{equation}

These are, respectively, the rate of parallel thermal conduction and the drift frequency associated with the electron-temperature gradient. Note that the dispersion relation (3.13) is quadratic in the frequency $\omega$ because the parallel velocity $u_{\parallel e}$ is determined instantaneously in terms of the other fields, by (3.3), unlike in collisionless ETG theory (Adkins et al. Reference Adkins, Schekochihin, Ivanov and Roach2022).

If we consider the limit of long parallel wavelengths, viz.,

(3.15)\begin{equation} {\omega_\parallel} \ll \omega \ll \omega_{*e}, \end{equation}

then the balance of the first and last terms in (3.13) gives us

(3.16)\begin{equation} \omega^2 = \left( 1 + \frac{c_2 }{c_1}\right) {\rm i} {\omega_\parallel} \omega_{*e} \bar{\tau} \quad \Rightarrow \quad \omega ={\pm} \frac{1+ {\rm i} \mathrm{sgn}(k_y)}{\sqrt{2}} \left( 1 + \frac{c_2 }{c_1}\right)^{1/2} \left( {\omega_\parallel}|\omega_{*e}| \bar{\tau} \right)^{1/2}. \end{equation}

We recognise this as the collisional slab ETG (sETG) instability (Adkins et al. Reference Adkins, Schekochihin, Ivanov and Roach2022), a cousin of the collisionless sETG (Lee et al. Reference Lee, Dong, Guzdar and Liu1987).

The minimal set of equations that elucidate this process physically can be obtained from (3.2)–(3.4) under the ordering (3.15): using (3.7) to express the density perturbations in terms of the electrostatic potential, we have:

(3.17)\begin{equation} \frac{\partial}{\partial t} \bar{\tau}^{{-}1} \varphi = \frac{\partial u_{{\parallel} e}}{\partial z} , \quad \frac{\nu_{ei}}{c_1} u_{{\parallel} e} ={-} \left(1 + \frac{c_2}{c_1} \right) \frac{v_{{\rm th}e}^2}{2} \frac{\partial}{\partial z} \frac{\delta T_e}{T_{0e}}, \quad \frac{\mathrm{d}}{\mathrm{d} t} \frac{\delta T_e}{T_{0e} } ={-} \frac{\rho_e v_{{\rm th}e}}{2 L_T} \frac{\partial \varphi}{\partial y}. \end{equation}

In this limit, the instability works as follows. Suppose that a small perturbation of the electron temperature is created with $k_y \neq 0$ and $k_\parallel \neq 0$, bringing the plasma from regions with higher $T_{0e}$ to those with lower $T_{0e}$ ($\delta T_{e} >0$), and vice versa ($\delta T_{e} < 0$). This temperature perturbation produces alternating hot and cold regions along the equilibrium magnetic field. The resulting perturbed temperature (and, therefore, pressure) gradients drive electron flows – determined instantaneously by the balance between the pressure gradient and collisional drag – from the hot regions to the cold regions [the second equation in (3.17)], giving rise to increased electron density in the cold regions [the first equation in (3.17)]. By quasineutrality, the electron density perturbation gives rise to an exactly equal ion density perturbation, and that, via the Boltzmann response (3.7), creates an electric field that produces an ${\boldsymbol {E}} \times {\boldsymbol {B}}$ drift that in turn pushes hotter particles further into the colder region, and vice versa [the third equation in (3.17)], reinforcing the initial temperature perturbation and thus completing the feedback loop required for the instability.

At short enough parallel wavelengths, the collisional sETG instability is quenched by rapid thermal conduction that leads to the damping of the associated temperature perturbation. To see this, we relax the assumption (3.15) and consider the exact stability boundary of (3.13), determined by the requirement that, assuming $\omega$ to be purely real, the real and imaginary parts of (3.13) must vanish individually. The resultant equations can be straightforwardly combined to yield

(3.18)\begin{equation} \left( \frac{ {\omega_\parallel}}{\omega_{*e}} \right)^2 = \frac{\displaystyle \frac{3}{2} \left(1 + \frac{c_2}{c_1} \right)^2 \bar{\tau}^2}{\displaystyle (1+ \bar{\tau}) \frac{c_3 }{c_1}\left[ 1+ \bar{\tau} + \frac{2}{3} \left(1 + \frac{c_2}{c_1}\right)^2 + \frac{2}{3} \frac{c_3}{c_1} \right]^2}. \end{equation}

This is a curve $k_y \propto k_\parallel ^2$ in wavenumber space, plotted as the grey dashed line in figure 1(a). Above this line, corresponding to the limit ${\omega _\parallel } \gg \omega _{*e}$, all modes are purely damped due to rapid thermal conduction, as in figure 1(b).

Figure 1. The growth rate of the collisional sETG instability: these are solutions to (3.13) for $\tau = Z = 1$, normalised to $\omega _{*e}$. Panel (a) is a contour plot of the positive growth rates ($\mathrm {Im} \omega >0)$ in the $(k_y, k_\parallel )$ plane; panel (b) shows cuts of the growth rate at constant $k_y \rho _\perp$, plotted as a function of $k_\parallel L_T/\sqrt {\sigma }$. The normalisations $\rho _\perp$ and $\sigma$ are defined in (3.11). The stability boundary (3.18) is indicated by grey dashed line in (a).

At any given $k_y$, the maximum growth rate of the collisional sETG is, therefore, reached when

(3.19)\begin{equation} {\omega_\parallel} \sim \omega \sim \omega_{*e}, \end{equation}

which is a balance between dissipation (through conduction) and energy injection due to the background temperature gradient. Indeed, maximising the growth rate from (3.13) with respect to ${\omega _\parallel }$, one finds $\gamma _\text {max} = C(\tau,Z) \omega _{*e}$, where $C(\tau,Z)$ is a constant formally of order unity, e.g. $C(1,1) \approx 0.094$ (cf. the maximum values in figure 1b). The increase of the maximum growth rate of the sETG instability with the perpendicular wavenumber, $\omega _{*e} \propto k_y$, can only be checked by the effects of the electron perpendicular thermal diffusion due to finite electron Larmor motion, which, as discussed in § 3.2, occurs outside the range of wavenumbers in which (3.8)–(3.9) are valid [see (3.11)], meaning that the instability grows fastest at the smallest perpendicular scales. The consequences of the intrinsic reliance of the collisional sETG on dissipative physics will be discussed in a nonlinear setting in § 5.1.

4 Numerical verification of scale invariance

4.1 Numerical setup

In what follows, the system (3.8)–(3.9) is solved numerically in a triply periodic box of size $L_x \times L_y \times L_\parallel$ using a pseudo-spectral algorithm. Numerical integration is done in Fourier space ($N_x$, $N_y$ and $N_\parallel$ are the number of Fourier harmonics in the respective directions) with the nonlinear term calculated in real space using the 2/3 rule for de-aliasing (Orszag Reference Orszag1971). We integrate the linear terms implicitly in time using the Crank–Nicolson method, while the nonlinear term is integrated explicitly using the Adams–Bashforth three-step method. This integration scheme is similar to the one implemented in the popular gyrokinetic code $\texttt{GS2}$ (Kotschenreuther, Rewoldt & Tang Reference Kotschenreuther, Rewoldt and Tang1995b; Dorland et al. Reference Dorland, Jenko, Kotschenreuther and Rogers2000).

Perpendicular hyperviscosity is introduced in order to provide an ultraviolet (large-wavenumber) cutoff for the instabilities, achieved by the replacement of the time derivative on the left-hand sides of (3.8) and (3.9) with

(4.1)\begin{equation} \frac{\partial}{\partial t} + ({-}1)^{N_\nu} \nu_\perp \left(\rho_\perp {\boldsymbol{\nabla}}_\perp \right)^{2 {N_\nu}}, \end{equation}

where $\nu _\perp$ is the ‘hypercollision’ frequency and ${N_\nu } \geqslant 2$. With this change, our equations now depend only on the following dimensionless parameters: the perpendicular and parallel box sizes $L_x/\rho _\perp$, $L_y/\rho _\perp$ and $L_\parallel \sqrt {\sigma }/L_T$, the hyper-collision frequency $2(\rho _\perp /\rho _e)^{2{N_\nu }}\nu _\perp /\nu _{ei}$, and the power of the hyperviscous diffusion operator ${N_\nu }$. Convergence scans in $N_x$, $N_y$, and the perpendicular box size $L_x = L_y = L_\perp$ were carried out on a baseline simulation (see table 1) to ensure that the chosen resolution adequately captured the dynamics, and to verify that $L_\perp$ was large enough so that it did not significantly affect the simulation results, as was required for the arguments of § 2.

Table 1. The parameters used in the ‘baseline’ and ‘higher-resolution’ simulations. Both simulations had $\tau = Z = 1$.

We have also found that our results do not depend on the specific details of the hyperviscosity, viz., on the values of $\nu _\perp$ and ${N_\nu }$. It can be viewed as a numerical tool that allows us to capture the dynamics of the system within a finite simulation domain and resolution, and is not intended to model a specific physical process. Ultimately, (4.1) is a stand-in for the physical sinks of energy that exist at higher perpendicular wavenumbers. The fact that our results end up being independent of hyperviscosity is, however, significant. The addition of (4.1) breaks the scale invariance associated with the transformation (2.3), similarly to the way in which FLR effects would break the drift-kinetic scale invariance in the context of gyrokinetics, a point that we shall revisit in § 6. One could thus question the inevitability of obtaining the scaling of the heat flux (2.6) in our system of equations with the modification (4.1). Furthermore, the fact that the growth rate of the sETG instability peaks at a perpendicular scale determined by the hyperviscosity – since (3.8) and (3.9) contain no intrinsic perpendicular wavenumber cutoff – may also be a cause for concern, as the most unstable perpendicular scale is often thought to play a central role in determining turbulent transport. Both of these concerns can be dispelled by the realisation that the arguments of § 2 did not rely on the details of the state of the system at small perpendicular scales; indeed, the behaviour of the heat flux is determined by the parallel system size $L_\parallel$, which is manifestly an equilibrium-scale quantity. In § 5.2, we will show that this is a consequence of the fact that the outer scale is central in (dynamically) determining the transport and that this outer scale turns out to be independent of hyperviscosity.

4.2 Scan in $L_\parallel /L_T$

In order to test the dependence of the turbulent heat flux on $L_\parallel$ predicted by (2.6), we performed a series of simulations in which $L_\parallel \sqrt {\sigma }/L_T$ was varied between 15 and 55 at fixed parallel resolution (viz., fixed ratio of $L_\parallel \sqrt {\sigma }/L_T$ to $N_\parallel$), while keeping all other parameters the same as in the baseline simulation (see table 1). Each simulation was run to long enough times for it to reach saturation and stay in a statistically stationary state for a while, as can be seen from the time traces of the instantaneous heat fluxes plotted in figure 2. That such a stationary state exists confirms one of the assumptions necessary for (2.6)Footnote 2 , the other assumption, also confirmed numerically, being that this state is independent of $L_x$ and $L_y$.

Figure 2. Time traces of the instantaneous heat flux from simulations in which $L_\parallel \sqrt {\sigma }/L_T$ was varied from 15 to 55, normalised to $(\rho _\perp /\rho _e) Q_{\text {gB}e}$.

In figure 3, we plot the time average of the turbulent heat flux – as defined in (2.1) for $s=e$, and normalised to $(\rho _\perp /\rho _e) Q_{\text {gB}e}$, where $Q_{\text {gB}e} = n_{0e} T_{0e} v_{{\rm th}e} (\rho _e/L_T)^2$ is the (electron) ‘gyro-Bohm’ flux. It is clear that the simulation data agrees extremely well with the theoretical scaling (2.6). This agreement, however, should not be a cause for complacency: though these results suggest that (2.6) correctly predicts the transport, we would like to understand how the system manages this, i.e. how it contrives to satisfy the assumptions underpinning the prediction (2.6). To explain this, we shall consider the dynamics in the inertial range. This is the subject of the following section.

Figure 3. The scaling of the turbulent heat flux with $L_\parallel /L_T$, normalised to $(\rho _\perp /\rho _e) Q_{\text {gB}e}$ and plotted against logarithmic axes. The points are the simulation data, while the theoretical prediction [see (2.6)] is shown by the dashed black line. A logarithmic fit to the data gives the slope of $2.02$.

5 Inertial-range dynamics

To ensure that we had sufficient numerical resolution to resolve adequately the dynamics of the inertial range, we conducted a ‘higher-resolution’ simulation (see table 1), on which we shall now focus. Due to the computational demands introduced by the higher resolution, this simulation was run only up to 5000 $(\rho _e/\rho _\perp )^2\nu _{ei} t/2$; this was sufficient to ensure that the heat flux had converged to a well-defined average value (see figure 4).

Figure 4. Turbulent heat flux in the higher-resolution simulation (see table 1), normalised to $(\rho _\perp /\rho _e) Q_{\text {gB}e}$. The upper and lower show, respectively, the instantaneous and (rolling) time-averaged heat fluxes in solid black. The dashed horizontal line in the lower panel is the average value – as calculated over the entire time interval – while the transparent grey region around this value shows the error bar associated with the mean, calculated by means of a moving window average. The time-averaged heat flux converges to within the final error bar by $(\rho _e/\rho _\perp )^2\nu _{ei} t/2 \sim 2000$.

5.1 Free-energy budget

Magnetised plasma systems containing small perturbations around a Maxwellian equilibrium nonlinearly conserve free energy, which is a quadratic norm of the magnetic perturbations and the perturbations of the distribution functions of both ions and electrons away from the Maxwellian (see, e.g. Abel et al. Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013). In the system of equations that we are considering, the (normalised) free energy reduces to the form

(5.1)\begin{equation} \frac{W}{n_{0e} T_{0e}} = \int \frac{\mathrm{d}^3 {\boldsymbol{r}}}{V} \left[\frac{1}{2\bar{\tau}}\left(1 + \frac{1}{\bar{\tau}} \right) \varphi^2 + \frac{3}{4} \frac{\delta T_{e}^2}{T_{0e}^2} \right]. \end{equation}

The free energy is a nonlinear invariant, i.e. it is conserved by nonlinear interactions, but can be injected into the system by equilibrium gradients and is dissipated by collisions. It is straightforward to show from (3.8) and (3.9) (with the hyperviscosity (4.1) appended) that the free-energy budget is

(5.2)\begin{equation} \frac{1}{n_{0e} T_{0e}} \frac{\mathrm{d} W}{\mathrm{d} t} = \varepsilon - D_\parallel{-} D_\perp, \end{equation}

where

(5.3)\begin{equation} \varepsilon = \frac{1}{L_T} \int \frac{\mathrm{d}^3 {\boldsymbol{r}}}{V} \frac{3}{2} \frac{\delta T_{e}}{T_{0e}} v_{Ex}, \quad v_{Ex} ={-}\frac{\rho_e v_{{\rm th}e}}{2} \frac{\partial \varphi }{\partial y}, \end{equation}

is the energy-injection rate from the equilibrium temperature gradient, and

(5.4)\begin{align} D_\parallel &= \frac{c_1v_{{\rm th}e}^2}{2 \nu_{ei}} \int \frac{\mathrm{d}^3 {\boldsymbol{r}}}{V} \left\{\left[\left(1 + \frac{1}{\bar{\tau}} \right) \frac{\partial \varphi}{\partial z} - \left(1 + \frac{c_2}{c_1} \right) \frac{\partial}{\partial z} \frac{\delta T_{e}}{T_{0e}} \right]^2 + \frac{c_3}{c_1} \left(\frac{\partial}{\partial z} \frac{\delta T_{e}}{T_{0e}} \right)^2\right\}, \end{align}
(5.5)\begin{align}D_\perp &= \nu_\perp \int \frac{\mathrm{d}^3 {\boldsymbol{r}}}{V} \left[\left( \rho_\perp^{{N_\nu}} {\boldsymbol{\nabla}}_\perp^{{N_\nu}} \varphi \right)^2 + \frac{3}{2} \left( \rho_\perp^{{N_\nu}} {\boldsymbol{\nabla}}_\perp^{{N_\nu}} \frac{\delta T_{e}}{T_{0e}}\right)^2 \right], \end{align}

are the dissipation rates due to (parallel) thermal conduction and (perpendicular) hyperviscosity, respectively. The corresponding one-dimensional (1D) perpendicular wavenumber spectrum of the energy injection is

(5.6)\begin{equation} \varepsilon_{{\boldsymbol{k}}} (k_\perp) = 2 {\rm \pi}k_\perp \int_{-\infty}^\infty \mathrm{d} k_\parallel \frac{3}{2} \mathrm{Re} \left\langle {\rm i} \omega_{*e} \varphi_{{\boldsymbol{k}}}^* \frac{{\delta T_{e}}_{{\boldsymbol{k}}}}{T_{0e}} \right\rangle, \end{equation}

while those of the parallel and perpendicular dissipation are

(5.7)\begin{align} {D_\parallel}_{{\boldsymbol{k}}}(k_\perp)& = 2 {\rm \pi}k_\perp \int_{-\infty}^\infty \mathrm{d} k_\parallel \left\langle \omega_\parallel\left[ c_{1}\left| \left(1+ \frac{1}{\bar{\tau}}\right) \varphi_{{\boldsymbol{k}}} - \left(1 + \frac{c_2}{c_1} \right) \frac{{\delta T_{e}}_{{\boldsymbol{k}}}}{T_{0e}}\right|^2 + c_3 \left|\frac{{\delta T_{e}}_{{\boldsymbol{k}}}}{T_{0e}}\right|^2 \right] \right\rangle, \end{align}
(5.8)\begin{align}{D_\perp}_{{\boldsymbol{k}}}(k_\perp) &= 2 {\rm \pi}k_\perp \int_{-\infty}^\infty \mathrm{d} k_\parallel \left\langle (k_\perp \rho_\perp)^{2{N_\nu}} \nu_\perp \left(\left|\varphi_{{\boldsymbol{k}}}\right|^2 + \frac{3}{2} \left|\frac{{\delta T_{e}}_{{\boldsymbol{k}}}}{T_{0e}}\right|^2 \right)\right\rangle . \end{align}

In (5.6), the asterisk denotes complex conjugation, and the angle brackets an ensemble average. Note that when analysing the output of simulations, we consider ensemble averages to be equal to time averages over a period following saturation and the establishment of a statistical steady state (e.g. after $(\rho _e/\rho _\perp )^2\nu _{ei} t/2 \sim 2000$ in figure 2).

Plotting the injection and dissipation spectra (5.6)–(5.8) in figure 5(a) allows us to make a series of important observations. The first, and unsurprising, one is that the perpendicular dissipation due to hyperviscosity is dominant only at the very smallest scales, where ${D_\perp }_{{\boldsymbol {k}}}$ peaks. This confirms the assertion made in § 4.1 that it can be viewed as a sink of energy that exists at higher perpendicular wavenumbers and has no significant effect on the dynamics. The outer scale – at which energy is primarily injected into the turbulence and which we define as corresponding to the perpendicular wavenumber where the maximum of (5.6) is achievedFootnote 3 – appears to be independent of hyperviscosity, being localised on much larger scales, where ${D_\perp }_{{\boldsymbol {k}}}$ is negligible. The arguments of § 3.2 leading to the heat-flux scaling (2.6) relied on the scale invariance of the drift-kinetic system, which, as we have discussed previously, is broken by the introduction of hyperviscosity. The fact that the energy injection is both independent of hyperviscosity and localised at the largest scales supports the prediction of (2.6) that the heat flux should be determined by the inviscid dynamics at scales where scale-invariant drift kinetics is valid.

Figure 5. (a) The 1D perpendicular spectra of the energy injection (5.6) (solid red), parallel dissipation (5.7) (dashed blue) and perpendicular dissipation (5.8) (dotted blue), normalised to $(\rho _e/L_T)^2 \nu _{ei}/2$. The location of the outer scale is shown by the black dot. The rate of parallel dissipation is significant at the largest scales, while perpendicular dissipation takes over at the smallest scales. (b) The cumulative perpendicular wavenumber integrals of the quantities plotted in (a), as well as the nonlinear energy flux (5.9) (solid black line). The latter is approximately constant in the inertial range, displaying only an order-unity variation, due to the finite simulation domain.

Considering scales that are larger than the injection scale, it is clear from figure 5 that the parallel dissipation ${D_\parallel }_{{\boldsymbol {k}}}$ is dominant there, peaking on scales comparable to the outer scale. This is because the existence of the collisional sETG instability depends intrinsically on the presence of thermal conduction (see § 3.3), which is a dissipative effect. Indeed, the maximum growth rate (3.19) occurs where the rates of thermal conduction and energy injection are comparable, ${\omega _\parallel } \sim \omega _{*e}$. Thus, in order to inject energy, the system has to dissipate a finite fraction of it. The energy that survives this dissipation then cascades to small scales through a constant-flux inertial range. This can be seen in figure 5(b), where we plot the cumulative perpendicular wavenumber integrals of (5.6), (5.7), and (5.8), as well as the nonlinear energy flux, which can be inferred from the difference between injection and dissipation:

(5.9)\begin{equation} \varGamma(k_\perp) = \int_0^{k_\perp} \mathrm{d} k_\perp' \left[\varepsilon_{{\boldsymbol{k}}'}(k_\perp') -{D_\parallel}_{{\boldsymbol{k}}'}(k_\perp') - {D_\perp}_{{\boldsymbol{k}}'}(k_\perp')\right]. \end{equation}

Both the injection and parallel-dissipation rates reach an approximate plateau at scales smaller than the outer scale and are much larger than the nonlinear energy flux, which is approximately constant in the inertial range, displaying an order-unity variation due to the finite width of the latter in our numerical simulations. The remainder of § 5 is devoted to characterising the dynamics in the inertial range in order to explain how the system organises itself to maintain a constant-flux cascade to small scales despite the presence of significant (parallel) dissipation.

5.2 Constant flux and critical balance in the inertial range

The results of the previous section suggest that our fully developed electrostatic turbulence organises itself into a state wherein there is a local cascade of the free energy (5.1) that carries the injected power from the outer scale, through an inertial range, to the (perpendicular) dissipation scale. This injected power is the (order-unity) fraction of $\varepsilon$ that survives the parallel dissipation at larger scales, viz., $\varepsilon - D_\parallel$, which, for brevity, we shall call $\varepsilon$ in the scaling arguments that follow.

The only nonlinearity in our equations is the advection of the temperature fluctuations by the fluctuating ${\boldsymbol {E}}\times {\boldsymbol {B}}$ flows in (3.9). Therefore, we take the nonlinear cascade time to be the nonlinear ${\boldsymbol {E}}\times {\boldsymbol {B}}$ advection time:

(5.10)\begin{equation} t_\text{nl}^{{-}1} \sim k_\perp v_E \sim \rho_e v_{{\rm th}e} k_\perp^2 \bar{\varphi} \sim \varOmega_e (k_\perp \rho_e)^2 \bar{\varphi}. \end{equation}

Here and in what follows, $\bar {\varphi }$ refers to the characteristic amplitude of the electrostatic potential at the scale $k_\perp ^{-1}$. Formally, $\bar {\varphi }$ can be defined by

(5.11)\begin{equation} \bar{\varphi}^2 = \int_{k_\perp}^\infty \mathrm{d} k_\perp' \, E_\perp^\varphi(k_\perp'),\quad E_\perp^\varphi(k_\perp) \equiv \int_{-\infty}^\infty \mathrm{d} k_\parallel \, 2{\rm \pi} k_\perp \left\langle|\varphi_{{\boldsymbol{k}}}|^2 \right\rangle, \end{equation}

where $E_\perp ^\varphi (k_\perp )$ is the 1D perpendicular spectrum of $\varphi$, $\varphi _{{\boldsymbol {k}}}$ is the spatial Fourier transform of the potential, and the angle brackets denote an ensemble average. The corresponding quantities for the temperature perturbations, $\delta \bar {T}_e$, $E_\perp ^T(k_\perp )$, and ${\delta T_{e}}_{{\boldsymbol {k}}}$, are defined analogously.

Assuming that any possible anisotropy in the perpendicular plane can be neglected (an assumption that will be verified in § 5.5), a Kolmogorov-style constant-flux argument leads to a scaling of the amplitudes in the inertial range:

(5.12)\begin{equation} t_\text{nl}^{{-}1} \frac{\delta \bar{T}_e^2}{T_{0e}^2 } \sim \varepsilon = \text{const} \quad \Rightarrow \quad \bar{\varphi} \frac{\delta \bar{T}_e^2}{T_{0e}^2} \sim \frac{\varepsilon}{\varOmega_e} (k_\perp \rho_e)^{{-}2}. \end{equation}

We are using the $\delta T_{e}/T_{0e}$ part of the free energy (5.1) because, as we noted earlier, in (3.8)–(3.9), $\delta T_{e}/T_{0e}$ is the only field that is advected nonlinearly whereas $\varphi$ is ‘sourced’ by the temperature perturbations through the second term on the left-hand side of (3.8).

To estimate the size of the electrostatic potential, we therefore balance the two terms in (3.8), yielding

(5.13)\begin{equation} \bar{\varphi} \sim \frac{{\omega_\parallel}}{\omega } \frac{\delta \bar{T}_{e}}{T_{0e}}, \end{equation}

which should hold at every scale. This implies that the potential and temperature perturbations will be comparable in magnitude and have the same wavenumber scaling throughout the inertial range if we posit, scale by scale, that

(5.14)\begin{equation} t_\text{nl}^{{-}1} \sim \omega \sim {\omega_\parallel}. \end{equation}

This is the conjecture of critical balance, whereby the characteristic time associated with parallel dynamics along the field lines is assumed comparable to the nonlinear advection rate $t_\text {nl}^{-1}$ at each perpendicular scale $k_\perp ^{-1}$, as in Barnes et al. (Reference Barnes, Parra and Schekochihin2011) and Adkins et al. (Reference Adkins, Schekochihin, Ivanov and Roach2022). The original rationale for this conjecture comes from the causality argument proposed in the context of MHD turbulence (Goldreich & Sridhar Reference Goldreich and Sridhar1995; Boldyrev Reference Boldyrev2005; Schekochihin Reference Schekochihin2022): two points along a field line can only remain correlated with one another if information can propagate between them faster than they are decorrelated by the (perpendicular) nonlinearity; in MHD, this information is carried by Alfvén waves (similarly, it can be carried by other waves in different plasma and hydro-dynamical systems; see Cho & Lazarian Reference Cho and Lazarian2004, Nazarenko & Schekochihin Reference Nazarenko and Schekochihin2011 and Adkins et al. Reference Adkins, Schekochihin, Ivanov and Roach2022). In our system, the parallel dynamics are dissipative, with the relevant timescale being set by the parallel conduction rate ${\omega _\parallel }$. Since there is no mechanism to preserve the parallel coherence of structures created by perpendicular mechanisms (via injection due to the sETG instability, or nonlinear cascade), one expects them to break up in the parallel direction to as fine scales as the system will allow, i.e. structures for which ${\omega _\parallel } \ll t_\text {nl}^{-1}$ should be immediately decorrelated by the nonlinearity and broken up into shorter pieces in the parallel direction. The limiting factor for this parallel refinement is that if structures reach parallel scales such that ${\omega _\parallel } \gg t_\text {nl}^{-1}$, they are wiped out by heat conduction. As a result, the ‘dissipation ridge’ (the line of critical balance) ${\omega _\parallel } \sim t_\text {nl}^{-1}$ will form a natural locus for turbulent structures. This is a version of critical balance that is appropriate for a system where parallel dissipation is present everywhere [which may also be true for collisionless plasmas, where ${\omega _\parallel }$ is instead the Landau (Reference Landau1946) damping rateFootnote 4]. In § 5.1, we saw that the actual amount of parallel dissipation that happens in the inertial range is small – free energy chooses to stay just shy of the dissipation region (${\omega _\parallel } > t_\text {nl}^{-1}$) and instead cascade, at an approximately constant rate, along the dissipation ridge (${\omega _\parallel } \sim t_\text {nl}^{-1}$).

Combining (5.12), (5.13) and (5.14), we find the following scaling of the amplitudes in the inertial range:

(5.15)\begin{equation} \bar{\varphi} \sim \frac{\delta \bar{T}_{e}}{T_{0e}} \sim \left( \frac{\varepsilon}{\varOmega_e} \right)^{1/3} (k_\perp \rho_e)^{{-}2/3}. \end{equation}

Then, recalling (5.11), the 1D perpendicular energy spectra in the inertial range are:

(5.16)\begin{equation} E_\perp^\varphi(k_\perp) \sim E_\perp^T(k_\perp) \sim \frac{\bar{\varphi}^2}{k_\perp} \propto k_\perp^{{-}7/3}. \end{equation}

Using (5.10) and (5.15), the critical balance (5.14) translates into the following relationship between parallel and perpendicular scales in the inertial range:

(5.17)\begin{equation} k_\parallel \lambda_{ei} \sim \frac{\varOmega_e^{1/3} \varepsilon^{1/6}}{\nu_{ei}^{1/2}} (k_\perp \rho_e)^{2/3}. \end{equation}

If we define the 1D parallel spectrum

(5.18)\begin{equation} E_\parallel^\varphi(k_\parallel) \equiv \int_0^\infty \mathrm{d} k_\perp \, 2{\rm \pi} k_\perp \left\langle|\varphi_{{\boldsymbol{k}}}|^2 \right\rangle, \end{equation}

and the corresponding temperature spectrum $E_\parallel ^T(k_\parallel )$ analogously, (5.17) and (5.15) imply the following inertial-range scaling of amplitudes with parallel wavenumbers:

(5.19)\begin{equation} \bar{\varphi} \propto k_\parallel^{{-}1} \quad \Rightarrow \quad E_\parallel^\varphi(k_\parallel) \sim E_\parallel^T(k_\parallel) \sim \frac{\bar{\varphi}^2}{k_\parallel} \propto k_\parallel^{{-}3}. \end{equation}

These simple scaling arguments are vindicated by simulation data. The 1D spectra (5.11) and (5.18) for both $\varphi$ and $\delta T_{e}$ are plotted in figure 6. They follow quite well the predicted scalings (5.16) and (5.19), respectively, below the outer scale and up to the wavenumbers at which the spectra begin to steepen due to perpendicular dissipation.

Figure 6. The 1D (a) perpendicular (5.11) and (b) parallel (5.18) spectra, normalised to their value at the outer scale. The spectra of the electrostatic potential are plotted in blue, those of the temperature perturbations are in red. The predicted inertial-range scalings (5.16) and (5.19) are shown by the dashed black lines. The location of the outer scale (see § 5.3) is indicated by the black dot. In (a), this is calculated from the maximum of (5.6), while in (b), it is calculated from the maximum of the 1D parallel spectrum of the energy injection, defined analogously to (5.6).

We shall return to these inertial-range scalings in § 5.4, where we will study the full 2D spectra of the turbulence and provide further support for the argument that the cascade follows the dissipation ridge (the line of critical balance), but first let us demonstrate how the simple scaling theory developed above allows one to recover – now on physically motivated dynamical grounds – the scaling of the heat flux (2.6) that was previously inferred from a formal scaling symmetry of our equations.

5.3 Outer scale and scaling of heat flux

From (3.19), we know that, for a given $k_y$, the most unstable collisional sETG modes satisfy

(5.20)\begin{equation} {\omega_\parallel} \sim \omega_{*e} \sim k_y \rho_e \frac{ v_{{\rm th}e}}{L_T}, \end{equation}

and thus grow at a rate $\sim \omega _{*e} \propto k_y$. This means that the linear instability will be overwhelmed by nonlinear interactions in the inertial range, because their characteristic rate increases more quickly with perpendicular wavenumber: from (5.10) and (5.15),

(5.21)\begin{equation} t_\text{nl}^{{-}1} \sim \varOmega_e \left( \frac{\varepsilon}{\varOmega_e} \right)^{1/3} (k_\perp \rho_e)^{4/3}. \end{equation}

The outer scale is then the scale at which these two rates are comparable: balancing (5.21) and (5.20), we get

(5.22)\begin{equation} \varOmega_e (k_\perp^o \rho_e )^2 \bar{\varphi}^o \sim \omega_\parallel^o \sim \omega_{*e} \quad \Rightarrow \quad \bar{\varphi}^o \sim (k_\perp^o L_T)^{{-}1},\quad k_y^o \rho_e \sim (k_\parallel^o)^2 L_T \lambda_{ei}, \end{equation}

where the superscript ‘$o$’ refers to outer-scale quantities. Thus, we have two relationships between $k_\perp ^o$, $\bar {\varphi }^o$ and $k_\parallel ^o$, but in order to determine the outer-scale quantities uniquely, we need a third constraint. Given that our system (3.8)–(3.9) is scale invariant, there is no special (microscopic) perpendicular scale that can be used to fix $k_y^o$. Then, assuming that the heat flux is independent of the perpendicular system size, the only remaining physically meaningful length scale that can set the outer scale is the parallel system size $L_\parallel$. The same should be true for more general systems described by electrostatic drift kinetics, as the scaling (2.6) would suggest and as we shall discuss shortly.

Assuming, then, that the outer scale is indeed set by the parallel system size, we find from (5.22) that

(5.23)\begin{equation} k_\parallel^o L_\parallel\sim 1 \quad \Rightarrow \quad \left( \frac{L_T}{\rho_\perp} \right) \bar{\varphi}^o \sim \left( k_\perp^o \rho_\perp \right)^{{-}1}, \quad k_\perp^o \rho_\perp \sim \left( \frac{L_T}{L_\parallel \sqrt{\sigma}} \right)^2,\end{equation}

where $\rho _\perp$ and $\sigma$ are defined in (3.11), and the magnitude of $\sigma$ only matters for the purposes of normalising amplitudes and wavenumbers in plots. Figure 7 shows that these theoretical predictions agree very well with the data from the scan in $L_\parallel /L_T$ that was presented in § 4.2.

Figure 7. (a) The scaling of the perpendicular outer scale $k_\perp ^o$ (defined as the peak wavenumber of the energy injection (5.6)) with $L_\parallel /L_T$. (b) The scaling of the amplitude of the electrostatic potential $\bar {\varphi }^o$ [defined as the amplitude of $\varphi$ at $k_\perp = k_\perp ^o$, via (5.11)] with the perpendicular outer scale. The black points are the simulation data, while the theoretical predictions (5.23) are shown by the black dashed lines. A logarithmic fit to the data gives the slopes of $-$1.99 and $-$0.95 in (a) and (b), respectively.

Let us now estimate the energy flux that is injected by the collisional sETG instability at the outer scale (5.23): using its definition (5.3), and ignoring any possibility of a non-order-unity contribution from phase factors, we have, from (5.13), (5.22) and (5.23),

(5.24)\begin{equation} \varepsilon \sim \omega_{*e}^o \bar{\varphi}^o \frac{\delta \bar{T}_{e}^o}{T_{0e}} \sim \frac{v_{{\rm th}e} \rho_e^2}{L_T^3 }(k_\perp^o \rho_e)^{{-}1} \sim \frac{v_{{\rm th}e} \rho_e^2 L_\parallel^2}{\lambda_{ei} L_T^4}. \end{equation}

Recalling (2.1), the combination of (5.24) and (5.22) yields the following expression for the turbulent heat flux:

(5.25)\begin{equation} Q_e \sim n_{0e}T_{0e} \varepsilon L_T \sim Q_{\text{gB}e} \frac{L_T}{\lambda_{ei}} \left( \frac{L_\parallel}{L_T} \right)^2 \propto \frac{L_\parallel^2}{L_T^3}, \end{equation}

where once again $Q_{\text {gB}e} = n_{0e} T_{0e} v_{{\rm th}e} (\rho _e/L_T)^2$ is the ‘gyro-Bohm’ flux. Unsurprisingly, this reproduces the scaling with $L_\parallel$ given by (2.6) for $\alpha = 2$ [cf., also, (2.7) for $G = L_T/\lambda _{ei}$]. Note that, apart from the inevitable dimensional factors, the $L_\parallel$ scaling determines (the nontrivial part of) the dependence of the turbulent heat flux on the temperature gradient, as we anticipated following (2.6). The fact that our equations (3.8) and (3.9) are invariant under the same transformation as drift kinetics (2.3) means that obtaining this scaling was, in a sense, a foregone conclusion. That being said, in arriving at (5.25) via this alternative route, we have been able to elucidate the dynamical origin of this scaling, viz., that it is consistent with a critically balanced, constant-flux nonlinear cascade of free energy to small perpendicular scales. This conclusion is not exclusive to the collisional model considered in this paper. Starting from (5.10), one can construct an entirely analogous theory for the turbulence driven by the collisionless sETG instability, obtaining a result equivalent to (5.25), which reproduces the scaling (2.6), this time for $\alpha = 1$ (see Adkins et al. Reference Adkins, Schekochihin, Ivanov and Roach2022).

Let us discuss the significance of our finding that the outer scale is fixed by the assumption that $k_\parallel ^o L_\parallel \sim 1$. Such a choice goes back to the work by Barnes et al. (Reference Barnes, Parra and Schekochihin2011), who conjectured, and numerically verified, that the outer scale of electrostatic, gyrokinetic ITG turbulence in tokamak geometry was set by the connection length $L_\parallel \sim qR$. While in their case, like ours, this was the only scale that could be reasonably viewed as the characteristic system size (the spatial inhomogeneity of the magnetic equilibrium), there was also another, seemingly more physically intuitive, justification available for its role in determining the large-scale cutoff for the ITG turbulence: one could assume that any turbulent structures correlated on parallel scales longer than the connection length would be damped in the stable (‘good-curvature’) region on the inboard side of the tokamak. Thus, one could believe that the operative reason for the significance of $L_\parallel \sim qR$ was the presence of large-scale dissipation, rather than, as we have now concluded, just the breaking of scale invariance – in our case, by the finiteness of a periodic box in the parallel directionFootnote 5 . A practical implication of this conclusion for more realistic systems appears to be that any long-scale parallel inhomogeneity should be sufficient to set $k_\parallel ^o$, without the need for it to be tied to an energy sink – this could matter for the analysis of turbulence in, e.g. edge plasmas (Parisi et al. Reference Parisi, Parra, Roach, Giroud, Dorland, Hatch, Barnes, Hillesheim, Aiba and Ball2020, Reference Parisi, Parra, Roach, Hardman, Schekochihin, Abel, Aiba, Ball, Barnes and Chapman-Oplopoiou2022) or in stellarators (Roberg-Clark, Plunk & Xanthopoulos Reference Roberg-Clark, Plunk and Xanthopoulos2022), where magnetic fields have parallel structure on scales shorter than the connection length.

5.4 Two-dimensional spectra

To provide a more detailed description of the critically balanced cascade (and to provide more evidence that it is indeed a critically balanced cascade), it is interesting to consider the 2D spectra:

(5.26)\begin{align} E^\varphi_{2\text{D}}(k_\perp,k_\parallel) &= 2{\rm \pi} k_\perp \left\langle|\varphi_{{\boldsymbol{k}}}|^2 \right\rangle , \end{align}
(5.27)\begin{align}E^T_{2\text{D}}(k_\perp,k_\parallel) &= 2{\rm \pi} k_\perp \left\langle|{\delta T_{e}}_{{\boldsymbol{k}}}/T_{0e}|^2 \right\rangle . \end{align}

Unlike in § 5.2, we can no longer assume that $E^\varphi _{2\text {D}} \sim E^T_{2\text {D}}$; this was true only for the ‘integrated’ 1D spectra dominated by the wavenumbers where the critical-balance conjecture (5.14) was assumed satisfied, and the two fields thus had the same scaling (5.13) for $\omega \sim \omega _\parallel$. With this in mind, we will first consider the spectrum of the temperature perturbations, from which the spectrum of the potential perturbations can then be inferred via (5.13).

We consider two wavenumber regions, above and below the ‘critical-balance line’ (5.17):

(5.28)\begin{equation} E_{2\text{D}}^T(k_\perp, k_\parallel) \sim\left\{\begin{array}{@{}ll} k_\parallel^{{-}a} k_\perp^b, & k_\parallel \gtrsim k_\perp^{2/3} , \\ k_\perp^{{-}c} k_\parallel^d, & k_\parallel \lesssim k_\perp^{2/3}, \end{array}\right. \end{equation}

where $a$, $b$, $c$, and $d$ are positive constants to be determined. Here, and in what follows, whenever our expressions appear to be dimensionally incorrect, this is because we have implicitly chosen to normalise our wavenumbers to the outer scale $k_\parallel /k_\parallel ^o \rightarrow k_\parallel$, $k_\perp /k_\perp ^o \rightarrow k_\perp$ so as to reduce notational clutter. To determine the scaling exponents, we follow the general scheme, which, for MHD turbulence, was laid out by Schekochihin (Reference Schekochihin2022) (see his appendix C).

Evidently, the scalings in the two regions in (5.28) must match along the boundary $k_\parallel \sim k_\perp ^{2/3}$, giving

(5.29)\begin{equation} a + d = \tfrac{3}{2}(b+c). \end{equation}

If $a>1$ and $d> -1$, $k_\parallel \sim k_\perp ^{2/3}$ will be the energy-containing parallel wavenumber at a given $k_\perp$. The 1D perpendicular spectrum is, therefore,

(5.30)\begin{equation} E_\perp^T(k_\perp) = \int \mathrm{d} k_\parallel \, E_{2\text{D}}^T(k_\perp, k_\parallel) \sim \int_0^{k_\perp^{2/3}} \mathrm{d} k_\parallel \, k_\perp^{{-}c}k_\parallel^d \sim k_\perp^{{-}c + 2(1+d)/3}. \end{equation}

This must match the scaling (5.16) of the 1D perpendicular spectrum derived from the constant-flux conjecture, implying that

(5.31)\begin{equation} c = \tfrac{2}{3}(1+d) + \tfrac{7}{3}. \end{equation}

Two further constraints follow from imposing boundary conditions as $k_\parallel$ or $k_\perp \rightarrow 0$ at constant $k_\perp$ or $k_\parallel$, respectively. The scaling of the spectrum as $k_\perp \rightarrow 0$ (in the region $k_\parallel \gg k_\perp ^{2/3}$) can be determined purely kinematically: the low-$k_\perp$ asymptotic behaviour of a homogenous 2D-isotropic field must be $k_\perp ^3$, implying that

(5.32)\begin{equation} b = 3. \end{equation}

This is a fairly standard resultFootnote 6 (see, e.g. appendix A of Schekochihin et al. Reference Schekochihin, Parker, Highcock, Dellar, Dorland and Hammett2016). Finally, the scaling as $k_\parallel \rightarrow 0$ (in the region $k_\parallel \ll k_\perp ^{2/3}$) follows from causality. Indeed, in § 5.2, we argued that fluctuations become decorrelated for $\omega _\parallel \lesssim t_\text {nl}^{-1}$ because they cannot communicate across parallel distances $\sim k_\parallel ^{-1}$, but such $k_\parallel$ are also too small for the fluctuations to be erased by thermal conduction. Therefore, the parallel spectrum at $k_\parallel \ll k_\perp ^{2/3}$ must be the spectrum of a 1D white noise:

(5.33)\begin{equation} d = 0. \end{equation}

Combining (5.32) and (5.33) with (5.29) and (5.31), we find

(5.34)\begin{equation} a = 9, \quad c = 3. \end{equation}

This gives us the following scalings for the 2D spectrum of the temperature perturbationsFootnote 7:

(5.35)\begin{equation} E_{2\text{D}}^T(k_\perp, k_\parallel) \sim \left\{\begin{array}{@{}ll} k_\parallel^{{-}9} k_\perp^3, & k_\parallel \gtrsim k_\perp^{2/3} , \\ k_\perp^{{-}3} k_\parallel^0, & k_\parallel \lesssim k_\perp^{2/3}. \end{array}\right. \end{equation}

Turning now to the 2D spectrum of the potential perturbations, analogously to (5.28), the conditions (5.29) and (5.31) are unmodified – the spectrum must still be continuous along $k_\parallel \sim k_\perp ^{2/3}$, and match the scaling of the 1D perpendicular spectrum that follows from the constant-flux conjecture, which is the same for both the potential and temperature perturbations. Similarly, the scaling of the spectrum as $k_\perp \rightarrow 0$ (in the region $k_\parallel \gtrsim k_\perp ^{2/3}$) will once again be $k_\perp ^3$ by the same kinematic argument, implying (5.32). From (5.29) and (5.31), we again have $a=9$. However, the causality argument that led to the white-noise scaling (5.33) at $k_\parallel \lesssim k_\perp ^{2/3}$ now no longer holds, because $\varphi$ is not directly decorrelated by the nonlinearity. Instead, it inherits its scaling from $\delta T_{e}/T_{0e}$ via the balance (5.13), viz.,

(5.36)\begin{equation} E_{2\text{D}}^\varphi \sim \frac{\omega_\parallel^2}{\omega^2} E_{2\text{D}}^T \sim \frac{k_\parallel^4 k_\perp^{{-}3}}{\omega^2}, \end{equation}

where we have used $\omega _\parallel \propto k_\parallel ^2$ and the second expression in (5.35). Now, in the region $k_\parallel \lesssim k_\perp ^{2/3}$, we expect thermal conductivity in the temperature equation to be subdominant to the nonlinear rate, and so estimating $\omega \sim t_\text {nl}^{-1} \propto k_\perp ^{4/3}$ in (5.36), we find that

(5.37)\begin{equation} d=4, \quad c = \tfrac{17}{3}. \end{equation}

This gives us the following scalings of the 2D spectrum of the potential fluctuations:

(5.38)\begin{equation} E_{2\text{D}}^\varphi(k_\perp, k_\parallel) \sim \left\{\begin{array}{@{}ll} k_\parallel^{{-}9} k_\perp^3, & k_\parallel \gtrsim k_\perp^{2/3} , \\ k_\perp^{{-}17/3} k_\parallel^4, & k_\parallel \lesssim k_\perp^{2/3}. \end{array}\right. \end{equation}

The full 2D spectrum of the temperature perturbations is plotted in figure 8. The organisation of the system about the critical-balance line is manifest here. Cuts of the 2D spectra (5.26) and (5.27) at constant $k_\parallel$ and $k_\perp$ are shown in figures 9 and 10, respectively, for both potential and temperature perturbations, showing good agreement with the theoretical scalings (5.35) and (5.38). The white-noise spectra at $k_\parallel \lesssim k_\perp ^{2/3}$, in particular, are another confirmation of the causal nature of the critical balance.

Figure 8. A contour plot of the logarithm of the 2D spectrum (5.27) of the temperature perturbations in the $(k_\perp,k_\parallel )$ plane, normalised to its value at the outer scale. The line of critical balance is shown as the dashed black line, while the outer scale is shown by the black dot. The horizontal dotted line shows the upper bound on the parallel-wavenumber cuts plotted in the right panels of figure 9. Similarly, the vertical dotted lines show the lower and upper bounds on the perpendicular-wavenumber cuts plotted in the right panels of figure 10.

Figure 9. Cuts of the 2D spectra of (a) the electrostatic potential and (b) the temperature perturbations at constant $k_\parallel$, normalised to $(\rho _\perp /L_T)^2$. The colours indicate the value of $k_\parallel L_T/\sqrt {\sigma }$ for a given cut. The left panels show the entire spectrum plotted as a function of $k_\perp \rho _\perp$. The right panels show selected cuts for $k_\parallel L_T$ within the inertial range, with $k_\perp$ rescaled according to the critical-balance relation (5.17). The black dashed lines show the theoretical scalings (5.38) and (5.35) in (a), and (b), respectively. The spectra show reasonable agreement with theory at both small and large perpendicular scales, despite the effects of hyperviscosity being present at the smallest scales.

Figure 10. Cuts of the 2D spectra of (a) the electrostatic potential and (b) the temperature perturbations at constant $k_\perp$, normalised to $(\rho _\perp /L_T)^2$. The colours indicate the value of $k_\perp \rho _\perp$ for a given cut. The left panels show the entire spectrum plotted as a function of $k_\parallel L_T$. The right panels show selected cuts of the spectrum for $k_\perp \rho _\perp$ within the inertial range, with $k_\parallel$ rescaled according to the critical-balance relation (5.17). The black dashed lines show the theoretical scalings (5.38) and (5.35) in (a) and (b), respectively. There is very good agreement with theory, especially at $k_\parallel \lesssim k_\perp ^{2/3}$, where the scalings extend well beyond the inertial range to higher $k_\perp \rho _\perp$, as can be seen from the left panels – this is because the causality argument is not sensitive to the precise details of the decorrelation physics.

The extraordinarily steep parallel-wavenumber scaling of the 2D spectra (5.35) and (5.38) in the region $k_\parallel \gtrsim k_\perp ^{2/3}$ can also be viewed as further evidence for the version of critical balance proposed following (5.14). In terms of timescales, this wavenumber constraint corresponds to $\omega _\parallel \gtrsim t_\text {nl}^{-1}$, and thus to a region of dominant thermal conduction that attempts to erase parallel structure created by the turbulence. The $k_\parallel$ scaling in this region proves to be so steep that the free-energy sink due to parallel dissipation is ineffective: the free energy cannot be nonlinearly transferred into this region in an efficient way, and instead cascades towards higher perpendicular wavenumbers along the critical-balance line, eventually encountering perpendicular dissipation, introduced, in our model, through hyperviscosity. Parallel dissipation thus acts not as a sink for the cascade, but instead creates the aforementioned ‘dissipation ridge’, constraining the cascade of energy in wavenumber space to be along the critical-balance line $k_\parallel \sim k_\perp ^{2/3}$.

One could dismiss this feature as being a peculiarity of our collisional model, given that the dissipative nature of collisional sETG instability (3.16) is hard-wired into it by construction. However, this picture might not be entirely dissimilar from what is observed in more generic systems of plasma turbulence: e.g. Hatch et al. (Reference Hatch, Terry, Jenko, Merz and Nevins2011) observed an overlap of the spatial scales of energy injection and dissipation in electrostatic, ion-scale toroidal gyrokinetic simulations, as did Told et al. (Reference Told, Jenko, TenBarge, Howes and Hammett2015) in the context of Alfvénic turbulence. The same behaviour could also be relevant in the context of kinetic ETG-driven turbulence. The growth rate of the collisionless sETG instability is limited by the parallel streaming rate ${\omega _\parallel } \sim k_\parallel v_{{\rm th}e}$ (see, e.g. Adkins et al. Reference Adkins, Schekochihin, Ivanov and Roach2022), which is also the rate of Landau damping; viewed in the context of the current discussion, this suggests, perhaps, that Landau damping could play a dissipative role similar to that of the thermal conduction in determining the way in which the system organises itself in order to support a constant-flux cascade of energy to small scales. Then, the rates of either parallel streaming or thermal conduction appearing in the critical balance ${\omega _\parallel } \sim \omega \sim t_\text {nl}^{-1}$ can also be interpreted as being there because they are the rates of parallel dissipation, rather than of the parallel propagation of information, limiting any further refinement of the parallel scale of the turbulent structures.

5.5 Perpendicular isotropy

Throughout §§ 5.2 to 5.4, our theoretical deductions were carried out under the assumption that $k_x \sim k_y \sim k_\perp$. This assumption of perpendicular isotropy is not obviously true and must be tested. Indeed, the maximum sETG growth rate (3.19) is at $k_x=0$, and so the outer-scale energy injection is predominantly into the so-called ‘streamers’: highly anisotropic ($k_x \ll k_y$) structures that can be identified in real space by their alternating pattern of horizontal bands stretched along the radial ($x$) direction (Cowley et al. Reference Cowley, Kulsrud and Sudan1991). In the context of ITG-driven turbulence, it is often assumed (and usually confirmed numerically) that these streamers are broken apart by zonal flows (see Barnes et al. (Reference Barnes, Parra and Schekochihin2011) and references therein), restoring isotropy at the outer scale; isotropy in the inertial range is then assumed as well. In ETG-driven turbulence, however, the role of zonal flows is less obvious (see, e.g. Dorland et al. Reference Dorland, Jenko, Kotschenreuther and Rogers2000; Jenko et al. Reference Jenko, Dorland, Kotschenreuther and Rogers2000; Colyer et al. Reference Colyer, Schekochihin, Parra, Roach, Barnes, Ghim and Dorland2017), and the existence of an isotropic state far from guaranteed – indeed, the real-space snapshots shown in figure 11 suggest that the system is in fact dominated by streamer-like structures on the largest scales, and there is little zonal-flow activity. Qualitatively, this is quite similar to what ETG turbulence has been reported to look like in gyrokinetic simulations (Joiner et al. Reference Joiner, Applegate, Cowley, Dorland and Roach2006; Candy et al. Reference Candy, Waltz, Fahey and Holland2007; Roach et al. Reference Roach, Abel, Akers, Arter, Barnes, Camenen, Casson, Colyer, Connor and Cowley2009; Guttenfelder & Candy Reference Guttenfelder and Candy2011).

Figure 11. Real-space snapshots of (a) the electrostatic potential and (b) the temperature perturbations from the higher-resolution simulation at $(\rho _e/\rho _\perp )^2\nu _{ei} t/2 = 3000$ (see table 1). The coordinate axes are as shown, while the red and blue colours correspond to regions of positive and negative fluctuation amplitudes. The turbulence does not appear to be isotropic on the large scales that are visible in these plots (streamers are manifest), but turns out to be isotropic in the inertial range (see figure 12).

To assess how isotropic the saturated state is, in particular in the inertial range, we plot the 2D spectra

(5.39)$$\begin{align} E^T(k_x, k_y) &= \int_{-\infty}^\infty \mathrm{d} k_\parallel \left\langle|{\delta T_{e}}_{{\boldsymbol{k}}}/T_{0e}|^2 \right\rangle, \end{align}$$
(5.40)$$\begin{align}E^T(k_\perp, \theta) &= k_\perp \int_{-\infty}^\infty \mathrm{d} k_\parallel \left\langle|{\delta T_{e}}_{{\boldsymbol{k}}}/T_{0e}|^2 \right\rangle, \end{align}$$

in figure 12. Here and in what follows, $\theta = \tan ^{-1}(k_y/k_x)$ is the polar angle in the perpendicular wavenumber plane. In both Cartesian and polar representations, we see that the spectrum is approximately isotropic with respect to the perpendicular wavevectors: at scales sufficiently smaller than the outer scale (viz., in the inertial range) contours of constant $E^T$ are either circles, in the case of (5.39), or horizontal lines, in the case of (5.40). The spectra of the potential perturbations, defined analogously to (5.39) and (5.40), display a similar isotropy. Thus, despite the fact that the largest scales are anisotropic due to the existence of streamers and the lack of vigorous zonal flows to break them apart, isotropy is restored in the inertial range.

Figure 12. Contour plots of the 2D spectra of the temperature perturbations, normalised to $(\rho _\perp /L_T)^2$: (a) in Cartesian coordinates, with the radial and poloidal wavenumbers plotted on the horizontal and vertical axes, respectively; contours of constant $E^T(k_x, k_y)$ (5.39) (black dashed lines) are approximately circular away from the origin, where injection is localised and the presence of streamers is manifested by the spectral power being shifted towards $k_y >k_x$; (b) in polar coordinates, with $\theta = \tan ^{-1}(k_y/k_x)$ and $k_\perp \rho _\perp$ plotted on the horizontal and vertical axes, respectively; contours of constant $E^T(k_\perp, \theta )$ (5.40) (black dashed lines) are approximately horizontal far away from $k_\perp \rho _\perp \lesssim 1$, where injection is localised.

The observed lack of zonal-flow activity is linked to the fact that there is nothing on electron scales to give zonal flows privileged status. This makes ETG turbulence quite unlike its ITG cousin on ion scales, where one encounters the modified adiabatic electron response (see, e.g. Hammett et al. Reference Hammett, Beer, Dorland, Cowley and Smith1993; Abel & Cowley Reference Abel and Cowley2013):

(5.41)\begin{equation} \frac{\delta n_i}{n_{0i}} = \frac{\delta n_{e}}{n_{0e}} = \varphi', \end{equation}

where $\varphi '$ is the non-zonal component of the (non-dimensionalised) electrostatic potential. Indeed, (5.41) has been found to be crucial for capturing essential zonal-flow physics (Rogers, Dorland & Kotschenreuther Reference Rogers, Dorland and Kotschenreuther2000; Ivanov et al. Reference Ivanov, Schekochihin, Dorland, Field and Parra2020, Reference Ivanov, Schekochihin and Dorland2022). Physically, this can be explained by the fact that (5.41) reserves a special status for zonal flows, in that it allows there to be non-trivial nonlinear interactions between the zonal flows and $\varphi '$ through the electrostatic ${\boldsymbol {E}} \times {\boldsymbol {B}}$ nonlinearity contained in the convective derivative (3.5) of the density perturbation. However, as we discussed in § 5.2, the adiabatic ion response (3.7) causes the nonlinearity in the electron continuity equation (3.8) to vanish identically. Crucially, this means that the system lacks any nonlinearity capable of generating 2D secondary instabilities that are responsible for the generation of zonal flows and destruction of streamer structures (Hasegawa & Mima Reference Hasegawa and Mima1978; Hasegawa & Wakatani Reference Hasegawa and Wakatani1983; Terry & Horton Reference Terry and Horton1983; Diamond et al. Reference Diamond, Itoh, Itoh and Hahm2005; Ivanov et al. Reference Ivanov, Schekochihin, Dorland, Field and Parra2020; Zhu, Zhou & Dodin Reference Zhu, Zhou and Dodin2020). A further consequence of this is that the model system (3.8)–(3.9) proves to be incapable of generating zonal flows even on longer timescales than those over which the observed isotropisation occurs, unlike what was observed in, e.g. Colyer et al. (Reference Colyer, Schekochihin, Parra, Roach, Barnes, Ghim and Dorland2017) and Tirkas et al. (Reference Tirkas, Chen, Merlo, Jenko and Parker2023).

6 Summary and discussion

We have considered the transport properties of electrostatic, drift-kinetic plasma turbulence, with a particular focus on the connection between its macroscopic transport properties and microscale (inertial-range) dynamics. In the presence of constant perpendicular equilibrium gradients, it has been observed that the equations of electrostatic drift kinetics possess a symmetry associated with their intrinsic scale invariance, in both the collisionless and collisional limits. Under the assumptions of spatial periodicity, stationarity, and locality, this symmetry has been shown to imply a particular scaling (2.6) of the heat flux $Q_s$ with the parallel system size $L_\parallel$, viz., $Q_s \propto L_\parallel$ (or $\propto L_\parallel ^2$, in the collisional limit), with the dependence on the equilibrium temperature gradient following from dimensional analysis (in the absence of other equilibrium gradients). This macroscopic transport prediction was then confirmed numerically in an electron-scale, collisional fluid model of electrostatic turbulence driven by the electron-temperature gradient. A critically balanced, constant-flux cascade of energy from some large, outer scale – at which energy is effectively injected by the (collisional) slab ETG instability – to small scales was then shown to be the microscale dynamics consistent with these macroscopic transport properties. This is one of only two extant numerical demonstrations of the existence of such a state in gradient-driven turbulence of this kind – the other being Barnes et al. (Reference Barnes, Parra and Schekochihin2011), for gyrokinetic ITG turbulence.

Two key observations can be made from our results: (i) the effects of dissipation associated with parallel thermal conduction play a key role in determining the saturated state of the turbulence, limiting the cascade of free energy in parallel wavenumbers by clamping it to the ‘dissipation ridge’, or line of ‘critical-balance’, in wavenumber space (Landau damping could play an analogous role in the collisionless limit); (ii) the outer scale of the turbulence is determined by the breaking of drift-kinetic scale invariance due to the existence of some long-scale parallel inhomogeneity – in our case, this was the finite size of our periodic box. In more realistic plasma systems like the tokamak, this could be the connection length $L_\parallel \sim qR$, or some shorter scale associated with inhomogeneities in the equilibrium magnetic field, such as in the tokamak edge or in stellarators.

These results demonstrate that the details of the (parallel) plasma equilibrium play a central role in determining the microscopic outer scale for the turbulence, and thus the saturated amplitudes to which turbulent fluctuations will grow, which in turn determine the observed macroscopic transport properties. The fact that turbulence in gradient-driven systems appears to behave similarly to those in which energy is injected by explicitly large-scale processes is encouraging from the perspective of theory, as it suggests that existing insights into, and experience of, the latter can be applied to the former, significantly less well-studied case.

Though the implications of drift-kinetic scale invariance were investigated here in a reduced model of ETG-driven turbulence, they nevertheless have implications for more realistic plasma systems due to strong constraints placed on the system by the resultant symmetry of the governing equations. Indeed, it must be stressed that while the physical regimes covered by the reduced model are limited in scope – lacking dynamics associated with, e.g. gradients of the plasma density or equilibrium magnetic field, kinetic effects, etc. – the scaling (2.6) of the heat flux suffers from no such limitations since it follows directly from the scale invariance of the electrostatic drift-kinetic system of equations. The existence of this scaling, however, is predicated on the adoption of the drift-kinetic limit. Restoring FLR effects by reverting to the (electrostatic) gyrokinetic equation will evidently break the scale invariance, as the scales $k_\perp \sim \rho _s^{-1}$ will now appear explicitly in the equations through the Bessel functions (see, e.g. Abel et al. Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013). Apart from some interesting exceptions (Parisi et al. Reference Parisi, Parra, Roach, Giroud, Dorland, Hatch, Barnes, Hillesheim, Aiba and Ball2020, Reference Parisi, Parra, Roach, Hardman, Schekochihin, Abel, Aiba, Ball, Barnes and Chapman-Oplopoiou2022), the general effect of these Bessel functions is to provide a cutoff for instabilities at large perpendicular wavenumbers, restricting the region of instability on the ultraviolet side and providing a sink of energy (dissipation) beyond the wavenumbers where the sETG growth rate peaks. The constant-flux arguments of § 5 assumed that there was sufficient separation between the outer scale and these dissipation regions in order to allow an inertial range to develop at the intermediate scales. Should such a separation exist, the system will effectively be drift-kinetic in the inertial range and, crucially, at the outer scale, where our results will continue to apply, despite the system being fully gyrokinetic. In other words, even if drift-kinetic scale invariance is broken at small scales, the assumption behind (2.6) is that the transport is set by the outer scale, which is in the drift-kinetic limit, and the relevant breaking of scale-invariance is done by $L_\parallel$. This is indeed what we observed in § 5: despite the breaking of scale invariance at small perpendicular scales due to the introduction of hyperviscosity, our simulation results still confirmed the scaling (2.6) as the well-defined outer scale was still set by $L_\parallel$. We acknowledge, however, that the scale separation required for such a state is far from guaranteed: non-zero magnetic shear, for example, can create long-wavelength modes with binormal wavenumbers $k_y \rho _i \sim 1$ but narrow radial structures near mode-rational surfaces on the scale $k_x \rho _e \sim 1$ (Hardman et al. Reference Hardman, Parra, Chong, Adkins, Anastopoulos-Tzanis, Barnes, Dickinson, Parisi and Wilson2022, Reference Hardman, Parra, Patel, Roach, Ruiz Ruiz, Barnes, Dickinson, Dorland, Parisi and St-Onge2023; Parisi et al. Reference Parisi, Parra, Roach, Hardman, Schekochihin, Abel, Aiba, Ball, Barnes and Chapman-Oplopoiou2022). Whether our results are robust to the effects of significant magnetic shear and other forms of equilibrium shaping that can amplify the importance of FLR – and thus undermine the possible separation between FLR effects and a putative outer scale – is a subject for future work.

A key assumption behind the scaling (2.6) is that the heat flux is able to reach a (statistical) steady state. The existence of such a state, however, is less assured than one might think: indeed, it has been known for some time that nonlinear saturation can fail to occur in simulations of electron-scale, gradient-driven turbulence due to the persistence of large-scale streamer structures in the absence of flow shear or a non-adiabatic ion response (Joiner et al. Reference Joiner, Applegate, Cowley, Dorland and Roach2006; Candy et al. Reference Candy, Waltz, Fahey and Holland2007; Roach et al. Reference Roach, Abel, Akers, Arter, Barnes, Camenen, Casson, Colyer, Connor and Cowley2009; Guttenfelder & Candy Reference Guttenfelder and Candy2011). In our simple fluid model, we too find that introducing magnetic drifts associated with an inhomogeneous equilibrium magnetic field is sufficient to reproduce this behaviour. In such simulations, the curvature-mediated ETG instability (Horton, Hong & Tang Reference Horton, Hong and Tang1988; Adkins et al. Reference Adkins, Schekochihin, Ivanov and Roach2022), absent in a straight magnetic field, gives rise to nonlinearly robust, large-scale streamer structures that cause unbounded growth of the heat flux with time. Further details of these simulations can be found in appendix C. This behaviour is consistent with the view that the adiabatic ion response (3.7) is insufficient to saturate ETG-scale turbulence in the presence of finite magnetic-field gradients (Hammett et al. Reference Hammett, Beer, Dorland, Cowley and Smith1993). It must be stressed that this lack of saturation does not break the drift-kinetic scale invariance (2.3), which is valid for any constant perpendicular equilibrium gradients, including those of the equilibrium magnetic field – it merely demonstrates that the steady state required to deduce (2.6) from (2.3) may not always be achievable in the regime of interest. Indeed, if we had been able to find a case of turbulence driven by the curvature-mediated ETG instability that saturated, we would have expected the scaling (2.6) for the corresponding heat flux (although not necessarily the same detailed inertial-range structure as for the sETG turbulence that we studied above), but in any event, no such saturated cases have so far presented themselves.

Another limiting assumption of our work was the electrostatic nature of the turbulence. The existence of finite electromagnetic perturbations also formally breaks the (electrostatic) drift-kinetic symmetry observed in § 2 (this is manifest on inspection of the equations of electromagnetic gyrokinetics). However, this does not necessarily imply that the heat-flux scaling (2.6) can never be realised in systems with finite beta. Indeed, we argued above that this scaling would still hold in the presence of FLR effects if the outer scale of the turbulence remained within the drift-kinetic limit, despite scale invariance being formally broken at the smallest spatial scales. A similar argument is applicable here. If the outer scale lies at scales sufficiently smaller than those on which electromagnetic effects are important (the ‘flux-freezing scale’, determined by the electron inertia in the collisionless limit, or resistivity in the collisional one; see Adkins et al. (Reference Adkins, Schekochihin, Ivanov and Roach2022)), then the scaling (2.6) will continue to hold as, once again, the assumption behind it is that the transport is set by the outer scale located in the electrostatic, drift-kinetic limit, and the relevant breaking of scale invariance is done by $L_\parallel$, rather than the flux-freezing scale. For example, Chapman-Oplopoiou et al. (Reference Chapman-Oplopoiou, Hatch, Field, Frassinetti, Hillesheim, Horvath, Maggi, Parisi, Roach and Saarelma2022) performed nonlinear, electromagnetic simulations of JET–ILW pedestals for $k_\perp \rho _i \gtrsim 1$, and observed the same scaling of the heat flux with $L_T$ as (2.7) at gradients sufficiently far above the linear threshold. However, if the outer scale lies on scales larger than the flux-freezing scale, i.e. if the turbulence is truly electromagnetic, then the constraints imposed by the scale invariance of electrostatic drift kinetics is lifted. Given that such regimes will likely be realised within tokamak-relevant reactor scenarios (see, e.g. Shimomura et al. Reference Shimomura, Murakami, Polevoi, Barabaschi, Mukhovatov and Shimada2001; Sips Reference Sips2005; Patel et al. Reference Patel, Dickinson, Roach and Wilson2021) due to higher experimental values of the plasma beta (the ratio of the thermal to magnetic pressures), a central focus of ongoing research is the extent to which any of the general physical conclusions of this paper carry over into truly electromagnetic systems of tokamak turbulence.

Acknowledgements

We are indebted to G. Acton, M. Barnes, W. Clarke, S. Cowley and W. Dorland for helpful discussions and suggestions at various stages of this project.

Editor N. Loureiro thanks the referees for their advice in evaluating this article.

Funding

This work was supported by the Engineering and Physical Sciences Research Council (EPSRC) [EP/R034737/1]. T.A. was previously supported by a UK EPSRC studentship. His work was carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programmes 2014–2018 and 2019–2020 under grant agreement no. 633053, and from the UKRI Energy Programme (EP/T012250/1). The views and opinions expressed herein do not necessarily reflect those of the European Commission. The work of A.A.S. was also supported in part by a grant from STFC (ST/W000903/1) and by the Simons Foundation via a Simons Investigator award.

Declaration of interests

The authors report no conflict of interest.

Appendix A. Derivation of scale invariance

In this appendix, we demonstrate explicitly that electrostatic drift kinetics remains invariant under the transformation (2.3), which leads to the heat-flux scaling (2.6).

We take as our starting point the electrostatic drift-kinetic system, in which the perturbed distribution function for species $s$ consists of a Boltzmann and gyrokinetic parts

(A1)\begin{equation} \delta \! f_s({\boldsymbol{r}}, {\boldsymbol{v}}, t) ={-} \frac{q_s \phi}{T_{0s}} f_{0s}({\boldsymbol{r}}, {\boldsymbol{v}}) + h_s ({\boldsymbol{r}},{\boldsymbol{v}}, t), \end{equation}

and $h_s$ evolves according to

(A2)\begin{equation} \frac{\partial }{\partial t} \left(h_s - \frac{q_s \phi }{T_{0s}} f_{0s} \right) + \left( v_\parallel {\boldsymbol{b}}_0 + {\boldsymbol{v}}_{ds} \right) {\boldsymbol{\cdot}} {\boldsymbol{\nabla}} h_s + \frac{c}{B_0} {\boldsymbol{b}}_0 {\boldsymbol{\cdot}} \left[ {\boldsymbol{\nabla}} \phi \times {\boldsymbol{\nabla}} \left( h_s + f_{0s} \right) \right] = \sum_{s'} C_{ss'}^{(l)}[h_s].\end{equation}

Here, and in what follows, $f_{0s}$ is the Maxwellian equilibrium distribution of species $s$ with density $n_{0s}$ and temperature $T_{0s}$, and $\phi$ is the perturbed electrostatic potential. The magnetic-drift velocity arising from the inhomogeneities in the equilibrium magnetic field is given by

(A3)\begin{equation} {\boldsymbol{v}}_{ds} = \frac{{\boldsymbol{b}}_0}{\varOmega_s} \times \left( v_\parallel^2 {\boldsymbol{b}}_0{\boldsymbol{\cdot}} \boldsymbol{\nabla}{\boldsymbol{b}}_0 + \frac{1}{2}v_\perp^2 \boldsymbol{\nabla}\log B_0 \right), \end{equation}

where ${\boldsymbol {b}}_0 = {\boldsymbol {B}}_0/B_0$ is the direction of the equilibrium magnetic field, $B_0 = |{\boldsymbol {B}}_0|$ is its magnitude, and $\varOmega _s = q_s B_0/m_s c$, $q_s$ and $m_s$ are the Larmor frequency, charge and mass of species $s$, respectively. The collision term on the right-hand side of (A2) is the linearised Landau collision operator

(A4) \begin{align} C_{s s'}^{(l)}\left[h_s \right] = \frac{\gamma_{s s'}}{m_s} {\boldsymbol{\nabla}}_v {\boldsymbol{\cdot}} \int \mathrm{d}^3 {\boldsymbol{v}}' \, & f_{0s}(v) f_{0s '}(v') ({\boldsymbol{\nabla}}_w {\boldsymbol{\nabla}}_w w ) \nonumber\\ & \boldsymbol{\cdot} \left[\frac{1}{m_s} {\boldsymbol{\nabla}}_v \frac{h_s({\boldsymbol{v}})}{f_{0s}(v)} -\frac{1}{m_{s '}}{\boldsymbol{\nabla}}_{v'} \frac{h_{s '}({\boldsymbol{v}}')}{f_{0s '}(v')} \right], \end{align}

where $w= |{\boldsymbol {w}}|$, ${\boldsymbol {w}} = {\boldsymbol {v}} - {\boldsymbol {v}}'$, $\gamma _{s s'} = 2{\rm \pi} q_s^2 q_{s '}^2 \log \varLambda$, and all velocity derivatives are evaluated at constant position ${\boldsymbol {r}}$. Finally, (A2) is closed by quasineutrality,

(A5)\begin{equation} 0 = \sum_{s} q_s \delta n_s = \sum_s q_s \left[ -\frac{q_s \phi}{T_{0s}} n_{0s} + \int \mathrm{d}^3 {\boldsymbol{v}} h_s \right]. \end{equation}

In what follows, it will be useful to decompose $h_s$ into parts that are even and odd in the parallel velocity $v_\parallel$, viz.,

(A6)\begin{gather} h_s^\text{even} ({\boldsymbol{r}}, v_\parallel, v_\perp, t) = \tfrac{1}{2 }\left[h_s({\boldsymbol{r}}, v_\parallel, v_\perp, t) + h_s({\boldsymbol{r}}, -v_\parallel, v_\perp, t)\right], \end{gather}
(A7)\begin{gather}h_s^\text{odd} ({\boldsymbol{r}}, v_\parallel, v_\perp, t) = \tfrac{1}{2 }\left[h_s({\boldsymbol{r}}, v_\parallel, v_\perp, t) - h_s({\boldsymbol{r}}, -v_\parallel, v_\perp, t)\right]. \end{gather}

It follows straightforwardly from (A2) and (A4) that $h_s^\text {even}$ and $h_s^\text {odd}$ satisfy, respectively,

(A8)\begin{align} & \frac{\partial }{\partial t} \left(h_s^\text{even} - \frac{q_s \phi }{T_{0s}} f_{0s} \right) +v_\parallel {\boldsymbol{b}}_0 {\boldsymbol{\cdot}} {\boldsymbol{\nabla}} h_s^\text{odd} + {\boldsymbol{v}}_{ds} {\boldsymbol{\cdot}} {\boldsymbol{\nabla}} h_s^\text{even} \nonumber\\ &\quad + \frac{c}{B_0} {\boldsymbol{b}}_0 {\boldsymbol{\cdot}} \left[ {\boldsymbol{\nabla}} \phi \times {\boldsymbol{\nabla}} h_s^\text{even}\right] + \frac{c}{B_0} {\boldsymbol{b}}_0 {\boldsymbol{\cdot}} \left[{\boldsymbol{\nabla}} \phi \times {\boldsymbol{\nabla}} f_{0s}\right] = \sum_{s'} C_{ss'}^{(\ell)} \left[h_s^\text{even} \right], \end{align}
(A9)\begin{align} &\frac{\partial h_s^\text{odd}}{\partial t} +v_\parallel {\boldsymbol{b}}_0 {\boldsymbol{\cdot}} {\boldsymbol{\nabla}} h_s^\text{even} + {\boldsymbol{v}}_{ds} {\boldsymbol{\cdot}} {\boldsymbol{\nabla}} h_s^\text{odd} + \frac{c}{B_0} {\boldsymbol{b}}_0 {\boldsymbol{\cdot}} \left[{\boldsymbol{\nabla}} \phi \times {\boldsymbol{\nabla}} h_s^\text{odd}\right] = \sum_{s'} C_{ss'}^{(\ell)} \left[h_s^\text{odd} \right]. \end{align}

The quasineutrality condition (A5) becomes

(A10)\begin{equation} 0 = \sum_s q_s \left[ -\frac{q_s \phi}{T_{0s}} n_{0s} + \int \mathrm{d}^3 {\boldsymbol{v}} \, h_s^\text{even} \right]. \end{equation}

Note that in (A8) and (A9), we have assumed that the (radial) gradient of the equilibrium distribution function ${\boldsymbol {\nabla }} f_{0s}$ is an even function of $v_\parallel$ – this is only the case in systems without any equilibrium flows.

We now wish to consider transformations of the system of equations (A8)–(A10) that can be made whilst preserving the size of perpendicular equilibrium gradients. It is obvious from considering, e.g. the magnetic-drift velocity (A3) that any rescaling of the velocity variables $v_\parallel$ and $v_\perp$ – at fixed equilibrium magnetic-field strength – would require a compensatory rescaling of ${\boldsymbol {\nabla }} \log B_0$ and $|{\boldsymbol {b}}_0 {\boldsymbol {\cdot }} {\boldsymbol {\nabla }} {\boldsymbol {b}}_0|$ in order to preserve the magnitude and direction of ${\boldsymbol {v}}_{ds}$. Therefore, we will henceforth restrict ourselves to transformations involving only the spatial and time coordinates. In a similar vein to Connor & Taylor (Reference Connor and Taylor1977), we consider the following one-parameter transformation:

(A11)\begin{align} \tilde{h}_s^\text{even} &= \lambda^{a_{e}} h_s^\text{even} (x/\lambda^{a_\perp}, y/\lambda^{a_\perp}, z/\lambda^{a_\parallel}, t/\lambda^{a_t}), \end{align}
(A12)\begin{align} \tilde{h}_s^\text{odd} &= \lambda^{a_{o}} h_s^\text{odd} (x/\lambda^{a_\perp}, y/\lambda^{a_\perp}, z/\lambda^{a_\parallel}, t/\lambda^{a_t}), \end{align}
(A13)\begin{align} \tilde{\phi} &= \lambda^{a_{e}} \phi (x/\lambda^{a_\perp}, y/\lambda^{a_\perp}, z/\lambda^{a_\parallel}, t/\lambda^{a_t}), \end{align}

where $x$, $y$ and $z$ are the radial, binormal and parallel (to the magnetic field) coordinates, respectively, the tildes indicate the transformed distribution functions and fields, and $a_i$ are real constants parametrising the transformation. Quasineutrality (A10) implies that the amplitudes of the ‘even’ fields must be rescaled in the same way, as in (A11) and (A13), while the rescaling of the amplitude of $h_s^\text {odd}$ remains unconstrained. The spatial and time coordinates can then be rescaled independently, with the caveat that the radial and binormal coordinates should be rescaled in the same way in order not to rule out perpendicular isotropy. The rescaling (A11)–(A13) is the most general one-parameter transformation of electrostatic drift kinetics that can be made while allowing (although not requiring) the spatial isotropy of structures in the perpendicular plane.

The constants $a_i$ can be fixed by demanding that the transformation leave (A8) and (A9) invariant. In the collisionless limit, the collision operator can be neglected and it is easy to show that $a_{e} = a_{o} = a_\perp = a_\parallel = a_t$ is the only choice that fulfils this condition. The collisional limit is somewhat more subtle. As we have done throughout this paper, we order $\omega \sim (k_\parallel v_{{\rm th}e})^2/\nu _{ss'} \ll \nu _{ss'}$ and $\omega h_s^{\text {even}} \sim k_\parallel v_{{\rm th}s} h_s^\text {odd}$. In the resultant collisional expansion, the collision operator is forced to vanish at leading order (see appendix B.3.1), and can only survive at higher order due to the presence of finite-Larmor-radius effects (see appendix B.3.3), which are neglected within the drift-kinetic approximation. At first order, one obtains, from (A9), a balance between the parallel streaming of $h_s^{\text {even}}$ and the collision operator acting on $h_s^\text {odd}$ (see appendix B.3.2). At second order, one evolves $h_s^\text {even}$ via (A8) with the collision operator neglected (see appendix B.3.3). One can then show that $a_{e} = 2 a_{o} = a_\perp = 2 a_\parallel = a_t$ is the only choice of parameters that leaves the drift-kinetic equations invariant. Any constraints on $a_i$ inferred from (A11)–(A10) in this way are only valid to second order within the collisional expansion, and not to any higher orders. However, given that the solvability conditions (B25) and (B36) guarantee that a closed system can be obtained solely from these two orders, this is not a problematic limitation.

Thus, it follows from the above discussion that electrostatic drift kinetics is invariant under the transformation

(A14)\begin{align} \tilde{h}_s^\text{even} &= \lambda^{2} h_s^\text{even} (x/\lambda^{2}, y/\lambda^{2}, z/\lambda^{2/\alpha}, t/\lambda^{2}), \end{align}
(A15)\begin{align}\tilde{h}_s^\text{odd} &= \lambda^{2/\alpha} h_s^\text{odd} (x/\lambda^{2}, y/\lambda^{2}, z/\lambda^{2/\alpha}, t/\lambda^{2}), \end{align}
(A16)\begin{align}\tilde{\phi} &= \lambda^{2} \phi (x/\lambda^{2}, y/\lambda^{2}, z/\lambda^{2/\alpha}, t/\lambda^{2}), \end{align}

where we have chosen $a_\mathrm {e} = 2$ without loss of generality, and $\alpha = 1,2$ in the collisionless and collisional limits, respectively. The transformation of the odd and even parts of the distribution function $h_s$ is inherited by its moments that are odd and even in $v_\parallel$; e.g. the temperature perturbation, being a velocity moment that is even in $v_\parallel$, viz.,

(A17)\begin{equation} \frac{\delta T_s}{T_{0s}} = \frac{2}{3n_{0s}} \int \mathrm{d}^3 {\boldsymbol{v}} \left(\frac{m_s v^2}{2T_{0s}} - \frac{3}{2} \right) h_s^{\text{even}}, \end{equation}

will transform according to (A14). This, when combined with the transformation (A16) of the electrostatic potential $\phi$ gives exactly (2.3), which is the starting point for the deductions presented in the main text.

Appendix B. Derivation of collisional fluid model

This appendix details a self-contained derivation of the electron-fluid equations (3.2)–(3.4). An alternative route to these via a subsidiary expansion of a more general system of (electromagnetic) equations can be found in Adkins (Reference Adkins2023) (see also appendix G of Adkins et al. Reference Adkins, Schekochihin, Ivanov and Roach2022). In what follows, appendix B.1 describes and physically motivates our electron-scale, collisional ordering, which is then implemented to derive equations describing our ion and electron dynamics in appendices B.2 and B.3, respectively. Although the magnetic geometry adopted throughout the majority of this paper is that of a conventional slab (see § 3), we shall here consider the more general case in which the equilibrium (mean) magnetic field ${\boldsymbol {B}}_0$ is assumed to have the scale length and radius of curvature

(B1)\begin{equation} L_B^{{-}1} ={-} \frac{1}{B_0} \frac{\mathrm{d} B_0}{\mathrm{d} x},\quad R^{{-}1} = \left| {\boldsymbol{b}}_0 {\boldsymbol{\cdot}} {\boldsymbol{\nabla}} {\boldsymbol{b}}_0 \right|, \end{equation}

assumed constant across the domain. Doing so will allow us to capture the effects of the magnetic drifts on our plasma while retaining most of the simplicity associated with conventional slab gyrokinetics (Howes et al. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2006; Newton et al. Reference Newton, Cowley and Loureiro2010; Ivanov et al. Reference Ivanov, Schekochihin, Dorland, Field and Parra2020, Reference Ivanov, Schekochihin and Dorland2022; Adkins et al. Reference Adkins, Schekochihin, Ivanov and Roach2022). Note that for a low-beta plasma, $R = L_B$.

B.1 Collisional, electron-scale ordering

In our model, we would like to be able to capture, at a minimum, the physics associated with drift waves, perpendicular advection by both magnetic drifts and ${\boldsymbol {E}}\times {\boldsymbol {B}}$ flows, and parallel heat conduction. Therefore, we postulate an asymptotic ordering in which the frequencies $\omega$ of the perturbations in the plasma are comparable to the characteristic frequencies associated with these phenomena, viz.,

(B2)\begin{equation} \nu_{ee} \sim \nu_{ei} \gg \omega \sim \omega_{*s} \sim \omega_{ds} k_\perp v_E \sim \kappa k_\parallel^2, \end{equation}

where

(B3)\begin{equation} \omega_{*s} = \frac{k_y \rho_s v_{{\rm th}s}}{2L_{T_s}}, \quad \omega_{ds} = \frac{k_y \rho_s v_{{\rm th}s}}{2L_B} \end{equation}

are the drift and magnetic-drift frequencies, respectively, ${\boldsymbol {v}}_E = c {\boldsymbol {E}} \times {\boldsymbol {B}}/B^2$ is the ${\boldsymbol {E}} \times {\boldsymbol {B}}$ drift velocity ($c$ is the speed of light), $\kappa \sim v_{{\rm th}e}^2/\nu _{ei}$ is the electron thermal diffusivity, and

(B4)\begin{equation} \nu_{ei} = \frac{4\sqrt{2{\rm \pi}}}{3} \frac{e^4 n_{0e} \log \varLambda}{m_e^{1/2} T_{0e}^{3/2}},\quad \nu_{ee}= \frac{\nu_{ei}}{Z} \end{equation}

are the electron–ion and electron–electron collision frequencies, respectively, $\log \varLambda$ being the Coulomb logarithm (Braginskii Reference Braginskii1965; Helander & Sigmar Reference Helander and Sigmar2005).

The ordering of the parallel conduction rate with respect to the drift frequency gives us a constraint relating parallel and perpendicular wavenumbers:

(B5)\begin{equation} \kappa k_\parallel^2 \sim \omega_{*s} \sim \omega_{ds} \sim k_\perp \rho_e \frac{v_{{\rm th}e}}{L} \quad \Rightarrow \quad (k_\parallel L)^2 \sim \frac{L}{\lambda_{ei}} k_\perp \rho_e, \end{equation}

where $\lambda _{ei} = v_{{\rm th}e}/\nu _{ei}$ is the electron–ion mean free path and $L$ is some (perpendicular) equilibrium length scale, $L\sim L_{T_s} \sim L_B \sim R$. The ordering of the parallel conduction rate with respect to the ${\boldsymbol {E}} \times {\boldsymbol {B}}$ drifts determines the size of perpendicular flows within our system:

(B6)\begin{equation} \kappa k_\parallel^2 \sim k_\perp v_E \quad \Rightarrow \quad \frac{v_E}{v_{{\rm th}e}} \sim \frac{k_\parallel}{k_\perp} k_\parallel \lambda_{ei} \sim \frac{\rho_e}{L} \equiv \epsilon, \end{equation}

where $\epsilon = \rho _e/L$ is the gyrokinetic small parameter (see, e.g. Abel et al. Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013), mandating small-amplitude, anisotropic perturbations. The frequency of these perturbations is small compared with the Larmor frequencies of both the electrons and ions:

(B7)\begin{equation} \frac{\omega}{\varOmega_e} \sim \frac{k_\perp v_E}{\varOmega_e} \sim k_\perp \rho_e \epsilon, \quad \frac{\omega}{\varOmega_i} = \frac{m_i}{Z m_e} \frac{\omega}{\varOmega_e}\sim k_\perp \rho_e \epsilon \frac{ m_i}{m_e}. \end{equation}

The ordering (B6) of $v_E$ relative to the electron thermal velocity allows us to order the amplitude of the perturbed scalar potential $\phi$:

(B8)\begin{equation} \frac{v_E}{v_{{\rm th}e}} \sim \frac{c}{B_0} \frac{k_\perp \phi}{v_{{\rm th}e}} \sim k_\perp\rho_e \frac{e\phi}{T_{0e}} \quad \Rightarrow \quad \frac{e\phi}{T_{0e}} \sim \frac{\epsilon}{k_\perp \rho_e}. \end{equation}

The density perturbations $\delta n_s$ are ordered anticipating a Boltzmann density response and the temperature perturbations $\delta T_s$ are assumed comparable to them:

(B9)\begin{equation} \frac{\delta T_{e}}{T_{0e}} \sim \frac{\delta T_i}{T_{0i}} \sim \frac{\delta n_{i}}{n_{0i}} = \frac{\delta n_{e}}{n_{0e}} \sim \frac{e\phi}{T_{0e}} \sim \frac{\epsilon}{k_\perp \rho_e}. \end{equation}

Finally, for the ordering of perpendicular magnetic-field perturbations, we demand that the effects of Lorentz tension (equivalently, of parallel compressions) must always be large enough to have an effect on the electron density perturbation, viz. [cf. (3.2)],

(B10)\begin{equation} \frac{\mathrm{d}}{\mathrm{d} t} \frac{\delta n_{e}}{n_{0e}} \sim \boldsymbol{\nabla}_\parallel u_{{\parallel} e} \sim \frac{c}{4{\rm \pi} e n_{0e}} \boldsymbol{\nabla}_\parallel \left[ {\boldsymbol{b}}_0 {\boldsymbol{\cdot}} \left({\boldsymbol{\nabla}}_\perp{\times} \delta \! {\boldsymbol{B}}_\perp \right) \right] \quad \Rightarrow \quad \frac{\delta \! {\boldsymbol{B}}_\perp}{B_0} \sim \beta_e\frac{k_\parallel \lambda_{ei}}{k_\perp \rho_e} \frac{e\phi}{T_{0e}}, \end{equation}

where $\beta _e = 8{\rm \pi} n_{0e} T_{0e}/B_0^2$ the electron plasma beta. The (compressive) parallel magnetic-field perturbations are ordered anticipating pressure balance:

(B11)\begin{equation} \frac{\delta \! B_\parallel}{B_0} = \frac{4{\rm \pi}}{B_0^2} \delta \left( \frac{B^2}{8{\rm \pi}} \right) \sim \frac{4{\rm \pi}}{B_0^2} \delta (n_s T_s) \sim \beta_e \frac{\delta T_{e}}{T_{0e}} \sim \frac{\epsilon \beta_e}{k_\perp \rho_e}. \end{equation}

The orderings (B5)–(B11) still allow for a choice of ordering for perpendicular wavenumbers $k_\perp$ with respect to the electron and ion Larmor radii. Given that we would like to obtain a set of electrostatic equations that exhibit the scale invariance discussed in § 2, we consider wavenumbers

(B12)\begin{equation} \beta_e\frac{\lambda_{ei}}{L} \ll k_\perp \rho_e \ll \frac{\lambda_{ei}}{L}, \end{equation}

for which the physical motivation is discussed in § 3.2. In terms of time scales, (B12) is equivalent to demanding that

(B13)\begin{equation} (k_\perp \rho_e)^2 \nu_{ee} \ll \omega \sim \omega_{*s} \sim k_\perp v_E \sim \kappa k_\parallel^2 \ll (k_\perp d_e )^2 \nu_{ei}, \end{equation}

where $d_e = \rho _e/\sqrt {\beta _e}$ is the electron inertial scale. We shall formalise (B12) by demanding that $k_\perp \rho _e \sim \sigma \lambda _{ei}/L$, where $\sigma$ is a placeholder constant satisfying $\beta _e \ll \sigma \ll 1$. In other words, $k_\perp \rho _\perp \sim 1$, where $\rho _\perp = \rho _e L/\lambda _{ei} \sigma$, an ‘intermediate’ spatial scale [cf. (3.11)].

To summarise, (B5)–(B12) imply the following ordering of frequencies:

(B14)\begin{equation} \frac{\omega}{\varOmega_e} \sim \sigma \frac{\lambda_{ei}}{L} \epsilon,\quad \frac{\omega}{\varOmega_i} \sim \frac{m_i}{m_e} \sigma \frac{\lambda_{ei}}{L} \epsilon; \end{equation}

length scales:

(B15)\begin{equation} k_\perp \rho_i \sim \sigma \frac{\lambda_{ei}}{L} \sqrt{\frac{m_i}{m_e}},\quad k_\perp \rho_e \sim \sigma \frac{\lambda_{ei}}{L}, \quad k_\parallel L \sim \sqrt{\sigma}, \quad \frac{k_\parallel}{k_\perp} \sim \frac{L}{\sqrt{\sigma}\lambda_{ei}} \epsilon; \end{equation}

and amplitudes:

(B16)\begin{equation} \frac{e\phi}{T_{0e}} \sim \frac{\delta n_{e}}{n_{0e}} \sim \frac{\delta n_{i}}{n_{0i}} \sim \frac{\delta T_{e}}{T_{0e}} \sim \frac{\delta T_i}{T_{0i}} \sim \frac{L}{\sigma\lambda_{ei}} \epsilon,\quad \frac{\delta \! {\boldsymbol{B}}_\perp}{B_0} \sim \frac{\beta_e}{\sqrt{\sigma}} \frac{e \phi}{T_{0e}}, \quad \frac{\delta \! B_\parallel}{B_0}\sim \beta_e \frac{e \phi}{T_{0e}}. \end{equation}

All relevant quantities are thus naturally ordered with respect to some combination of $m_e/m_i$, $\sigma$, $\lambda _{ei}/L$, and the gyrokinetic small parameter $\epsilon = \rho _e/L$. The above ordering of frequencies, length scales and amplitudes with respect to $\epsilon$ is the standard gyrokinetic ordering (see, e.g. Abel et al. Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013). We choose to treat the ordering in $\lambda _{ei}/L$ – the fact that this should be formally small following straightforwardly from, e.g. $\nu _{ei} \gg \omega _{*s}$ – as subsidiary to both the orderings in $\epsilon$ and in the mass ratio [see the first expression in (B15)], meaning that the formal hierarchy of our expansions is

(B17)\begin{equation} \epsilon \ll \sqrt{\frac{m_e}{m_i}} \ll \sigma \frac{\lambda_{ei}}{L} \ll 1, \end{equation}

with all other dimensionless parameters treated as finite.

B.2 Ion kinetics

Given that the ordering of perpendicular wavenumbers (B15) implies that $k_\perp \rho _i \gg 1$ under the expansion in the mass ratio, the ion distribution function $h_i$ will satisfy the gyrokinetic equation (see, e.g. Abel et al. Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013), rather than the drift-kinetic one (A2). It is straightforward to show (by, e.g. expanding the Bessel functions therein for $k_\perp \rho _i \gg 1$) that, to leading order in the mass-ratio expansion, the gyrokinetic equation is solved by

(B18)\begin{equation} h_i = 0. \end{equation}

The contributions to quasineutrality (A5) arising from the next-order solution will be of the size

(B19)\begin{equation} \frac{\left\langle h_i \right\rangle_{{\boldsymbol{r}}}}{f_{0i}} \sim \left\langle\left\langle\varphi \right\rangle_{{\boldsymbol{R}}_i} \right\rangle_{{\boldsymbol{r}}} \sim \frac{\varphi}{k_\perp \rho_i}, \end{equation}

which can be safely neglected. Thus, the ion dynamics do not enter anywhere into our equations, which is the approximation of ‘adiabatic ions’. Given that no further reference will be made to the ion temperature gradient $L_{T_i}$, we henceforth denote the electron temperature gradient $L_{T_e} = L_T$.

B.3 Electron fluid equations

We now proceed with our derivation of the electron fluid equations. It will turn out that the ordering (B15) of perpendicular length scales means that no FLR effects need be retained within our equations – these can only enter at second order within our expansion, but they are negligible even at this order (see appendix B.3.3). Furthermore, the ordering of the perpendicular and parallel magnetic field perturbations (B16) implies that both $\delta \! {\boldsymbol {B}}_\perp$ and $\delta \! B_\parallel$ can be neglected at all orders in our expansion. We thus adopt the drift-kinetic equation (A2) for $s=e$ as the starting point, expanding our distribution function $h_e$ in $\sigma \lambda _{ei}/L \ll 1$ as

(B20)\begin{equation} h_e = \sum_{n=0}^{\infty} h_e^{(n)}, \quad h_e^{(n)} \sim \left(\sigma \frac{\lambda_{ei}}{L} \right)^{n} \frac{e \phi}{T_{0e}} f_{0e}. \end{equation}
B.3.1 Zeroth order: perturbed Maxwellian

Given the ordering of timescales (B2), the collision operator on the right-hand side of (A2) is dominant to leading order:

(B21)\begin{equation} C_{ee}^{(l)}\left[ h_e^{(0)} \right] + \mathcal{L}_{ei}\left[ h_e^{(0)} \right] = 0, \end{equation}

where $C_{ee}^{(l)}$ is given by (A4) for $s = s' = e$, and

(B22)\begin{equation} \mathcal{L}_{ei} \left[ h_e \right] = \frac{\gamma_{ei} n_{0e}}{m_e^2} {\boldsymbol{\nabla}}_v \left[ f_{0e} {\boldsymbol{\cdot}} \left({\boldsymbol{\nabla}}_v {\boldsymbol{\nabla}}_v v \right) {\boldsymbol{\cdot}} {\boldsymbol{\nabla}}_v \frac{h_e}{f_{0e}} \right]\end{equation}

is the pitch-angle scattering (Lorentz) collision operator, valid to leading order in the mass ratio. We multiply (B21) by $h_e^{(0)}/f_{0e}$ and integrate over the entire phase space, yielding

(B23)\begin{equation} \int \frac{\mathrm{d}^3 {\boldsymbol{r}}}{V} \int \mathrm{d}^3 {\boldsymbol{v}} \,C_{ee}^{(l)}\left[ h_e^{(0)} \right] + \int \frac{\mathrm{d}^3 {\boldsymbol{r}}}{V} \int \mathrm{d}^3 {\boldsymbol{v}} \,\mathcal{L}_{ei}\left[ h_e^{(0)} \right] = 0. \end{equation}

Both terms in (B23) are negative definite and must vanish individually, meaning that the solution is constrained to be a perturbed Maxwellian with no mean flow (Helander & Sigmar Reference Helander and Sigmar2005), viz.,

(B24)\begin{equation} h_e^{(0)} = \left[\frac{\delta n_{e}}{n_{0e}} - \varphi + \frac{\delta T_{e}}{T_{0e}} \left(\frac{v^2}{v_{{\rm th}s}^2} - \frac{3}{2} \right) \right] f_{0e}, \end{equation}

where $\varphi = e\phi /T_{0e}$, and we have imposed the solvability conditions

(B25)\begin{equation} \int \mathrm{d}^3 {\boldsymbol{v}} \, h_e^{(n)} = \int \mathrm{d}^3 {\boldsymbol{v}} \, v^2 h_e^{(n)} = 0, \quad n \geqslant 1, \end{equation}

in order to determine uniquely the density $\delta n_{e}$ and temperature $\delta T_{e}$ moments in (B24). Note that, in general, the Lorentz collision operator constrains the electron distribution function to be isotropic in the frame moving with the parallel ion velocity. However, the parallel ion velocity is zero to all orders within our expansion in $\sigma \lambda _{ei}/L$ (given the adiabatic ion solution (B18)), meaning that the electron distribution function will have no parallel velocity moment to leading order.

We are now in a position to simplify the quasineutrality constraint (A5). Using the solutions (B18) and (B24), it straightforwardly becomes (3.7).

B.3.2 First order: parallel flows

The parallel flows are determined self-consistently from the leading-order perturbations at the next order in our expansion, viz., $h_e^{(1)}$ is determined by the solution of the Spitzer–Härm problem (Spitzer & Härm Reference Spitzer and Härm1953; Braginskii Reference Braginskii1965; Helander & Sigmar Reference Helander and Sigmar2005):

(B26)\begin{equation} v_\parallel \frac{\partial}{\partial z} \left[\left(\frac{\delta n_{e}}{n_{0e}} -\varphi + \frac{\delta T_{e}}{T_{0e}} \right) + \left( \frac{v^2}{v_{{\rm th}e}^2} - \frac{5}{2} \right) \frac{\delta T_{e}}{T_{0e}} \right]f_{0e} = C_{ee}^{(l)}\left[ h_e^{(1)} \right] + \mathcal{L}_{ei}\left[ h_e^{(1)} \right]. \end{equation}

This can be inverted for $h_e^{(1)}$ by means of a standard variational method. We define the functional:

(B27)\begin{align} \varSigma[h_e] &={-} \left\langle h_e, C_{ee}^{(l)} \left[ h_e\right] \right\rangle - \left\langle h_e, \mathcal{L}_{ei}\left[ h_e\right] \right\rangle \nonumber\\ &\quad + 2 \left\langle h_e, v_\parallel \frac{\partial}{\partial z} \left[\left(\frac{\delta n_{e}}{n_{0e}} -\varphi + \frac{\delta T_{e}}{T_{0e}} \right) + \left( \frac{v^2}{v_{{\rm th}e}^2} - \frac{5}{2} \right) \frac{\delta T_{e}}{T_{0e}} \right] f_{0e} \right\rangle, \end{align}

where $\langle \cdots, \cdots \rangle$ denotes an inner product in velocity space weighted by the inverse of the electron (Maxwellian) equilibrium $f_{0e}$. Then, considering small variations $h_e = h_\text {min} + \delta h$ and using the self-adjointness of the linearised collision operator, it is straightforward to show that the functional $\varSigma [h_e]$ has a minimum at $h_\text {min} = h_e^{(1)}$, for any variation $\delta h$ (see, e.g. Helander & Sigmar Reference Helander and Sigmar2005). Given that the spherical harmonics are eigenfunctions of the linearised collision operator, we choose to expand our distribution function in terms of spherical coordinates in velocity space $(x, \alpha, \beta )$, with $x = v^2 /v_{{\rm th}e}^2$, as

(B28)\begin{equation} h_e^{(1)} = \sum_{p=0}^\infty a_p L_p^{(3/2)}(x) v_\parallel f_{0e}(v) = \sum_{p=0}^\infty a_p L_p^{(3/2)}(x) v \cos \alpha f_{0e}(v), \end{equation}

where $L_p^{(3/2)}(x)$ are the generalised Laguerre polynomials and $a_p$ coefficients to be determined. Using this in (B27), one obtains

(B29) \begin{align} \varSigma \left[h_e^{(1)} \right] = n_{0e} v_{{\rm th}e}^2 \left[\sum_{p=0}^\infty \sum_{q=0}^\infty \frac{a_p a_q}{2}\right. & \left. \left( \nu_{ee} K_{pq}^{ee} + \nu_{ei} K_{pq}^{ei} \right) \right. \nonumber\\ & \left.\vphantom{\sum_{p=0}^\infty \sum_{q=0}^\infty} + a_0 \frac{\partial}{\partial z}\left(\frac{\delta n_{e}}{n_{0e}} -\varphi + \frac{\delta T_{e}}{T_{0e}} \right) - \frac{5}{2} a_1 \frac{\partial}{\partial z} \frac{\delta T_{e}}{T_{0e}}\right], \end{align}

where

(B30)\begin{align} K_{pq}^{ee} &={-} \frac{2}{n_e \nu_{ee}} \left\langle x^{1/2} L_p^{(3/2)}(x) f_{0e}(v) \cos \alpha, C_{ee}^{(l)} \left[x^{1/2} L_q^{(3/2)}(x) f_{0e}(v) \cos \alpha \right] \right\rangle, \end{align}
(B31)\begin{align}K_{pq}^{ei} &={-}\frac{2}{n_e \nu_{ei}} \left\langle x^{1/2} L_p^{(3/2)}(x) f_{0e}(v) \cos \alpha, \mathcal{L}_{ei} \left[x^{1/2} L_q^{(3/2)}(x) f_{0e}(v) \cos \alpha \right] \right\rangle , \end{align}

are the coefficients calculated in, e.g. Hardman et al. (Reference Hardman, Parra, Chong, Adkins, Anastopoulos-Tzanis, Barnes, Dickinson, Parisi and Wilson2022) (and references therein). Truncating (B28) at $p = 3$, and demanding that the functional (B29) be stationary with respect to variations in the coefficients $a_p$, we find that

(B32)\begin{equation} h_e^{(1)} = \left[ a_0 + a_1 L_1^{(3/2)}(x) + a_2 L_2^{(3/2)}(x) \right] v_\parallel f_{0e}, \end{equation}

where the coefficients are given by

(B33)\begin{align} \nu_{ei} a_0 &={-} \frac{\dfrac{217}{64} + \dfrac{151}{8\sqrt{2}Z} + \dfrac{9}{2Z^2}}{1 + \dfrac{61}{8 \sqrt{2}Z} + \dfrac{9}{2Z^2}}\frac{\partial}{\partial z}\left(\frac{\delta n_{e}}{n_{0e}} -\varphi + \frac{\delta T_{e}}{T_{0e}} \right) - \frac{\dfrac{5}{2}\left( \dfrac{33}{16} + \dfrac{45}{8 \sqrt{2}Z} \right)}{1 + \dfrac{61}{8 \sqrt{2}Z} + \dfrac{9}{2Z^2}} \frac{\partial}{\partial z} \frac{\delta T_{e}}{T_{0e}}, \end{align}
(B34)\begin{align}\nu_{ei} a_1 &= \frac{\dfrac{33}{16} + \dfrac{45}{8 \sqrt{2}Z} }{1 + \dfrac{61}{8 \sqrt{2}Z} + \dfrac{9}{2Z^2}}\frac{\partial}{\partial z} \left(\frac{\delta n_{e}}{n_{0e}} -\varphi + \frac{\delta T_{e}}{T_{0e}} \right) + \frac{ \dfrac{5}{2} \left( \dfrac{13}{4} + \dfrac{45}{8\sqrt{2}Z} \right) }{1 + \dfrac{61}{8 \sqrt{2}Z} + \dfrac{9}{2Z^2}} \frac{\partial}{\partial z} \frac{\delta T_{e}}{T_{0e}}, \end{align}
(B35)\begin{align}\nu_{ei} a_2 &={-} \frac{\dfrac{3}{8} - \dfrac{3}{2\sqrt{2}Z}}{1 + \dfrac{61}{8 \sqrt{2}Z} + \dfrac{9}{2Z^2}}\frac{\partial}{\partial z} \left(\frac{\delta n_{e}}{n_{0e}} -\varphi + \frac{\delta T_{e}}{T_{0e}} \right)- \frac{\dfrac{5}{2 }\left(\dfrac{3}{2} + \dfrac{3}{2\sqrt{2}Z} \right)}{1 + \dfrac{61}{8 \sqrt{2}Z} + \dfrac{9}{2Z^2}} \frac{\partial}{\partial z} \frac{\delta T_{e}}{T_{0e}}, \end{align}

which can easily be shown to satisfy the Onsager (Reference Onsager1931) relations.

The solution (B32) allows us to determine, subject to the solvability condition

(B36)\begin{equation} \int \mathrm{d}^3 {\boldsymbol{v}} \, v_\parallel h_{e}^{(n)} = 0, \quad n \geqslant 2, \end{equation}

the parallel electron flow:

(B37)\begin{equation} u_{{\parallel} e} = \frac{1}{n_{0e}} \int \mathrm{d}^3 {\boldsymbol{v}} \,v_\parallel h_e^{(1)}. \end{equation}

Using (B32) for $h_e^{(1)}$ in (B37) and defining the (ion-charge-dependent) coefficients [cf. for $Z=1$, (C16) and (C17) in Hardman et al. Reference Hardman, Parra, Chong, Adkins, Anastopoulos-Tzanis, Barnes, Dickinson, Parisi and Wilson2022]

(B38)\begin{equation} c_1 \!=\! \frac{\dfrac{217}{64} + \dfrac{151}{8\sqrt{2}Z} + \dfrac{9}{2Z^2}}{1 + \dfrac{61}{8 \sqrt{2}Z} + \dfrac{9}{2Z^2}},\quad c_2 \!=\! \frac{\dfrac{5}{2}\left( \dfrac{33}{16} + \dfrac{45}{8 \sqrt{2}Z} \right)}{1 + \dfrac{61}{8 \sqrt{2}Z} + \dfrac{9}{2Z^2}},\quad c_3 \!=\! \frac{\dfrac{25}{4} \left( \dfrac{13}{4} + \dfrac{45}{8 \sqrt{2} Z} \right)}{1 + \dfrac{61}{8 \sqrt{2}Z} + \dfrac{9}{2Z^2}} - \frac{c_2^2}{c_1}, \end{equation}

we obtain (3.3).

B.3.3 Second order: density and temperature evolution

At second order, the electron drift-kinetic equation

(B39)\begin{align} & \frac{\mathrm{d}}{\mathrm{d} t} \left( h_e^{(0)} + \varphi f_{0e} \right) + v_\parallel \frac{\partial h_e^{(1)}}{\partial z} + {\boldsymbol{v}}_{de} {\boldsymbol{\cdot}} {\boldsymbol{\nabla}}_\perp h_e^{(0)} + \frac{\rho_e v_{{\rm th}e}}{2 L_T} \frac{\partial \varphi}{\partial y} \left( \frac{v^2}{v_{{\rm th}e}^2} - \frac{3}{2} \right) f_{0e} \nonumber\\ &\quad = C^{(l)}_{e e }\left[ h_e^{(2)} \right] + \mathcal{L}_{e i } \left[h_e^{(2)} \right], \end{align}

describes the evolution of the density and temperature perturbations in (B24). When taking the density and temperature moments of (B39), the contributions from the collision operator on the right-hand side vanish, as the electron–electron and Lorentz collision operators conserve particle number and energy to this order in our expansion. An observant reader may have noticed, however, that in starting our expansion from the drift-kinetic equation (A2), we ruled out the possibility of retaining higher-order collisional terms due to FLR motions of the electrons. Indeed, if we had instead considered the gyrokinetic equation (see, e.g. Abel et al. Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013) and expanded the oscillatory exponential factors $\exp ({\pm {\rm i} {\boldsymbol {k}} {\boldsymbol {\cdot }} {\boldsymbol {\rho }}_e})$ arising from the presence of the gyroaverages of the collision operator on its right-hand side, we would have obtained terms of the form ${\sim }\nu _{ee} \rho _e^2 {\boldsymbol {\nabla }}_\perp ^2 h_e^{(0)}$ at order $(k_\perp \rho _e)^2$. These represent electron thermal diffusion (cf. Newton et al. Reference Newton, Cowley and Loureiro2010; Ivanov et al. Reference Ivanov, Schekochihin, Dorland, Field and Parra2020; Hardman et al. Reference Hardman, Parra, Chong, Adkins, Anastopoulos-Tzanis, Barnes, Dickinson, Parisi and Wilson2022). Recalling the ordering of timescales (B13), however, it is clear that these terms are negligible in comparison with those on the left-hand side of (B39) – this justifies post factum the choice to perform our expansion starting from the drift-kinetic equation (A2), rather than from the gyrokinetic one.

Thus, taking the density moment of (B39), employing the identity [see (A3) and (B1)]

(B40)\begin{equation} {\boldsymbol{v}}_{de} {\boldsymbol{\cdot}} {\boldsymbol{\nabla}}_\perp h_e^{(0)} = \frac{\rho_e v_{{\rm th}e}}{2} \left( \frac{2}{R} \frac{v_\parallel^2}{v_{{\rm th}s}^2} + \frac{1}{L_B} \frac{v_\perp^2}{v_{{\rm th}s}^2} \right) \frac{\partial h_e^{(0)}}{\partial y}, \end{equation}

and making use of the fact that $R = L_B$ for a low-beta plasma, we find:

(B41)\begin{equation} \frac{\mathrm{d}}{\mathrm{d} t} \frac{\delta n_{e}}{n_{0e}}+ \frac{\partial u_{{\parallel} e}}{\partial z} + \frac{\rho_e v_{{\rm th}e}}{L_B} \frac{\partial}{\partial y} \left( \frac{\delta n_{e}}{n_{0e}} - \varphi + \frac{\delta T_{e}}{T_{0e}} \right) =0. \end{equation}

For the temperature moment, we first note that, from (B32),

(B42)\begin{equation} \frac{1}{n_{0e}} \int \mathrm{d}^3 {\boldsymbol{v}} \,v_\parallel \left( \frac{v^2}{v_{{\rm th}e}^2} - \frac{3}{2}\right) h_e^{(1)} = \left( 1 + \frac{c_2}{c_1} \right) u_{{\parallel} e} + \frac{\delta q_e}{n_{0e} T_{0e}}, \end{equation}

where $\delta q_e$ is defined in (3.6). Therefore, taking the temperature moment of (B39) and dividing throughout by $3/2$ yields

(B43)\begin{align} & \frac{\mathrm{d}}{\mathrm{d} t} \frac{\delta T_{e}}{T_{0e}} + \frac{2}{3} \frac{\partial}{\partial z} \frac{\delta q_e}{n_{0e} T_{0e}} + \frac{2}{3} \left( 1+ \frac{c_2}{c_1} \right) \frac{\partial u_{{\parallel} e}}{\partial z} + \frac{2}{3} \frac{\rho_e v_{{\rm th}e}}{L_B} \frac{\partial}{\partial y} \left( \frac{\delta n_{e}}{n_{0e}} - \varphi + \frac{7}{2} \frac{\delta T_{e}}{T_{0e}} \right) \nonumber\\ &\quad ={-}\frac{\rho_e v_{{\rm th}e}}{2 L_T} \frac{\partial \varphi}{\partial y}. \end{align}

Neglecting the magnetic drifts, one straightforwardly obtains (3.2) and (3.4) from (B41) and (B43), respectively.

Appendix C. Case with finite magnetic-field gradients

This appendix details the behaviour of our model system in the presence of magnetic drifts associated with an inhomogeneous equilibrium magnetic field. Assuming its scale length $L_B$ [see (B1)] to be constant across the domain, our evolution equations for the electrostatic potential and temperature perturbations are now [see (B41) and (B43)]:

(C1)\begin{align} \frac{\partial}{\partial t} \bar{\tau}^{{-}1} \varphi & - \frac{c_1 v_{{\rm th}e}^2}{2 \nu_{ei}} \frac{\partial^2 }{\partial z^2}\left[\left(1 + \frac{1}{\bar{\tau}} \right)\varphi - \left( 1 + \frac{c_2}{c_1} \right) \frac{\delta T_{e}}{T_{0e}} \right] \nonumber\\ & + \frac{\rho_e v_{{\rm th}e}}{L_B} \frac{\partial}{\partial y} \left[ \left( 1 + \frac{1}{\bar{\tau}} \right) \varphi - \frac{\delta T_{e}}{T_{0e}} \right]= 0, \end{align}
(C2)\begin{align} \frac{\mathrm{d}}{\mathrm{d} t} \frac{\delta T_{e}}{T_{0e}} &+ \frac{2}{3 } \frac{c_1 v_{{\rm th}e}^2}{2 \nu_{ei}} \frac{\partial^2}{\partial z^2} \left\{\left( 1 + \frac{1}{\bar{\tau}} \right) \left( 1 + \frac{c_2}{c_1} \right) \varphi - \left[\frac{c_3}{c_1} + \left( 1 + \frac{c_2}{c_1} \right)^2 \right] \frac{\delta T_{e}}{T_{0e}} \right\} \nonumber\\ & - \frac{2}{3} \frac{\rho_e v_{{\rm th}e}}{L_B} \frac{\partial}{\partial y} \left[ \left( 1 + \frac{1}{\bar{\tau}} \right) \varphi - \frac{7}{2}\frac{\delta T_{e}}{T_{0e}}\right]={-} \frac{\rho_e v_{{\rm th}e}}{2 L_T}\frac{\partial \varphi}{\partial y}. \end{align}

The reduction of these equations to (3.8)–(3.9) occurs for very steep electron-temperature gradients, in the limit $L_B/L_T \rightarrow \infty$.

The presence of the magnetic-drift terms in (C1)–(C2) introduces another instability into the system, the curvature-mediated ETG (cETG) instability (Horton et al. Reference Horton, Hong and Tang1988; Adkins et al. Reference Adkins, Schekochihin, Ivanov and Roach2022), which can, and generally does, modify its turbulent-transport properties. In particular, the turbulence theory of § 5 assumed that the sETG instability was the dominant source of energy injection; this is only the case at sufficiently large $L_B/L_T$, meaning that we would expect departures from the behaviour observed in § 5 to be most significant for $L_B/L_T$ of order unity. A series of simulations were conducted in which $L_B/L_T$ was varied, with all other parameters being kept the same as in the baseline simulation (see table 1); the heat flux from these simulations is plotted in figure 13. It is readily apparent that the introduction of finite magnetic-field gradients leads to a failure of saturation for all simulations where $L_B/L_T$ is above the linear critical gradient for the cETG instability:

(C3)\begin{equation} \frac{L_B}{L_T} > \frac{1}{2} \left( \bar{\tau} + \frac{40}{9} \frac{1}{\bar{\tau}^2} \right), \end{equation}

a threshold that can be derived straightforwardly from (C1) and (C2) in the 2D limit. This lack of saturation appears to persist irrespective of changes in box size, aspect ratio, and resolution in any (or all) of the coordinate directions.

Figure 13. Time traces of the instantaneous heat flux from simulations with finite $L_B/L_T$, with the limit of $L_B/L_T \rightarrow \infty$ shown for comparison. All parameters are the same as the baseline simulation (see table 1), and the heat flux is normalised to $(\rho _\perp /\rho _e) Q_{\text {gB}e}$. The heat flux grows without bound in all simulations with (finite) $L_B/L_T$ above the linear critical gradient (C3) (${\approx }2.72$ for $\bar {\tau }= 1$), with the rate of divergence decreasing as $L_B/L_T$ is increased.

As discussed in § 5.5, the adiabatic ion response (3.7) causes the nonlinearity in the electron-scale continuity equation (3.8) to vanish identically, a property that is shared by (C1), meaning that the system lacks any nonlinearity capable of generating 2D secondary instabilities that are responsible for the generation of zonal flows and destruction of streamer structures (see references in § 5.5). Indeed, the lack of saturation in the case of our ETG simulations appears to be due to the inability of the system to break apart the streamers created by the cETG instability; the existence of such streamers causes the heat flux to diverge as they ‘short circuit’ the heat transport across the radial domain. Even if the simulation initially appears to saturate after the linear phase, it eventually forms these large-scale streamers, which appear to be immune to all types of nonlinear shearing, as seen clearly in the real-space snapshots of cETG turbulence shown in figures 14 and 15.

Figure 14. Real-space snapshots of the (a) electrostatic potential and (b) temperature perturbations from the simulation with $L_B/L_T = 1000$ from figure 13, taken at $(\rho _e/\rho _\perp )^2\nu _{ei} t/2 = 200$. The coordinate axes are as shown, while the red and blue colours correspond to regions of positive and negative fluctuation amplitude. At these early times, the turbulence appears similar to that of saturated sETG turbulence for $L_B/L_T \rightarrow \infty$ (cf. figure 11), despite the eventual lack of saturation (see figure 15).

Figure 15. The same as figure 14, except taken at $(\rho _e/\rho _\perp )^2\nu _{ei} t/2 = 1000$. The unbounded growth of the heat flux is associated with the formation of large-scale, approximately 2D streamer structures that appear to be immune to all types of nonlinear shearing.

This perhaps confirms the view of Hammett et al. (Reference Hammett, Beer, Dorland, Cowley and Smith1993) that the adiabatic ion response (3.7) is insufficient to saturate ETG-scale turbulence in the presence of finite magnetic-field gradients, and one may have to resort to more inclusive closures for the ions. One such closure including scales comparable to the ion-Larmor radius is (see, e.g. Adkins et al. Reference Adkins, Schekochihin, Ivanov and Roach2022)

(C4)\begin{equation} \frac{\delta n_{e}}{n_{0e}} ={-} \bar{\tau}^{{-}1} \varphi + \frac{1}{n_{0i}} \int \mathrm{d}^3 {\boldsymbol{v}} \left\langle g_i \right\rangle_{{\boldsymbol{r}}}, \end{equation}

where $\bar {\tau }^{-1}$ is now an operator defined as follows:

(C5)\begin{equation} - \bar{\tau}^{{-}1} \varphi ={-} \frac{Z}{\tau} (1 - \hat{\varGamma}_0) \varphi \approx \left\{\begin{array}{@{}cc} \displaystyle\dfrac{Z}{2\tau} \rho_i^2 {\boldsymbol{\nabla}}_\perp^2 \varphi, & k_\perp \rho_i \ll 1,\\ \displaystyle - \dfrac{Z}{\tau} \varphi, & k_\perp \rho_i \gg 1, \end{array}\right. \end{equation}

and the operator $\hat {\varGamma }_0$ can be expressed in Fourier space in terms of the modified Bessel function of the first kind: $\varGamma _0 = I_0(\alpha _i)\,{\rm e}^{-\alpha _i}$, where $\alpha _i = (k_\perp \rho _i)^2/2$. The presence of the non-adiabatic ion distribution function $g_i$ in (C4), however, means that one would have also to include a self-consistent treatment of ions in order to make use of this closure ($g_i =0$ is not a solution to the ion gyrokinetic equation in the presence of finite magnetic drifts). Should this, or other similar closures, allow for saturation, this would imply that one must always appeal to (elements of) ion-scale physics for saturation of electrostatic cETG-driven turbulence. The extent to which such considerations are practically relevant, however, depends on whether or not the system being considered contains any electromagnetic physics, and thus on the value of the (electron) plasma beta $\beta _e$. Indeed, for $\beta _e \gtrsim m_e/m_i$, the ‘flux-freezing scale’ $d_e = \rho _e/\sqrt {\beta _e}$ is encountered before (i.e. is smaller than) the ion Larmor radius $\rho _i$ when moving towards larger perpendicular scales. Provided that the wavenumber interval between $d_e$ and $\rho _i$ is sufficiently wide to allow for the presence of electron-scale, electromagnetic instabilities, the mechanisms of saturation in such a system could be vastly different than in the electrostatic regime. This is a subject of ongoing research.

Footnotes

1 Implicit in (3.10) is the assumption that the electron inertial scale $d_e = \rho _e/\sqrt {\beta _e}$ is smaller than the ion Larmor radius $\rho _i$; this is only the case if $\beta _e \gtrsim m_e/m_i$, which is well-satisfied in most systems of interest. Should $d_e$ lie on the large-scale side of $\rho _i$ (i.e. $\beta _e \lesssim m_e/m_i$), the lower bound in (3.12) must be replaced with $\sqrt {m_e/m_i}$.

2 Refining our consideration beyond this assumption of stationarity, we observe that the characteristic timescale of the fluctuations of the instantaneous heat flux increases with the parallel system size – this is manifest in figure 2, where the simulations with larger $L_\parallel \sqrt {\sigma }/L_T$ exhibit higher-amplitude, longer-timescale fluctuations. The origin of this trend can be understood as follows. Relaxing the assumption of stationarity, instead of (2.6), we have, from (2.5), $\tilde {Q}_s(\lambda ^{2/\alpha } L_\parallel, t/\lambda ^2) = \lambda ^2 Q_s (L_\parallel, t)$. If $Q_s$ exhibits fluctuations on some characteristic timescale $\tau$, then, if we assume that that both solutions must be periodic with the same period, the corresponding timescale for the transformed heat flux will be $\tilde {\tau } = \lambda ^2 \tau$. Given that the parallel system sizes for both solutions are related by $\tilde {L}_{\parallel } = \lambda ^{2/\alpha }L_\parallel$, it follows that $\tau \propto L_\parallel ^\alpha$. This dependence was confirmed numerically for the set of simulations shown in figure 2.

3 In standard turbulence literature, the outer scale is often defined to be the integral scale of the 1D perpendicular energy spectrum (5.11), viz., $k_\perp ^o \equiv \int _0^\infty \mathrm {d} k_\perp \, E_\perp ^\varphi (k_\perp )/ \int _0^\infty \mathrm {d} k_\perp \,k_\perp ^{-1} E_\perp ^\varphi (k_\perp )$. However, given that, physically, we are interested in the outer scale as the scale at which the free energy is predominantly injected, the choice to maximise (5.6) seems to be better motivated physically.

4 Although it remains to be seen whether the dominant effect in enforcing (5.14) is plain linear dissipation or its suppression via stochastic echos (Schekochihin et al. Reference Schekochihin, Parker, Highcock, Dellar, Dorland and Hammett2016; Adkins & Schekochihin Reference Adkins and Schekochihin2018).

5 Our system does of course also have parallel dissipation via heat conduction, at the rate $ \sim \omega _\parallel$, but $\omega _\parallel$ decreases with increasing parallel scale and, at any rate, does not break scale invariance, so cannot set $k_\parallel ^o$.

6 Though not one that can be taken for granted. For example, Hosking & Schekochihin (Reference Hosking and Schekochihin2022) (see also appendix C of Schekochihin Reference Schekochihin2022) showed that a $k_\perp ^1$ scaling could emerge instead through a balance between turbulent diffusion at large scales and the nonlinear ‘source’ that would otherwise give rise to the $k_\perp ^3$ scaling. Let us estimate the rate of turbulent diffusion in our system. The dominant contribution to the turbulent-diffusion coefficient $D$ will be from $k_\perp \sim k_\parallel ^{3/2}$, which, at any given $k_\parallel$, plays the role of the energy-containing scale. Then $D \sim v_E^2 t_\text {nl} \sim \tilde {\omega }_\parallel /{\tilde {k}_\perp }^2$, where have used the critical-balance condition (5.14), and the tildes denote quantities evaluated at $\tilde {k}_\perp \sim k_\parallel ^{3/2} \gg k_\perp$, where $k_\perp$ is the wavenumber at which turbulent diffusion is acting. The rate of turbulent diffusion at this wavenumber will thus be $k_\perp ^2 D \sim \tilde {\omega }_\parallel (k_\perp /\tilde {k}_\perp )^2 \ll \omega _\parallel$. Turbulent diffusion is, therefore, negligible, and so we are justified in adopting the $k_\perp ^3$ scaling. Note that the survival of the $k_\perp ^3$ scaling is a noteworthy feature of our system, where the dynamics at $k_\parallel \gg k_\perp ^{2/3}$ are dominated by parallel dissipation due to thermal conductivity and do not produce significant turbulent diffusion – unlike waves, which, e.g. in reduced hydrodynamics (RMHD), do (Schekochihin Reference Schekochihin2022).

7 Schekochihin et al. (Reference Schekochihin, Parker, Highcock, Dellar, Dorland and Hammett2016) obtained, by a similar method, an analogous result for long-wavelength electrostatic ITG turbulence (which, in this approach, is no different for ETG). Specifically, they found that $a = 5$ and $c = 11/3$ – this was a consequence of the fact that they considered collisionless turbulence, for which the critical-balance condition is ${\omega _\parallel } \sim k_\parallel v_{{\rm th}e} \sim t_\text {nl}^{-1}$ implying that $k_\parallel \sim k_\perp ^{4/3}$.

References

Abel, I.G. & Cowley, S.C. 2013 Multiscale gyrokinetics for rotating tokamak plasmas: II. Reduced models for electron dynamics. New J. Phys. 15, 023041.CrossRefGoogle Scholar
Abel, I.G., Plunk, G.G., Wang, E., Barnes, M., Cowley, S.C., Dorland, W. & Schekochihin, A.A. 2013 Multiscale gyrokinetics for rotating tokamak plasmas: fluctuations, transport and energy flows. Rep. Prog. Phys. 76, 116201.CrossRefGoogle ScholarPubMed
Adkins, T., Schekochihin, A.A., Ivanov, P.G. & Roach, C.M. 2022 Electromagnetic instabilities and plasma turbulence driven by electron-temperature gradient. J. Plasma Phys. 88, 905880410.CrossRefGoogle Scholar
Adkins, T.G. 2023 Electromagnetic instabilities and plasma turbulence driven by the electron-temperature gradient. PhD thesis, University of Oxford.CrossRefGoogle Scholar
Adkins, T.G. & Schekochihin, A.A. 2018 A solvable model of Vlasov-kinetic plasma turbulence in Fourier-Hermite phase space. J. Plasma Phys. 84, 905840107.CrossRefGoogle Scholar
Barnes, M., Parra, F.I. & Schekochihin, A.A. 2011 Critically balanced ion temperature gradient turbulence in fusion plasmas. Phys. Rev. Lett. 107, 115003.CrossRefGoogle ScholarPubMed
Beer, M.A. & Hammett, G.W. 1996 Toroidal gyrofluid equations for simulations of tokamak turbulence. Phys. Plasmas 3, 4046.CrossRefGoogle Scholar
Boldyrev, S. 2005 On the spectrum of magnetohydrodynamic turbulence. Astrophys. J. 626, L37.CrossRefGoogle Scholar
Braginskii, S.I. 1965 Transport processes in a plasma. Rev. Plasma Phys. 1, 205.Google Scholar
Candy, J., Waltz, R.E., Fahey, M.R. & Holland, C. 2007 The effect of ion-scale dynamics on electron-temperature-gradient turbulence. Plasma Phys. Control. Fusion 49, 1209.CrossRefGoogle Scholar
Chapman-Oplopoiou, B., Hatch, D.R., Field, A.R., Frassinetti, L., Hillesheim, J., Horvath, L., Maggi, C.F., Parisi, J., Roach, C.M., Saarelma, S., et al. 2022 The role of ETG modes in JET–ILW pedestals with varying levels of power and fuelling. Nucl. Fusion 62, 086028.CrossRefGoogle Scholar
Cho, J. & Lazarian, A. 2004 The anisotropy of electron magnetohydrodynamic turbulence. Astrophys. J. 615, L41.CrossRefGoogle Scholar
Colyer, G.J., Schekochihin, A.A., Parra, F.I., Roach, C.M., Barnes, M.A., Ghim, Y.-c. & Dorland, W. 2017 Collisionality scaling of the electron heat flux in ETG turbulence. Plasma Phys. Control. Fusion 59, 055002.CrossRefGoogle Scholar
Connor, J.W. & Taylor, J.B. 1977 Scaling laws for plasma confinement. Nucl. Fusion 17, 1047.CrossRefGoogle Scholar
Cowley, S.C., Kulsrud, R.M. & Sudan, R. 1991 Considerations of ion-temperature-gradient-driven turbulence. Phys. Fluids B 3, 2767.CrossRefGoogle Scholar
Diamond, P.H., Itoh, S. -I., Itoh, K. & Hahm, T.S. 2005 Zonal flows in plasma—a review. Plasma Phys. Control. Fusion 47, 35.CrossRefGoogle Scholar
Dorland, W. & Hammett, G.W. 1993 Gyrofluid turbulence models with kinetic effects. Phys. Fluids B 5, 812.CrossRefGoogle Scholar
Dorland, W., Jenko, F., Kotschenreuther, M. & Rogers, B.N. 2000 Electron temperature gradient turbulence. Phys. Rev. Lett. 85, 5579.CrossRefGoogle ScholarPubMed
Field, A., Chapman-Oplopoiou, B., Connor, J., Frassinetti, L., Hatch, D., Roach, C., Saarelma, S., JET Contributors 2023 Comparing pedestal structure in JET-ILW H-mode plasmas with a model for stiff ETG turbulent heat transport. Phil. Trans. R. Soc. Lond. A 381, 20210228.Google Scholar
Goldreich, P. & Sridhar, S. 1995 Toward a theory of interstellar turbulence. 2: Strong Alfvénic turbulence. Astrophys. J. 438, 763.CrossRefGoogle Scholar
Guttenfelder, W., Battaglia, D.J., Belova, E., Bertelli, N., Boyer, M.D., Chang, C.S., Diallo, A., Duarte, V.N., Ebrahimi, F., Emdee, E.D., et al. 2022 NSTX-U theory, modeling and analysis results. Nucl. Fusion 62, 042023.CrossRefGoogle Scholar
Guttenfelder, W. & Candy, J. 2011 Resolving electron scale turbulence in spherical tokamaks with flow shear. Phys. Plasmas 18, 022506.CrossRefGoogle Scholar
Guttenfelder, W., Groebner, R.J., Canik, J.M., Grierson, B.A., Belli, E.A. & Candy, J. 2021 Testing predictions of electron scale turbulent pedestal transport in two DIII-D ELMy H-modes. Nucl. Fusion 61, 056005.CrossRefGoogle Scholar
Hammett, G.W., Beer, M.A., Dorland, W., Cowley, S.C. & Smith, S.A. 1993 Developments in the gyrofluid approach to tokamak turbulence simulations. Plasma Phys. Control. Fusion 35, 973.CrossRefGoogle Scholar
Hammett, G.W., Dorland, W. & Perkins, F.W. 1992 Fluid models of phase mixing, Landau damping, and nonlinear gyrokinetic dynamics. Phys. Fluids B 4, 2052.CrossRefGoogle Scholar
Hammett, G.W. & Perkins, F.W. 1990 Fluid moment models for Landau damping with application to the ion-temperature-gradient instability. Phys. Rev. Lett. 64, 3019.CrossRefGoogle Scholar
Hardman, M.R., Parra, F.I., Chong, C., Adkins, T., Anastopoulos-Tzanis, M.C., Barnes, M., Dickinson, D., Parisi, J.F. & Wilson, H.R. 2022 Extended electron tails in electrostatic microinstabilities and the nonadiabatic response of passing electrons. Plasma Phys. Control. Fusion 64, 055004.CrossRefGoogle Scholar
Hardman, M.R., Parra, F.I., Patel, B.S., Roach, C.M., Ruiz Ruiz, J., Barnes, M., Dickinson, D., Dorland, W., Parisi, J.F., St-Onge, D., et al. 2023 New stability parameter to describe low-$\beta$ electromagnetic microinstabilities driven by passing electrons in axisymmetric toroidal geometry. Plasma Phys. Control. Fusion 65, 045011.CrossRefGoogle Scholar
Hasegawa, A. & Mima, K. 1978 Pseudo-three-dimensional turbulence in magnetized nonuniform plasma. Phys. Fluids 21, 87.CrossRefGoogle Scholar
Hasegawa, A. & Wakatani, M 1983 Plasma edge turbulence. Phys. Rev. Lett. 50, 682.CrossRefGoogle Scholar
Hatch, D.R., Kotschenreuther, M., Mahajan, S.M., Merlo, G., Field, A.R., Giroud, C., Hillesheim, J.C., Maggi, C.F., Perez von Thun, C., Roach, C.M., et al. 2019 Direct gyrokinetic comparison of pedestal transport in JET with carbon and ITER-like walls. Nucl. Fusion 59, 086056.CrossRefGoogle Scholar
Hatch, D.R., Michoski, C., Kuang, D., Chapman-Oplopoiou, B., Curie, M., Halfmoon, M., Hassan, E., Kotschenreuther, M., Mahajan, S.M., Merlo, G., et al. 2022 Reduced models for ETG transport in the tokamak pedestal. Phys. Plasmas 29, 062501.CrossRefGoogle Scholar
Hatch, D.R., Terry, P.W., Jenko, F., Merz, F. & Nevins, W.M. 2011 Saturation of gyrokinetic turbulence through damped eigenmodes. Phys. Rev. Lett. 106, 115003.CrossRefGoogle ScholarPubMed
Helander, P. & Sigmar, D.J. 2005 Collisional Transport in Magnetized Plasmas. Cambridge University Press.Google Scholar
Horton, W., Hong, B.G. & Tang, W.M. 1988 Toroidal electron temperature gradient driven drift modes. Phys. Fluids 31, 2971.CrossRefGoogle Scholar
Hosking, D.N. & Schekochihin, A.A. 2022 Emergence of long-range correlations and thermal spectra in forced turbulence. arXiv:2202.00462.Google Scholar
Howes, G.G., Cowley, S.C., Dorland, W., Hammett, G.W., Quataert, E. & Schekochihin, A.A. 2006 Astrophysical gyrokinetics: basic equations and linear theory. Astrophys. J. 651, 590.CrossRefGoogle Scholar
Ivanov, P.G., Schekochihin, A.A. & Dorland, W. 2022 Dimits transition in three-dimensional ion-temperature-gradient turbulence. J. Plasma Phys. 88, 905880506.CrossRefGoogle Scholar
Ivanov, P.G., Schekochihin, A.A., Dorland, W., Field, A.R. & Parra, F.I. 2020 Zonally dominated dynamics and Dimits threshold in curvature-driven ITG turbulence. J. Plasma Phys. 86, 855860502.CrossRefGoogle Scholar
Jenko, F., Dorland, W., Kotschenreuther, M. & Rogers, B.N. 2000 Electron temperature gradient driven turbulence. Phys. Plasmas 7, 1904.CrossRefGoogle Scholar
Joiner, N., Applegate, D., Cowley, S.C., Dorland, W. & Roach, C.M. 2006 Electron temperature gradient driven transport in a MAST H-mode plasma. Plasma Phys. Control. Fusion 48, 685.CrossRefGoogle Scholar
Kolmogorov, A.N. 1941 Local structure of turbulence in incompressible viscous fluid at very large Reynolds numbers. Dokl. Acad. Nauk SSSR 30, 299.Google Scholar
Kotschenreuther, M., Dorland, W., Beer, M.A. & Hammett, G.W. 1995 a Quantitative predictions of tokamak energy confinement from first-principles simulations with kinetic effects. Phys. Plasmas 2, 2381.CrossRefGoogle Scholar
Kotschenreuther, M., Rewoldt, G. & Tang, W.M. 1995 b Comparison of initial value and eigenvalue codes for kinetic toroidal plasma instabilities. Comput. Phys. Commun. 88, 128.CrossRefGoogle Scholar
Landau, L. 1946 On the vibrations of the electronic plasma. Zh. Eksp. Teor. Fiz. 16, 574.Google Scholar
Lee, Y.C., Dong, J.Q., Guzdar, P.N. & Liu, C.S. 1987 Collisionless electron temperature gradient instability. Phys. Fluids 30, 1331.CrossRefGoogle Scholar
Liu, C.S. 1971 Instabilities in a magnetoplasma with skin current. Phys. Rev. Lett. 27, 1637.CrossRefGoogle Scholar
Nazarenko, S.V. & Schekochihin, A.A. 2011 Critical balance in magnetohydrodynamic, rotating and stratified turbulence: towards a universal scaling conjecture. J. Fluid Mech. 677, 134.CrossRefGoogle Scholar
Newton, S.L., Cowley, S.C. & Loureiro, N.F. 2010 Understanding the effect of sheared flow on microinstabilities. Plasma Phys. Control. Fusion 52, 125001.CrossRefGoogle Scholar
Onsager, L. 1931 Reciprocal relations in irreversible processes. I. Phys. Rev. 37, 405.CrossRefGoogle Scholar
Orszag, S.A. 1971 On the elimination of aliasing in finite-difference schemes by filtering high-wavenumber components. J. Atmos. Sci. 28, 1074.2.0.CO;2>CrossRefGoogle Scholar
Parisi, J.F., Parra, F.I., Roach, C.M., Giroud, C., Dorland, W., Hatch, D.R., Barnes, M., Hillesheim, J.C., Aiba, N., Ball, J., et al. 2020 Toroidal and slab ETG instability dominance in the linear spectrum of JET-ILW pedestals. Nucl. Fusion 60, 126045.CrossRefGoogle Scholar
Parisi, J.F., Parra, F.I., Roach, C.M., Hardman, M.R., Schekochihin, A.A., Abel, I.G., Aiba, N., Ball, J., Barnes, M., Chapman-Oplopoiou, B., et al. 2022 Three-dimensional inhomogeneity of electron-temperature-gradient turbulence in the edge of tokamak plasmas. Nucl. Fusion 62, 086045.CrossRefGoogle Scholar
Patel, B.S., Dickinson, D., Roach, C.M. & Wilson, H.R. 2021 Linear gyrokinetic stability of a high $\beta$ non-inductive spherical tokamak. Nucl. Fusion 62, 016009.CrossRefGoogle Scholar
Ren, Y., Belova, E., Gorelenkov, N., Guttenfelder, W., Kaye, S.M., Mazzucato, E., Peterson, J.L., Smith, D.R., Stutman, D., Tritz, K., et al. 2017 Recent progress in understanding electron thermal transport in NSTX. Nucl. Fusion 57, 072002.CrossRefGoogle Scholar
Roach, C.M., Abel, I.G., Akers, R.J., Arter, W., Barnes, M., Camenen, Y., Casson, F.J., Colyer, G., Connor, J.W., Cowley, S.C., et al. 2009 Gyrokinetic simulations of spherical tokamaks. Plasma Phys. Control. Fusion 51, 124020.CrossRefGoogle Scholar
Roberg-Clark, G.T., Plunk, G.G. & Xanthopoulos, P. 2022 Coarse-grained gyrokinetics for the critical ion temperature gradient in stellarators. Phys. Rev. Res. 4, L032028.CrossRefGoogle Scholar
Rogers, B.N., Dorland, W. & Kotschenreuther, M. 2000 Generation and stability of zonal flows in ion-temperature-gradient mode turbulence. Phys. Rev. Lett. 85, 5336.CrossRefGoogle ScholarPubMed
Schekochihin, A.A. 2022 MHD turbulence: a biased review. J. Plasma Phys. 88, 155880501.CrossRefGoogle Scholar
Schekochihin, A.A., Cowley, S.C., Dorland, W., Hammett, G.W., Howes, G.G., Quataert, E. & Tatsuno, T. 2009 Astrophysical gyrokinetics: kinetic and fluid turbulent cascades in magnetized weakly collisional plasmas. Astrophys. J. Suppl. 182, 310.CrossRefGoogle Scholar
Schekochihin, A.A., Parker, J.T., Highcock, E.G., Dellar, P.J., Dorland, W. & Hammett, G.W. 2016 Phase mixing versus nonlinear advection in drift-kinetic plasma turbulence. J. Plasma Phys. 82, 905820212.CrossRefGoogle Scholar
Shimomura, Y., Murakami, Y., Polevoi, A.R., Barabaschi, P., Mukhovatov, V. & Shimada, M. 2001 ITER: opportunity of burning plasma studies. Plasma Phys. Control. Fusion 43, 385.CrossRefGoogle Scholar
Sips, A.C.C. 2005 Advanced scenarios for ITER operation. Plasma Phys. Control. Fusion 47, A19.CrossRefGoogle Scholar
Snyder, P.B., Hammett, G.W. & Dorland, W. 1997 Landau fluid models of collisionless magnetohydrodynamics. Phys. Plasmas 4, 3974.CrossRefGoogle Scholar
Spitzer, L. & Härm, R. 1953 Transport phenomena in a completely ionized gas. Phys. Rev. 89, 977.CrossRefGoogle Scholar
Terry, P.W. & Horton, W. 1983 Drift wave turbulence in a low-order $k$ space. Phys. Fluids 26, 106.CrossRefGoogle Scholar
Tirkas, S., Chen, H., Merlo, G., Jenko, F. & Parker, P. 2023 Zonal flow excitation in electron-scale tokamak turbulence. Nucl. Fusion 63, 026015.CrossRefGoogle Scholar
Told, D., Jenko, F., TenBarge, J.M., Howes, G.G. & Hammett, G.W. 2015 Multiscale nature of the dissipation range in gyrokinetic simulations of Alfvénic turbulence. Phys. Rev. Lett. 115, 025003.CrossRefGoogle ScholarPubMed
Waltz, R.E. 1988 Three-dimensional global numerical simulation of ion temperature gradient mode turbulence. Phys. Fluids 31, 1962.CrossRefGoogle Scholar
Zhu, H., Zhou, Y. & Dodin, I.Y. 2020 Theory of the tertiary instability and the dimits shift within a scalar model. J. Plasma Phys. 86, 905860405.CrossRefGoogle Scholar
Figure 0

Figure 1. The growth rate of the collisional sETG instability: these are solutions to (3.13) for $\tau = Z = 1$, normalised to $\omega _{*e}$. Panel (a) is a contour plot of the positive growth rates ($\mathrm {Im} \omega >0)$ in the $(k_y, k_\parallel )$ plane; panel (b) shows cuts of the growth rate at constant $k_y \rho _\perp$, plotted as a function of $k_\parallel L_T/\sqrt {\sigma }$. The normalisations $\rho _\perp$ and $\sigma$ are defined in (3.11). The stability boundary (3.18) is indicated by grey dashed line in (a).

Figure 1

Table 1. The parameters used in the ‘baseline’ and ‘higher-resolution’ simulations. Both simulations had $\tau = Z = 1$.

Figure 2

Figure 2. Time traces of the instantaneous heat flux from simulations in which $L_\parallel \sqrt {\sigma }/L_T$ was varied from 15 to 55, normalised to $(\rho _\perp /\rho _e) Q_{\text {gB}e}$.

Figure 3

Figure 3. The scaling of the turbulent heat flux with $L_\parallel /L_T$, normalised to $(\rho _\perp /\rho _e) Q_{\text {gB}e}$ and plotted against logarithmic axes. The points are the simulation data, while the theoretical prediction [see (2.6)] is shown by the dashed black line. A logarithmic fit to the data gives the slope of $2.02$.

Figure 4

Figure 4. Turbulent heat flux in the higher-resolution simulation (see table 1), normalised to $(\rho _\perp /\rho _e) Q_{\text {gB}e}$. The upper and lower show, respectively, the instantaneous and (rolling) time-averaged heat fluxes in solid black. The dashed horizontal line in the lower panel is the average value – as calculated over the entire time interval – while the transparent grey region around this value shows the error bar associated with the mean, calculated by means of a moving window average. The time-averaged heat flux converges to within the final error bar by $(\rho _e/\rho _\perp )^2\nu _{ei} t/2 \sim 2000$.

Figure 5

Figure 5. (a) The 1D perpendicular spectra of the energy injection (5.6) (solid red), parallel dissipation (5.7) (dashed blue) and perpendicular dissipation (5.8) (dotted blue), normalised to $(\rho _e/L_T)^2 \nu _{ei}/2$. The location of the outer scale is shown by the black dot. The rate of parallel dissipation is significant at the largest scales, while perpendicular dissipation takes over at the smallest scales. (b) The cumulative perpendicular wavenumber integrals of the quantities plotted in (a), as well as the nonlinear energy flux (5.9) (solid black line). The latter is approximately constant in the inertial range, displaying only an order-unity variation, due to the finite simulation domain.

Figure 6

Figure 6. The 1D (a) perpendicular (5.11) and (b) parallel (5.18) spectra, normalised to their value at the outer scale. The spectra of the electrostatic potential are plotted in blue, those of the temperature perturbations are in red. The predicted inertial-range scalings (5.16) and (5.19) are shown by the dashed black lines. The location of the outer scale (see § 5.3) is indicated by the black dot. In (a), this is calculated from the maximum of (5.6), while in (b), it is calculated from the maximum of the 1D parallel spectrum of the energy injection, defined analogously to (5.6).

Figure 7

Figure 7. (a) The scaling of the perpendicular outer scale $k_\perp ^o$ (defined as the peak wavenumber of the energy injection (5.6)) with $L_\parallel /L_T$. (b) The scaling of the amplitude of the electrostatic potential $\bar {\varphi }^o$ [defined as the amplitude of $\varphi$ at $k_\perp = k_\perp ^o$, via (5.11)] with the perpendicular outer scale. The black points are the simulation data, while the theoretical predictions (5.23) are shown by the black dashed lines. A logarithmic fit to the data gives the slopes of $-$1.99 and $-$0.95 in (a) and (b), respectively.

Figure 8

Figure 8. A contour plot of the logarithm of the 2D spectrum (5.27) of the temperature perturbations in the $(k_\perp,k_\parallel )$ plane, normalised to its value at the outer scale. The line of critical balance is shown as the dashed black line, while the outer scale is shown by the black dot. The horizontal dotted line shows the upper bound on the parallel-wavenumber cuts plotted in the right panels of figure 9. Similarly, the vertical dotted lines show the lower and upper bounds on the perpendicular-wavenumber cuts plotted in the right panels of figure 10.

Figure 9

Figure 9. Cuts of the 2D spectra of (a) the electrostatic potential and (b) the temperature perturbations at constant $k_\parallel$, normalised to $(\rho _\perp /L_T)^2$. The colours indicate the value of $k_\parallel L_T/\sqrt {\sigma }$ for a given cut. The left panels show the entire spectrum plotted as a function of $k_\perp \rho _\perp$. The right panels show selected cuts for $k_\parallel L_T$ within the inertial range, with $k_\perp$ rescaled according to the critical-balance relation (5.17). The black dashed lines show the theoretical scalings (5.38) and (5.35) in (a), and (b), respectively. The spectra show reasonable agreement with theory at both small and large perpendicular scales, despite the effects of hyperviscosity being present at the smallest scales.

Figure 10

Figure 10. Cuts of the 2D spectra of (a) the electrostatic potential and (b) the temperature perturbations at constant $k_\perp$, normalised to $(\rho _\perp /L_T)^2$. The colours indicate the value of $k_\perp \rho _\perp$ for a given cut. The left panels show the entire spectrum plotted as a function of $k_\parallel L_T$. The right panels show selected cuts of the spectrum for $k_\perp \rho _\perp$ within the inertial range, with $k_\parallel$ rescaled according to the critical-balance relation (5.17). The black dashed lines show the theoretical scalings (5.38) and (5.35) in (a) and (b), respectively. There is very good agreement with theory, especially at $k_\parallel \lesssim k_\perp ^{2/3}$, where the scalings extend well beyond the inertial range to higher $k_\perp \rho _\perp$, as can be seen from the left panels – this is because the causality argument is not sensitive to the precise details of the decorrelation physics.

Figure 11

Figure 11. Real-space snapshots of (a) the electrostatic potential and (b) the temperature perturbations from the higher-resolution simulation at $(\rho _e/\rho _\perp )^2\nu _{ei} t/2 = 3000$ (see table 1). The coordinate axes are as shown, while the red and blue colours correspond to regions of positive and negative fluctuation amplitudes. The turbulence does not appear to be isotropic on the large scales that are visible in these plots (streamers are manifest), but turns out to be isotropic in the inertial range (see figure 12).

Figure 12

Figure 12. Contour plots of the 2D spectra of the temperature perturbations, normalised to $(\rho _\perp /L_T)^2$: (a) in Cartesian coordinates, with the radial and poloidal wavenumbers plotted on the horizontal and vertical axes, respectively; contours of constant $E^T(k_x, k_y)$ (5.39) (black dashed lines) are approximately circular away from the origin, where injection is localised and the presence of streamers is manifested by the spectral power being shifted towards $k_y >k_x$; (b) in polar coordinates, with $\theta = \tan ^{-1}(k_y/k_x)$ and $k_\perp \rho _\perp$ plotted on the horizontal and vertical axes, respectively; contours of constant $E^T(k_\perp, \theta )$ (5.40) (black dashed lines) are approximately horizontal far away from $k_\perp \rho _\perp \lesssim 1$, where injection is localised.

Figure 13

Figure 13. Time traces of the instantaneous heat flux from simulations with finite $L_B/L_T$, with the limit of $L_B/L_T \rightarrow \infty$ shown for comparison. All parameters are the same as the baseline simulation (see table 1), and the heat flux is normalised to $(\rho _\perp /\rho _e) Q_{\text {gB}e}$. The heat flux grows without bound in all simulations with (finite) $L_B/L_T$ above the linear critical gradient (C3) (${\approx }2.72$ for $\bar {\tau }= 1$), with the rate of divergence decreasing as $L_B/L_T$ is increased.

Figure 14

Figure 14. Real-space snapshots of the (a) electrostatic potential and (b) temperature perturbations from the simulation with $L_B/L_T = 1000$ from figure 13, taken at $(\rho _e/\rho _\perp )^2\nu _{ei} t/2 = 200$. The coordinate axes are as shown, while the red and blue colours correspond to regions of positive and negative fluctuation amplitude. At these early times, the turbulence appears similar to that of saturated sETG turbulence for $L_B/L_T \rightarrow \infty$ (cf. figure 11), despite the eventual lack of saturation (see figure 15).

Figure 15

Figure 15. The same as figure 14, except taken at $(\rho _e/\rho _\perp )^2\nu _{ei} t/2 = 1000$. The unbounded growth of the heat flux is associated with the formation of large-scale, approximately 2D streamer structures that appear to be immune to all types of nonlinear shearing.