2. Gyrokinetics

We consider turbulent fluctuations in magnetised plasmas that satisfy the GK ordering $k_\perp \rho _s \sim k_\parallel L \sim 1$ and $\omega / \Omega _s \sim \rho _s / L \ll 1$ , where $k_\perp$ and $k_\parallel$ are the typical perpendicular and parallel (to the mean magnetic field) wavenumbers, $\rho _s$ and $\Omega _s$ are the Larmor radius and frequency of the charged particles of species $s$ , $\omega$ is the inverse time scale associated with the turbulent fluctuations, and $L$ is the length scale of variation of the plasma equilibrium. Under this ordering, we expand the distribution function for each species into equilibrium and fluctuating parts $f_{s} = F_{s} + \delta f_{s}$ to obtain the GK equation that governs the dynamics of the fluctuations.Footnote ¹ With the additional assumption that the plasma beta $\beta _s \equiv 8\pi n_{s} T_{s} / B^2$ is small, $n_{s}$ and $T_{s}$ being the equilibrium density and temperature of species $s$ , respectively, we can neglect the fluctuations of the magnetic field, leading to

(2.1) \begin{align} \left (\frac {\partial }{\partial t} + \boldsymbol{u}\boldsymbol {\cdot }\frac {\partial }{\partial \boldsymbol{R}_s}\right ) \left (h_s - \frac {q_s \langle \phi \rangle _{\boldsymbol{R}_s}}{T_{s}} F_{s}\right ) + w_\parallel \hat {\boldsymbol{b}} \boldsymbol {\cdot } \frac {\partial h_s}{\partial \boldsymbol{R}_s} + \boldsymbol{v}_{ds}& \boldsymbol {\cdot } \frac {\partial h_s}{\partial \boldsymbol{R}_s} + \boldsymbol{v}_E \boldsymbol {\cdot } \frac {\partial }{\partial \boldsymbol{R}_s} \left (F_{s} + h_s\right ) \nonumber \\ & = \sum _{s^\prime }\left \langle C_{s s^\prime } \right \rangle _{\boldsymbol{R}_s} \end{align}

where the perturbed distribution function of species $s$ is

(2.2) \begin{align} \delta f_{s}(\boldsymbol{r}, \boldsymbol{w}) = h_s(\boldsymbol{R}_s, \varepsilon _s, \mu _s) - \frac {q_s \phi (\boldsymbol{r})}{T_{s}}F_{s}, \end{align}

$\boldsymbol{R}_s = \boldsymbol{r} - \hat {\boldsymbol{b}}\times \boldsymbol{w}/\Omega _s$ is the gyrocentre, $\boldsymbol{w} = \boldsymbol{v} - \boldsymbol{u}$ is the peculiar velocity, $\varepsilon _s = m_s w^2/2$ , $\mu _s = m_s w_\perp ^2 / 2B$ , $\phi$ is the perturbed electrostatic potential, $F_{s}$ is the equilibrium Maxwellian distribution with density $n_{s}$ and temperature $T_{s}$ , $\boldsymbol{u}$ is the equilibrium plasma flow (same for all species, see Abel et al. Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013), the magnetic drifts are $\boldsymbol{v}_{ds} = (\hat {\boldsymbol{b}}/2\Omega _s)\times (2w_\parallel ^2 \hat {\boldsymbol{b}}\boldsymbol {\cdot }{\boldsymbol {\nabla }}\hat {\boldsymbol{b}} + w_\perp ^2 {\boldsymbol {\nabla }} \ln B)$ , the perturbed $\boldsymbol{E}\times \boldsymbol{B}$ drift is $\boldsymbol{v}_E = (c/ B) \hat {\boldsymbol{b}}\times {\boldsymbol {\nabla }} \! \left \langle \phi \right \rangle _{\boldsymbol{R}_s}$ , $\hat {\boldsymbol{b}}$ is the unit vector parallel to the mean magnetic field, $q_s$ is the charge of species $s$ , $C_{s s^\prime }$ is the linearised Fokker–Planck operator for collisions between particles of species $s$ and $s^{\prime}$ , and $\left \langle \dots \right \rangle _{\boldsymbol{R}_s}$ denotes the standard gyroaverage. A comprehensive derivation of the GK equations can be found from, e.g. Abel et al. (Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013) or Catto (Reference Catto2019). Note that the theoretical analysis presented in § 3 does not depend on a particular coordinate system, i.e. the precise choice of radial, poloidal and parallel coordinates, labelled $x$ , $y$ and $z$ , respectively, will be irrelevant.

The nonlinear-advection and linear-drive terms in (2.1) are

(2.3) \begin{align} \boldsymbol{v}_E \boldsymbol {\cdot } \frac {\partial h_s}{\partial \boldsymbol{R}_s} &= \frac {c}{B} \hat {\boldsymbol{b}} \boldsymbol {\cdot } ({\boldsymbol {\nabla }} x \times {\boldsymbol {\nabla }} y) \left \lbrace \langle \phi \rangle _{\boldsymbol{R}_s}, h_s \right \rbrace, \end{align}

(2.4) \begin{align} \boldsymbol{v}_E \boldsymbol {\cdot } \frac {\partial F_{s}}{\partial \boldsymbol{R}_s} &= -\frac {c}{B} \hat {\boldsymbol{b}} \boldsymbol {\cdot } ({\boldsymbol {\nabla }} x \times {\boldsymbol {\nabla }} y) \frac {\partial \langle \phi \rangle _{\boldsymbol{R}_s}}{\partial y} \frac {\partial F_{s}}{\partial x}, \\[6pt] \nonumber \end{align}

respectively, where $\lbrace f, g \rbrace = (\partial _x f) (\partial _y g) - (\partial _y f) (\partial _x g)$ . The nonlinear term (2.3) expresses the advection of the perturbed distribution function by the perturbed $\boldsymbol{E}\times \boldsymbol{B}$ flow, while the linear term (2.4) represents the injection of free energy by the radial gradients of the equilibrium (via the advection of that equilibrium by the perturbed flows). The electrostatic GK equation is closed by the quasineutrality condition:

(2.5) \begin{align} \sum _s q_s \int \text {d}^{3} \boldsymbol{w} \ \delta f_{s} = 0, \end{align}

where the velocity integral is evaluated at fixed $\boldsymbol{r}$ .

Finally, fluctuations that evolve according to (2.1) can be shown to satisfy a free-energy conservation law (Abel et al. Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013) of the form

(2.6) \begin{align} \frac {\textrm {d} W}{\textrm {d} t} = I - D, \end{align}

where the free-energy density $W$ in a plasma of volume $V$ is given by

(2.7) \begin{align} W = \sum _s \int \frac {\textrm {d}^3{\boldsymbol{r}}}{V} \int \text {d}^{3} \boldsymbol{w} \ \frac {T_{s} \delta f_{s}^2}{2F_{s}}. \end{align}

The dissipation $D$ in (2.6) arises due to particle collisions and is a sink of free energy. Its precise form will not be needed here. The free-energy injection rate $I$ depends on the gradients of the equilibrium distribution $F_{s}$ and can be written asFootnote ²

(2.8) \begin{align} I = -\sum _s \left [\varGamma _s T_{s} \left (\frac {\partial \ln n_{s}}{\partial x} - \frac {3}{2}\frac {\partial \ln T_{s}}{\partial x}\right ) + Q_s \frac {\partial \ln T_{s}}{\partial x}\right ], \end{align}

where we have defined the flux of particles $\varGamma _s$ and the heat (or energy) flux $Q_s$ due to species $s$ as

(2.9) \begin{align} \varGamma _s \equiv &\int \frac {\textrm {d}^3 \boldsymbol{r}}{V} \int \text {d}^{3} \boldsymbol{w} \ (\boldsymbol{v}_E \boldsymbol {\cdot } {\boldsymbol {\nabla }} x) \delta f_{s}, \end{align}

(2.10) \begin{align} Q_s \equiv &\int \frac {\textrm {d}^3 \boldsymbol{r}}{V} \int \text {d}^{3} \boldsymbol{w} \ (\boldsymbol{v}_E \boldsymbol {\cdot } {\boldsymbol {\nabla }} x) \frac {m_s v^2}{2} \delta f_{s} . \\[6pt] \nonumber \end{align}

In the most general case, $I$ depends on both fluxes and can be estimated as

(2.11) \begin{align} I \sim \frac {\varGamma _s T_{s}}{L_{n_s}}\sim \frac {Q_s}{L_{T_s}}, \end{align}

where no additional orderings have been imposed on the density $L_{n_s}^{-1} \equiv - \partial _x \ln n_{s}$ and the temperature $L_{T_s}^{-1} \equiv -\partial _x \ln T_{s}$ gradients, viz. $L_{n_s} \sim L_{T_s} \sim L$ . In this work, we consider only temperature-gradient-driven instabilities in cases where $\varGamma _s=0$ and so our main focus will be on $Q_s$ . Nevertheless, the arguments presented in § 3 are readily generalisable to cases where the injection of energy is dominated by the particle flux rather than the heat flux.

3. Nonlinear saturation

Magnetised-plasma turbulence exhibits a broad range of different saturation mechanisms and so a universal theory of turbulent saturation under the influence of flow shear is not feasible. Instead, here we focus on one particular type of saturated state, viz. that for which the following two assumptions hold: (i) there is a scale separation between the energy-injection scale (the outer scale), where linear instabilities inject free energy into the fluctuations, and the dissipation scale, where fluctuations lose energy to dissipative effects (viz., ultimately, particle collisions); and (ii) the transfer of energy between these scales is realised by a local (in scale) energy cascade. This allows us to define the inertial range as the range of scales between the outer and dissipation scales where the local energy cascade takes place. Let us revisit the current understanding of how such an energy cascade determines the properties of the saturated turbulence.

3.1. Outer scale, free-energy cascade and turbulent heat flux

The form of the GK nonlinearity (2.3) implies that the nonlinear time $\tau_{ \mathrm{nl}}$ at poloidal scale $k_y$ satisfies

(3.1) \begin{align} \tau _{\mathrm{nl}}^{-1} \sim \Omega _s \rho _s^2 k_x k_y \overline {\varphi }, \end{align}

where $\overline {\varphi }$ is a measure of the characteristic amplitude of the normalised electrostatic potential $\varphi \equiv q_s \phi / T_{s}$ at scale $k_y$ and $s$ is any reference particle species. In (3.1), the radial scale $k_x$ is an implicit function of $k_y$ , viz. at each $k_y$ , the turbulent fluctuations have a typical radial scale $k_x$ that depends on $k_y$ . With this in mind, (3.1) can also be written as

(3.2) \begin{align} \tau _{\mathrm{nl}}^{-1} \sim {\mathcal {A}} \Omega _s \rho _s^2 k_y^2 \overline {\varphi }, \end{align}

where we have defined the fluctuation aspect ratio at scale $k_y$ as ${\mathcal {A}} \equiv k_x/k_y$ . This aspect ratio will play a critical role in the theory of sheared turbulence laid out in § 3.2. Note that the precise definition of $\overline {\varphi }$ is not important because the phenomenological theory that is to follow predicts only scalings; nevertheless, to make things specific, one possible such definition is

