Scale invariance and critical balance in electrostatic drift-kinetic turbulence

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).

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. (2009) 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 1988; Cowley, Kulsrud & Sudan 1991; Kotschenreuther et al. 1995a) and electron-temperature gradient (ETG) ones (see, e.g. Liu 1971; Lee et al. 1987; Dorland et al. 2000; Jenko et al. 2000). 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. 2017; Hatch et al. 2019; Parisi et al. 2020; Guttenfelder et al. 2021, 2022; Chapman-Oplopoiou et al. 2022; Parisi et al. 2022; Field et al. 2023), 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 (1941) 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 (2011) 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 (1965) (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. 2013) – 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 (1941)-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. 1991; Newton, Cowley & Loureiro 2010; Ivanov et al. 2020; Ivanov, Schekochihin & Dorland 2022). 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 1990; Hammett, Dorland & Perkins 1992; Dorland & Hammett 1993; Hammett et al. 1993; Beer & Hammett 1996; Snyder, Hammett & Dorland 1997). 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. 2006). 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. 2021, Hatch et al. 2022, Chapman-Oplopoiou et al. 2022 and Field et al. 2023, 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 1965; Helander & Sigmar 2005), 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 (1965).

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. 2022), 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. 2022; Adkins 2023) – both of these bring in a special perpendicular scale that would break the drift-kinetic scale invariance¹ . 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. 2022).

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. 2022), a cousin of the collisionless sETG (Lee et al. 1987).

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 1971). 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 1995b; Dorland et al. 2000).

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)² , 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. 2013). 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 achieved³ – 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. (2011) and Adkins et al. (2022). The original rationale for this conjecture comes from the causality argument proposed in the context of MHD turbulence (Goldreich & Sridhar 1995; Boldyrev 2005; Schekochihin 2022): 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 2004, Nazarenko & Schekochihin 2011 and Adkins et al. 2022). 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 (1946) damping rate⁴]. 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. 2022).

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. (2011), 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 direction⁵ . 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. 2020, 2022) or in stellarators (Roberg-Clark, Plunk & Xanthopoulos 2022), 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 (2022) (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 result⁶ (see, e.g. appendix A of Schekochihin et al. 2016). 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 perturbations⁷:

(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. (2011) observed an overlap of the spatial scales of energy injection and dissipation in electrostatic, ion-scale toroidal gyrokinetic simulations, as did Told et al. (2015) 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. 2022), 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. 1991). 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. (2011) 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. 2000; Jenko et al. 2000; Colyer et al. 2017), 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. 2006; Candy et al. 2007; Roach et al. 2009; Guttenfelder & Candy 2011).

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. 1993; Abel & Cowley 2013):

(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 2000; Ivanov et al. 2020, 2022). 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 1978; Hasegawa & Wakatani 1983; Terry & Horton 1983; Diamond et al. 2005; Ivanov et al. 2020; Zhu, Zhou & Dodin 2020). 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. (2017) and Tirkas et al. (2023).

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. (2011), 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. 2013). Apart from some interesting exceptions (Parisi et al. 2020, 2022), 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. 2022, 2023; Parisi et al. 2022). 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. 2006; Candy et al. 2007; Roach et al. 2009; Guttenfelder & Candy 2011). 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 1988; Adkins et al. 2022), 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. 1993). 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. (2022)), 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. (2022) 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. 2001; Sips 2005; Patel et al. 2021) 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.