Scale invariance and critical balance in electrostatic drift-kinetic turbulence

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


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 et al. 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. 2021Guttenfelder et al. , 2022Field 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 et al. (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 (k ⊥ , k ∥ ) 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.

Electrostatic drift-kinetic scale invariance
For systems adequately described by electrostatic drift kinetics, the heat flux through some volume V is given by where ∇x is the (radial) direction of the equilibrium gradients, n 0s and T 0s are the equilibrium density and temperature, respectively, of species s, δT s is the corresponding temperature perturbation [see (A 17)], and is the E × B drift velocity due to the perturbed electrostatic potential ϕ, B 0 and b 0 being the magnitude and direction of the equilibrium magnetic field, respectively. In what follows, L ns and L Ts 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 ths = 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 oneparameter transformation. Under this transformation, the perturbed temperature and electrostatic potential transform as, for any λ, δT s = λ 2 δT s (x/λ 2 , y/λ 2 , z/λ 2/α , t/λ 2 ),φ = λ 2 ϕ(x/λ 2 , y/λ 2 , z/λ 2/α , t/λ 2 ). (2.3) Here, x, y and z are the radial, binormal and parallel coordinates, respectively, the tildes indicate the transformed fields, and α = 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 ω is comparable to rate of thermal conduction, but much smaller than ν ss ′ , the characteristic collision frequency between species s and s ′ , viz., ω ∼ (k ∥ v ths ) 2 /ν ss ′ ≪ ν ss ′ , as in Braginskii (1965) (where k ∥ 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 ρ s -the Larmor radius of species s, manifest in the gyroaverages and the resultant Bessel functions appearing in gyrokinetics (see, e.g., )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 asQ s = λ 2 Q s .
(2.4) Now suppose that our original solutions for δT s and ϕ were periodic in x, y and z with domain sizes L x , L y , and L ∥ , respectively. Then, the transformed solutions δT s andφ are still periodic in x, y and z, except with domain sizes λ 2 L x , λ 2 L y , and λ 2/α L ∥ , implying thatQ s (λ 2 L x , λ 2 L y , λ 2/α L ∥ , t/λ 2 ) = λ 2 Q s (L x , L y , L ∥ , t). (2.5) 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 −1 ∥ ∼ L ∥ ∼ L Ts ), but remain microscopic in the perpendicular direction (k −1 ⊥ ∼ ρ s ). This means that, as the perpendicular domain size L ⊥ (ordered as L ⊥ ∼ L x ∼ L y ∼ ρ 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 ⊥ /ρ s → ∞, 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 λ can be chosen arbitrarily, (2.5) directly implies that where once again α = 1, 2 in the collisionless and collisional limits, respectively. Physically, L ∥ can be thought of either as a measure of a quantity analogous to the connection length 2πqR 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, where Q gBs = n 0s T 0s v ths (ρ s /L Ts ) 2 is the "gyro-Bohm" heat flux, G is an unknown function, ν * s = (L T /v ths ) s ′ ν ss ′ is the normalised collisionality, and ". . . " 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 Ts . 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 Ts follows directly from its dependence on L ∥ . 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 ⊥ (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.

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 et al. 2010;Ivanov et al. 2020Ivanov et al. , 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 et al. 1992Hammett et al. , 1993Beer & Hammett 1996;Snyder et al. 1997). Here we consider an electron-scale, collisional (ν * e ≫ 1) fluid model of electrostatic turbulence driven by the electrontemperature 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 b 0 =ẑ 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 ϕ by E = −∇ϕ, and has no mean part. The equilibrium profile of the electron temperature T 0e varies radially, with the scale length 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 −1 T ; 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;Field et al. 2023, where a strong dependence of the heat flux on the equilibrium density gradient was identified).

Moment equations
With this local equilibrium, we derive, in appendix B, evolution equations for the density (δn e ), parallel velocity (u ∥e ), and temperature (δT e ) perturbations of the electrons: d dt 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 E × B motion (2.2) of the electrons: and their compression or rarefaction due to the perturbed electron flow u ∥e b 0 parallel to the equilibrium magnetic-field direction. This flow velocity is determined (instantaneously) from a balance between the electronion frictional force -proportional to the electron-ion collision frequency ν ei [see (B 4)] 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 (A 4)]. 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 (B 38)]: e.g., for Z = 1, c 1 ≈ 1.94, c 2 ≈ 1.39 and c 3 ≈ 3.16, in agreement with Braginskii (1965).
The temperature perturbation in (3.4) is advected by the local E × B flow (2.2), again according to (3.5), and is locally increased (or decreased) by compressional heating (or rarefaction cooling) due to u ∥e , as well as by the perturbed parallel collisional heat flux 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 E × B flow) responsible for extracting free energy from the equilibrium temperature gradient L −1 T , defined in (3.1). Finally, the electron-density perturbation is related to the non-dimensionalised potential φ via quasineutrality: where τ = 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 ρ 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: Note that, due to the Boltzmann density response (3.7), the advection term in (3.8) has vanished, leaving a purely linear relationship between φ and δT e . The only nonlinearity left in the system is thus the advection of δT e by the E × 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.

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 α = 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), 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 1 . In other words, (3.8) and (3.9) describe physics on scales where σ is, formally, some arbitrary constant satisfying β e ≪ σ ≪ 1. (3.12) 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 σ.

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: where we have introduced the characteristic parallel and perpendicular frequencies (3.14) 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 ω because the parallel velocity u ∥e is determined instantaneously in terms of the other fields, by (3.3), unlike in collisionless ETG theory ).
If we consider the limit of long parallel wavelengths, viz., 1 Implicit in (3.10) is the assumption that the electron inertial scale de = ρe/ √ βe is smaller than the ion Larmor radius ρi; this is only the case if βe ≳ me/mi, which is well-satisfied in most systems of interest. Should de lie on the large-scale side of ρi (i.e., βe ≲ me/mi), the lower bound in (3.12) must be replaced with me/mi. then the balance of the first and last terms in (3.13) gives us We recognise this as the collisional slab ETG (sETG) instability ), 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: In this limit, the instability works as follows. Suppose that a small perturbation of the electron temperature is created with k y ̸ = 0 and k ∥ ̸ = 0, bringing the plasma from regions with higher T 0e to those with lower T 0e (δT e > 0), and vice-versa (δ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 E × 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 ω to be purely real, the real and imaginary parts of (3.13) must vanish individually. The resultant equations can be straightforwardly combined to yield This is a curve k y ∝ k 2 ∥ in wavenumber space, plotted as the grey dashed line in figure 1(a). Above this line, corresponding to the limit ω ∥ ≫ ω * e , all modes are purely damped due to rapid thermal conduction, as in figure 1 At any given k y , the maximum growth rate of the collisional sETG is, therefore, reached when (3.19) 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 ω ∥ , one finds γ max = C(τ, Z)ω * e , where C(τ, Z) is a constant formally of order unity, e.g., C(1, 1) ≈ 0.094 [cf. the maximum values in figure 1(b)]. The increase of the maximum growth rate of the sETG instability with the perpendicular wavenumber, ω * e ∝ 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.

Numerical setup
In what follows, the system (3.8)-(3.9) is solved numerically in a triply periodic box of size L x × L y × L ∥ using a pseudo-spectral algorithm. Numerical integration is done in Fourier space (N x , N y and N ∥ are the number of Fourier harmonics in the respective directions) with the nonlinear term calculated in real space using the 2/3 rule for dealiasing (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 GS2 (Kotschenreuther et al. 1995b;Dorland et al. 2000).
Perpendicular hyperviscosity is introduced in order to provide an ultraviolet (largewavenumber) cutoff for the instabilities, achieved by the replacement of the time derivative on the left-hand sides of (3.8) and (3.9) with where ν ⊥ is the "hypercollision" frequency and N ν ⩾ 2. With this change, our equations now depend only on the following dimensionless parameters: the perpendicular and parallel box sizes L x /ρ ⊥ , L y /ρ ⊥ and L ∥ √ σ/L T , the hyper-collision frequency 2σ(ρ ⊥ /ρ e ) 2 ν ⊥ /ν ei , and the power of the hyperviscous diffusion operator N ν . Convergence Baseline  40  40  20  191 191 31  0.00050  2  Higher-resolution  40  40  20  383 383 63  0.00015  2   Table 1: The parameters used in the "baseline" and "higher-resolution" simulations. Both simulations had τ = Z = 1. scans in N x , N y , and the perpendicular box size L x = L y = L ⊥ 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 ⊥ was large enough so that it did not significantly affect the simulation results, as was required for the arguments of §2.
We have also found that our results do not depend on the specific details of the hyperviscosity, viz., on the values of ν ⊥ and N ν . 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 driftkinetic 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 ∥ , 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.

Scan in L ∥ /L T
In order to test the dependence of the turbulent heat flux on L ∥ predicted by (2.6), we performed a series of simulations in which L ∥ √ σ/L T was varied between 15 and 55 at fixed parallel resolution (viz., fixed ratio of L ∥ √ σ/L T to N ∥ ), 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) 2 , the other assumption, also confirmed numerically, being that this state is independent of L x and L y .
In figure 3, we plot the time average of the turbulent heat flux -as defined in (2.1) for s = e, and normalised to (ρ ⊥ /ρ e )Q gBe , where Q gBe = n 0e T 0e v the (ρ e /L T ) 2 is the 2 Refining our consideration beyond this assumption of stationarity, we observe that the characteristic timescale of the fluctuations of the instantaneous heat flux increases with the parallel system size -this is manifest in figure 2, where the simulations with larger L ∥ √ σ/LT exhibit higher-amplitude, longer-timescale fluctuations. The origin of this trend can be understood as follows. Relaxing the assumption of stationarity, instead of (2.6), we have, from (2.5),Qs(λ 2/α L ∥ , t/λ 2 ) = λ 2 Qs(L ∥ , t). If Qs exhibits fluctuations on some characteristic timescale τ , then, if we assume that that both solutions must be periodic with the same period, the corresponding timescale for the transformed heat flux will beτ = λ 2 τ . Given that the parallel system sizes for both solutions are related byL ∥ = λ 2/α L ∥ , it follows that τ ∝ L α ∥ . This dependence was confirmed numerically for the set of simulations shown in figure 2. (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.

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 (ρ e /ρ ⊥ ) 2 ν ei t/2σ; this was sufficient to ensure that the heat flux had converged to a well-defined average value (see figure 4).

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., ). In the system of equations that we are considering, the (normalised) free energy reduces to the form . (5.1) 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 is the energy-injection rate from the equilibrium temperature gradient, and are the dissipation rates due to (parallel) thermal conduction and (perpendicular) hyperviscosity, respectively. The corresponding 1D perpendicular wavenumber spectrum of the energy injection is while those of the parallel and perpendicular dissipation are (5.8) 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 (ρ e /ρ ⊥ ) 2 ν ei t/2σ ∼ 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 ⊥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 3 -appears to be independent of hyperviscosity, being 3 In standard turbulence literature, the outer scale is often defined to be the integral scale of the 1D perpendicular energy spectrum (5.11), viz., . However, given that, physically, we are interested in the outer scale as the scale at which the free energy is predominantly injected, the choice to maximise (5.6) seems to be better motivated physically. localised on much larger scales, where D ⊥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 hypervisocisty. The fact that the energy injection is both independent of hypervisocisty 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. Considering scales that are larger than the injection scale, it is clear from figure 5 that the parallel dissipation D ∥ 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, ω ∥ ∼ ω * 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: 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.

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 ε that survives the parallel dissipation at larger scales, viz., ε − D ∥ , which, for brevity, we shall call ε in the scaling arguments that follow.
The only nonlinearity in our equations is the advection of the temperature fluctuations by the fluctuating E × B flows in (3.9). Therefore, we take the nonlinear cascade time to be the nonlinear E × B advection time: (5.10) Here and in what follows,φ refers to the characteristic amplitude of the electrostatic potential at the scale k −1 ⊥ . Formally,φ can be defined bȳ is the 1D perpendicular spectrum of φ, φ k is the spatial Fourier transform of the potential, and the angle brackets denote an ensemble average. The corresponding quantities for the temperature perturbations, δT e , E T ⊥ (k ⊥ ), and δT ek , 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: We are using the δT e /T 0e part of the free energy (5.1) because, as we noted earlier, in (3.8)-(3.9), δT e /T 0e is the only field that is advected nonlinearly whereas φ 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φ 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 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 −1 nl at each perpendicular scale k −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;): 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;Adkins et al. 2022). In our system, the parallel dynamics are dissipative, with the relevant timescale being set by the parallel conduction rate ω ∥ . 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 ω ∥ ≪ t −1 nl 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 ω ∥ ≫ t −1 nl , they are wiped out by heat conduction. As a result, the "dissipation ridge" (the line of critical balance) ω ∥ ∼ t −1 nl 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 ω ∥ is instead the Landau (1946) damping rate 4 ]. In §5.1, we saw that the actual amount of parallel dissipation that happens in the inertial range is small -free energy chooses to stay just shy of the dissipation region (ω ∥ > t −1 nl ) and instead cascade, at an approximately constant rate, along the dissipation ridge (ω ∥ ∼ t −1 nl ). Combining (5.12), (5.13) and (5.14), we find the following scaling of the amplitudes in the inertial rangeφ Then, recalling (5.11), the 1D perpendicular energy spectra in the inertial range are: (5.16) 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: If we define the 1D parallel spectrum 18) and the corresponding temperature spectrum E T ∥ (k ∥ ) analogously, (5.17) and (5.15) imply the following inertial-range scaling of amplitudes with parallel wavenumbers: These simple scaling arguments are vindicated by simulation data. The 1D spectra (5.11) and (5.18) for both φ and δ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. 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.

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 20) and thus grow at a rate ∼ ω * e ∝ 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), The outer scale is then the scale at which these two rates are comparable: balancing (5.21) and (5.20), we get where the superscript 'o' refers to outer-scale quantities. Thus, we have two relationships between k o ⊥ ,φ o and k 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 o y . 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 ∥ . 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): where ρ ⊥ and σ are defined in (3.11), and the magnitude of σ 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 ∥ /L T that was presented in §4.2.
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) Recalling (2.1), the combination of (5.24) and (5.22) yields the following expression for the turbulent heat flux: where once again Q gBe = n 0e T 0e v the (ρ e /L T ) 2 is the "gyro-Bohm" flux. Unsurprisingly, this reproduces the scaling with L ∥ given by (2.6) for α = 2 [cf., also, (2.7) for G = L T /λ ei ]. Note that, apart from the inevitable dimensional factors, the L ∥ 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 α = 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 o ∥ L ∥ ∼ 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 ∥ ∼ 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 ∥ ∼ 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 5 . A practical implication of this conclusion for more realistic systems appears to be that any long-scale parallel inhomogeneity should be sufficient to set k 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 or in stellarators (Roberg-Clark et al. 2022), where magnetic fields have parallel structure on scales shorter than the connection length.

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: Unlike in §5.2, we can no longer assume that E φ 2D ∼ E T 2D ; 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 ω ∼ ω ∥ . 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): 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 ∥ /k o ∥ → k ∥ , k ⊥ /k o ⊥ → k ⊥ 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 If a > 1 and d > −1, k ∥ ∼ k 2/3 ⊥ will be the energy-containing parallel wavenumber at a given k ⊥ . The 1D perpendicular spectrum is, therefore, This must match the scaling (5.16) of the 1D perpendicular spectrum derived from the constant-flux conjecture, implying that c = 2 3 (1 + d) + 7 3 .
(5.31) Two further constraints follow from imposing boundary conditions as k ∥ or k ⊥ → 0 at constant k ⊥ or k ∥ , respectively. The scaling of the spectrum as k ⊥ → 0 (in the region k ∥ ≫ k 2/3 ⊥ ) can be determined purely kinematically: the low-k ⊥ asymptotic behaviour of a homogenous 2D-isotropic field must be k 3 ⊥ , implying that b = 3.
(5.32) This is a fairly standard result 6 (see, e.g., appendix A of Schekochihin et al. 2016). 6 Though not one that can be taken for granted. For example, Hosking & Schekochihin (2022) (see also appendix C of Schekochihin 2022) showed that a k 1 ⊥ scaling could emerge instead through a balance between turbulent diffusion at large scales and the nonlinear 'source' that would otherwise give rise to the k 3 ⊥ scaling. Let us estimate the rate of turbulent diffusion in our system. The dominant contribution to the turbulent-diffusion coefficient D will be from k ⊥ ∼ k 3/2 ∥ , which, at any given k ∥ , plays the role of the energy-containing scale. Then D ∼ v 2 E t nl ∼ω ∥ /k 2 ⊥ , where have used the critical-balance condition (5.14), and the tildes denote quantities evaluated atk ⊥ ∼ k 3/2 ∥ ≫ k ⊥ , where k ⊥ is the wavenumber at which turbulent diffusion is acting. The Finally, the scaling as k ∥ → 0 (in the region k ∥ ≪ k 2/3 ⊥ ) follows from causality. Indeed, in §5.2, we argued that fluctuations become decorrelated for ω ∥ ≲ t −1 nl because they cannot communicate across parallel distances ∼ k −1 ∥ , but such k ∥ are also too small for the fluctuations to be erased by thermal conduction. Therefore, the parallel spectrum at k ∥ ≪ k 2/3 ⊥ must be the spectrum of a 1D white noise: (5.33) Combining (5.32) and (5.33) with (5.29) and (5.31), we find a = 9, c = 3. (5.34) This gives us the following scalings for the 2D spectrum of the temperature perturbations 7 : (5.35) 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 ∥ ∼ k 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 ⊥ → 0 (in the region k ∥ ≳ k 2/3 ⊥ ) will once again be k 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 ∥ ≲ k 2/3 ⊥ now no longer holds, because φ is not directly decorrelated by the nonlinearity. Instead, it inherits its scaling from δT e /T 0e via the balance (5.13), viz., where we have used ω ∥ ∝ k 2 ∥ and the second expression in (5.35). Now, in the region k ∥ ≲ k 2/3 ⊥ , we expect thermal conductivity in the temperature equation to be subdominant to the nonlinear rate, and so estimating ω ∼ t −1 nl ∝ k 4/3 ⊥ in (5.36), we find that d = 4, c = 17 3 .
(5.37) rate of turbulent diffusion at this wavenumber will thus be k 2 ⊥ D ∼ω ∥ (k ⊥ /k ⊥ ) 2 ≪ ω ∥ . Turbulent diffusion is, therefore, negligible, and so we are justified in adopting the k 3 ⊥ scaling. Note that the survival of the k 3 ⊥ scaling is a noteworthy feature of our system, where the dynamics at k ∥ ≫ k 2/3 ⊥ are dominated by parallel dissipation due to thermal conductivity and do not produce significant turbulent diffusion -unlike waves, which, e.g., in RMHD, do . 7 Schekochihin et al. (2016) obtained, by a similar method, an analogous result for long-wavelength electrostatic ITG turbulence (which, in this approach, is no different for ETG). Specifically, they found that a = 5 and c = 11/3 -this was a consequence of the fact that they considered collisionless turbulence, for which the critical-balance condition is This gives us the following scalings of the 2D spectrum of the potential fluctuations: (5.38) 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 ∥ and k ⊥ 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 ∥ ≲ k 2/3 ⊥ , in particular, are another confirmation of the causal nature of the critical balance.
The extraordinarily steep parallel-wavenumber scaling of the 2D spectra (5.35) and (5.38) in the region k ∥ ≳ k 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 ω ∥ ≳ t −1 nl , and thus to a region of dominant thermal conduction that attempts to erase parallel structure created by the turbulence. The k ∥ 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 ∥ ∼ k 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 Figure 9: Cuts of the 2D spectra of (a) the electrostatic potential and (b) the temperature perturbations at constant k ∥ , normalised to (ρ ⊥ /LT ) 2 . The colours indicate the value of k ∥ LT / √ σ for a given cut. The left panels show the entire spectrum plotted as a function of k ⊥ ρ ⊥ . The right panels show selected cuts for k ∥ LT within the inertial range, with k ⊥ rescaled according to the critical-balance relation (5.17). The black dashed lines show the theoretical scalings (5.38) and (5.35) in panels (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.
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, ionscale 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- Figure 10: Cuts of the 2D spectra of (a) the electrostatic potential and (b) the temperature perturbations at constant k ⊥ , normalised to (ρ ⊥ /LT ) 2 . The colours indicate the value of k ⊥ ρ ⊥ for a given cut. The left panels show the entire spectrum plotted as a function of k ∥ LT . The right panels show selected cuts of the spectrum for k ⊥ ρ ⊥ within the inertial range, with k ∥ rescaled according to the critical-balance relation (5.17). The black dashed lines show the theoretical scalings (5.38) and (5.35) in panels (a) and (b), respectively. There is very good agreement with theory, especially at k ∥ ≲ k 2/3 ⊥ , where the scalings extend well beyond the inertial range to higher k ⊥ ρ ⊥ , 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. driven turbulence. The growth rate of the collisionless sETG instability is limited by the parallel streaming rate ω ∥ ∼ k ∥ v the (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 Figure 11: Real-space snapshots of (a) the electrostatic potential and (b) the temperature perturbations from the higher-resolution simulation at (ρe/ρ ⊥ ) 2 νeit/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). 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 ω ∥ ∼ ω ∼ t −1 nl 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.

Perpendicular isotropy
Throughout sections 5.2 to 5.4, our theoretical deductions were carried out under the assumption that k x ∼ k y ∼ k ⊥ . 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 ≪ 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).
To assess how isotropic the saturated state is, in particular in the inertial range, we plot the two-dimensional spectra Here and in what follows, θ = 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. 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;: where φ ′ 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 Ivanov et al. 2020Ivanov et al. , 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 φ ′ through the electrostatic E × 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 twodimensional 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 et al. 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).

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 ∥ , viz., Q s ∝ L ∥ (or ∝ L 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 ∥ ∼ 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 gradientdriven 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 ⊥ ∼ ρ −1 s will now appear explicitly in the equations through the Bessel functions (see, e.g., . Apart from some interesting exceptions (Parisi et al. 2020, 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 constantflux 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 driftkinetic 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 ∥ . 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 ∥ . 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 ρ i ∼ 1 but narrow radial structures near mode-rational surfaces on the scale k x ρ e ∼ 1 (Hardman et al. 2022(Hardman et al. , 2023. 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 et al. 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 ETGscale turbulence in the presence of finite magnetic-field gradients . 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 inertialrange 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 ∥ , rather than the flux-freezing scale. For example, Chapman-Oplopoiou et al. (2022) performed nonlinear, electromagnetic simulations of JET-ILW pedestals for k ⊥ ρ i ≳ 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.
(EP/T012250/1). The views and opinions expressed herein do not necessarily reflect those of the European Commission. The work of AAS was also supported in part by a grant from STFC (ST/W000903/1) and by the Simons Foundation via a Simons Investigator award.

Declaration of interests
The authors report no conflicts of interest.

Appendix A. Derivation of scale invariance
In this appendix, we demonstrate explicitly that electrostatic drift kinetics remains invariant under the transformation (2.3), which leads to the heat-flux scaling (2.6).
We take as our starting point the electrostatic drift-kinetic system, in which the perturbed distribution function for species s consists of a Boltzmann and gyrokinetic parts and h s evolves according to Here, and in what follows, f 0s is the Maxwellian equilibrium distribution of species s with density n 0s and temperature T 0s , and ϕ is the perturbed electrostatic potential. The magnetic-drift velocity arising from the inhomogeneities in the equilibrium magnetic field is given by where b 0 = B 0 /B 0 is the direction of the equilibrium magnetic field, B 0 = |B 0 | is its magnitude, and Ω s = q s B 0 /m s c, q s and m s are the Larmor frequency, charge and mass of species s, respectively. The collision term on the right-hand side of (A 2) is the linearised Landau collision operator In what follows, it will be useful to decompose h s into parts that are even and odd in the parallel velocity v ∥ , viz., It follows straightforwardly from (A 2) and (A 4) that h even s and h odd s satisfy, respectively, The quasineutrality condition (A 5) becomes Note that in (A 8) and (A 9), we have assumed that the (radial) gradient of the equilibrium distribution function ∇f 0s is an even function of v ∥ -this is only the case in systems without any equilibrium flows. We now wish to consider transformations of the system of equations (A 8)-(A 10) that can be made whilst preserving the size of perpendicular equilibrium gradients. It is obvious from considering, e.g., the magnetic-drift velocity (A 3) that any rescaling of the velocity variables v ∥ and v ⊥ -at fixed equilibrium magnetic-field strengthwould require a compensatory rescaling of ∇ log B 0 and |b 0 · ∇b 0 | in order to preserve the magnitude and direction of v ds . Therefore, we will henceforth restrict ourselves to transformations involving only the spatial and time coordinates. In a similar vein to Connor & Taylor (1977), we consider the following one-parameter transformation: h even s = λ ae h even s (x/λ a ⊥ , y/λ a ⊥ , z/λ a ∥ , t/λ at ), (A 11) h odd s = λ ao h odd s (x/λ a ⊥ , y/λ a ⊥ , z/λ a ∥ , t/λ at ), (A 12) ϕ = λ ae ϕ(x/λ a ⊥ , y/λ a ⊥ , z/λ a ∥ , t/λ at ), (A 13) where x, y and z are the radial, binormal and parallel (to the magnetic field) coordinates, respectively, the tildes indicate the transformed distribution functions and fields, and a i are real constants parametrising the transformation. Quasineutrality (A 10) implies that the amplitudes of the 'even' fields must be rescaled in the same way, as in (A 11) and (A 13), while the rescaling of the amplitude of h odd s remains unconstrained. The spatial and time coordinates can then be rescaled independently, with the caveat that the radial and binormal coordinates should be rescaled in the same way in order not to rule out perpendicular isotropy. The rescaling (A 11)-(A 13) is the most general one-parameter transformation of electrostatic drift kinetics that can be made while allowing (although not requiring) the spatial isotropy of structures in the perpendicular plane.
The constants a i can be fixed by demanding that the transformation leave (A 8) and (A 9) invariant. In the collisionless limit, the collision operator can be neglected and it is easy to show that a e = a o = a ⊥ = a ∥ = a t is the only choice that fulfils this condition. The collisional limit is somewhat more subtle. As we have done throughout this paper, we order ω ∼ (k ∥ v the ) 2 /ν ss ′ ≪ ν ss ′ and ωh even s ∼ k ∥ v ths h odd s . In the resultant collisional expansion, the collision operator is forced to vanish at leading order (see appendix B.3.1), and can only survive at higher order due to the presence of finite-Larmor-radius effects (see appendix B.3.3), which are neglected within the drift-kinetic approximation. At first order, one obtains, from (A 9), a balance between the parallel streaming of h even s and the collision operator acting on h odd s (see appendix B.3.2). At second order, one evolves h even s via (A 8) with the collision operator neglected (see appendix B.3.3). One can then show that a e = 2a o = a ⊥ = 2a ∥ = a t is the only choice of parameters that leaves the drift-kinetic equations invariant. Any constraints on a i inferred from (A 11)-(A 10) in this way are only valid to second order within the collisional expansion, and not to any higher orders. However, given that the solvability conditions (B 25) and (B 36) guarantee that a closed system can be obtained solely from these two orders, this is not a problematic limitation.
Thus, it follows from the above discussion that electrostatic drift kinetics is invariant under the transformatioñ where we have chosen a e = 2 without loss of generality, and α = 1, 2 in the collisionless and collisional limits, respectively. The transformation of the odd and even parts of the distribution function h s is inherited by its moments that are odd and even in v ∥ ; e.g., the temperature perturbation, being a velocity moment that is even in v ∥ , viz., will transform according to (A 14). This, when combined with the transformation (A 16) of the electrostatic potential ϕ gives exactly (2.3), which is the starting point for the deductions presented in the main text.

Appendix B. Derivation of collisional fluid model
This appendix details a self-contained derivation of the electron-fluid equations (3.2)-(3.4). An alternative route to these via a subsidiary expansion of a more general system of (electromagnetic) equations can be found in Adkins (2023) (see also appendix G of Adkins et al. 2022). In what follows, appendix B.1 describes and physically motivates our electron-scale, collisional ordering, which is then implemented to derive equations describing our ion and electron dynamics in appendices B.2 and B.3, respectively. Although the magnetic geometry adopted throughout the majority of this paper is that of a conventional slab (see §3), we shall here consider the more general case in which the equilibrium (mean) magnetic field B 0 is assumed to have the scale length and radius of curvature assumed constant across the domain. Doing so will allow us to capture the effects of the magnetic drifts on our plasma while retaining most of the simplicity associated with conventional slab gyrokinetics (Howes et al. 2006;Newton et al. 2010;Ivanov et al. 2020Ivanov et al. , 2022Adkins et al. 2022). Note that for a low-beta plasma, R = L B .

B.1. Collisional, electron-scale ordering
In our model, we would like to be able to capture, at a minimum, the physics associated with drift waves, perpendicular advection by both magnetic drifts and E × B flows, and parallel heat conduction. Therefore, we postulate an asymptotic ordering in which the frequencies ω of the perturbations in the plasma are comparable to the characteristic frequencies associated with these phenomena, viz., where are the drift and magnetic-drift frequencies, respectively, v E = cE × B/B 2 is the E × B drift velocity (c is the speed of light), κ ∼ v 2 the /ν ei is the electron thermal diffusivity, and are the electron-ion and electron-electron collision frequencies, respectively, log Λ being the Coulomb logarithm (Braginskii 1965;Helander & Sigmar 2005). The ordering of the parallel conduction rate with respect to the drift frequency gives us a constraint relating parallel and perpendicular wavenumbers: where λ ei = v the /ν ei is the electron-ion mean free path and L is some (perpendicular) equilibrium length scale, L ∼ L Ts ∼ L B ∼ R. The ordering of the parallel conduction rate with respect to the E × B drifts determines the size of perpendicular flows within our system: where ϵ = ρ e /L is the gyrokinetic small parameter (see, e.g., ), mandating small-amplitude, anisotropic perturbations. The frequency of these perturbations is small compared to the Larmor frequencies of both the electrons and ions: ω The ordering (B 6) of v E relative to the electron thermal velocity allows us to order the amplitude of the perturbed scalar potential ϕ: The density perturbations δn s are ordered anticipating a Boltzmann density response and the temperature perturbations δT s are assumed comparable to them: , (B 10) where β e = 8πn 0e T 0e /B 2 0 the electron plasma beta. The (compressive) parallel magneticfield perturbations are ordered anticipating pressure balance: The orderings (B 5)-(B 11) still allow for a choice of ordering for perpendicular wavenumbers k ⊥ with respect to the electron and ion Larmor radii. Given that we would like to obtain a set of electrostatic equations that exhibit the scale invariance discussed in §2, we consider wavenumbers for which the physical motivation is discussed in §3.2. In terms of time scales, (B 12) is equivalent to demanding that where d e = ρ e / √ β e is the electron inertial scale. We shall formalise (B 12) by demanding that k ⊥ ρ e ∼ σλ ei /L, where σ is a placeholder constant satisfying β e ≪ σ ≪ 1. In other words, k ⊥ ρ ⊥ ∼ 1, where ρ ⊥ = ρ e L/λ ei σ, an "intermediate" spatial scale [cf. (3.11)].
To summarise, (B 5)-(B 12) imply the following ordering of frequencies: length scales: and amplitudes: All relevant quantities are thus naturally ordered with respect to some combination of m e /m i , σ, λ ei /L, and the gyrokinetic small parameter ϵ = ρ e /L. The above ordering of frequencies, length scales and amplitudes with respect to ϵ is the standard gyrokinetic ordering (see, e.g., . We choose to treat the ordering in λ ei /L -the fact that this should be formally small following straightforwardly from, e.g., ν ei ≫ ω * s -as subsidiary to both the orderings in ϵ and in the mass ratio [see the first expression in (B 15)], meaning that the formal hierarchy of our expansions is with all other dimensionless parameters treated as finite.

B.2. Ion kinetics
Given that the ordering of perpendicular wavenumbers (B 15) implies that k ⊥ ρ i ≫ 1 under the expansion in the mass ratio, the ion distribution function h i will satisfy the gyrokinetic equation (see, e.g., , rather than the drift-kinetic one (A 2). It is straightforward to show [by, e.g., expanding the Bessel functions therein for k ⊥ ρ i ≫ 1] that, to leading order in the mass-ratio expansion, the gyrokinetic equation is solved by The contributions to quasineutrality (A 5) arising from the next-order solution will be of the size which can be safely neglected. Thus, the ion dynamics do not enter anywhere into our equations, which is the approximation of 'adiabatic ions'. Given that no further reference will be made to the ion temperature gradient L Ti , we henceforth denote the electron temperature gradient L Te = L T .

B.3. Electron fluid equations
We now proceed with our derivation of the electron fluid equations. It will turn out that the ordering (B 15) of perpendicular length scales means that no finite-Larmorradius (FLR) effects need be retained within our equations -these can only enter at second order within our expansion, but they are negligible even at this order (see appendix B.3.3). Furthermore, the ordering of the perpendicular and parallel magnetic field perturbations (B 16) implies that both δB ⊥ and δB ∥ can be neglected at all orders in our expansion. We thus adopt the drift-kinetic equation ( Given the ordering of timescales (B 2), the collision operator on the right-hand side of (A 2) is dominant to leading order: where C (l) ee is given by (A 4) for s = s ′ = e, and is the pitch-angle scattering (Lorentz) collision operator, valid to leading order in the mass ratio. We multiply (B 21) by h (0) e /f 0e and integrate over the entire phase space, yielding Both terms in (B 23) are negative definite and must vanish individually, meaning that the solution is constrained to be a perturbed Maxwellian with no mean flow (Helander & Sigmar 2005), viz., where φ = eϕ/T 0e , and we have imposed the solvability conditions (B 25) in order to determine uniquely the density δn e and temperature δT e moments in (B 24). Note that, in general, the Lorentz collision operator constrains the electron distribution function to be isotropic in the frame moving with the parallel ion velocity. However, the parallel ion velocity is zero to all orders within our expansion in σλ ei /L [given the adiabatic ion solution (B 18)], meaning that the electron distribution function will have no parallel velocity moment to leading order.
We are now in a position to simplify the quasineutrality constraint (A 5). Using the solutions (B 18) and (B 24), it straightforwardly becomes (3.7).

B.3.2. First order: parallel flows
The parallel flows are determined self-consistently from the leading-order perturbations at the next order in our expansion, viz., h (1) e is determined by the solution of the Spitzer-Härm problem (Spitzer & Härm 1953;Braginskii 1965;Helander & Sigmar 2005): This can be inverted for h (1) e by means of a standard variational method. We define the functional: where ⟨. . . , . . . ⟩ denotes an inner product in velocity space weighted by the inverse of the electron (Maxwellian) equilibrium f 0e . Then, considering small variations h e = h min + δh and using the self-adjointness of the linearised collision operator, it is straightforward to show that the functional Σ[h e ] has a minimum at h min = h (1) e , for any variation δh (see, e.g., Helander & Sigmar 2005). Given that the spherical harmonics are eigenfunctions of the linearised collision operator, we choose to expand our distribution function in terms of spherical coordinates in velocity space (x, α, β), with x = v 2 /v 2 the , as Truncating (B 28) at p = 3, and demanding that the functional (B 29) be stationary with respect to variations in the coefficients a p , we find that h (1) e = a 0 + a 1 L (3/2) 1 The reduction of these equations to (3.8)-(3.9) occurs for very steep electron-temperature gradients, in the limit L B /L T → ∞.
The presence of the magnetic-drift terms in (C 1)-(C 2) introduces another instability into the system, the curvature-mediated ETG (cETG) instability (Horton et al. 1988;Adkins et al. 2022), which can, and generally does, modify its turbulent-transport properties. In particular, the turbulence theory of §5 assumed that the sETG instability was the dominant source of energy injection; this is only the case at sufficiently large L B /L T , meaning that we would expect departures from the behaviour observed in §5 to be most significant for L B /L T of order unity. A series of simulations were conducted in which L B /L T was varied, with all other parameters being kept the same as in the baseline simulation (see table 1); the heat flux from these simulations is plotted in figure 13. It is readily apparent that the introduction of finite magnetic-field gradients leads to a failure of saturation for all simulations where L B /L T is above the linear critical gradient for the cETG instability: L B L T > 1 2 τ + 40 9 1 τ 2 , (C 3) a threshold that can be derived straightforwardly from (C 1) and (C 2) in the twodimensional limit. This lack of saturation appears to persist irrespective of changes in box size, aspect ratio, and resolution in any (or all) of the coordinate directions. As discussed in §5.5, the adiabatic ion response (3.7) causes the nonlinearity in the electron-scale continuity equation (3.8) to vanish identically, a property that is shared by (C 1), meaning that the system lacks any nonlinearity capable of generating two-dimensional secondary instabilities that are responsible for the generation of zonal flows and destruction of streamer structures (see references in §5.5). Indeed, the lack of saturation in the case of our ETG simulations appears to be due to the inability of the system to break apart the streamers created by the cETG instability; the existence of such streamers causes the heat flux to diverge as they 'short circuit' the heat transport across the radial domain. Even if the simulation initially appears to saturate after the linear phase, it eventually forms these large-scale streamers, which appear to be immune to all types of nonlinear shearing, as seen clearly in the real-space snapshots of cETG turbulence shown in figures 14 and 15.
This perhaps confirms the view of Hammett et al. (1993) that the adiabatic ion response (3.7) is insufficient to saturate ETG-scale turbulence in the presence of finite magnetic-field gradients, and one may have to resort to more inclusive closures for the ions. One such closure including scales comparable to the ion-Larmor radius is (see, e.g., Figure 13: Time traces of the instantaneous heat flux from simulations with finite LB/LT , with the limit of LB/LT → ∞ shown for comparison. All parameters are the same as the baseline simulation (see table 1), and the heat flux is normalised to (ρ ⊥ /ρe)QgBe. The heat flux grows without bound in all simulations with (finite) LB/LT above the linear critical gradient (C 3) (≈ 2.72 forτ = 1), with the rate of divergence decreasing as LB/LT is increased.

Adkins et al. 2022)
δn e n 0e = −τ −1 φ + 1 n 0i d 3 v ⟨g i ⟩ r , (C 4) whereτ −1 is now an operator defined as follows: and the operatorΓ 0 can be expressed in Fourier space in terms of the modified Bessel function of the first kind: Γ 0 = I 0 (α i )e −αi , where α i = (k ⊥ ρ i ) 2 /2. The presence of the non-adiabatic ion distribution function g i in (C 4), however, means that one would have also to include a self-consistent treatment of ions in order to make use of this closure (g i = 0 is not a solution to the ion gyrokinetic equation in the presence of finite magnetic drifts). Should this, or other similar closures, allow for saturation, this would imply that one must always appeal to (elements of) ion-scale physics for saturation of electrostatic cETG-driven turbulence. The extent to which such considerations are practically relevant, however, depends on whether or not the system being considered contains any electromagnetic physics, and thus on the value of the (electron) plasma beta β e . Indeed, for β e ≳ m e /m i , the "flux-freezing scale" d e = ρ e / √ β e is encountered before (i.e., is smaller than) the ion Larmor radius ρ i when moving towards larger perpendicular scales. Provided that the wavenumber interval between d e and ρ i is sufficiently wide to allow for the presence of electron-scale, electromagnetic instabilities, the mechanisms of saturation in such a system could be vastly different than in the electrostatic regime. This is a subject of ongoing research. At these early times, the turbulence appears similar to that of saturated sETG turbulence for LB/LT → ∞ (cf. figure 11), despite the eventual lack of saturation (see figure 15).