(3.3) \begin{align} \overline {\varphi }^2 \equiv \int _{|k_y^{\prime}|\gt k_y} \text {d} k_y^{\prime} \int _{-\infty }^{+\infty } \text {d} k_x^{\prime} \int \frac {\textrm {d} z}{L_\parallel } \: |\varphi _{\boldsymbol{k}_\perp ^{\prime}}|^2, \end{align}

where $\varphi _{{\boldsymbol{k}}_\perp }$ is the two-dimensional spatial Fourier transform of $\varphi$ in $x$ and $y$ , i.e. in the plane perpendicular to the equilibrium magnetic field, and $L_\parallel$ is the parallel size of the integration domain of (2.1) (assumed finite). By Parseval’s theorem, the contributions to the free energy (2.7) at each scale are proportional to the squared amplitude of the fluctuations and so, in view of (2.2) and (2.5), to $\overline {\varphi }^2$ .

If $k_x \sim k_y \sim k_\perp$ (i.e. ${\mathcal {A}} \sim 1$ ) is satisfied throughout the inertial range, (3.1) implies that the free-energy flux $\varepsilon$ through scale $k_\perp$ satisfiesFootnote ³

(3.4) \begin{align} \varepsilon \sim \tau _{\mathrm{nl}}^{-1} n_{s} T_{s} \overline {\varphi }^2 \propto k_\perp ^2 \overline {\varphi }^3, \end{align}

where $\text {d}/\text {d} t \sim \tau _{\mathrm {nl}}^{-1}$ in the inertial range because the dynamics there is dominated by the nonlinear effects. Assuming a local, constant-flux cascade, viz. that $\varepsilon$ is constant throughout the inertial range, and using (3.1), we conclude that

(3.5) \begin{align} \overline {\varphi } \propto k_\perp ^{-2/3} \implies \tau _{ \mathrm{nl}}^{-1} \propto k_\perp ^{4/3} \end{align}

in the inertial range. The assumption that ${\mathcal {A}} \sim 1$ in the inertial range is motivated by the fact that the nonlinearity in (2.1) is isotropic in the perpendicular plane. However, this is not a sufficient condition for ${\mathcal {A}} \sim 1$ . For example, reduced magnetohydrodynamic (RMHD) turbulence, whose nonlinearity is also isotropic in the plane perpendicular to the mean magnetic field, is known to have fluctuations that are anisotropic in the perpendicular plane and whose anisotropy depends on the scale, thus introducing a nontrivial scale-dependent factor into the RMHD version of (3.1) (Boldyrev Reference Boldyrev2006; Boldyrev, Mason & Cattaneo Reference Boldyrev, Mason and Cattaneo2009; Mallet et al. Reference Mallet, Schekochihin, Chandran, Chen, Horbury, Wicks and Greenan2016; Mallet & Schekochihin Reference Mallet and Schekochihin2017; Schekochihin Reference Schekochihin2022). Nevertheless, assuming ${\mathcal {A}} \sim 1$ in the inertial range is not unreasonable and it agrees with our numerical observations reported in § 4.

Assuming that the rate of energy injection is determined by the linear-instability growth rate $\gamma _{\boldsymbol{k}}$ and that the latter satisfies $\gamma _{\boldsymbol{k}} \propto k_y$ ,Footnote ⁴ the inertial-range nonlinear rate (3.5) has a steeper dependence on $k_y$ than the injection rate. The outer scale will then be the scale at which the rates of nonlinear mixing and linear injection balance (Barnes et al. Reference Barnes, Parra, Highcock, Schekochihin, Cowley and Roach2011; Adkins et al. Reference Adkins, Schekochihin, Ivanov and Roach2022, Reference Adkins, Ivanov and Schekochihin2023):

(3.6) \begin{align} \left (\tau _{\mathrm{nl}}^{\textrm {o}}\right )^{-1} \sim \gamma ^{\textrm {o}}. \end{align}

Here and in what follows, the superscript ‘ $^{\textrm {o}}$ ’ denotes quantities associated with the outer scale. The inertial range is thus located at $k_y \gt k_y^{\textrm {o}}$ , where the nonlinear interactions dominate the linear injection rate (see figure 1).

Figure 1. An illustration of the relationship between the nonlinear mixing rate $\tau _{\mathrm{nl}}^{-1}$ , the energy-injection rate $\gamma _{\boldsymbol{k}}$ and the location of the outer scale, where $\tau _{\mathrm{nl}}^{-1} \sim \gamma _{\boldsymbol{k}}$ . The scaling $\tau _{\mathrm{nl}}^{-1}\propto k_y^{4/3}$ is a consequence of the local energy cascade and is thus valid only in the inertial range $k_y \gt k_y^{\textrm {o}}$ (see the discussion in § 3.1).

3.1.1. Heat flux

Assuming that the heat flux $Q_s$ is dominated by contributions from the outer scale, we can estimate it, in view of its definition (2.10), as follows:

(3.7) \begin{align} Q_s \sim Q^{\textrm {o}}_s \sim n_{s} T_{s} v_{\textrm {th}s} k_y^{\textrm {o}} \rho _s \left (\overline {\varphi }^{\textrm {o}}\right )^2. \end{align}

This is justified as long as the spectrum of the fluctuations decays sufficiently quickly in the inertial range; specifically, we require the fluctuation amplitudes to decay faster than $k_\perp ^{-1/2}$ , which is readily satisfied by (3.5). We also assume that the phase relationship between $\varphi$ and $\delta f_{s}$ does not introduce any nontrivial factors – technically, (3.7) is an upper bound for (2.10). Using (3.2) and (3.6), we can rewrite (3.7) as

(3.8) \begin{align} \frac {Q_s}{n_{s} T_{s} v_{\textrm {th}s}} \sim \left (\frac {\gamma ^{\textrm {o}}}{\Omega _s}\right )^2 \frac {1}{(k_y^{\textrm {o}} \rho _s)^3 ({\mathcal {A}}^{\textrm {o}})^2}. \end{align}

Therefore, to determine $Q_s$ , we need to know the energy-injection rate $\gamma ^{\textrm {o}}$ (or equivalently $\tau _{\mathrm{nl}}^{\textrm {o}}$ ), the poloidal wavenumber $k_y^{\textrm {o}}$ and the fluctuation aspect ratio ${\mathcal {A}}^{\textrm {o}}$ at the outer scale. If $\gamma ^{\textrm {o}}\sim \gamma _{{\boldsymbol{k}}^{\textrm {o}}}$ , where $\gamma _{{\boldsymbol{k}}^{\textrm {o}}}$ is the growth rate at the outer scale, then only two of $\gamma ^{\textrm {o}}$ , $k_y^{\textrm {o}}$ and ${\mathcal {A}}^{\textrm {o}}$ are independent. Thus, we require additional assumptions. There are multiple ways to proceed.

In the absence of flow shear, Barnes et al. (Reference Barnes, Parra, Highcock, Schekochihin, Cowley and Roach2011) posit: (i) that the outer scale is governed by the ‘critical balance’ of $\gamma ^{\textrm {o}}$ and $(\tau _{\mathrm{nl}}^{\textrm {o}})^{-1}$ with the parallel-streaming rate across the plasma connection length, $\omega _\parallel ^{\textrm {o}} \sim v_{\textrm {th}s} / qR$ , where $q$ and $R$ are the safety factor and major radius, respectively; (ii) that the outer-scale fluctuations are isotropic, ${\mathcal {A}}^{\textrm {o}} \sim 1$ ; and (iii) that the energy-injection rate is given by a simple estimate of the growth rate of temperature-gradient-driven instabilities, $\gamma ^{\textrm {o}} \sim k_y^{\textrm {o}} \rho _s v_{\textrm {th}s}/L_{T_s}$ . Combined with (3.8), assumptions (i)–(iii) imply

(3.9) \begin{align} \frac {Q_s}{n_{s} T_{s} v_{\textrm {th}s}} \sim \left (\frac {\rho _s}{R}\right )^2 \Big (\frac {R}{L_{T_s}}\Big )^3 q. \end{align}

Note that Barnes et al. (Reference Barnes, Parra, Highcock, Schekochihin, Cowley and Roach2011) studied ion-scale turbulence, which amounts to setting $s = i$ in the above arguments.

A modification of these results, backed by experimental (Ghim et al. Reference Ghim, Schekochihin, Field, Abel, Barnes, Colyer, Cowley, Parra, Dunai and Zoletnik2013) and theoretical (Nies et al. Reference Nies, Parra, Barnes, Mandell and Dorland2024) evidence, is to replace assumption (ii) in the arguments by Barnes et al. (Reference Barnes, Parra, Highcock, Schekochihin, Cowley and Roach2011) by a ‘grand critical balance’

(3.10) \begin{align} \gamma ^{\textrm {o}} \sim (\tau _{\mathrm{nl}}^{\textrm {o}})^{-1} \sim \omega _\parallel ^{\textrm {o}} \sim \omega ^{\textrm {o}}_{d,x} \end{align}

between all the aforementioned rates and the radial magnetic-drift frequency $\omega ^{\textrm {o}}_{d,x} \sim k_x^{\textrm {o}} \rho _s v_{\textrm {th}s}/R$ at the outer scale. This implies

(3.11) \begin{align} {\mathcal {A}}^{\textrm {o}} \sim \frac {R}{L_{T_s}}, \end{align}

which, together with (3.8), results in the following scaling for the heat flux:

(3.12) \begin{align} \frac {Q_s}{n_{s} T_{s} v_{\textrm {th}s}} \sim \left (\frac {\rho _s}{R}\right )^2 \frac {R}{L_{T_s}} q. \end{align}

In the rest of this paper, we consider the influence of mean flow shear on the saturated state. We will not discuss the details of how the outer scale is determined in the case of zero imposed flow shear, but assume that the system does indeed have a well-defined zero-shear saturated state, that the outer-scale nonlinear rate is governed by (3.2) and (3.6), and that (3.8) is a good estimate for the heat flux. Thus, our arguments will hold regardless of whether the zero-shear outer scale is chosen à la Barnes et al. (Reference Barnes, Parra, Highcock, Schekochihin, Cowley and Roach2011), through a grand critical balance, or otherwise.

3.2. Perpendicular flow shear

For the remainder of this article, we assume an equilibrium shear flow in the direction perpendicular to the mean magnetic field and with a linear profile: $\boldsymbol{u} = \gamma _{E} x \hat {\boldsymbol{y}}$ , where $\gamma _{E}$ is the shearing rate.Footnote ⁵ In the presence of such a flow, the GK equation (2.1) is no longer homogeneous in $x$ . For brevity, we henceforth drop the species subscript from the heat flux $Q$ .

To understand the effect of the flow shear on the fluctuations, it is instructive to consider a patch of turbulence in which the magnetic field can be considered locally constant and oriented along the $z$ -direction; i.e. this patch is approximated as a ‘slab’. One can then perform a coordinate transformation from the original (laboratory) frame to the so-called shearing frame (Newton, Cowley & Loureiro Reference Newton, Cowley and Loureiro2010; Schekochihin, Highcock & Cowley Reference Schekochihin, Highcock and Cowley2012):

