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 ¹ . 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.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 ² , 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.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 ³ – 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 ⁴]. 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 ⁵ . 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 ⁶ (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 ⁷:

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