(3.13) \begin{align} t^{\prime} = t, \ x^{\prime} = x, \ y^{\prime} = y - x\gamma _{E} t, \ z^{\prime} = z. \end{align}

The substitution of (3.13) into the GK equation (2.1) eliminates the radially inhomogeneous advection term $\boldsymbol{u}\boldsymbol {\cdot }{\boldsymbol {\nabla }}$ at the cost of introducing an inhomogeneity in time via the $\partial _x$ derivatives. Consequently, the laboratory-frame radial wavenumber $k_x$ of a fluctuation advected by the mean flow, i.e. of a fluctuation with a given fixed wavenumber ${\boldsymbol{k}}^{\prime}$ in the shearing frame, satisfies

(3.14) \begin{align} k_x = k_x^{\prime} - k_y^{\prime} \gamma _{E} t^{\prime}. \end{align}

Crucially, the nonlinear interactions (2.3) and the linear drive (2.4) have the same form in both the laboratory frame and the shearing frame; therefore, (3.14) captures completely the effects of flow shear in the shearing frame. Equation (3.14) tells us that the shearing action of the perpendicular flow, which results in a ‘tilting’ of the eddies (Fox et al. Reference Fox, van Wyk, Field, Ghim, Parra and Schekochihin2017), is equivalent to a ‘drift’ in Fourier space of the radial wavenumber $k_x$ of the turbulent fluctuations.

Let us consider introducing flow shear into a system that, in its absence, would reach a saturated state by establishing a local energy cascade. As discussed in § 3.1, the transport properties (e.g. the radial heat flux $Q$ ) of such a system are dominated by the fluctuations at the outer scale. The lifetime of these fluctuations is given by the outer-scale nonlinear time, which, according to (3.6), is $\tau _{\mathrm{nl}}^{\textrm {o}}(0) \sim \gamma ^{\textrm {o}}(0)^{-1}$ , where we will use the notation $\gamma ^{\textrm {o}}(\gamma _{E})$ to denote the dependence of outer-scale quantities on the flow shear, so $\gamma ^{\textrm {o}}(0)$ is the outer-scale injection rate in the absence of it. We shall distinguish two different regimes of sheared turbulence: a weak-shear regime with $\gamma _{E} \lt \gamma ^{\textrm {o}}(0)$ and a strong-shear regime with $\gamma _{E} \gt \gamma ^{\textrm {o}}(0)$ . This distinction is motivated by the so-called ‘quench’ rule (Waltz, Kerbel & Milovich Reference Waltz, Kerbel and Milovich1994, Reference Waltz, Dewar and Garbet1998; Kobayashi & Rogers Reference Kobayashi and Rogers2012; Ivanov et al. Reference Ivanov, Schekochihin, Dorland, Field and Parra2020), according to which flow shear is able to suppress the energy injection associated with some linearly unstable modes only if the shearing rate is comparable to the growth rate of those modes. If this is true, the flow shear should be unable to stifle energy injection at the outer scale if $\gamma _{E} \lt \gamma ^{\textrm {o}}(0)$ and so the outer-scale injection rate should remain independent of $\gamma _{E}$ in the weak-shear regime, i.e. $\gamma ^{\textrm {o}}(\gamma _{E}) \approx \gamma ^{\textrm {o}}(0)$ . In contrast, in the strong-shear regime, we expect that the injection rate will be modified by the presence of the flow shear.

Let us analyse the physics of both regimes, starting with the weak-shear one.

3.2.1. Weak-shear regime

Let us consider more carefully the influence of flow shear on the outer scale in the case $\gamma _{E} \lt \gamma ^{\textrm {o}}(0)$ . As just discussed, in this regime, the outer-scale-injection and nonlinear-mixing rates should remain approximately the same as those at $\gamma _{E} = 0$ , viz.

(3.15) \begin{align} \tau _{\mathrm{nl}}^{\textrm {o}}(\gamma _{E})^{-1} \sim \gamma ^{\textrm {o}}(\gamma _{E}) \sim \gamma ^{\textrm {o}}(0) \sim \tau _{\mathrm{nl}}^{\textrm {o}}(0)^{-1}. \end{align}

Assuming that the injection rate $\gamma ^{\textrm {o}}$ is determined by the linear growth rate at the outer scaleFootnote ⁶ and that the latter is (at least approximately) only a function of $k_y$ ,Footnote ⁷ we conclude that the poloidal wavenumber of the outer-scale eddies is also set by its value at $\gamma _{E}=0$ and independent of $\gamma _{E}$ in the weak-shear regime, viz.

(3.16) \begin{align} k_y^{\textrm {o}}(\gamma _{E}) \sim k_{y}^{\textrm {o}}(0). \end{align}

However, the assumption that $\gamma _{\boldsymbol{k}}$ is only a weak function of $k_x$ means that one cannot make a similar statement about the radial wavenumber $k_x^{\textrm {o}}(\gamma _{E})$ . Indeed, approximating the lifetime of the outer-scale fluctuations as equal to the nonlinear mixing time $\tau _{\mathrm{nl}}^{\textrm {o}}$ , the wavenumber drift (3.14), together with (3.15) and (3.16), suggests that

(3.17) \begin{align} k_x^{\textrm {o}}(\gamma _{E}) &\sim k_{x}^{\textrm {o}}(0) + k_{y}^{\textrm {o}}(0)\tau _{\mathrm{nl}}^{\textrm {o}}(0)\gamma _{E} \nonumber \\ &\sim k_{x}^{\textrm {o}}(0) \left [1 + \frac {\gamma _{E}}{{\mathcal {A}}^{\textrm {o}}(0) \gamma ^{\textrm {o}}(0)}\right ], \\[6pt] \nonumber \end{align}

where ${\mathcal {A}}^{\textrm {o}}(0)=k_{x}^{\textrm {o}}(0)/k_{y}^{\textrm {o}}(0)$ is the fluctuation aspect ratio at the outer scale at $\gamma _{E} = 0$ . Therefore,

(3.18) \begin{align} \frac {k_x^{\textrm {o}}(\gamma _{E})}{k_{x}^{\textrm {o}}(0)}\sim \frac {{\mathcal {A}}^{\textrm {o}}(\gamma _{E})}{{\mathcal {A}}^{\textrm {o}}(0)} \sim 1 + \frac {\gamma _{E}}{\gamma _{\textrm {c}}}, \end{align}

where we have introduced the critical shearing rate

(3.19) \begin{align} \gamma _{\textrm {c}} \equiv {\mathcal {A}}^{\textrm {o}}(0) \gamma ^{\textrm {o}}(0). \end{align}

Then, (3.8) implies that the radial turbulent heat flux satisfies

(3.20) \begin{align} \frac {Q(\gamma _{E})}{Q(0)}\sim \left [\frac {{\mathcal {A}}^{\textrm {o}}(0)}{{\mathcal {A}}^{\textrm {o}}(\gamma _{E})}\right ]^2 \sim \frac {1}{(1+\gamma _{E}/\gamma _{\textrm {c}})^2}, \end{align}

where $Q(0)$ is the heat flux at $\gamma _{E} = 0$ . Note that at no step leading to (3.20) did we use any formulae from § 3.1 that relied on isotropy, which would otherwise have restricted us to ${\mathcal {A}}^{\textrm {o}} \sim 1$ . Expressions (3.19) and (3.20) predict that the transport properties in the weak-shear regime are determined by the ratio of the radial and poloidal wavenumbers of the outer-scale eddies, ${\mathcal {A}}^{\textrm {o}}(0) = k_{x}^{\textrm {o}}(0)/k_{y}^{\textrm {o}}(0)$ . If the unsheared fluctuations have ${\mathcal {A}}^{\textrm {o}}(0) \sim 1$ , then (3.19) implies that $\gamma _{\textrm {c}} \sim \gamma ^{\textrm {o}}(0)$ . As the weak-shear regime is characterised by $\gamma _{E} \lt \gamma ^{\textrm {o}}(0)$ , (3.18) implies that ${\mathcal {A}}^{\textrm {o}}(\gamma _{E}) \sim {\mathcal {A}}^{\textrm {o}}(0) \sim 1$ and thus $Q(\gamma _{E}) \sim Q(0)$ throughout the weak-shear regime. In other words, if the unsheared turbulence has ${\mathcal {A}}^{\textrm {o}}(0) \sim 1$ at the outer scale, shearing it with any $\gamma _{E} \lt \gamma ^{\textrm {o}}(0)$ will not reduce the turbulent transport by more than an order-unity amount – an unsurprising outcome.

However, due to the nature of the underlying linear instabilities, it is, in fact, often the case that the outer-scale eddies in temperature-gradient-driven turbulence are radially elongated. Such eddies, often called ‘streamers’, are a well-documented feature of this type of turbulence, especially in its electron-scale variety (Drake, Guzdar & Hassam Reference Drake, Guzdar and Hassam1988; Jenko et al. Reference Jenko, Dorland, Kotschenreuther and Rogers2000; Dorland et al. Reference Dorland, Jenko, Kotschenreuther and Rogers2000; Jenko Reference Jenko2006; Roach et al. Reference Roach2009; Colyer et al. Reference Colyer, Schekochihin, Parra, Roach, Barnes, Ghim and Dorland2017). A turbulent state dominated by streamers satisfies ${\mathcal {A}}^{\textrm {o}}(0) \ll 1$ and so $\gamma _{\textrm {c}} \ll \gamma ^{\textrm {o}}(0)$ . In this case, (3.18) implies that the outer-scale aspect ratio increases linearly with the flow shear due to the tilting of the eddies, viz.

(3.21) \begin{align} {\mathcal {A}}^{\textrm {o}}(\gamma _{E}) \sim {\mathcal {A}}^{\textrm {o}}(0) + \frac {\gamma _{E}}{\gamma ^{\textrm {o}}(0)}. \end{align}

Furthermore, (3.20) predicts that the heat flux will be suppressed if $\gamma _{\textrm {c}} \lesssim \gamma _{E} \ll \gamma ^{\textrm {o}}(0)$ . In particular, for intermediate values of the shearing rate that satisfy $\gamma _{\textrm {c}} \ll \gamma _{E} \ll \gamma ^{\textrm {o}}(0)$ , (3.21) becomes

(3.22) \begin{align} {\mathcal {A}}^{\textrm {o}}(\gamma _{E}) \sim \frac {\gamma _{E}}{\gamma ^{\textrm {o}}(0)}, \end{align}

and so, by (3.20),

(3.23) \begin{align} Q(\gamma _{E}) \propto \gamma _{E}^{-2}. \end{align}

If the shear is increased further, (3.21) implies that, at the transition from the weak- to the strong-shear regime, where $\gamma _{E} \sim \gamma ^{\textrm {o}}(0)$ , the outer-scale aspect ratio is ${\mathcal {A}}^{\textrm {o}}(\gamma _{E}) \sim 1$ , and so, by (3.20), the heat flux has been reduced by a large factor:

(3.24) \begin{align} \frac {Q[\gamma ^{\textrm {o}}(0)]}{Q(0)} \sim \left [{\mathcal {A}}^{\textrm {o}}(0)\right ]^2 \ll 1. \end{align}

A cautious reader may have spotted a potential clash between having ${\mathcal {A}}^{\textrm {o}}(0) \ll 1$ at the outer scale and the theory of the energy cascade laid out in § 3.1: there, we assumed that the fluctuations in the inertial range had ${\mathcal {A}} \sim 1$ , yet the inertial range must connect to the outer scale, where now, ${\mathcal {A}} \ll 1$ . There are two possible resolutions to this problem: (i) the scaling arguments presented in § 3.1 are, in fact, unchanged if the inertial range inherits the aspect ratio at the outer scale, i.e. if $k_x/k_y$ is a scale-independent, even if numerically small, number below the outer scale; or (ii) there exists a transition region below the outer scale, wherein the dependence of $k_x$ on $k_y$ is faster than linear so that $k_x$ gradually increases to match $k_y$ at some smaller scale, below which the scalings from § 3.1 become valid. Our numerical results, presented in § 4.1, are consistent with option (ii). In Appendix A, we develop a simple theory for the transition region.

Finally, let us mention that, while here we shall consider only cases where ${\mathcal {A}} \lesssim 1$ , this is not necessarily satisfied in all instances of fusion-relevant turbulence. For example, the large-temperature-gradient environment of the pedestal has been shown numerically to give rise to poloidally elongated turbulent fluctuations with ${\mathcal {A}} \gg 1$ (Parisi et al. Reference Parisi2020, Reference Parisi2022). As discussed in § 3.1.1, the ‘grand critical balance’ (3.10) leads to poloidally elongated eddies at large temperature gradients, as per (3.11). Recent numerical and analytical work by Nies et al. (Reference Nies, Parra, Barnes, Mandell and Dorland2024) suggests that such behaviour may indeed be consistent with strongly driven ion-temperature-gradient turbulence in axisymmetric toroidal geometry. Here, leaving the reader cognisant of these recent developments, we shall nevertheless focus on the case ${\mathcal {A}} \lesssim 1$ .

3.2.2. Strong-shear regime

In the strong-shear regime, defined by $\gamma _{E} \gt \gamma ^{\textrm {o}}(0)$ , the flow shear is strong enough to affect energy injection at the outer scale. In particular, it is no longer possible to excite fluctuations at wavenumbers where the growth rate is $\gamma _{\boldsymbol{k}} \lesssim \gamma _{E}$ (Waltz et al. Reference Waltz, Kerbel and Milovich1994, Reference Waltz, Dewar and Garbet1998). To compensate for this, the outer scale must adjust to match the shearing rate. Thus, we propose that, for $\gamma _{E} \gg \gamma ^{\textrm {o}}(0)$ , the outer scale will be governed by the balance of nonlinear, injection and shearing rates:

(3.25) \begin{align} \tau _{\textrm {nl}}^{\textrm {o}}(\gamma _{E})^{-1} \sim \gamma ^{\textrm {o}}(\gamma _{E}) \sim \gamma _{E}, \end{align}

as illustrated in figure 2. As always, the lifetime of turbulent eddies at the outer scale is set by the nonlinear time; (3.14) then implies that ${\mathcal {A}}^{\textrm {o}}(\gamma _{E}) \sim 1$ throughout the strong-shear regime. Assuming that the linear growth rate is $\gamma _{\boldsymbol{k}} \propto k_y$ , we expect that

(3.26) \begin{align} k_y^{\textrm {o}} \propto \gamma _{E}. \end{align}

This is intuitively clear: stronger flow shear pushes turbulence towards smaller (and thus faster) scales since the larger (and slower) eddies are more strongly affected by the shear. Consequently, (3.8), together with (3.25) and (3.26), implies

(3.27) \begin{align} Q(\gamma _{E}) \propto \gamma _{E}^{-1} \end{align}

in the strong-shear regime $\gamma _{E} \gt \gamma ^{\textrm {o}}(0)$ .

Figure 2. A qualitative illustration, analogous to figure 1, of the effect of strong flow shear $\gamma _{E} \gg \gamma ^{\textrm {o}}(0)$ , leading to the time scale balance (3.25) determining the outer scale. In this regime, $\gamma ^{\textrm {o}}(\gamma _{E}) \sim \gamma _{E}$ .

The outer-scale balance (3.25), and thus the scaling (3.27), cannot be satisfied for arbitrarily large values of flow shear because the linear growth rate $\gamma _{\boldsymbol{k}}$ is bounded by some $\gamma _{{\rm max}}$ , normally found at much larger wavenumbers than those associated with the dominant energy injection.Footnote ⁸ For $\gamma _{E} \gtrsim \gamma _{{\rm max}}$ , the system is no longer able to sustain the turbulent fluctuations because the shearing rate $\gamma _{E}$ cannot be matched by the rate of energy injection at any scale. Therefore, we expect a sharp cut-off in the fluctuations’ amplitude, and thus in the heat flux, as $\gamma _{E}$ becomes comparable to $\gamma _{{\rm max}}$ . Figure 3 summarises the expected dependence of the heat flux on $\gamma _{E}$ in both regimes.

4. Numerical results

To test the validity of the theory presented in § 3.2, we consider two different models of turbulence. The first (§ 4.1) is a two-fluid model that captures the dynamics of electrostatic fluctuations of density and electron temperature in a straight magnetic field. This turbulence is driven by the collisional slab ETG (sETG) instability (Adkins et al. Reference Adkins, Schekochihin, Ivanov and Roach2022) on scales between the ion and electron gyroradii. While this model is extremely simple, even simplistic, the benefit of using it is that its saturation mechanism has already been investigated in great detail and has been shown to conform to the picture of a local energy cascade outlined in § 3.1 (Adkins, Ivanov & Schekochihin Reference Adkins, Ivanov and Schekochihin2023). Therefore, it is a prime candidate for confirming the validity of the theory laid out in this paper. For our second set of simulations (§ 4.2), we employ the GK code GENE (Jenko et al. Reference Jenko, Dorland, Kotschenreuther and Rogers2000; Jenko Reference Jenko2000) to perform gyrokinetic flux-tube simulations of ITG-driven turbulence. This is a much more realistic model of plasma turbulence, and the ‘Cyclone base case’ used here is a setup that has been extensively studied in the literature (Lin et al. Reference Lin, Hahm, Lee, Tang and Diamond1999; Dimits et al. Reference Dimits2000; Barnes et al. Reference Barnes, Abel, Dorland, Ernst, Hammett, Ricci, Rogers, Schekochihin and Tatsuno2009; Highcock et al. Reference Highcock, Schekochihin, Cowley, Barnes, Parra, Roach and Dorland2012; Peeters et al. Reference Peeters, Rath, Buchholz, Camenen, Candy, Casson, Grosshauser, Hornsby, Strintzi and Weikl2016; Li et al. Reference Li, Terry, Whelan and Pueschel2021; C. J. et al. et al., Reference Ajay, Brunner and Ball2021; Volčokas et al. Reference Volčokas, Ball and Brunner2022; Hoffmann, Frei & Ricci Reference Hoffmann, Frei and Ricci2023; Lippert, Rath & Peeters Reference Lippert, Rath and Peeters2023; Tirkas et al. Reference Tirkas, Chen, Merlo, Jenko and Parker2023 constitute a small sample) – it is thus a natural testbed for any theory aspiring to tokamak relevance.

Figure 3. A qualitative diagram of the heat flux $Q$ as a function of the flow shear $\gamma _{E}$ in the case of (a) ${\mathcal {A}}^{\textrm {o}}(0) \ll 1$ and (b) ${\mathcal {A}}^{\textrm {o}}(0) \sim 1$ . In each case, there are two distinct regimes. (i) For $\gamma _{E} \lt \gamma ^{\textrm {o}}(0)$ , we have the weak-shear regime (§ 3.2.1), where, in the ${\mathcal {A}}^{\textrm {o}}(0) \ll 1$ case, we find $Q(\gamma _{E}) \propto (1 + \gamma _{E}/\gamma _{\textrm {c}})^{-2}$ [see (3.20)]. In contrast, if ${\mathcal {A}}^{\textrm {o}}(0) \sim 1$ , the flow shear is unable to affect significantly the fluctuations at the outer scale and, consequently, the heat flux is approximately independent of $\gamma _{E}$ . (ii) For $\gamma ^{\textrm {o}}(0) \lt \gamma _{E} \lt \gamma _{{\rm max}}$ , the system is in the strong-shear regime (§ 3.2.2), where the outer-scale injection rate is determined by the flow shear, viz. $\gamma ^{\textrm {o}}(\gamma _{E}) \sim \gamma _{E}$ . Here, ${\mathcal {A}}^{\textrm {o}}(\gamma _{E}) \sim 1$ at the outer scale, regardless of ${\mathcal {A}}^{\textrm {o}}(0)$ , and $Q(\gamma _{E})\propto \gamma _{E}^{-1}$ . Finally, the fluctuations, and hence the heat flux, are completely suppressed at $\gamma _{E} \gtrsim \gamma _{{\rm max}}$ .

4.1. Fluid ETG turbulence

In this section, we report numerical simulations in a triply periodic domain of size $L_x$ , $L_y$ and $L_\parallel$ in $x$ , $y$ and $z$ , respectively, of the following collisional slab ETG model (Adkins et al. Reference Adkins, Ivanov and Schekochihin2023):

(4.1) \begin{align} &\frac {\textrm {d}}{\textrm {d} t} \frac {\delta n_e}{n_{e}} + \frac {\partial u_{\parallel e}}{\partial z} = 0, \end{align}

(4.2) \begin{align} &\frac {\nu _{ei}}{c_1}\frac {u_{\parallel e}}{v_{\textrm {th}e}} = -\frac {v_{\textrm {th}e}}{2} \frac {\partial }{\partial z}\left [ \frac {\delta n_e}{n_{e}} - \varphi + \left (1 + \frac {c_2}{c_1}\right )\frac {\delta T_e}{T_{e}} \right ], \end{align}

(4.3) \begin{align} &\frac {\textrm {d}}{\textrm {d} t}\frac {\delta T_e}{T_{e}} - \frac {c_3v_{\textrm {th}e}^2}{3\nu _{ei}}\frac {\partial ^2}{\partial z^2} \frac {\delta T_e}{T_{e}} + \frac {2}{3} \left (1 + \frac {c_2}{c_1}\right )\frac {\partial u_{\parallel e}}{\partial z} = -\frac {\rho _e v_{\textrm {th}e}}{2 L_T}\frac {\partial \varphi }{\partial y}, \\[6pt] \nonumber \end{align}

where the ‘convective derivative’

(4.4) \begin{align} \frac {\textrm {d}}{\textrm {d} t} = \frac {\partial }{\partial t} + \gamma _{E} x \frac {\partial }{\partial y} + \frac {\rho _ev_{\textrm {th}e}}{2} \left (\hat {\boldsymbol{z}} \times {\boldsymbol {\nabla }} \varphi \right )\boldsymbol {\cdot }{\boldsymbol {\nabla }} + \nu _\perp \rho _e^4 \nabla _\perp ^4 \end{align}

includes the mean flow shear, the nonlinear advection by the perturbed $\boldsymbol{E}\times \boldsymbol{B}$ drift and hyperviscous dissipation. Appendix B describes some important details of the numerical implementation of the flow-shearing term $\gamma _{E} x \partial _y$ . The electron density is related to the electrostatic potential $\varphi \equiv e \phi / T_{e}$ by quasineutrality (2.5) combined with the assumption of adiabatic ions:

(4.5) \begin{align} \frac {\delta n_e}{n_{e}} = - \frac {Z T_{e}}{T_{i}} \varphi, \end{align}

where $Z = q_i / e$ , $q_i$ being the ion charge. The numerical coefficients $c_1$ , $c_2$ and $c_3$ arise from the physics of collisions and depend on $Z$ : 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). We used $Z = 1$ and $T_{i} = T_{e}$ for all simulations reported here. Finally, the electron heat flux (2.10) can be expressed using the Fourier amplitudes of the fluctuations as follows:

(4.6) \begin{align} Q = \frac {3}{2} n_{e}T_{e}v_{\textrm {th}e} \sum _{\boldsymbol{k}} i k_y\rho _e \varphi _{\boldsymbol{k}}^* \frac {\delta T_{e,{\boldsymbol{k}}}}{T_{e}}. \end{align}

Together, (4.1)–(4.3) and (4.5) form an asymptotic model derived in an electrostatic, collisional limit of GK in a straight and uniform magnetic field with $\hat {\boldsymbol{b}} = \hat {\boldsymbol{z}}$ (Adkins et al. Reference Adkins, Ivanov and Schekochihin2023). They describe the electrostatic dynamics of fluctuations on perpendicular and parallel scales that satisfy

(4.7) \begin{align} k_\parallel L_T \sim \sqrt {\sigma }, \quad k_\perp \rho _\perp \sim 1, \quad \rho _\perp \equiv \frac {\rho _e}{\sigma } \frac {L_T}{\lambda _{ei}}, \end{align}

where $L_T^{-1} \equiv -\partial \ln T_{e}/\partial x$ is the electron-temperature gradient, $\lambda _{ei}$ is the electron-ion mean free path and $\sigma$ is a formal scaling parameter that is arbitrary provided it satisfies $\beta _e \ll \sigma \ll 1$ . The fact that this parameter is arbitrary is a consequence of the scale invariance of the model (Adkins et al. Reference Adkins, Ivanov and Schekochihin2023). This implies a particular scaling of the heat flux with the square of the normalised parallel system size, viz.

(4.8) \begin{align} Q \propto \left (\frac {L_\parallel \sqrt {\sigma }}{L_T}\right )^2. \end{align}

Similarly, any intrinsic time scales in (4.1)–(4.3) (e.g. the outer-scale injection rate $\gamma ^{\textrm {o}}$ ) can be shown to be proportional to the inverse square of the normalised parallel box size. The numerical results for the (electron) heat flux $Q$ and any relevant rate $\gamma$ (e.g. $\gamma _{E}, \gamma _{\boldsymbol{k}}, \gamma ^{\textrm {o}}$ , etc.) can, therefore, be presented in terms of the following normalised quantities:

(4.9) \begin{align} \hat {Q} &\equiv \left (\frac {L_T}{L_\parallel \sqrt {\sigma }} \right )^2 \frac {Q}{(\rho _\perp /\rho _e) Q_{\textrm {gB}e}}, \end{align}

(4.10) \begin{align} \hat {\gamma } & \equiv \left (\frac {L_T}{L_\parallel \sqrt {\sigma }} \right )^{-2} \frac {\gamma }{\omega _\perp }, \\[6pt] \nonumber \end{align}

where $Q_{\textrm {gB}e} = n_{e} T_{e} v_{\textrm {th}e} (\rho _e/L_T)^2$ is the gyro-Bohm heat flux and

(4.11) \begin{align} \omega _\perp = \rho _e v_{\textrm {th}e} / 2\rho _\perp L_T \end{align}

is the value of the electron drift frequency at $k_y \rho _\perp = 1$ . A direct consequence of the scale invariance of (4.1)–(4.3) is that $\hat {Q}$ must be independent of $L_T$ and the perpendicular box size (in any direction, provided $L_x$ and $L_y$ are sufficiently large); therefore, it is a function of $\hat {\gamma }_{E}$ only. Note that all of the aforementioned scalings are valid only when the hyperviscous cut-off is far from the outer scale, i.e. when $\nu _\perp$ is small enough and so does not upset the scale invariance of the outer-scale quantities.

Table 1. A summary of the simulation parameters used in § 4.1. The simulation domain is taken to be ‘square’ with $L_x = L_y = L_\perp$ and $n_x = n_y = n_\perp$ , where $n_x$ , $n_y$ and $n_\parallel$ are the number of resolved (i.e. after dealiasing – see Appendix B) Fourier modes in the $x$ , $y$ and $z$ coordinates, respectively. The last column shows the maximum growth rate $\gamma _{ \rm{max}}$ normalised according to (4.10).

In the absence of flow shear, the nonlinear saturated state of (4.1)–(4.3) has been investigated extensively and exhibits a critically balanced local energy cascade (Adkins et al. Reference Adkins, Ivanov and Schekochihin2023). Therefore, (4.1)–(4.3) with $\gamma _{E}\neq 0$ should give rise to the kind of turbulence that is described by the theory laid out in § 3. Here, we show data from four sets of simulations where we varied $\gamma _{E}$ while keeping all other parameters fixed, as detailed in table 1.

Figure 4(a) shows the dependence of $\hat {Q}$ on $\hat {\gamma }_{E}$ for these simulations. The agreement with figure 3 is evident. The predicted scaling of the heat flux in the weak-shear regime (3.20) holds very well up to $\hat {\gamma }_{E} \approx 100$ . This is followed by a swift transition to the strong-shear scaling (3.27). Recall that the theory of § 3.2.1 predicts that the transition between the two regimes should occur at $\gamma _{E} \sim \gamma ^{\textrm {o}}(0)$ , where the energy-injection rate at the outer scale $\gamma ^{\textrm {o}}(0)$ is approximately the linear-instability growth rate $\gamma _{\boldsymbol{k}}$ at the outer scale. Using the outer-scale estimates presented in figure 4(b), we find that $k_{y}^{\textrm {o}}(0)\rho _\perp \approx 0.3$ . Solving the linear dispersion relation for (4.1)–(4.3) for Sim1 at $k_x = 0$ , $k_y\rho _\perp =0.35$ and substituting $k_\parallel L_T/\sqrt {\sigma } = 2\pi /50$ for the box-sized mode in the parallel direction, we obtain an approximation of the normalised outer-scale injection rate of $\hat {\gamma }^{\textrm {o}}(0) \sim \hat {\gamma }_{\boldsymbol{k}} \approx 70$ , consistent with the numerically observed transition at $\hat {\gamma }_{E} \approx 100$ .Footnote ⁹

Figure 4. (a) Time-averaged, saturated radial turbulent heat flux, normalised to its value at zero flow shear, as a function of normalised flow shear $\hat {\gamma }_{E}$ [normalised per (4.10)] for the sets of simulations detailed in table 1. The data from all four sets overlay due to the scale invariance of (4.1)–(4.3). The black dashed and dash-dotted lines show the theoretical predictions (3.20) and (3.27), respectively, where, for the former, the curve is plotted using $\hat {\gamma }_{\textrm {c}} \approx 39$ , found by fitting to the data presented here. The vertical black dotted line marks the approximate shearing rate $\hat {\gamma }_{E}\approx 100$ where the system transitions from the weak- to the strong-shear regime. The values of $\gamma _{E} \approx \gamma _{{\rm max}}$ are shown using vertical dotted lines of the same colour as the data points for each respective set of simulations. (b) Outer-scale wavenumbers $k_x^{\textrm {o}}(\gamma _{E})$ and $k_y^{\textrm {o}}(\gamma _{E})$ , defined as those that maximise (4.12) and (4.13), respectively, for the Sim1 set of simulations. The dashed line indicates a linear dependence on the flow shear, $k \propto \gamma _{E}$ . The left vertical dotted line is the same as in panel (a) and marks the location $\gamma _{E} \approx 100$ where the system transitions from the weak- to the strong-shear regime. In the former, $k_y^{\textrm {o}}$ is (approximately) pinned to $k_{y}^{\textrm {o}}(0)$ , but $k_x^{\textrm {o}}$ increases linearly with $\gamma _{E}$ . In the strong-shear regime, $k_x^{\textrm {o}} \sim k_y^{\textrm {o}} \propto \gamma _{E}$ . The right vertical dotted line indicates the value of flow shear that is equal to the largest growth rate $\gamma _{ \rm{max}}$ , where the outer scale $k_y^{\textrm {o}}$ reaches, at least approximately, the scale of the most unstable mode $k_{y, \rm{max}}\rho _\perp \approx 3.7$ . Note that, at low $\gamma _{E}$ , (4.12) is sometimes maximised at $k_x = 0$ . In those cases, represented by the hollow triangles, we have set $k_x^{\textrm {o}}\rho _\perp$ to the box scale $2\pi \rho _\perp /L_x \approx 0.063$ .

Figure 4(b) shows that $k_x^{\textrm {o}}(\gamma _{E})$ and $k_y^{\textrm {o}}(\gamma _{E})$ are consistent with the behaviour of the outer scale predicted in § 3.2. As $\gamma _{E}$ is increased, $k_y^{\textrm {o}}$ stays roughly constant until $k_x^{\textrm {o}}$ catches up with it, whereafter the system transitions into the strong-shear regime where $k_x^{\textrm {o}}$ and $k_y^{\textrm {o}}$ remain comparable and both increase linearly with $\gamma _{E}$ , as predicted in § 3.2.2. This behaviour persists until $k_y^{\textrm {o}}$ reaches (approximately) the poloidal wavenumber at which the linear growth rate is maximal. For $\gamma _{E}$ larger than $\gamma _{ \rm {max}}$ , the turbulence is completely suppressed. This behaviour is also visually confirmed by figure 5, where we show real-space snapshots from simulations with different flow shear.

The outer-scale wavenumbers shown in figure 4 are defined as those that maximise the (steady-state) poloidally and radially averaged heat fluxes, defined respectively as

(4.12) \begin{align} \langle Q\rangle _{y} (k_x) &\equiv \frac {3}{2} n_{e} T_{e} v_{\textrm {th}e} \sum _{k_y, k_\parallel } i k_y \rho _e{\varphi }^*_{\boldsymbol{k}} \frac {\delta {T}_{e, {\boldsymbol{k}}}}{T_{e}}, \end{align}

(4.13) \begin{align} \langle Q\rangle _{x} (k_y) &\equiv \frac {3}{2} n_{e} T_{e} v_{\textrm {th}e} \sum _{k_x, k_\parallel } i k_y \rho _e{\varphi }^*_{\boldsymbol{k}} \frac {\delta {T}_{e, {\boldsymbol{k}}}}{T_{e}}, \\[6pt] \nonumber \end{align}

where $\varphi _{\boldsymbol{k}}$ and $\delta {T}_{e, {\boldsymbol{k}}}$ are the three-dimensional Fourier amplitudes of the fields. The ratio of the heat flux at the transition to that at zero flow shear is found to be approximately $Q(0)/Q(\hat {\gamma }_{E}=100)\approx 13$ , which, according to (3.20), would require ${\mathcal {A}}^{\textrm {o}}(0) \approx 0.4$ , consistent with $k_x^{\textrm {o}}/k_y^{\textrm {o}} \approx 0.3$ , as seen in figure 4 at low values of flow shear. Given the remarkably good fit for $Q$ as a function of $\gamma _{E}$ , the small discrepancy is likely due to our estimates of $k_x^{\textrm {o}}$ and $k_y^{\textrm {o}}$ via (4.12) and (4.13) being an imperfect measure of the outer scale. For instance, (4.12) fails to produce a nonzero estimate for $k_x$ if $\gamma _{E}$ is too low or zero. Note that the scale invariance of the fluid system implies that the streamer aspect ratio ${\mathcal {A}}^{\textrm {o}}(0)$ of (4.1)–(4.3) is not a function of any of the parameters of the system; i.e. it is just an order-unity (albeit measurably, and consequentially, smaller than unity) constant, so it is not possible to perform a scan in ${\mathcal {A}}^{\textrm {o}}(0)$ .

Figure 5. Snapshots of $\varphi$ (top row) and $\delta T_e/T_{e}$ (bottom row) in the $(x, y)$ plane for Sim1 simulations with four different values of $\gamma _{E}$ , as specified above each column. For each snapshot, the amplitudes are normalised to lie in the range $[{-}1, 1]$ , with the values in this interval corresponding to colours between dark blue and dark red, respectively. The second column corresponds to the weak-shear regime (i) from figure 3(a), where the flow shear is too weak to influence the saturated state significantly. The third column also corresponds to the weak-shear regime, with $k_y^{\textrm {o}}(\gamma _{E})$ pinned to $k_{y}^{\textrm {o}}(0)$ but with $k_x^{\textrm {o}}(\gamma _{E})$ increased by the influence of the flow shear, which here clearly manifests itself as the tilting of the eddies. In this case, the structures have a similar size in $y$ to those in the first- and second-row panels, but a shorter length scale in $x$ due to being sheared. The last column shows the saturated state in the strong-shear regime (ii) of figure 3(a), where the flow shear has manifestly pushed the outer scale to much shorter wavelengths.

Figure 6. Radial localisation of turbulent perturbations at very large values of flow shear. Taken from a Sim4 simulation with $\hat {\gamma }_{E} = 540$ , which is just over the largest growth rate $\hat {\gamma }_{ \rm{max}} \approx 493$ . The simulation has achieved a steady state with time-averaged normalised heat flux $\hat {Q}(\gamma _{E})/\hat {Q}(0) \approx 4 \times 10^{-7}$ , which is why it is not visible in figure 4.

For values of $\gamma _{E}$ comparable to, or larger than, the maximum growth rate $\gamma _{ \rm{max}}$ , the turbulence is strongly suppressed, as expected. However, this occurs in a surprising and nontrivial way – turbulence becomes radially localised into disjoint turbulent patches (see figure 6). This localisation is reminiscent of the formation of coherent structures called ‘ferdinons’ in sheared ITG turbulence (van Wyk et al. Reference van Wyk, Highcock, Schekochihin, Roach, Field and Dorland2016, Reference van Wyk, Highcock, Field, Roach, Schekochihin, Parra and Dorland2017; Ivanov et al. Reference Ivanov, Schekochihin, Dorland, Field and Parra2020). Despite their qualitative similarities, the localised ETG structures reported here differ in a number of ways from the ITG ferdinons observed in similar fluid simulations by Ivanov et al. (Reference Ivanov, Schekochihin, Dorland, Field and Parra2020): the former are three-dimensional structures (i.e. require a nonzero $k_\parallel$ ) that drift radially in the absence of magnetic drifts, while the latter exist in two dimensions and depend on the poloidal magnetic drift for their radial motion. The radial localisation of sheared turbulence through the formation of coherent structures appears to be a universal phenomenon, whose detailed investigation is the subject of our ongoing work that falls outside the scope of this paper. Here, we note that, apart from the poloidally localised structures shown in figure 6 and reported by van Wyk et al. (Reference van Wyk, Highcock, Schekochihin, Roach, Field and Dorland2016) and Ivanov et al. (Reference Ivanov, Schekochihin, Dorland, Field and Parra2020), solitary travelling structures that are localised only radially and consist of a single poloidal harmonic (i.e. a single $k_y$ ) have been seen in various models of sheared plasma turbulence (Pringle, Mcmillan & Teaca Reference Pringle, Mcmillan and Teaca2017; McMillan et al. Reference McMillan, Jolliet, Tran, Villard, Bottino and Angelino2009, Reference McMillan, Pringle and Teaca2018; Zhou et al. Reference Zhou, Zhu and Dodin2019, Reference Zhou, Loureiro and Uzdensky2020). We find no such structures and it remains an open question exactly what is the determining factor for their appearance.

4.2. Gyrokinetic ITG turbulence

We now explore the validity of our theory for plasma turbulence in axisymmetric toroidal geometry. Using the GK code GENE, we performed numerical flux-tube simulations of ITG-driven turbulence with modified adiabatic electrons (Hammett et al. Reference Hammett, Beer, Dorland, Cowley and Smith1993) in a Cyclone-base-case (CBC) (Lin et al. Reference Lin, Hahm, Lee, Tang and Diamond1999; Dimits et al. Reference Dimits2000) geometry with normalised magnetic shear $\hat {s}=0.796$ , safety factor $q = 1.4$ and inverse aspect ratio $\epsilon = 0.18$ . We focus on two values of the ion-temperature gradient, $R/L_{T_i} = 10$ and $R/L_{T_i} = 14$ , box sizes corresponding to the smallest wavenumbers $k_{x, \rm{min}}\rho _i = 1.6\times 10^{-2}$ and $k_{y, \rm{min}}\rho _i = 6.25\times 10^{-3}$ , $n_x = 288$ radial modes, $n_y=256$ poloidal modes for $R/L_{T_i} = 10$ and $n_y=512$ for $R/L_{T_i} = 14$ , $n_z = 16$ parallel grid points, and velocity-space resolution $n_v = 32$ (parallel velocity), $n_\mu = 8$ (magnetic moment). We performed scans of the equilibrium perpendicular flow shear $\gamma _{E}$ , with parallel flow shear turned off. Note that the temperature gradients that we consider are significantly above the nominal CBC value of $R/L_{T_i} = 6.92$ . This is done to ensure that the system is strongly driven and far from any marginal states where other physics, e.g. zonal flows, could be setting the saturated fluctuation levels.Footnote ¹⁰

Figure 7. (a) Time-averaged, saturated-state, radial turbulent heat flux, normalised to its value at $\gamma _{E} = 0$ and (b) outer-scale poloidal wavenumber $k_y^{\textrm {o}}\rho _i$ as a function of the flow shear for two different values of the ion-temperature gradient (see § 4.2 for other relevant numerical parameters). The black dashed line corresponds to the trend $Q \propto \gamma _{E}^{-1}$ , while the blue and red dashed lines are linear fits for $k_y^{\textrm {o}}$ as a function of $\gamma _{E}$ for each temperature gradient. The vertical dotted lines correspond to $1.5\gamma _{\rm {max}}$ for each of the simulations. The flow shear is normalised to $a/c_s$ , where $a$ is the minor radius and $c_s$ is the ion sound speed.

Figure 7 shows the dependence of the heat flux and poloidal outer-scale wavenumber on the shearing rate. The strong-shear scaling $Q \propto \gamma _{E}^{-1}$ (3.27) is followed reasonably well for both values of the ion-temperature gradient. Also, the dependence of $k_y^{\textrm {o}}$ on $\gamma _{E}$ is approximately linear, as expected from (3.26). Note that (3.26) is only an asymptotic scaling and for $\gamma _{E}$ close to $\gamma ^{\textrm {o}}(0)$ , we have $k_y^{\textrm {o}}(\gamma _{E}) - k_{y}^{\textrm {o}}(0) \propto \gamma _{E}$ . To calculate $k_y^{\textrm {o}}$ for figure 7(b), we used the GK equivalent of (4.13), i.e. we approximated $k_y^{\textrm {o}}$ by the $k_y$ that has the largest contribution to the heat flux.

According to the theory presented in § 3.2.2, the strong-shear regime should end at $\gamma _{E} \sim \gamma _{{\rm max}}$ . Indeed, for both values of $R/L_{T_i}$ , the $Q \propto \gamma _{E}^{-1}$ dependence lasts roughly up to $\gamma _{E}/\gamma _{{\rm max}} \approx 1.5$ , after which the heat flux is sharply suppressed. Note the absence of the weak-shear scaling (3.20) – this is expected because the unsheared ( $\gamma _{E} = 0$ ) ITG turbulence in this case has ${\mathcal {A}}^{\textrm {o}}(0) \sim 1$ . Consequently, the heat flux is (approximately) constant in the weak-shear regime and the dependence of $Q$ on $\gamma _{E}$ resembles figure 3(b).

4.2.1. Bistability of high-shear states

Figure 8. Radial turbulent heat flux versus time for $R/L_{T_i} = 14$ and four different values of $\gamma _{E}$ , as labelled in the title of each panel. The blue lines are time traces from simulations initialised with small-amplitude noise, while the red ones represent simulations restarted from a saturated $\gamma _{E}=0$ run.

Figure 9. Snapshots of the electrostatic potential in the perpendicular plane from the (a) saturated high-transport and (b) low-transport states with $R/L_{T_i} = 14$ and $a\gamma _{E}/c_s = 1.3$ , whose heat-flux time traces are shown in blue and red, respectively, in the bottom right panel of figure 8. The perpendicular coordinates are normalised to the ion sound radius $\rho _s$ . Note that the aspect ratio of the panels corresponds to that of the simulation domain.

At large values of flow shear $\gamma _{E} \gt \gamma _{ \rm {max}}$ , i.e. beyond the strong-shear regime, we find that the system can saturate at (at least) two different levels of heat transport. Figure 8 shows the time traces of the turbulent heat flux for $R/L_{T_i} = 14$ and four different values of $\gamma _{E}$ , where, for each value of the flow shear, the simulations were initialised either with small-amplitude noise or with data from a saturated $\gamma _{E} = 0$ simulation. For $a\gamma _{E}/c_s \leqslant 1$ , the saturated state is found to be independent of the initial conditions. In contrast, a simulation with $a\gamma _{E}/c_s = 1.3$ initialised with a small-amplitude noise (corresponding to the rightmost point in figure 7) saturates with a time-averaged radial heat flux that is nearly two orders of magnitude smaller than that obtained by restarting it from an already saturated high-amplitude state. Similar bistability in gyrokinetic turbulence with mean flow shear has been reported by Christen et al. (Reference Christen, Barnes, Hardman and Schekochihin2022). The physics of this phenomenon, as observed in the simulations presented in figure 8, falls outside of the range of validity of the theory presented in § 3 because, at least in the case investigated here, it happens only at $\gamma _{E} \gt \gamma _{ \rm {max}}$ , where the assumption of a balance between the rates of shearing and energy injection (3.25) is not expected to hold.

Finally, we note that the appearance of coherent structures (whether of the poloidally localised kind or not, see the discussion at the end of § 4.1) can naturally lead to nonunique saturation. Indeed, if the saturated state is a collection of localised structures, then volume-averaged turbulent quantities, like the heat flux, are proportional to the number of structures in the simulation domain (van Wyk et al. Reference van Wyk, Highcock, Schekochihin, Roach, Field and Dorland2016). If these structures happen to be well localised and noninteracting (or only weakly interacting), then their number in the simulation domain is not necessarily uniquely determined, but can be a function of the initial conditions and/or of the perpendicular box size. The turbulent fluctuations in the low-transport saturated state displayed in figure 9 show signs of spatial localisation (as opposed to those in the high-transport state), even if not quite as obviously as they do in the ETG turbulence shown in figure 6.

5. Momentum transport

We have thus far focused on the influence of imposed flow shear on the turbulent heat transport. In sheared systems, another important quantity of interest is the turbulent transport of momentum, which is crucial for driving and maintaining equilibrium differential rotation in fusion experiments. Such rotation is typically achieved using neutral beams that deposit both energy and momentum into the plasma. The profiles of temperature and equilibrium flow are then determined by the turbulent heat diffusivity and turbulent viscosity, which are usually larger than their collisional counterparts. The simplifying choice of a purely perpendicular flow made in § 3.2 means that we are unable to describe all of the relevant physics: e.g. we are missing the effect of the parallel-velocity-gradient instability, which can be the main driver of transport at high $\gamma _{E}$ (Highcock et al. Reference Highcock, Barnes, Schekochihin, Parra, Roach and Cowley2010; Barnes et al. Reference Barnes, Parra, Highcock, Schekochihin, Cowley and Roach2011).Footnote ⁵ Nevertheless, the theory derived in § 3.2 makes numerically falsifiable predictions for the momentum transport of turbulence with imposed flow shear. Let us investigate them here.

Consider the fluid ETG model described in § 4.1, wherein we have imposed a mean flow in the poloidal $y$ direction that varies along the radial $x$ direction. This flow then allows the plasma to have nonzero radial transport of poloidal momentum, defined, similarly to (2.10), as

(5.1) \begin{align} \varPi \equiv \sum _s m_s \int \frac {\textrm {d}^3 \boldsymbol{r}}{V} \int \text {d}^{3} \boldsymbol{w} \ (\boldsymbol{v}_E \boldsymbol {\cdot } {\boldsymbol {\nabla }} x) (\boldsymbol{w} \boldsymbol {\cdot } {\boldsymbol {\nabla }} y) \delta f_{s}. \end{align}

It can be shown (see Appendix C) that, in the fluid model described in § 4.1, (5.1) becomes

(5.2) \begin{align} \varPi = -\frac {n_{e} T_{e} \rho _e^2 }{2} \int \frac {\textrm {d}^3 \boldsymbol{r}}{V} \frac {\partial \varphi }{\partial y} \frac {\partial }{\partial x} \left [\left (1 + \frac {Z T_{e}}{T_{i}}\right ) \varphi - \frac {\delta T_e}{T_{e}}\right ]. \end{align}

As expected, $\varPi$ is written as the sum of the Reynolds stress of the $\boldsymbol{E}\times \boldsymbol{B}$ flow and a diamagnetic stress (Smolyakov, Diamond & Medvedev Reference Smolyakov, Diamond and Medvedev2000; Ivanov et al. Reference Ivanov, Schekochihin, Dorland, Field and Parra2020; Sarazin et al. Reference Sarazin2021). Note that $\varPi$ itself is of no relevance to the dynamics of the fluctuations as it does not enter (4.1)–(4.3) at all, unlike $Q$ , which is responsible for the injection of energy, as per (2.7).

Just like we did for the heat flux $Q$ in § 3.1.1, we can estimate $\varPi$ as

(5.3) \begin{align} \varPi \sim n_{e} T_{e} k_{x, \rm{tilt}}^{\textrm {o}} k_y^{\textrm {o}} \rho _e^2 (\overline {\varphi }^{\textrm {o}})^2, \end{align}

where

(5.4) \begin{align} k_{x,\rm {tilt}}^{\textrm {o}}(\gamma _{E}) \equiv k_x^{\textrm {o}}(\gamma _{E}) - k_{x}^{\textrm {o}}(0) \sim k_y^{\textrm {o}}(\gamma _{E})\tau _{\mathrm{nl}}^{\textrm {o}}(\gamma _{E})\gamma _{E}. \end{align}

Using $k_{x, \rm{tilt}}^{\textrm {o}}$ instead of $k_x^{\textrm {o}}$ in (5.3) is necessary because, in the absence of flow shear, the GK equation (2.1) obeys the symmetry

(5.5) \begin{align} (x, y, z) \mapsto ({-}x, y, -z), \ w_\parallel \mapsto -w_\parallel, \ \phi \mapsto -\phi, \ h_s \mapsto -h_s, \end{align}

under which $\varPi \mapsto -\varPi$ and so the time-averaged $\varPi$ must vanish if $\gamma _{E} = 0$ (Parra et al. Reference Parra, Barnes and Peeters2011b; Sugama et al. Reference Sugama, Watanabe, Nunami and Nishimura2011; Fox et al. Reference Fox, van Wyk, Field, Ghim, Parra and Schekochihin2017). Thus, only the part of $k_x$ associated with the eddy tilting contributes to the momentum flux.Footnote ¹¹

Figure 10. A qualitative diagram of the momentum flux $\varPi$ versus flow shear $\gamma _{E}$ in the fluid ETG model (see figure 4 for a similar diagram for the heat flux $Q$ ). The indicated ratio between the plateau in the strong-shear regime and the peak of $\varPi$ at $\gamma _{E}=\gamma _{\textrm {c}}$ formally holds when ${\mathcal {A}}^{\textrm {o}}(0) \ll 1$ .

Using (3.2) and (3.6), we write (5.3) as

(5.6) \begin{align} \frac {\varPi }{n_{e} T_{e}} \sim \frac {\gamma _{E} \gamma ^{\textrm {o}}}{\Omega _e^2} \frac {1}{(k_x^{\textrm {o}} \rho _e)^2}\sim \frac {\gamma _{E} \gamma ^{\textrm {o}}}{\Omega _e^2} \frac {1}{({\mathcal {A}}^{\textrm {o}} k_y^{\textrm {o}} \rho _e)^2}. \end{align}

In the weak-shear regime (§ 3.2.1), where $\gamma _{E} \lt \gamma ^{\textrm {o}}(0)$ , the linear and nonlinear times at the outer scale are given by (3.15), and the fluctuation aspect ratio at the outer scale satisfies (3.18), we find

(5.7) \begin{align} \frac {\varPi (\gamma _{E})}{\varPi (0)} \propto \frac {\gamma _{E}}{(1 + \gamma _{E}/\gamma _{\textrm {c}})^2}. \end{align}

Note that if ${\mathcal {A}}^{\textrm {o}}(0)\sim 1$ , i.e. if the turbulence is isotropic in the absence of flow shear, then $\gamma _{\textrm {c}}\sim \gamma ^{\textrm {o}}(0)$ and $\varPi$ is (approximately) proportional to $\gamma _{E}$ in the weak-shear regime, unlike $Q$ , which would be constant in this case (see figure 3 b). In the strong-shear regime (§ 3.2.2), where $\gamma _{E} \gt \gamma ^{\textrm {o}}(0)$ , the fluctuations satisfy ${\mathcal {A}}^{\textrm {o}}(\gamma _{E}) \sim 1$ , and (3.25) and (3.26) hold, we find

(5.8) \begin{align} \varPi (\gamma _{E}) \propto \gamma _{E}^0, \end{align}

i.e. the momentum flux is independent of the imposed flow shear (see figure 10). Crucially, in both regimes,

(5.9) \begin{align} \frac {\varPi }{Q} \propto \gamma _{E}. \end{align}

This can also be seen directly from (3.8) and (5.6), which imply

(5.10) \begin{align} \frac {\varPi v_{\textrm {th}e}}{Q} \sim \frac {\gamma _{E}}{\gamma ^{\textrm {o}}} k_y^{\textrm {o}} \rho _e, \end{align}

and so (5.9) follows immediately if $\gamma ^{\textrm {o}} \propto k_y^{\textrm {o}}$ , as we have assumed throughout § 3.2. If we define the normalised turbulent viscosity $\nu$ and diffusivity $\chi$ as

(5.11) \begin{align} \nu \equiv \frac {\varPi }{n_{e} m_e \gamma _{E}}, \quad \chi \equiv \frac {Q}{n_{e} T_{e}L_T^{-1}}, \end{align}

then the turbulent Prandtl number is

(5.12) \begin{align} \text {Pr} \ = \frac {\nu }{\chi } = \frac {\varPi }{Q} \frac {v_{\textrm {th}e}^2}{2\gamma _{E} L_T}. \end{align}

According to (5.9), this is independent of $\gamma _{E}$ . This prediction is consistent with previous GK studies of turbulent transport in sheared turbulence (Highcock et al. Reference Highcock, Barnes, Schekochihin, Parra, Roach and Cowley2010; Barnes et al. Reference Barnes, Parra, Highcock, Schekochihin, Cowley and Roach2011).

Figure 11. (a) Radial flux of poloidal momentum in the fluid ETG model (5.2) as a function of $\gamma _{E}$ from simulation sets Sim1 and Sim2 (see table 1). The black dashed line is the best-fit line of the form (5.7) to the data up to $\gamma _{E}/\omega _\perp = 0.04$ (denoted by the vertical dotted line), where the system transitions from the weak- to the strong-shear regime (see also figure 4). (b) Prandtl number, defined as (5.12) for the same simulations as in panel (a).

Figure 11(a) shows the values of turbulent heat and momentum flux from the Sim1 and Sim2 sets of simulations, detailed in table 1. There is good agreement with the weak-shear scaling (5.7). According to (5.7), $\varPi$ should peak in the weak-shear regime at $\gamma _{E} = \gamma _{\textrm {c}} = {\mathcal {A}}^{\textrm {o}}(0)\gamma ^{\textrm {o}}(0)$ . For the Sim1 and Sim2 simulations, whose data are presented in figure 11, the transition to the strong-shear regime happens at $\gamma ^{\textrm {o}}(0)/\omega _\perp \approx 0.04$ , where $\omega _\perp$ is given by (4.11) (see also figure 4). As discussed in § 4.1, we find that, numerically, ${\mathcal {A}}^{\textrm {o}}(0) \approx 0.4$ , which gives $\gamma _{\textrm {c}}/\omega _\perp \approx 0.016$ , consistent with the observed peak in $\varPi$ at $\gamma _{E}/\omega _\perp \approx 0.01$ . However, the prediction that $\varPi$ should be constant in the strong-shear regime is not observed, likely due to finite-hyperviscosity effects and the finite extent of the inertial range. Our theory does not take into account the finite value of hyperviscosity necessary to achieve a saturated state numerically. Since $\varPi$ is a higher-order Fourier-space moment of the fluctuation spectrum than $Q$ , it is more sensitive to the spectrum at high $k$ , where hyperviscous effects matter. These effects become more important as $\gamma _{E}$ increases and the inertial range shortens (see also figure 12 in Appendix A). As figure 11(a) shows, the simulations with lower hyperviscosity (Sim2) have a weaker dependence of $\varPi$ on $\gamma _{E}$ in the strong-shear regime, consistent with the hypothesis that the hyperviscous cut-off is responsible for the discrepancy with the theoretical prediction (5.8). Figure 11(b) shows that the numerically measured Prandtl number is weakly dependent on $\gamma _{E}$ , varying only by approximately a factor of two over a range for $\gamma _{E}$ of nearly three orders of magnitude, before decreasing at large $\gamma _{E}$ where hyperviscous effects become important.

The nonmonotonic weak-shear dependence (5.7), which has a local maximum at $\gamma _{E} = \gamma _{\textrm {c}}$ , implies that, for the same value of $\varPi$ , two distinct values of $\gamma _{E}$ , and so of $Q$ , are possible. This suggests that sheared anisotropic turbulence may be prone to transport bifurcations similar to those discussed by Highcock et al. (Reference Highcock, Barnes, Schekochihin, Parra, Roach and Cowley2010, Reference Highcock, Barnes, Parra, Schekochihin, Roach and Cowley2011) and Parra et al. (Reference Parra, Barnes, Highcock, Schekochihin and Cowley2011a). Of course, our oversimplified model of electron-scale turbulence cannot be applied directly to any experimental studies. Nevertheless, our theory suggests that if experimentally relevant turbulence is dominated by streamers, transport bifurcations might exist.

6. Summary and discussion

Starting from the standard picture of turbulent saturation via a local energy cascade (§ 3.1), we have developed a theory for the effect of mean perpendicular flow shear on temperature-gradient-driven turbulence in fusion plasmas. As argued in § 3.2, it is meaningful to distinguish two different regimes depending on the ratio of the shearing rate $\gamma _{E}$ to the rate $\gamma ^{\textrm {o}}(0)$ of energy injection in the corresponding system with zero flow shear.

In the weak-shear regime, defined by $\gamma _{E} \lt \gamma ^{\textrm {o}}(0)$ , the poloidal outer scale and the energy-injection rate remain approximately the same as when $\gamma _{E} = 0$ , but the radial outer scale decreases with increasing $\gamma _{E}$ (3.18) due to the tilting of turbulent eddies by the shear. The extent to which the flow shear is able to suppress the turbulent transport in this regime is linked to the aspect ratio of the outer-scale fluctuations at zero flow shear, ${\mathcal {A}}^{\textrm {o}}(0) = k_{x}^{\textrm {o}}(0)/k_{y}^{\textrm {o}}(0)$ . We find that turbulence with ${\mathcal {A}}^{\textrm {o}}(0) \sim 1$ is largely unaffected by flow shear unless the shear is comparable to, or larger than, $\gamma ^{\textrm {o}}(0)$ . In contrast, heat transport in streamer-dominated turbulence with ${\mathcal {A}}^{\textrm {o}}(0) \ll 1$ , which is often encountered in fusion-relevant contexts, is shown to be strongly suppressed by flow shear even at $\gamma _{E} \ll \gamma ^{\textrm {o}}(0)$ as long as $\gamma _{E} \gtrsim \gamma _{\textrm {c}}$ , where the critical shearing rate $\gamma _{\textrm {c}}$ (3.19) is smaller than $\gamma ^{\textrm {o}}(0)$ by a factor of ${\mathcal {A}}^{\textrm {o}}(0)$ . This reflects the intuitive notion that radially elongated fluctuations should be more susceptible to sheared poloidal flows.

At $\gamma _{E} \gt \gamma ^{\textrm {o}}(0)$ , the system is in the strong-shear regime, where the outer scale is determined by the balance among the shearing, energy-injection and nonlinear mixing rates (3.25), and the outer-scale perpendicular wavenumbers grows proportionally to $\gamma _{E}$ . Turbulence in the strong-shear regime is found always to have $k_x^{\textrm {o}} \sim k_y^{\textrm {o}}$ , even if it was dominated by streamers at $\gamma _{E}=0$ . This is due to the shear-induced tilting of the turbulent fluctuations, which forces them to have ${\mathcal {A}}^{\textrm {o}}(\gamma _{E}) \sim 1$ in this regime.

Our theoretical predictions for the dependence of the radial turbulent heat flux on the rate of perpendicular flow shear, (3.20) and (3.27), are confirmed to hold over a range of four orders of magnitude for the flow shear in idealised fluid ETG simulations (§ 4.1). Additionally, GK flux-tube simulations of ITG turbulence have demonstrated that our theory applies to more realistic models of plasma turbulence as well (§ 4.2). These two models are paradigmatic cases of turbulence with ${\mathcal {A}}^{\textrm {o}}(0) \ll 1$ and ${\mathcal {A}}^{\textrm {o}}(0) \sim 1$ , respectively.

In addition to the heat-flux scalings, in § 5, we use our theory to predict the dependence of momentum flux on flow shear in the fluid ETG simulations. While not directly applicable to more realistic numerical simulations due to the restriction of purely perpendicular flow shear, our results suggest that streamer-dominated turbulent transport might exhibit transport bifurcations.

There exists a body of work on the suppression of turbulence by flow shear that is based on ‘decorrelation theories’ (Shaing, Crume & Houlberg Reference Shaing, Crume and Houlberg1990; Biglari, Diamond & Terry Reference Biglari, Diamond and Terry1990; Zhang & Mahajan Reference Zhang and Mahajan1992, Reference Zhang and Mahajan1993; Hatch et al. Reference Hatch, Hazeltine, Kotschenreuther and Mahajan2018) stemming from Dupree (Reference Dupree1972) or the equivalent derivation by Zhang & Mahajan (Reference Zhang and Mahajan1993). There are some parallels that can be drawn between these theories and our approach: e.g. the definition of the normalised flow shear $W$ by Zhang & Mahajan (Reference Zhang and Mahajan1992) and Hatch et al. (Reference Hatch, Hazeltine, Kotschenreuther and Mahajan2018) includes the crucial influence of the fluctuation aspect ratio and can be mapped to $\gamma _{E} / {\mathcal {A}}^{\textrm {o}}(\gamma _{E}) \gamma ^{\textrm {o}}(0)$ in our notation. However, the derivation of our theory, as presented in § 3.2, is not related to these decorrelation theories and produces different scalings for the dependence of the heat flux on $\gamma _{E}$ . For example, the functional dependence (3.20), borne out by the numerical results presented in § 4.1, is not reproduced by these theories. While the decorrelation theory by Zhang & Mahajan (Reference Zhang and Mahajan1992) has been claimed to be applicable to some GK simulations (Hatch et al. Reference Hatch, Hazeltine, Kotschenreuther and Mahajan2018), we are (at this stage) unable to draw any conclusions about the more general applicability either of decorrelation theories or of the theory presented in § 3.2.

One possible application of the theory of sheared streamer-dominated turbulence described in § 3.2.1 is in the description of multiscale interactions in magnetised-plasma turbulence. Gyrokinetic theory and numerical analysis point to the existence of linear instabilities and turbulent fluctuations at two disparate perpendicular scales, viz. the gyroradii of the main ion species $\rho _i$ and of the electrons $\rho _e$ , often referred to as the ‘ion’ and ‘electron’ (perpendicular) scales, respectively. By virtue of the difference in the ion and electron masses, these two scales are well separated, viz. $\rho _e/\rho _i \sim \sqrt {m_e/m_i} \ll 1$ . Furthermore, the GK ordering implies that the growth rate of electron-scale instabilities satisfies $\gamma _e \sim v_{\textrm {th}e}/L$ , while that of the ion-scale instabilities is $\gamma _i \sim v_{\textrm {th}i}/L$ , where $L$ is some appropriate measure of the size of the device. Therefore, the fluctuations at these two spatial scales also occur on well-separated time scales: $\gamma _i/\gamma _e \sim \sqrt{m_e / m_i} \ll 1$ . Nevertheless, numerically expensive GK simulations with enough resolution to span both ion and electron scales have shown that ion-scale turbulence can suppress electron-scale fluctuations (Candy et al. Reference Candy, Waltz, Fahey and Holland2007; Waltz, Candy & Fahey Reference Waltz, Candy and Fahey2007; Maeyama et al. Reference Maeyama, Idomura, Watanabe, Nakata, Yagi, Miyato, Ishizawa and Nunami2015; Howard et al., Reference Howard, Holland, White, Greenwald and Candy2016 Reference Howard, Holland, White, Greenwald and Candya ) and experimental results in support of this hypothesis have been reported (Howard et al. Reference Howard, Holland, White, Greenwald, Candy and Creely2016 Reference Howard, Holland, White, Greenwald, Candy and Creelyb ).

One possible mechanism for this suppression is that the turbulent ion-scale $\boldsymbol{E}\times \boldsymbol{B}$ flows shear away the electron-scale perturbations. However, if we assume that the electron-scale fluctuations satisfy ${\mathcal {A}}_{e} \sim 1$ , this could not happen: the time scale associated with the ion-scale $\boldsymbol{E}\times \boldsymbol{B}$ shear is $\gamma _i$ , while the growth rate of the electron-scale instabilities is $\gamma _e \gg \gamma _i$ , so, by the quench rule, $\gamma _i$ is too small to suppress the electron-scale turbulence.Footnote ¹²

The theory of weakly sheared turbulence proposed in § 3.2.1 offers a possible alternative explanation. Simulations of electron-scale turbulence indicate that turbulence at those scales is dominated by radially elongated streamers, and, according to our theory, the ion-scale $\boldsymbol{E}\times \boldsymbol{B}$ shear will be relevant for shearing these streamers if they have ${\mathcal {A}}_{e} \ll 1$ . In particular, if $k_{y, e} \rho _e \sim 1$ but $k_{x, e} \rho _i \sim 1$ , i.e. if the ETG streamers have a radial size comparable to the scale of ion fluctuations, then the critical shear rate $\gamma _{\textrm {c},e}$ for ETG turbulence (3.19) would be of the order of the ion-scale fluctuating $\boldsymbol{E}\times \boldsymbol{B}$ shear. Such an ordering would put the ion-scale shearing rate in the weak-shear regime for the ETG streamers, where (3.20) holds, with $\gamma _{E}$ now representing the value of the flow shear that electron-scale fluctuations experience from the ion-scale ones. Therefore, the influence of this flow shear would be non-negligible, despite the naïve argument that might lead one to presume otherwise. Numerical experiments have indeed shown that the inclusion of ion-scale zonal flows (flux-surface-constant, poloidal $\boldsymbol{E}\times \boldsymbol{B}$ shear flows) in ETG simulations can drastically reduce the levels of electron-scale turbulence and transport, even though the ion-scale zonal-flow shear is much smaller than the growth rate of electron-scale instabilities (Waltz et al. Reference Waltz, Candy and Fahey2007). Whether this mechanism for ITG suppression of electron-scale fluctuations is indeed realised in multiscale plasma turbulence appears to be a promising subject for future work.