1. Introduction
Quantifying vertical transport (of heat, chemical tracers and momentum) in the stably stratified regions of stars and gas giant planets is paramount to a better understanding of their structure and evolution, as the transport rates required to match observations are largely inconsistent with purely diffusive processes. The same challenge arises in the modelling of stably stratified regions of the Earth's atmosphere and oceans. The solution in both cases is to include turbulent transport in the models, but this in turn requires reliable parameterisations of small-scale fluxes as functions of the local properties of the fluid and of its large-scale flow (Munk Reference Munk1966; Pinsonneault Reference Pinsonneault1997; Ivey, Winters & Koseff Reference Ivey, Winters and Koseff2008; Gregg et al. Reference Gregg, D'Asaro, Riley and Kunze2018; Aerts, Mathis & Rogers Reference Aerts, Mathis and Rogers2019; Caulfield Reference Caulfield2021).
Relatively small-scale turbulent vertical mixing in a stably stratified region typically incurs an energetic cost to raise dense fluid parcels irreversibly and correspondingly lower buoyant parcels and thus can only be sustained on long time scales by mechanisms that tap into some larger-scale energy reservoir. Fluid instabilities comprise one broad class of such mechanisms, with different fluid instabilities accessing distinct sources of energy: for instance, baroclinic and double-diffusive instabilities draw from the available potential energy of the fluid, magnetic instabilities rely on the magnetic energy of a finite-amplitude large-scale field and shear instabilities of a large-scale mean flow, or of the motions associated with internal waves, can tap into the kinetic energy reservoirs of those flows. The nature of the turbulence driven by each of these instabilities differs substantially, so that it ‘still remains extremely difficult to say anything generic about [stratified] mixing’ (Caulfield Reference Caulfield2021).
In this work, we restrict our attention to modelling stratified turbulence in stellar and planetary interiors driven exclusively by shear. Astrophysical fluids generally have a low Prandtl number 
$Pr$, where
\begin{equation} Pr= \frac{\nu^*}{\kappa^*}, \end{equation}
where 
$\nu ^*$ is the (dimensional) kinematic viscosity of the fluid, and 
$\kappa ^*$ is the corresponding (dimensional) buoyancy diffusivity (e.g. thermal or compositional diffusivity, depending on the main source of stratification); see figure 7 of Garaud et al. (Reference Garaud, Medrano, Brown, Mankovich and Moore2015b). Here, and throughout this work, dimensional quantities are starred, while non-dimensional quantities are not. We ignore rotation, magnetic fields, externally driven large-scale internal waves and the presence of multiple components contributing to the flow buoyancy. Although these effects undoubtedly are important in most situations, the dynamics of purely shear-driven stratified turbulence in low 
$Pr$ fluids is still far from well understood, justifying the apparently narrow scope of this study.
Under these simplifying assumptions, the principal remaining source of turbulent vertical transport in our model is a vertical shear, expected to lead to ‘overturning’ horizontally aligned vortices. Balancing the potential energy cost and the kinetic energy gain of adiabatic turbulent eddies in a vertically sheared flow, Richardson (Reference Richardson1920) concluded that turbulence can only be sustained provided
\begin{equation} J = \frac{N^{*2}}{S^{*2}} \le O(1), \end{equation}
where 
$N^*$ is the typical value of the buoyancy frequency of the stratification, and 
$S^*$ is the typical vertical shearing rate of the flow. Condition (1.2) is known today as the Richardson criterion (where the non-dimensional parameter 
$J$ is commonly referred to as a gradient Richardson number, which can vary with the vertical position 
$z^*$). For linear normal mode disturbances in an inviscid, non-diffusive parallel shear flow, Miles (Reference Miles1961) and Howard (Reference Howard1961) formalised this argument to establish rigorously the necessary condition for linear instability to be that 
$J(z^*) < 1/4$ somewhere within the flow.
The Richardson criterion can be relaxed when the flow disturbances are not purely adiabatic, because the energy cost of vertical motions is reduced if the advected fluid parcels radiatively (Townsend Reference Townsend1958; Dudis Reference Dudis1973) or diffusively (Zahn Reference Zahn1974; Jones Reference Jones1977; Lignières, Califano & Mangeney Reference Lignières, Califano and Mangeney1999) adjust their buoyancy to that of the background stratification on a time scale comparable to, or shorter than, their turnover time scale. In the diffusive case, this adjustment is sufficiently fast when the eddy Péclet number
\begin{equation} Pe_{\ell} \equiv \frac{u_{\ell}^* \ell^*}{\kappa^*} \le O(1), \end{equation}
where 
$u_{\ell }^*$ is the characteristic velocity of an eddy of size 
$\ell ^*$. Zahn (Reference Zahn1974) heuristically argued that an appropriate criterion for the instability of diffusive eddies of size 
$\ell ^*$ in a vertical shear flow is
\begin{equation} J Pe_{\ell} \le O(1), \end{equation}
provided (1.3) holds. According to his criterion, it ought to be possible to maintain turbulence in strongly stratified shear flows (
$J \gg 1$) provided eddies are small enough to ensure that 
$Pe_{\ell } \le O(J^{-1})$.
The criterion (1.4) for these so-called ‘diffusive’ shear instabilities should not be used in geophysical fluids where 
$Pr \ge O(1)$ because the condition 
$Pe_{\ell } \le O(1)$ equivalently implies that 
$Re_{\ell } \equiv u_{\ell }^* \ell ^*/\nu ^* \le O(1)$, where 
$Re_{\ell }$ is the corresponding eddy Reynolds number. As shear instabilities must have a relatively large Reynolds number to develop (otherwise viscous energy losses are too great), the requirement that 
$Re_{\ell } \gg 1$ is incompatible with 
$Pe_{\ell } = Pr Re_{\ell } \le O(1)$ when 
$Pr \ge O(1)$. The situation is quite different in stellar and planetary interiors, as first noted by Spiegel & Zahn (Reference Spiegel and Zahn1970) and Zahn (Reference Zahn1974), because of their intrinsically small Prandtl number. With 
$Pr \ll 1$, and regardless of the size of the outer-scale Reynolds and Péclet numbers
\begin{equation} Re \equiv \frac{U^*L^* }{\nu^*} \quad \mbox{and}\quad Pe \equiv \frac{U^*L^*}{\kappa^*}, \end{equation}
(where 
$U^*$ and 
$L^*$ are the system-scale characteristic flow velocity and length scale, respectively), there is always an intermediate range of scales where turbulent eddies satisfy 
$Pe_{\ell } \ll 1 \ll Re_{\ell }$, namely, where diffusion dominates their dynamics while viscous forces remain negligible. This property naturally allows for the maintenance of diffusive stratified turbulence provided the shear is large enough for (1.4) to hold (Zahn Reference Zahn1974).
An important question is the source of the vertical shear. Naively, one might expect it to arise from the spatial gradients of a mean flow on large vertical scales, which are commonly observed in geophysical and astrophysical fluids. However, these large-scale ‘mean’ gradients generally have corresponding Richardson numbers that are much larger than unity, especially in stellar or planetary interiors (Garaud Reference Garaud2021). This naive expectation would also fail to explain observations of stratified turbulence in laboratory experiments (Ruddick, McDougall & Turner Reference Ruddick, McDougall and Turner1989; Park, Whitehead & Gnanadeskian Reference Park, Whitehead and Gnanadeskian1994; Holford & Linden Reference Holford and Linden1999; Oglethorpe, Caulfield & Woods Reference Oglethorpe, Caulfield and Woods2013) and numerical experiments (Jacobitz & Sarkar Reference Jacobitz and Sarkar2000; Basak & Sarkar Reference Basak and Sarkar2006; Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007; Maffioli & Davidson Reference Maffioli and Davidson2016; Lucas, Caulfield & Kerswell Reference Lucas, Caulfield and Kerswell2017; Zhou, Taylor & Caulfield Reference Zhou, Taylor and Caulfield2017; Zhou & Diamessis Reference Zhou and Diamessis2019; Cope, Garaud & Caulfield Reference Cope, Garaud and Caulfield2020; Garaud Reference Garaud2020; Yi & Koseff Reference Yi and Koseff2023), where the mean flow has no vertical shear.
Instead, it is the emergent, instantaneous, local vertical shear that is responsible for much of the vertical mixing in the flow (Riley & DeBruynkops Reference Riley and DeBruynkops2003). The emergent shear intermittently appears and disappears as different horizontal layers of the fluid move past each other, driven by primarily horizontal and highly anisotropic large-scale flows. Crucially, it has been shown that these vertically decoupled flow structures, which are sometimes called ‘pancake eddies’, emerge even when the horizontal flow is forced in a vertically invariant manner, as in the above-mentioned studies. Furthermore, in the stably stratified regions of the ocean, there is increasing evidence that the emergent shear is often marginally stable to vertical shear instabilities, with the minimum value of 
$J(z^*)\simeq 1/4$ (Smyth & Moum Reference Smyth and Moum2013), suggestive of ‘self-organised criticality’ (Salehipour, Peltier & Caulfield Reference Salehipour, Peltier and Caulfield2018; Smyth, Nash & Moum Reference Smyth, Nash and Moum2019; Mashayek, Caulfield & Alford Reference Mashayek, Caulfield and Alford2022). Similar results have recently been reported by Garaud, Khan & Brown (Reference Garaud, Khan and Brown2024b) in numerical simulations of strongly stratified turbulence at 
$Pr = 0.1$. As demonstrated by Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) at 
$Pr = O(1)$, and in this paper at low 
$Pr$, this self-organised criticality seems to be an inevitable outcome of the highly anisotropic nature of the layerwise large-scale flows in strongly stratified fluids.
In the spirit of these previous studies, we further narrow our focus to studying vertical mixing in primarily horizontal flows driven by a vertically invariant horizontal force, first at 
$Pr \ge O(1)$, where most of the research on this topic has been focused, and then at 
$Pr \ll 1$. First, however, we review the existing literature in some detail, in order to place our work in context.
1.1. Stratified turbulence forced by horizontal flows at 
$Pr \ge O(1)$
In studies of stratified turbulence forced primarily by horizontal flows, the effects of stratification are generally quantified using the outer-scale Froude number, defined as
\begin{equation} Fr = \frac{U^*}{N^*L^*}. \end{equation}
In this expression 
$U^*$ and 
$L^*$ are now specifically assumed to be the characteristic outer scales for the horizontal velocity and horizontal length scale and over that outer scale there is also a well-defined (implicitly constant) ‘background’ buoyancy frequency 
$N^*$. Ruddick et al. (Reference Ruddick, McDougall and Turner1989), Park et al. (Reference Park, Whitehead and Gnanadeskian1994) and Holford & Linden (Reference Holford and Linden1999) demonstrated experimentally that regular horizontal motions of vertical rods in a salt-stratified fluid (having 
$Pr = O(1000)$) with an initially constant stratification, characterised by a buoyancy frequency 
$N^*$, can cause substantial vertical mixing even when the flow is ‘strongly’ stratified, in the sense that 
$Fr \ll 1$. Furthermore, they often observed the formation of steps in the density profile, with a characteristic vertical height 
$H^* \propto U^* / N^*$. Oglethorpe et al. (Reference Oglethorpe, Caulfield and Woods2013) similarly found substantial mixing and the spontaneous formation and long-time survival of a ‘staircase’ of relatively well-mixed layers with depth again proportional to 
$U^* / N^*$, separated by significantly thinner interfaces of locally enhanced density gradient in stratified Taylor–Couette flow experiments.
Direct numerical simulations (DNS) are arguably more practical for investigation of 
$Pr = \mathit {O}(1)$ fluids (such as thermally stratified air and water, due to the inherent difficulties associated with conductive heat losses at the boundaries), and over the last 20 years significant progress in quantifying stratified turbulence in this parameter regime has been enabled by advances in supercomputing. In particular, access to the full three-dimensional structure of the velocity and buoyancy fields has provided new insights into the nature of turbulence at very large Reynolds numbers and very strong stratification, a regime characteristic of stratified turbulence in the ocean and atmosphere.
Pioneering work by Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007), for instance, has demonstrated that the eddy field is highly anisotropic in this regime, with a vertical eddy scale once again proportional to 
$U^*/N^*$. This emergent ‘layer’ length scale thus seems to be a universal property of stratified turbulence at both moderate and high 
$Pr$ at sufficiently high Reynolds number, as theoretically predicted and experimentally demonstrated by Billant & Chomaz (Reference Billant and Chomaz2000, Reference Billant and Chomaz2001) (see also Caulfield (Reference Caulfield2021) for a review).
Using the insight gained from their DNS data, Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007) proposed an anisotropic rescaling of the governing equations, introducing the small vertical scale 
$H^* = Fr L^* =U^*/N^*$ and taking the asymptotic limit 
$Fr\to 0$. Crucially, the horizontal scales are assumed to remain 
$O(L^*)$ in their work. Inspection of the rescaled equations reveals the fundamental role of the so-called buoyancy Reynolds number
\begin{equation} Re_b \equiv {Fr^2 Re =} \left(\frac{H^*}{L^*}\right)^2 Re, \end{equation}
which needs to be substantially greater than one for viscous effects to be negligible (see Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007), Bartello & Tobias (Reference Bartello and Tobias2013) and § 2 for further details). In that case, the characteristic vertical velocities and buoyancy fluctuations are predicted to scale as 
$W^* \propto Fr U^*$ and 
$B^* \propto Fr L^* N^{*2}=H^* N^{*2}$, respectively. This flow regime where 
$Re \gg 1$ and 
$Fr \ll 1$ such that 
$Re_b \equiv Fr^2 Re \gg 1$ was referred to as the ‘strongly stratified regime, sometimes abbreviated to stratified turbulence’ by Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007). The nomenclature ‘strongly stratified turbulence’ (SST) has become quite common in the fluid dynamics literature, although ‘stratified turbulence’ is often (particularly in the geosciences literature) used more broadly to refer to any disordered motions in a stratified environment, as opposed to this particular asymptotic regime.
Interestingly and importantly, Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007) presented an alternative expression for 
$Re_b$ in terms of the kinetic energy dissipation rate 
$\varepsilon ^*$, by assuming an effective definition for the horizontal outer scale 
$L^*$ such that
\begin{equation} \varepsilon^* \equiv \frac{U^{*3}}{L^*} \rightarrow Re_b \equiv \frac{\varepsilon^*}{\nu^* N^{*2}} . \end{equation}
This expression reveals 
$Re_b$ to be an ‘intensity’ or ‘activity’ parameter for the turbulence (Gibson Reference Gibson1980; Gargett, Osborn & Nasmyth Reference Gargett, Osborn and Nasmyth1984). Considered in this fashion, a particular attraction of the SST nomenclature becomes apparent, as it succinctly describes the necessary properties of the bulk flow, i.e. that the flow is ‘strongly stratified’ in the sense that 
$Fr \ll 1$ when using global or bulk measures for the outer length and velocity scales and the buoyancy frequency and yet the flow is also ‘turbulent’ in the sense that the disorder is sufficiently ‘intense’ or ‘active’ to be considered turbulent.
However, there is emerging evidence that this classical inertial scaling may only be applicable when the effects of buoyancy are sufficiently weak (Mashayek et al. Reference Mashayek, Caulfield and Alford2022), also perhaps calling into question whether it is appropriate to think of such flows as ‘strongly stratified’.
Fundamental to the arguments presented in Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007) is the notion that there is a single (large, outer) horizontal scale of significance (indeed they referred to the assumption leading to (1.8) as ‘the perhaps most important assumption’ for the forward cascade to small scales essential to their theory), with the implicit consequence that spatio-temporal variability on smaller horizontal scales does not exert a controlling influence on the dynamics. The supposition of a single dominant horizontal scale, however, is not entirely self-consistent. Indeed, Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007) demonstrated that an inevitable consequence of their scaling is that 
$J \simeq O(1)$ at least somewhere in the flow, thus implying that they are locally weakly stratified, and potentially unstable to fast, small-scale isotropic perturbations. The coexistence of strongly stratified (
$J\gg 1$) quiescent regions and weakly stratified (
$J \sim O(1)$) turbulent regions in 
$Fr \ll 1$ flows has led Falder, White & Caulfield (Reference Falder, White and Caulfield2016) to propose that this strongly stratified turbulent regime rather be referred to as the ‘layered anisotropic stratified turbulence’ (LAST) regime. This nomenclature has the twin advantage that it is specifically descriptive, as the flow velocities are inherently anisotropic and layered, and also avoids ambiguity concerning the actual sense in which the flow is ‘strongly’ stratified, in that a small bulk or global value of 
$Fr$ may well mask local regions in space and/or time where the stratification is sufficiently ‘weak’ to allow such unstable, fast small-scale isotropic perturbations to develop (see also Portwood et al. (Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016), for further discussion). That said, we acknowledge that ‘strongly stratified turbulence’ is more concise terminology and does not refer to flow properties which are emergent, such as the emergent anisotropy and inevitable spatio-temporal intermittency. Nevertheless, DNS results in Garaud et al. (Reference Garaud, Chini, Cope, Shah and Caulfield2024a) clearly do show two such emergent distinct flow behaviours as stratification increases (see their figure 3), i.e. turbulent patches dominated by small-scale isotropic motions and layer-like, quiescent flow outside these regions. It is this regime characterised by the coexistence of large and small scales that we refer to as the LAST regime.
Later, Riley & Lindborg (Reference Riley and Lindborg2012) employed local but anisotropic energy cascade arguments to deduce that 
$W^*$ actually peaks at small scales and is proportional to 
$Fr^{1/2} U^*$ rather than 
$Fr U^*$. Direct numerical simulations of SST by Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007), Augier, Billant & Chomaz (Reference Augier, Billant and Chomaz2015) and Maffioli & Davidson (Reference Maffioli and Davidson2016) do, indeed, exhibit spatio-temporally intermittent motions on small scales. Furthermore, the 
$W^* \propto Fr^{1/2} U^*$ scaling law was tentatively verified by Maffioli & Davidson (Reference Maffioli and Davidson2016).
More recently, Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) leveraged the scale separation observed in these DNS to propose a new multiscale asymptotic model of stratified turbulence at 
$Pr = O(1)$, using the concept of marginal stability to constrain the representative minimum values of the gradient Richardson number to be 
$O(1)$, thus determining key properties of ‘fast’ motions associated with the presumed breakdown of these local shear instabilities. Importantly, this model recovers the usual (and empirically observed) scaling law 
$H^* \propto Fr L^*$ for the vertical eddy scale, but predicts a characteristic vertical velocity scale 
$W^* \propto Fr^{1/2} U^*$ and a characteristic buoyancy scale 
$B^* \propto {H^* N^{*2}}$. This scaling for 
$W^*$ differs from the corresponding predictions of Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007) as a direct consequence of relaxing the assumption of a single horizontal length scale to allow for the idealised modelling of certain important aspects of the spatio-temporally intermittent shear instabilities. Moreover, the analysis by Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) also differs from that of Riley & Lindborg (Reference Riley and Lindborg2012) in three key aspects. First, in Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022), the scaling of 
$W^*$ arises from a self-consistent asymptotic theory rather than by invoking heuristic Kolmogorov turbulence arguments. Second, the scaling behaviour is associated with a spectrally non-local energy cascade from large to small horizontal scales, whereas Riley & Lindborg (Reference Riley and Lindborg2012) rely on spectrally local cascade arguments that presume 
$\alpha$ increases from its minimum value at the large scale, 
$\alpha \sim Fr$, to unity at the Ozmidov scale, 
$\alpha \sim 1$, where isotropy is recovered. Third, the vertical velocity scaling in Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) is realised at horizontal scales of 
$\mathit {O}(H^*) = \mathit {O}(Fr L^*)$, which is asymptotically larger than the Ozmidov scale 
$l_O^* \propto Fr^{3/2} L^*$ at which Riley & Lindborg (Reference Riley and Lindborg2012) conclude it arises. Data from recent DNS by Garaud et al. (Reference Garaud, Chini, Cope, Shah and Caulfield2024a) tentatively support the argument of Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) that the vertical kinetic energy spectrum of stratified turbulence peaks at a scale proportional to 
$Fr$.
1.2. Stratified turbulence forced by horizontal flows at $Pr \ll 1$
By contrast with geophysical applications, the study of stratified turbulence driven by primarily horizontal flows in stellar and planetary interiors is still in its infancy. Most of the early work on the topic is summarised in the seminal paper of Zahn (Reference Zahn1992), who had used arguments similar to those of Riley & Lindborg (Reference Riley and Lindborg2012) to estimate the vertical turbulent diffusivity 
$D^*$ of momentum, or of a passive scalar, in terms of the kinetic energy dissipation rate 
$\varepsilon ^*$ and of the local fluid properties 
$\kappa ^*$ and 
$N^*$:
\begin{equation} D^* \propto \left(\frac{\varepsilon^* \kappa^*}{N^{*2}}\right)^{1/2}. \end{equation}
To arrive at this conclusion, he made two assumptions: (i) the energy spectrum of the anisotropic horizontal scales of the flow follows a classical inertial scaling, analogously to (1.8) such that 
$\varepsilon ^* \propto u^{*3}_{\ell }/\ell ^*$; and (ii) the vertical diffusivity is primarily due to the largest eddies whose emergent vertical shear 
$S^*$ is unstable according to (1.4), where 
$S^* = u^*_{\ell }/\ell ^*$. Lignières (Reference Lignières2020) recently revisited Zahn's argument, interpreting the relevant scale 
$\ell ^*$ as a ‘diffusive’ or ‘modified’ Ozmidov scale and explicitly writing it as
\begin{equation} \ell^*_{OM} = \left(\frac{\kappa^* \varepsilon^{*1/3}}{N^{*2}}\right)^{3/8}. \end{equation}
Using critical balance theory, Skoutnev (Reference Skoutnev2023) then argued that it is possible to relate 
$\ell ^*_{OM}$ and 
$u^*(\ell ^*_{OM})$ explicitly to properties of the larger-scale flow (
$U^*$, 
$W^*$, 
$H^*$ and 
$L^*$), arriving at the conclusion that
\begin{equation} H^* \propto \left(\frac{U^*\kappa^*}{N^{*2} L^{*3}}\right)^{1/4} L^* = \left(\frac{Fr^2}{Pe}\right)^{1/4} L^*\quad \mbox{and}\quad W^* \propto \left(\frac{Fr^2}{Pe}\right)^{1/4} U^*, \end{equation}so
\begin{equation} D^* \propto H^* W^* \propto \left(\frac{Fr^2}{Pe}\right)^{1/2} U^* L^* = \left(\frac{U^{*3}\kappa^*}{L^*N^{*2}}\right)^{1/2}, \end{equation}
which is consistent with Zahn's estimate (1.9) provided one further assumes that 
$\varepsilon ^* \propto U^{*3}/L^*$.
As turbulence in stars cannot be directly observed, and experimentation with 
$Pr \ll 1$ fluids is particularly difficult, numerical simulations are best suited to provide insight into the low 
$Pr$ regime. Cope et al. (Reference Cope, Garaud and Caulfield2020) and Garaud (Reference Garaud2020) recently presented a series of DNS of stratified turbulence driven by horizontal shear at low 
$Pr$, focusing on flows where the outer-scale Péclet number 
$Pe$ is low and high, respectively. They found, as in the 
$Pr = O(1)$ simulations of Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007), that the turbulence becomes highly anisotropic as stratification increases (
$Fr$ decreases), and identified three distinct anisotropic regimes: (i) a fully turbulent regime, (ii) a regime where the turbulence is spatially and temporally intermittent (as viscosity begins to affect the flow) and (iii) a fully viscous regime (where viscosity completely controls the flow dynamics). In the fully turbulent and intermittent regimes, the horizontal motions contain well-separated large and small scales.
When 
$Pe \ll 1$, diffusion is important at all scales, regardless of the stratification. Cope et al. (Reference Cope, Garaud and Caulfield2020) found empirically from their DNS that the vertical size and characteristic velocity of eddies in the fully turbulent regime have the following apparent dependence on non-dimensional parameters:
\begin{equation} H^* \propto \left(\frac{Fr^2}{Pe}\right)^{1/3} L^*,\quad W^* \propto \left(\frac{Fr^2}{Pe}\right)^{1/6} U^*, \end{equation}
for which a rigorous self-consistent theory has yet to be developed. In particular, (1.12a,b) notably differs from the predictions of Skoutnev (Reference Skoutnev2023). This discrepancy between the DNS results of Cope et al. (Reference Cope, Garaud and Caulfield2020) and Skoutnev's theory needs to be explained and is one of the principal motivations for this paper. Interestingly, the prediction for the turbulent diffusivity 
$D^* \propto H^*W^*$ yields the same expression (1.11b) as in Zahn (Reference Zahn1992), Lignières (Reference Lignières2020) and Skoutnev (Reference Skoutnev2023), despite 
$H^*$ and 
$W^*$ satisfying different scaling laws.
It is important to note, however, that Zahn (Reference Zahn1992) did not distinguish between viscous and non-viscous regimes or diffusive and non-diffusive regimes, assuming instead that the turbulence can always develop on scales where diffusion is important but viscous effects are negligible. Yet, the results of Cope et al. (Reference Cope, Garaud and Caulfield2020) and Garaud (Reference Garaud2020) demonstrate that viscosity always eventually begins to affect the turbulence as the stratification increases. By analogy with Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007) and Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022), one may expect the buoyancy Reynolds number 
$Re_b$ to be the relevant bifurcation parameter describing the impact of viscosity, a result that will be formally established in this paper. Garaud (Reference Garaud2020) also demonstrated numerically the existence of a regime of stratified turbulence at 
$Pe \gg 1$ and 
$Pr \ll 1$ where diffusion is negligible. This is not surprising in hindsight, but was not anticipated by Zahn (Reference Zahn1992). She tentatively proposed that her data are consistent with 
$H^* \propto Fr^{2/3} L^*$, and 
$W^* \propto Fr^{2/3} U^*$ as long as viscosity remains negligible, which is perhaps surprising given that one might naively have expected to recover the geophysical regime scalings in that case. She acknowledged, however, that her simulations may not have been performed at a sufficiently high Reynolds number to be in a meaningful asymptotic regime yet. Recently, Garaud et al. (Reference Garaud, Chini, Cope, Shah and Caulfield2024a) revisited this dataset and demonstrated that the inferred 
$W^* \propto Fr^{2/3} U^*$ scaling is an artefact of the coexistence of turbulent patches where 
$W^* \propto Fr^{1/2}U^*$, and more quiescent regions where 
$W^* \propto Fr U^*$.
The results presented above illustrate the complexity of stratified turbulence at low 
$Pr$, and motivate the need for additional work to determine how many regimes exist, how salient flow properties scale with input parameters in each regime and, finally, where the regime boundaries lie in parameter space. In this paper, we therefore approach the problem systematically by adapting the multiscale asymptotic methodology developed by Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) for the 
$Pr = O(1)$ regime to the low 
$Pr$ regime. We begin in § 2 by laying out the model equations and boundary conditions and performing the anisotropic scaling analysis in Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007) on the 
$Pr \ll 1$ case, explicitly comparing it with the 
$Pr = O(1)$ case. This analysis reveals the crucial role of the buoyancy Péclet number (thus named by analogy with 
$Re_b$), viz. 
$Pe_b = Pr Re_b$. We then perform a slow–fast decomposition of the governing equations in § 3. At high 
$Pe_b$, we recover the results of Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) and highlight the reason why this multiscale approach leads to a different vertical velocity scaling from the single-scale one derived by Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007). We then present new results at low 
$Pe_b$. Our findings are discussed in § 4, where we use this analysis to partition the parameter space into various regimes of stratified turbulence at low 
$Pr$ and explicitly provide scaling relationships for the vertical length scale and velocity in each case. Implications and potential applications of our results are presented in § 5.
2. Governing equations and anisotropic scalings
Consider a three-dimensional, non-rotating, incompressible flow expressed in a Cartesian coordinate system where the vertical coordinate 
$z^*$ (with unit vector 
$\hat {\boldsymbol {e}}_z$) is anti-aligned with gravity 
$(\boldsymbol {g}^* = -g^* \hat {\boldsymbol {e}}_z)$. The fluid is stably stratified, with mean density 
$\rho _0^*$ and constant background buoyancy frequency 
$N^*$. Buoyancy perturbations away from this mean state are incorporated in accordance with the Boussinesq approximation (Spiegel & Veronis Reference Spiegel and Veronis1960). A vertically invariant divergence-free horizontal body force 
$\boldsymbol {F}^*_h$, which only varies slowly with time and space, is applied to drive a mean horizontally sheared flow in a domain of size 
$(L^*_x,L^*_y,L^*_z)$. The dimensional governing equations are
$$\begin{gather} \frac{\partial \boldsymbol{u^*}}{\partial t^*} + \boldsymbol{u^*} \boldsymbol{\cdot} \boldsymbol{\nabla}^* \boldsymbol{u^*} ={-} \frac{1}{\rho_0^*} \boldsymbol{\nabla}^* p^* + b^* \hat{\boldsymbol{e}}_z + \nu^* \nabla^{*2} \boldsymbol{u^*} + \frac{{\boldsymbol{F}^*_h}}{\rho_0^*}, \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla}^* \boldsymbol{\cdot} \boldsymbol{u^*} = 0, \end{gather}$$
$$\begin{gather}\frac{\partial b^*}{\partial t^*} + \boldsymbol{u^*} \boldsymbol{\cdot} \boldsymbol{\nabla}^* b^* + N^{*2} w^* = \kappa^* \nabla^{*2} b^*, \end{gather}$$
where 
$t^*$ is time, 
$\boldsymbol {\nabla ^*} = (\partial /\partial x^*,\partial /\partial y^*,\partial /\partial z^*)$ where 
$x^*$ and 
$y^*$ are horizontal spatial coordinates, 
$\boldsymbol {u^*} = (u^*,v^*,w^*)$ denotes the velocity field, 
$p^*$ the pressure and 
$b^*$ the buoyancy perturbation with respect to the background stratification. In accord with the Boussinesq approximation, the fluid has a constant kinematic viscosity 
$\nu ^*$ and constant diffusivity 
$\kappa ^*$.
2.1. Anisotropic single-scale equations for $Pr=\mathit {O}(1)$
In the limit of strong stratification, vertical displacements are energetically costly, and hence, as discussed in § 1, fluid motions become strongly anisotropic (with 
$u^*, v^* \gg w^*$). Therefore, it seems natural to non-dimensionalise the governing equations (2.1) anisotropically. Billant & Chomaz (Reference Billant and Chomaz2001) and Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007) have argued that the (dimensional) horizontal and vertical length scales of turbulent eddies ought to be 
$L^*$ and 
$H^*$, respectively, with an aspect ratio
\begin{equation} \alpha \equiv \frac{H^*}{L^*} \ll 1, \end{equation}
where the dependence of 
$\alpha$ on the stratification (and other governing parameters) is determined from the following asymptotic analysis. We then introduce a new vertical coordinate 
$\zeta ^*$ such that
\begin{equation} z^* = \alpha \zeta^*. \end{equation}
Consequently, if the horizontal velocity scale is 
$U^*$, then the vertical velocity scale must be 
$\alpha U^*$ to respect the divergence free condition without overly restricting the allowable types of flows. Time should be scaled by the turnover time of the horizontal eddies 
$L^*/U^*$ and pressure by 
$\rho _0^* U^{*2}$. To ensure that the nonlinear terms balance the forcing, 
$U^*= (F^*_0 L^*/\rho _0^*)^{1/2}$, where 
$F_0^*$ is the characteristic forcing amplitude. Finally, the buoyancy scale is chosen to be 
$H^* N^{*2}$ to enable the vertical advection of the background stratification to enter the buoyancy equation at leading order and to balance the horizontal advection of the buoyancy fluctuations.
Denoting the horizontal components of the velocity as 
$\boldsymbol {u}_{\perp } = (u,v)$, and similarly the horizontal gradient as 
$\boldsymbol {\nabla }_{\perp } = (\partial /\partial x,\partial /\partial y)$, the dimensionless system is given by
$$\begin{gather} \frac{\partial \boldsymbol{u}_{{\perp}}}{\partial t} + \left( \boldsymbol{u}_{{\perp}} \boldsymbol{\cdot}\boldsymbol{\nabla}_\perp \right) \boldsymbol{u}_{{\perp}} + w \frac{\partial \boldsymbol{u}_{{\perp}}}{\partial \zeta} ={-} \boldsymbol{\nabla}_{{\perp}} p + \frac{1}{Re \alpha^2} \left( \alpha^2 \nabla_{{\perp}}^2 \boldsymbol{u}_{{\perp}} + \frac{\partial ^2 \boldsymbol{u}_{{\perp}}}{\partial \zeta^2} \right) + {\boldsymbol{F}_h}, \end{gather}$$
$$\begin{gather}\frac{\partial w}{\partial t} + \left( \boldsymbol{u}_{{\perp}} \boldsymbol{\cdot}\boldsymbol{\nabla}_\perp \right) w + w \frac{\partial w}{\partial \zeta} ={-}\frac{1}{\alpha^2} \frac{\partial p}{\partial \zeta} + \frac{b}{Fr^2} + \frac{1}{Re\alpha^2} \left(\alpha^2 \nabla_{{\perp}}^2 w + \frac{\partial ^2 w}{\partial \zeta^2} \right), \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla}_{{\perp}} \boldsymbol{\cdot} \boldsymbol{u}_{{\perp}} + \frac{\partial w}{\partial \zeta} = 0, \end{gather}$$
$$\begin{gather}\frac{\partial b}{\partial t} + \left( \boldsymbol{u}_{{\perp}} \boldsymbol{\cdot}\boldsymbol{\nabla}_\perp \right) b + w \frac{\partial b}{\partial \zeta} + w = \frac{1}{Pe \alpha^2} \left( \alpha^2 \nabla_{{\perp}}^2b + \frac{\partial ^2 b}{\partial \zeta^2} \right), \end{gather}$$where all variables are now non-dimensional and
\begin{equation} Re = \frac{U^*L^*}{\nu^*},\quad Fr = \frac{U^*}{N^*L^*},\quad Pe = \frac{U^*L^*}{\kappa^*} = Pr Re , \end{equation}are the usual Reynolds, Froude and Péclet numbers based on the characteristic horizontal scales of the flow.
When the stratification is strong, 
$Fr \ll 1$, and the buoyancy term in the vertical component of the momentum equation is unbalanced unless it is compensated by the vertical pressure gradient. In other words, the flow anisotropy implies that hydrostatic balance must be satisfied at lowest order, which then requires
\begin{equation} \alpha = Fr, \end{equation}
implying that the characteristic vertical velocity 
$W^* = FrU^*$ while the characteristic vertical scale of the flow 
$H^* = U^*/N^*$. As discussed in § 1, this scaling for 
$H^*$ is well established, and has been observed in several laboratory and numerical experiments (Holford & Linden Reference Holford and Linden1999; Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007; Oglethorpe et al. Reference Oglethorpe, Caulfield and Woods2013).
Keeping only the lowest-order terms in an asymptotic expansion of (2.4) in 
$\alpha = Fr \ll 1$ yields
$$\begin{gather} \frac{\partial \boldsymbol{u}_{{\perp}}}{\partial t} + (\boldsymbol{u}_{{\perp}} \boldsymbol{\cdot}\boldsymbol{\nabla}_\perp) \boldsymbol{u}_{{\perp}} + w \frac{\partial \boldsymbol{u}_{{\perp}}}{\partial \zeta} ={-} \boldsymbol{\nabla}_{{\perp}} p + \frac{1}{Re_b} \frac{\partial ^2 \boldsymbol{u}_{{\perp}}}{\partial \zeta^2} + {\boldsymbol{F}_h}, \end{gather}$$
$$\begin{gather}\frac{\partial p}{\partial \zeta} = b,\quad \boldsymbol{\nabla}_{{\perp}} \boldsymbol{\cdot} \boldsymbol{u}_{{\perp}} + \frac{\partial w}{\partial \zeta} = 0, \end{gather}$$
$$\begin{gather}\frac{\partial b}{\partial t} + (\boldsymbol{u}_{{\perp}} \boldsymbol{\cdot}\boldsymbol{\nabla}_\perp) b + w \frac{\partial b}{\partial \zeta} + w = \frac{1}{Pe_b } \frac{\partial ^2 b}{\partial \zeta^2}, \end{gather}$$where
\begin{equation} Re_b = \alpha^2 Re , \end{equation}is the usually defined buoyancy Reynolds number and
\begin{equation} Pe_b = Pr Re_b = \alpha^2 Pe , \end{equation}
is the corresponding buoyancy Péclet number. Note that 
$Pr = O(1)$, the condition 
$Re_b \ge O(1)$, which is necessary for viscous effects to be small or negligible, implies that 
$Pe_b = Pr Re_b \ge O(1)$. As such, the effects of buoyancy diffusion are also small or negligible. The set of equations (2.7), which will be referred to as the ‘anisotropic single-scale high-
$Pe_b$’ equations hereafter (‘single-scale high-
$Pe_b$’ equations for short), recovers the results of Billant & Chomaz (Reference Billant and Chomaz2001) in the inviscid and non-diffusive limit and of Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007) when viscous and diffusive effects are incorporated.
2.2. Anisotropic single-scale equations for $Pr \ll 1$
When the Prandtl number is asymptotically small, it is possible to have a regime where 
$Pe_b \ll 1 \le Re_b$. In this extreme limit, the diffusion term on the right-hand side of (2.4d) becomes unbalanced, unless the buoyancy field itself is much smaller than anticipated by the scaling 
$H^*N^{*2}$ used in the previous section. The strongly diffusive limit has in fact already been studied by Lignières (Reference Lignières1999), who showed that the buoyancy equation asymptotically reduces to a balance between the vertical advection of the background stratification and the diffusion of the buoyancy fluctuations. In dimensional terms, we therefore expect 
$N^{*2} w^* \simeq \kappa ^* \nabla ^2 b^*$. With this in mind, we let 
$b = Pe_b \hat b$ in (2.4a)–(2.4d), and anticipate that 
$\hat b =O(1)$ (this is equivalent to scaling the dimensional buoyancy by 
$\alpha ^3 N^{*2} U^*L^{*2}/\kappa ^*$ instead of 
$H^*N^{*2}$). The resulting dimensionless system becomes
$$\begin{gather} \frac{\partial \boldsymbol{u}_{{\perp}}}{\partial t} + (\boldsymbol{u}_{{\perp}} \boldsymbol{\cdot}\boldsymbol{\nabla}_\perp) \boldsymbol{u}_{{\perp}} + w \frac{\partial \boldsymbol{u}_{{\perp}}}{\partial \zeta} ={-} \boldsymbol{\nabla}_{{\perp}} p + \frac{1}{Re_b} \left( \alpha^2 \nabla_{{\perp}}^2 \boldsymbol{u}_{{\perp}} + \frac{\partial ^2 \boldsymbol{u}_{{\perp}}}{\partial \zeta^2} \right) + {\boldsymbol{F}_h}, \end{gather}$$
$$\begin{gather}\frac{\partial w}{\partial t} + (\boldsymbol{u}_{{\perp}} \boldsymbol{\cdot}\boldsymbol{\nabla}_\perp) w + w \frac{\partial w}{\partial \zeta} ={-}\frac{1}{\alpha^2} \frac{\partial p}{\partial \zeta} + \frac{Pe_b}{Fr^2} \hat b + \frac{1}{Re_b} \left( \alpha^2 \nabla_{{\perp}}^2 w + \frac{\partial ^2 w}{\partial \zeta^2} \right), \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla}_{{\perp}} \boldsymbol{\cdot} \boldsymbol{u}_{{\perp}} + \frac{\partial w}{\partial \zeta} = 0, \end{gather}$$
$$\begin{gather}\frac{\partial \hat b}{\partial t} + (\boldsymbol{u}_{{\perp}} \boldsymbol{\cdot}\boldsymbol{\nabla}_\perp) \hat b + w \frac{\partial \hat b}{\partial \zeta} + \frac{1}{Pe_b} w = \frac{1}{Pe_b} \left(\alpha^2 \nabla_{{\perp}}^2 \hat b + \frac{\partial ^2 \hat b}{\partial \zeta^2}\right). \end{gather}$$
For sufficiently strong stratification (i.e. 
$Fr^2 \ll Pe_b$), the vertical component of the momentum equation must again be in hydrostatic balance, which requires
\begin{equation} \alpha^2 = \frac{Fr^2}{Pe_b} = \frac{Fr^2}{\alpha^2 Pe} \rightarrow \alpha = Fr_M,\quad \mbox{where } Fr_M = \left(\frac{Fr^2}{Pe}\right)^{1/4} \end{equation}
is a modified Froude number (see Lignières Reference Lignières2020; Skoutnev Reference Skoutnev2023). We therefore find that the characteristic vertical length scale should be 
$H^* = (Fr^2/Pe)^{1/4} L^*$ and the characteristic vertical velocity scale should be 
$W^* = (Fr^2/Pe)^{1/4} U^*$, recovering the results of Skoutnev (Reference Skoutnev2023) albeit using a different argument.
In the limit 
$Fr_M \ll 1$ and 
$Pe_b \ll 1$, keeping only the lowest-order terms in (2.10a)–(2.10d) yields
$$\begin{gather} \frac{\partial \boldsymbol{u}_{{\perp}}}{\partial t} + (\boldsymbol{u}_{{\perp}} \boldsymbol{\cdot}\boldsymbol{\nabla}_\perp) \boldsymbol{u}_{{\perp}} + w \frac{\partial \boldsymbol{u}_{{\perp}}}{\partial \zeta} ={-} \boldsymbol{\nabla}_{{\perp}} p + \frac{1}{Re_b} \frac{\partial ^2 \boldsymbol{u}_{{\perp}}}{\partial \zeta^2} + {\boldsymbol{F}_h}, \end{gather}$$
$$\begin{gather}\frac{\partial p}{\partial \zeta} = \hat b ,\quad \boldsymbol{\nabla}_{{\perp}} \boldsymbol{\cdot} \boldsymbol{u}_{{\perp}} + \frac{\partial w}{\partial \zeta} = 0,\quad w = \frac{\partial ^2 \hat b}{\partial \zeta^2} . \end{gather}$$
These scaling laws and governing equations are the low 
$Pe_b$ analogues of (2.7). In what follows, we therefore refer to them as the ‘anisotropic single-scale low-
$Pe_b$’ equations (‘single-scale low-
$Pe_b$’ equations for short).
2.3. Motivation for considering multiscale dynamics
For strongly stratified flows, the reduced, anisotropic single-scale equations given by (2.7) for 
$Pe_b \ge O(1)$ and (2.12) for 
$Pe_b \ll 1$ assume, by construction, that horizontal scales are large, while the vertical scale is small (e.g. in the 
$Pe_b \ge O(1)$ case, for strongly stratified flows, 
$\alpha \sim Fr$ necessarily implies 
$\alpha \ll 1$). As such, these equations describe the dynamics of weakly coupled ‘pancake’ vortices or horizontal meanders of the mean flow. They cannot, however, capture the small-scale turbulence that is expected to develop from shear instabilities between these layerwise horizontal motions (Chini et al. Reference Chini, Michel, Julien, Rocha and Caulfield2022). Yet, these instabilities are ubiquitous when 
$Re_b$ is sufficiently large, and have been observed in laboratory experiments (Augier et al. Reference Augier, Billant, Negretti and Chomaz2014) and inferred from oceanographic in situ measurements (Falder et al. Reference Falder, White and Caulfield2016) as well as in DNS at 
$Pr \gtrsim O(1)$ (Riley & DeBruynkops Reference Riley and DeBruynkops2003; Waite Reference Waite2011; Augier, Chomaz & Billant Reference Augier, Chomaz and Billant2012).
Indeed, recent DNS of stratified turbulence has confirmed the presence and importance of flows on small horizontal scales. Small-scale features are clearly visible in figure 16 of Maffioli & Davidson (Reference Maffioli and Davidson2016), figure 6 of Cope et al. (Reference Cope, Garaud and Caulfield2020), figure 1 of Garaud (Reference Garaud2020) and figure 3 of Garaud et al. (Reference Garaud, Chini, Cope, Shah and Caulfield2024a). Furthermore, Cope et al. (Reference Cope, Garaud and Caulfield2020) empirically find that 
$H^* \propto (Fr^2/Pe)^{1/3} L^*$ instead of 
$H^* \propto (Fr^2/Pe)^{1/4} L^*$ and that 
$W^* \propto (Fr^2/Pe)^{1/6} U^*$ instead of 
$W^* \propto (Fr^2/Pe)^{1/4} U^*$. These results thus call for further work to clarify the conditions (if any) where the single-scale low-
$Pe_b$ equations (2.12) are relevant for modelling stratified turbulence, at least at large buoyancy Reynolds number.
As discussed in § 1, Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) have recently argued in the context of geophysical stratified turbulence that one must take into account the fast, small horizontal scales and study their (marginally stable) interaction with the slow, larger-scale anisotropic flow to obtain a more complete and more accurate model of stratified turbulence. Accordingly, we now hypothesise that a multiscale dynamics is the missing physics required to obtain a self-consistent model of low 
$Pr$ stratified turbulence and propose to extend their work to the low 
$Pr$ limit. The next section first outlines the work of Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) for pedagogical clarity, then extends the analysis to the low 
$Pr$ regime.
3. Multiscale models for stratified turbulence
We consider the same model set-up as introduced in § 2. Here, however, we make no assumption about the amplitude of the vertical motions when non-dimensionalising the governing equations, and instead allow the appropriate scaling to emerge naturally from the analysis. Accordingly, we non-dimensionalise the vertical velocity by 
$U^*$ and correspondingly choose the buoyancy scale to be 
$L^*N^{*2}$. Then, the dimensionless system is
$$\begin{gather} \frac{\partial \boldsymbol{u}_{{\perp}}}{\partial t} + (\boldsymbol{u}_{{\perp}} \boldsymbol{\cdot}\boldsymbol{\nabla}_\perp) \boldsymbol{u}_{{\perp}} + \frac{w}{\alpha} \frac{\partial \boldsymbol{u}_{{\perp}}}{\partial \zeta} ={-} \boldsymbol{\nabla}_{{\perp}} p + \frac{1}{Re_b} \left( \alpha^2 \nabla_{{\perp}}^2 \boldsymbol{u}_{{\perp}} + \frac{\partial ^2\boldsymbol{u}_{{\perp}}}{\partial \zeta^2} \right) + {\boldsymbol{F}_h}, \end{gather}$$
$$\begin{gather}\frac{\partial w}{\partial t} + (\boldsymbol{u}_{{\perp}} \boldsymbol{\cdot}\boldsymbol{\nabla}_\perp) w + \frac{w}{\alpha} \frac{\partial w}{\partial \zeta} ={-}\frac{1}{\alpha} \frac{\partial p}{\partial \zeta} + \frac{b}{Fr^2} + \frac{1}{Re_b} \left(\alpha^2 \nabla_{{\perp}}^2 w + \frac{\partial ^2 w}{\partial \zeta^2} \right), \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla}_{{\perp}} \boldsymbol{\cdot} \boldsymbol{u}_{{\perp}} + \frac{1}{\alpha} \frac{\partial w}{\partial \zeta} = 0, \end{gather}$$
$$\begin{gather}\frac{\partial b}{\partial t} + (\boldsymbol{u}_{{\perp}} \boldsymbol{\cdot}\boldsymbol{\nabla}_\perp) b + \frac{w}{\alpha} \frac{\partial b}{\partial \zeta} + w = \frac{1}{Pe_b} \left( \alpha^2 \nabla_{{\perp}}^2b + \frac{\partial ^2 b}{\partial \zeta^2} \right). \end{gather}$$
These equations are the starting point for our analysis. We assume that 
$Re_b \ge O(1)$, and that 
$Fr$ is sufficiently small to ensure that 
$\alpha \ll 1$, but make no other a priori assumption on the size of 
$Pe_b$ at this stage.
3.1. Slow–fast decomposition
We now perform a multiscale expansion of the system in the limit of small aspect ratio 
$\alpha$. Following Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022), we assume the existence of two distinct sets of horizontal length scales: the original large scales that are 
$O(1)$ in the chosen non-dimensionalisation, as well as much smaller horizontal scales that are 
$O(\alpha )$. With that choice, small-scale fluid motions are isotropic by construction. We further assume that the flow has two distinct time scales: a slow time scale associated with the turnover of the large horizontal eddies, as before, and a fast time scale inversely related to the emergent vertical shearing rate of the large-scale mean flow, 
$U^*/H^*$, that develops from the vertically decoupled horizontal flows. In practice, we thus define the slow and fast horizontal coordinates as 
$\boldsymbol {x}_{s} = \boldsymbol {x}_{\perp }$ and 
$\boldsymbol {x}_{f} = \boldsymbol {x}_{s}/\alpha$, respectively (henceforth, the subscript 
$f$ denotes fast and 
$s$ denotes slow). Similarly, we split time into slow and fast variables, such that 
$t_f = t_s/\alpha$ where 
$t_s = t$. Consequently, the partial derivatives with respect to time and to the horizontal variables become
\begin{equation} \frac{\partial }{\partial t} = \frac{1}{\alpha} \frac{\partial }{\partial t_f} + \frac{\partial }{\partial t_s},\quad \boldsymbol{\nabla}_{{\perp}} = \frac{1}{\alpha} \boldsymbol{\nabla}_{f} + \boldsymbol{\nabla}_{s}. \end{equation}
All dependent variables (collectively referred to as 
$q$) are now assumed to be functions of both fast and slow length and time scales: 
$q = q(\boldsymbol {x}_{f},\boldsymbol {x}_{s},\zeta,t_f,t_s;\alpha )$.
Assuming the fast and slow scales are sufficiently separated, Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) then define a fast-averaging operator 
$\overline {({\cdot })}$, such that
\begin{equation} \bar{q}(\boldsymbol{x}_{s},\zeta,t_s;\alpha) = \frac{1}{l_x l_y T}\int_\mathcal{T} \int_\mathcal{D} q (\boldsymbol{x}_{f},\boldsymbol{x}_{s},\zeta,t_f,t_s;\alpha)\,\mathrm{d}\kern0.07em \boldsymbol{x}_{f}\,\mathrm{d} t_f, \end{equation}
where 
$\mathcal {D}$ is a horizontal domain of size 
$(l_x \times l_y)$ centred on 
$\boldsymbol {x}_s$ where 
$\alpha \ll l_x,l_y \ll 1$, and 
$\mathcal {T}$ is a temporal domain of size 
$\alpha \ll T \ll 1$ centred on 
$t_s$. With this definition, 
$\bar {q}$ depends on slow variables only. Each quantity 
$q$ can then be split into a slowly varying field 
$\bar q$ and a rapidly fluctuating component 
$q' = q-\bar {q}$, which implies that the fast average of the fluctuation field must vanish, i.e. 
$\overline {q'} = 0$. Note that 
$q'$ itself can still depend on the slow length and time scales.
We first substitute the expressions (3.2a,b) for 
$\boldsymbol {\nabla }_{\perp }$ and 
$\partial /\partial t$ into (3.1), and split each field 
$q$ into 
$\bar q + q'$. We then take the fast average of each of the four governing equations to obtain the evolution equations for the mean flow, then subtract the mean from the total to obtain the evolution equations for the fluctuations.
Starting with the continuity equation, we have
\begin{equation}
\frac{1}{\alpha} \boldsymbol{\nabla}_{f} \boldsymbol{\cdot}
\boldsymbol{u'}_{\!\!{\perp}} + \boldsymbol{\nabla}_{s}
\boldsymbol{\cdot} \boldsymbol{\bar u}_{{\perp}} +
\boldsymbol{\nabla}_{s} \boldsymbol{\cdot}
\boldsymbol{u'}_{\!\!{\perp}} + \frac{1}{\alpha} \frac{\partial
\bar w}{\partial \zeta} + \frac{1}{\alpha} \frac{\partial
w'}{\partial \zeta} = 0,
\end{equation}
whose fast average reveals that
\begin{equation} \boldsymbol{\nabla}_{s} \boldsymbol{\cdot} \boldsymbol{\bar u}_{{\perp}} + \frac{1}{\alpha} \frac{\partial \bar w}{\partial \zeta} = 0 \end{equation}for the mean flow and
\begin{equation} \frac{1}{\alpha} \boldsymbol{\nabla}_{f} \boldsymbol{\cdot} \boldsymbol{u'}_{\!\!{\perp}} + \boldsymbol{\nabla}_{s}\boldsymbol{\cdot} \boldsymbol{u'}_{\!\!{\perp}} + \frac{1}{\alpha} \frac{\partial w'}{\partial \zeta} = 0 \end{equation}for the perturbations.
A similar procedure for the horizontal momentum equation yields
\begin{align} & \frac{\partial \bar{\boldsymbol{u}}_\perp}{\partial t_s} + \bar{\boldsymbol{u}}_\perp \boldsymbol{\cdot}\boldsymbol{\nabla}_s \bar{\boldsymbol{u}}_\perp{+} \frac{\bar{w}}{\alpha} \frac{\partial \bar{\boldsymbol{u}}_\perp}{\partial \zeta} + \frac{1}{\alpha} \left( \overline{\boldsymbol{u}_\perp' \boldsymbol{\cdot}\boldsymbol{\nabla}_f \boldsymbol{u}_\perp' + w' \frac{\partial \boldsymbol{u}_\perp'}{\partial \zeta}}\right) + \overline{\boldsymbol{u}_\perp' \boldsymbol{\cdot}\boldsymbol{\nabla}_s \boldsymbol{u}_\perp'} \nonumber\\ &\quad ={-} \boldsymbol{\nabla}_s \bar{p} + \frac{1}{Re_b} \left(\frac{\partial ^2 \bar{\boldsymbol{u}}_\perp}{\partial \zeta^2} + \alpha^2 \nabla_s^2 \bar{\boldsymbol{u}}_\perp \right) + {\bar{\boldsymbol{F}}_h}, \end{align}
for the mean flow (where the magnitude of the horizontal force is 
$O(1)$ by construction), and
\begin{align} & \frac{1}{\alpha}
\frac{\partial \boldsymbol{u}_\perp'}{\partial t_f}
+\frac{\partial \boldsymbol{u}_\perp'}{\partial t_s} +
\frac{1}{\alpha} \bar{\boldsymbol{u}}_\perp
\boldsymbol{\cdot}\boldsymbol{\nabla}_f
\boldsymbol{u}_\perp' + \bar{\boldsymbol{u}}_\perp
\boldsymbol{\cdot}\boldsymbol{\nabla}_s
\boldsymbol{u}_\perp' + \boldsymbol{u}_\perp'
\boldsymbol{\cdot}\boldsymbol{\nabla}_s
\bar{\boldsymbol{u}}_\perp \nonumber\\ &\qquad
+\frac{1}{\alpha} \boldsymbol{u}_\perp'
\boldsymbol{\cdot}\boldsymbol{\nabla}_f
\boldsymbol{u}_\perp' + \boldsymbol{u}_\perp'
\boldsymbol{\cdot}\boldsymbol{\nabla}_s
\boldsymbol{u}_\perp' + \frac{1}{\alpha} \left( \bar{w}
\frac{\partial \boldsymbol{u}'_\perp}{\partial \zeta} + w'
\frac{\partial \bar{\boldsymbol{u}}_\perp}{\partial \zeta}
+ w' \frac{\partial \boldsymbol{u}'_\perp}{\partial
\zeta}\right) \nonumber\\ &\quad ={-}\frac{1}{\alpha}
\boldsymbol{\nabla}_f p' - \boldsymbol{\nabla}_s p' +
\frac{1}{Re_b} \left(\nabla^2_f \boldsymbol{u}'_\perp + \frac{\partial ^2 \boldsymbol{u}'_\perp }{\partial
\zeta^2}\right) \nonumber\\ &\qquad + \frac{\alpha^2}{Re_b}
\nabla^2_s \boldsymbol{u}_\perp' + \frac{2 \alpha}{Re_b}
\boldsymbol{\nabla}_s
\boldsymbol{\cdot}\boldsymbol{\nabla}_f
\boldsymbol{u}_\perp' + \frac{1}{\alpha}
\left(\overline{\boldsymbol{u}_\perp'
\boldsymbol{\cdot}\boldsymbol{\nabla}_f
\boldsymbol{u}_\perp' + w' \frac{\partial
\boldsymbol{u}_\perp'}{\partial \zeta}}\right) +
\overline{\boldsymbol{u}_\perp'
\boldsymbol{\cdot}\boldsymbol{\nabla}_s
\boldsymbol{u}_\perp'},
\end{align}
for the fluctuations. Note that 
$\bar {\boldsymbol F}_h = {\boldsymbol F}_h$ because the forcing is slowly varying in time and space by assumption.
The mean buoyancy equation is
\begin{align} & \frac{\partial \bar{b}}{\partial t_s} + \bar{\boldsymbol{u}}_\perp \boldsymbol{\cdot}\boldsymbol{\nabla}_s \bar{b} + \frac{\bar{w}}{\alpha} \frac{\partial \bar{b}}{\partial \zeta} + \frac{1}{\alpha} \left(\overline{\boldsymbol{u}_\perp' \boldsymbol{\cdot}\boldsymbol{\nabla}_f b' + w' \frac{\partial b'}{\partial \zeta} } \right) + \overline{\boldsymbol{u}_\perp' \boldsymbol{\cdot}\boldsymbol{\nabla}_s b' } + \bar{w} \nonumber\\ &\quad = \frac{1}{Pe_b}\left(\frac{\partial ^2 \bar{b}}{\partial \zeta^2} + \alpha^2 \nabla_{s}^2 \bar{b} \right), \end{align}while the corresponding fluctuation equation becomes
\begin{align} & \frac{1}{\alpha} \frac{\partial b'}{\partial t_f} + \frac{\partial b'}{\partial t_s} + \frac{1}{\alpha} \bar{\boldsymbol{u}}_\perp \boldsymbol{\cdot}\boldsymbol{\nabla}_f b' + \boldsymbol{u}_\perp' \boldsymbol{\cdot}\boldsymbol{\nabla}_s \bar{b} + \frac{1}{\alpha} \boldsymbol{u}_\perp' \boldsymbol{\cdot}\boldsymbol{\nabla}_f b' + \boldsymbol{u}_\perp' \boldsymbol{\cdot}\boldsymbol{\nabla}_s b' \nonumber\\ &\quad + \frac{1}{\alpha} \left( w' \frac{\partial \bar{b}}{\partial \zeta} + \bar{w} \frac{\partial b'}{\partial \zeta} + w' \frac{\partial b'}{\partial \zeta} \right) + w' = \frac{1}{Pe_b} \left( \nabla^2_f b' + \frac{\partial ^2 b' }{\partial \zeta^2} \right) \nonumber\\ &\quad + \frac{\alpha^2}{Pe_b} \nabla^2_s b' + \frac{2 \alpha}{Pe_b} \boldsymbol{\nabla}_s \boldsymbol{\cdot} \boldsymbol{\nabla}_f b' + \frac{1}{\alpha} \left(\overline{\boldsymbol{u}_\perp' \boldsymbol{\cdot}\boldsymbol{\nabla}_f b' + w' \frac{\partial b'}{\partial \zeta}}\right) + \overline{\boldsymbol{u}_\perp' \boldsymbol{\cdot}\boldsymbol{\nabla}_s b'}. \end{align}Finally, the mean vertical momentum equation is
\begin{align} & \frac{\partial \bar{w}}{\partial t_s} + \bar{\boldsymbol{u}}_\perp \boldsymbol{\cdot}\boldsymbol{\nabla}_s \bar{w} + \frac{\bar{w} }{\alpha} \frac{\partial \bar{w}}{\partial \zeta} + \frac{1}{\alpha} \left(\overline{\boldsymbol{u}_\perp' \boldsymbol{\cdot} \boldsymbol{\nabla}_f w' + w' \frac{\partial w'}{\partial \zeta}}\right) + \overline{\boldsymbol{u}_\perp' \boldsymbol{\cdot}\boldsymbol{\nabla}_s w' } \nonumber\\ &\quad ={-} \frac{1}{\alpha} \frac{\partial \bar p}{\partial \zeta} + \frac{\bar b}{Fr^2} + \frac{1}{Re_b} \left(\frac{\partial ^2 \bar{w}}{\partial \zeta^2} + \alpha^2 \nabla_s^2 \bar{w}\right), \end{align}while the fluctuations satisfy
\begin{align} & \frac{1}{\alpha} \frac{\partial w'}{\partial t_f} + \frac{\partial w'}{\partial t_s} + \frac{1}{\alpha} \bar{\boldsymbol{u}}_\perp \boldsymbol{\cdot}\boldsymbol{\nabla}_f w' + \boldsymbol{u}_\perp' \boldsymbol{\cdot}\boldsymbol{\nabla}_s \bar{w} + \frac{1}{\alpha} \boldsymbol{u}_\perp' \boldsymbol{\cdot}\boldsymbol{\nabla}_f w' + \boldsymbol{u}_\perp' \boldsymbol{\cdot}\boldsymbol{\nabla}_s w' \nonumber\\ &\quad + \frac{1}{\alpha} \left( w' \frac{\partial \bar{w}}{\partial \zeta} + \bar{w} \frac{\partial w'}{\partial \zeta} + w' \frac{\partial w'}{\partial \zeta} \right) ={-}\frac{1}{\alpha} \frac{\partial p'}{\partial \zeta} + \frac{b'}{Fr^2} + \frac{1}{Re_b} \left(\nabla_f^2 w' + \frac{\partial ^2 w'}{\partial \zeta^2} \right)\nonumber\\ &\quad + \frac{\alpha^2}{Re_b} \nabla^2_s w' + \frac{2 \alpha}{Re_b} \boldsymbol{\nabla}_s \boldsymbol{\cdot} \boldsymbol{\nabla}_f w' + \frac{1}{\alpha} \left(\overline{\boldsymbol{u}_\perp' \boldsymbol{\cdot}\boldsymbol{\nabla}_f w' + w' \frac{\partial w'}{\partial \zeta}} \right) + \overline{\boldsymbol{u}_\perp' \boldsymbol{\cdot}\boldsymbol{\nabla}_s w'}. \end{align} We see, as noted by Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022), that the effective Reynolds and Péclet numbers of the fluctuation equations are 
$Re_b/\alpha$ and 
$Pe_b / \alpha$, respectively, which implies that the fluctuations are formally much less viscous and less diffusive than the mean. This perhaps counterintuitive conclusion is a direct consequence of the flow anisotropy, which enables the fluctuations to evolve on a faster time scale than the mean.
3.2. Multiscale model at $Pe_b \ge O(1)$
We begin by summarising the steps taken by Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) to derive a reduced multiscale model for stratified turbulence at 
$Re_b,Pe_b \ge O(1)$, as much of the analysis proves to be similar at low Prandtl number. As in that work, we posit the following asymptotic expansions:
\begin{align} [b,p,\boldsymbol{u}_{{\perp}},w] \sim [b_0, p_0, \boldsymbol{u}_{{\perp} 0}, w_0] + \alpha^{1/2}[b_1, p_1, \boldsymbol{u}_{{\perp} 1}, w_1] + \alpha [b_2, p_2, \boldsymbol{u}_{{\perp} 2}, w_2] + \cdots . \end{align}
The expansions start at 
$O(1)$ to reflect the expectation that the dominant contributions to the pressure and the horizontal velocity arise on large horizontal scales. Although the expansions for 
$b$ and 
$w$ also start at 
$\mathit {O}(1)$, we show below that 
$w_0$ and 
$b_0$ both vanish when 
$\alpha \rightarrow 0$. Finally, the expansions proceed as asymptotic series in 
$\alpha ^{1/2}$ following the results of Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022). They demonstrated that, because the small-scale fluctuations are isotropic, 
$\boldsymbol {u'}$ and 
$w'$ are necessarily of the same order. Inspection of the mean horizontal momentum equation then immediately reveals that both fields need to be 
$O(\alpha ^{1/2})$ to ensure that the Reynolds stresses feed back on 
$\boldsymbol {u}_{\perp 0}$ at leading order (see below for further details.)
With these choices, we substitute (3.9) into the equations with multiscale derivatives presented above, analysing in turn the continuity equation (3.5), the horizontal component of the momentum equation (3.6), the buoyancy equation (3.7) and finally, the vertical component of the momentum equation (3.8). At each step, we match terms at leading order to infer their sizes and respective evolution equations, and thus derive a reduced model for the flow.
Considering first the mean continuity equation (3.5a), we see that 
$\partial _\zeta \bar {w}_0=\partial _\zeta \bar {w}_1 = 0$, implying 
$\bar w_0 = \bar w_1 = 0$ to suppress unphysical ‘elevator modes’ from our model. Consequently,
\begin{equation} \boldsymbol{\nabla}_{s} \boldsymbol{\cdot} \boldsymbol{\bar u}_{{\perp} 0} + \frac{\partial \bar w_2}{\partial \zeta} = 0. \end{equation}In the fluctuation continuity equation (3.5b), the second term is clearly much smaller than the first and can therefore be neglected from the leading-order set of dominant terms. Substituting the asymptotic series (3.9) we then see that
\begin{equation} \boldsymbol{\nabla}_{f} \boldsymbol{\cdot} \boldsymbol{u'}_{\!\!{\perp} i} + \frac{\partial w'_i}{\partial \zeta} = 0, \end{equation}
for 
$i = 0, 1$ (for larger values of 
$i$, the slow derivative of 
$\boldsymbol {u'}_{\!\!\perp i-2}$ should be taken into account).
We now turn to the mean horizontal component of the momentum equation (3.6a). The Reynolds stresses must be 
$O(1)$ to feed back on the mean flow, so 
$\boldsymbol {u'}$ and 
$w'$ must both be 
$O(\alpha ^{1/2})$. Hence, 
$\boldsymbol {u}_{\perp 0}' = \boldsymbol {0}$ and 
$w'_{0} = 0$, yielding 
$\boldsymbol {u}_{\perp 0} = \bar {\boldsymbol {u}}_{\perp 0}$ and that 
$w_{0} = 0$, since 
$\bar {w}_0=0$, too. (See Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022), for a more detailed discussion of why 
$\boldsymbol {u}_{\perp 0}' = \boldsymbol {0}$.) In the horizontal momentum equation for the fluctuations, the fast dynamics takes place at 
$O(\alpha ^{-1/2}$) since 
$\boldsymbol {u'} = O(\alpha ^{1/2})$. Ensuring that pressure is a leading-order effect implies that 
$p'_0 = 0$, so 
$p_0 = \bar p_0$, as expected. Many of the remaining terms are formally higher order, including all fluctuation–fluctuation interactions, which are 
$O(1)$ or smaller. Of the nonlinear terms, the only ones that contribute at leading order are quasilinear: 
$\alpha ^{-1} \bar {\boldsymbol {u}}_\perp \boldsymbol {\cdot }\boldsymbol {\nabla }_f \boldsymbol {u}'_\perp$ and 
$\alpha ^{-1} w' \partial _z \bar {\boldsymbol {u}}_\perp$. Therefore, after substituting (3.9) into (3.6) and retaining only the leading terms, we obtain
\begin{align} & \frac{\partial \bar{\boldsymbol{u}}_{{\perp} 0} }{\partial t_s} + \bar{\boldsymbol{u}}_{{\perp} 0} \boldsymbol{\cdot} \boldsymbol{\nabla}_s \bar{\boldsymbol{u}}_{{\perp} 0} + \bar{w}_2 \frac{\partial \bar{\boldsymbol{u}}_{{\perp} 0}}{\partial \zeta} ={-} \boldsymbol{\nabla}_{s} \bar{p}_0 - \frac{\partial }{\partial \zeta} (\overline{w_{1}'\boldsymbol{u}_{{\perp} 1}'}) + \frac{1}{Re_b}\frac{\partial ^2 \bar{\boldsymbol{u}}_{{\perp} 0}}{\partial \zeta^2} + {\bar{\boldsymbol{F}}_h}, \end{align}
\begin{align} &\frac{\partial \boldsymbol{u}_{{\perp} 1}'}{\partial t_f} + \bar{\boldsymbol{u}}_{{\perp} 0} \boldsymbol{\cdot}\boldsymbol{\nabla}_{f} \boldsymbol{u}_{{\perp} 1}' + w_{1}' \frac{\partial \bar{\boldsymbol{u}}_{{\perp} 0}}{\partial \zeta} ={-} \boldsymbol{\nabla}_{f} p_1' + \frac{\alpha}{Re_b} \left(\nabla^2_f \boldsymbol{u}'_{{\perp} 1} + \frac{\partial ^2 \boldsymbol{u}'_{{\perp} 1}}{\partial \zeta^2} \right), \end{align}where the formally higher-order Laplacian term has been retained to regularise the fluctuation equation.
Next, we examine the buoyancy equation, which reveals the sizes of 
$\bar b$ and 
$b'$. In the mean equation (3.7a), 
$\bar w = O(\alpha )$, and the buoyancy flux term is formally also of that order (see Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022), and also below). Thus, 
$\bar b = O(\alpha )$ as well, as long as 
$Pe_b \ge O(1)$, which is implicit in the regime considered in this section, implying 
$\bar b_0 = \bar b_1 = 0$. The size of the buoyancy flux can be confirmed by inspection of the buoyancy fluctuation equation (3.7b). To be in a regime in which the stratification impacts the turbulent motions, the 
$O(b'/\alpha )$ fast dynamics must be of the same order as the advection of the background stratification, which is 
$O(w') = O(\alpha ^{1/2})$. This ordering implies that 
$b' = O(\alpha ^{3/2})$, so 
$b'_0 = b'_1 = b'_2 = 0$. Consequently, 
$b_0 = b_1 =0$, and 
$b_2 = \bar b_2$. Fluctuation-fluctuation interactions in this equation are again formally higher order, and the only remaining nonlinearities are quasilinear. At leading order, using (3.9) in (3.7) shows that the mean and fluctuation buoyancy equations are
$$\begin{gather} \frac{\partial \bar{b}_2}{\partial t_s} + \bar{\boldsymbol{u}}_{{\perp} 0} \boldsymbol{\cdot}\boldsymbol{\nabla}_s \bar{b}_2 + \bar{w}_2 \frac{\partial \bar{b}_2}{\partial \zeta} + \bar{w}_2 ={-} \frac{\partial }{\partial \zeta} (\overline{w_{1}' b_3'}) + \frac{1}{Pe_b} \frac{\partial ^2\bar{b}_2}{\partial \zeta^2}, \end{gather}$$
$$\begin{gather}\frac{\partial b_{3}'}{\partial t_f} + \bar{\boldsymbol{u}}_{{\perp} 0} \boldsymbol{\cdot}\boldsymbol{\nabla}_{f} b_{3}' + w_{1}' \frac{\partial \bar{b}_2}{\partial \zeta} + w_{1}' = \frac{\alpha}{Pe_b} \left(\nabla_f^2 b_{3}' + \frac{\partial ^2 b'_3}{\partial \zeta^2}\right), \end{gather}$$where, as before, the diffusion term in the fluctuation equation has been retained for regularisation.
Finally, we identify the leading-order terms in the vertical component of the momentum equation. In the mean equation (3.8a), a hydrostatic leading-order balance requires 
$\alpha = Fr$, as in the anisotropic single-scale high-
$Pe_b$ model described in § 2.1. This specification, when applied to the fluctuation equation (3.8b), is quite satisfactory as it implies that the fluctuating buoyancy force 
$b'/Fr^2$ arises at leading order. One implication of the 
$\alpha = Fr$ scaling is that the largest horizontal scale to exhibit an isotropic dynamics is predicted to be 
$\mathit {O}(U^*/N^*)$, which is much larger than the Ozmidov scale (since 
$\ell _O=\mathit {O}(Fr^{1/2}U^*/N^*)$). Physically, this prediction is consistent with the preferential excitation of stratified shear instabilities having horizontal length scales comparable to the layer thickness, as demonstrated theoretically in Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) and numerically in Augier et al. (Reference Augier, Billant and Chomaz2015). Nevertheless, the Ozmidov scale retains its significance as the largest horizontal scale to be (largely) unaffected by buoyancy. It can also easily be shown that, once again, fluctuation–fluctuation interactions are negligible at leading order, so the fast dynamics is quasilinear. Substituting (3.9) into (3.8) and using all of the information available, we obtain the leading-order mean and fluctuating components of the vertical momentum equation:
$$\begin{gather} \frac{\partial \bar{p}_0}{\partial \zeta} = \bar b_2, \end{gather}$$
$$\begin{gather}\frac{\partial w_{1}'}{\partial t_f} + \bar{\boldsymbol{u}}_{{\perp} 0} \boldsymbol{\cdot}\boldsymbol{\nabla}_{f} w_{1}' ={-} \frac{\partial p_1'}{\partial \zeta} + b_{3}' + \frac{\alpha}{Re_b} \left(\nabla_f ^2 w_{1}' + \frac{\partial ^2 w'_1}{\partial \zeta^2} \right), \end{gather}$$where the formally higher-order viscous term has been retained to regularise the fluctuation equation.
The system of equations formed by (3.10) (with 
$i=1$), (3.11), (3.12) and (3.13), is equivalent to (2.28)–(2.35) in Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022), the only differences arising from a different choice of non-dimensionalisation. Furthermore, we see that the mean flow equations (3.11a), (3.12a) and (3.13a) recover the single-scale high-
$Pe_b$ equations in (2.7) in the absence of fluctuations. The multiscale equations are fully closed and therefore self-consistent. Crucially, the resulting system is quasilinear, and can be solved efficiently by appealing to (plausible and empirically supported) marginal stability arguments for the fluctuations as discussed by Michel & Chini (Reference Michel and Chini2019) and Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022).
These equations are valid whenever 
$\alpha \ll 1$ (equivalently, 
$Fr \ll 1$ since 
$\alpha = Fr$), and the buoyancy Reynolds and Péclet numbers are both 
$O(1)$ or larger, thus allowing for the growth and saturation of the fluctuations via interactions with the mean. When 
$Pr = \mathit {O}(1)$, 
$Re_b \ge \mathit {O}(1)$ necessarily implies 
$Pe_b \ge \mathit {O}(1)$, so the two conditions are equivalent. Since 
$Re_b = \alpha ^2 Re$, this condition is equivalent to 
$Re,Pe \ge O(Fr^{-2})$.
For small 
$Pr$, these reduced equations remain valid as long as 
$Re_b, Pe_b \geq O(1)$. However, for sufficiently small 
$Pr$ it is possible to have an intermediate regime where 
$Re_b \ge O(1)$ while 
$Pe_b = Pr Re_b \ll 1$, and that regime is not captured by the reduced equations derived here and in Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022). Yet, this scenario is likely to be relevant in stellar interiors, where 
$Pr$ is asymptotically small (Garaud et al. Reference Garaud, Medrano, Brown, Mankovich and Moore2015b). We now perform a similar analysis to develop a multiscale model that is valid for low-
$Pe_b$ flows.
3.3. Multiscale model for $Pe_b \ll 1$
In the limit of small 
$Pe_b$, the arguments presented above continue to apply for the continuity equation and the horizontal component of the momentum equation, neither of which involve buoyancy terms. However, the procedure subsequently fails because diffusive effects dominate in the buoyancy equation and are unbalanced unless 
$b$ is much smaller than expected, mirroring the argument made in § 2.2. To capture this correctly within the context of our proposed multiscale model, we must introduce a two-parameter expansion for small 
$\alpha$ and small 
$Pe_b$ in lieu of (3.9).
As in Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022), we assume that the asymptotic expansion in 
$\alpha$ proceeds in powers of 
$\alpha ^{1/2}$ and show the self-consistency of this choice below. Inspired by Lignières (Reference Lignières1999), we now also assume that each field can be expanded as a series in powers of 
$Pe_b$ as well. We therefore have
\begin{equation} q = q_{00} + \alpha^{1/2} q_{01} + \alpha q_{02} + \cdots + Pe_b (q_{10} + \alpha^{1/2} q_{11} + \alpha q_{12} + \cdots) + O(Pe_b^2), \end{equation}
for 
$q \in \{ \boldsymbol {u}_\perp, p, w, b\}$, and assume that 
$\boldsymbol {\bar u}_{\perp 00} = O(1)$ and 
$\bar p_{00} = O(1)$ to balance the forcing. We then proceed exactly as before, substituting (3.14) in turn into the multiscale continuity equation (3.5), the horizontal component of the momentum equation (3.6), the buoyancy equation (3.7) and the vertical component of the momentum equation (3.8), to extract the relevant reduced equations at leading order.
Starting with the mean continuity equation (3.5a), we find again that 
$\bar {w} = O(\alpha )$ and hence 
$\bar {w}_{00} = \bar {w}_{01} = 0$. Equation (3.5b) further implies that 
$O(\boldsymbol {u}_\perp ') = O(w')$, so
$$\begin{gather} \boldsymbol{\nabla}_{s} \boldsymbol{\cdot} \bar{\boldsymbol{u}}_{{\perp} 00} + \frac{\partial \bar{w}_{02}}{\partial \zeta} = 0 , \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla}_{f} \boldsymbol{\cdot} \boldsymbol{u}_{{\perp} 01}' + \frac{\partial w_{01}'}{\partial \zeta} = 0. \end{gather}$$ From the mean horizontal component of the momentum equation (3.6a), 
$\boldsymbol {u'}_{\!\!\perp }$ and 
$w'$ must both be 
$O(\alpha ^{1/2})$ as before, hence 
$\boldsymbol {u}_{\perp 00}' = 0$ and 
$w'_{00} = 0$, so 
$\boldsymbol {u}_{\perp 00} = \bar {\boldsymbol {u}}_{\perp 00}$ and 
$w_{00} = 0$. The inferred sizes of the velocity fluctuations again shows that fluctuation–fluctuation interactions are formally higher order, and the only remaining nonlinearities in the corresponding fluctuation equation (3.6b) are quasilinear. Finally, for the fluctuating horizontal pressure gradient to influence the fluctuations of horizontal velocity at leading order, 
$p' = O(\boldsymbol {u}_\perp ') = O(\alpha ^{1/2})$, hence 
$p_{00}' = 0$. Substituting (3.14) into (3.6a) and (3.6b), using these deductions and retaining only the lowest-order terms, shows that
\begin{align} & \frac{\partial \bar{\boldsymbol{u}}_{{\perp} 00}}{\partial t_s} + ( \bar{\boldsymbol{u}}_{{\perp} 00} \boldsymbol{\cdot} \boldsymbol{\nabla}_{s}) \bar{\boldsymbol{u}}_{{\perp} 00} + \bar{w}_{02} \frac{\partial \bar{\boldsymbol{u}}_{{\perp} 00} }{\partial \zeta} \nonumber\\ &\quad ={-} \boldsymbol{\nabla}_{s} \bar{p}_{00} - \frac{\partial }{\partial \zeta} (\overline{w_{01}' \boldsymbol{u}_{{\perp} 01}'}) + \frac{1}{Re_b} \frac{\partial ^2 \bar{\boldsymbol{u}}_{{\perp} 00}}{\partial \zeta^2} + {\bar{\boldsymbol{F}}_h}, \end{align}
\begin{align} &\frac{\partial \boldsymbol{u}_{{\perp} 01}'}{\partial t_f} + (\bar{\boldsymbol{u}}_{{\perp} 00} \boldsymbol{\cdot} \boldsymbol{\nabla}_{f}) \boldsymbol{u}_{{\perp} 01}' + w_{01}' \frac{\partial \bar{\boldsymbol{u}}_{{\perp} 00}}{\partial \zeta} ={-} \boldsymbol{\nabla}_{f} p_{01}' + \frac{\alpha}{Re_b} \left(\nabla_f^2 \boldsymbol{u}_{{\perp} 01}' + \frac{\partial ^2 \boldsymbol{u}_{{\perp} 01}' }{\partial \zeta^2}\right), \end{align}where the viscous term in the fluctuation equation has been retained for regularisation.
Thus far, each step in the analysis has been identical to that taken in the previous section. Rapid diffusion, however, affects the size of the mean and fluctuating buoyancy fields 
$\bar b$ and 
$b'$, and does so in different ways because the effective Péclet number of the fluctuations is larger than that of the mean flow. More specifically, having assumed in this section that 
$Pe_b \ll 1$, we see that two possibilities arise when 
$\alpha \ll 1$: either 
$\alpha \ll Pe_b \ll 1$, in which case diffusion is dominant in the mean buoyancy equation but negligible in the fluctuation buoyancy equation, or 
$Pe_b \ll \alpha$ in which case diffusion is dominant at all scales. In what follows, we investigate both cases in turn.
Before doing so, however, we note that the mean buoyancy equation (3.7a) in both cases is unbalanced unless 
$\bar {b} = O(\alpha Pe_b)$ to match the 
$\bar w$ term (again mirroring the arguments given in § 2.2). This implies 
$\bar {b}_{0i} = 0$, 
$\forall i$, and 
$\bar b_{10} = \bar b_{11} = 0$. At lowest order in 
$Pe_b$, the only surviving terms in (3.7a) are therefore
\begin{equation} \frac{1}{\alpha} \left(\overline{\boldsymbol{u}_\perp' \boldsymbol{\cdot}\boldsymbol{\nabla}_f b' + w' \frac{\partial b'}{\partial \zeta}}\right) + \bar{w} = \frac{1}{Pe_b} \frac{\partial ^2 \bar{b}}{\partial \zeta^2}, \end{equation}
where the size of 
$b'$ is yet to be determined and differs depending on the relative sizes of 
$Pe_b$ and 
$\alpha$. For this reason, we have retained the turbulent buoyancy flux for now.
3.3.1. Case 1: $\alpha \ll Pe_b \ll 1$ (the intermediate regime)
We first consider a scenario in which 
$\alpha \ll Pe_b \ll 1$ and henceforth refer to this part of parameter space as the ‘intermediate regime’. While diffusion dominates the mean buoyancy equation, the fact that 
$\alpha /Pe_b \ll 1$ implies that it only formally enters the buoyancy fluctuation equation (3.7b) at higher order. Because of this, the evolution of 
$b'$ is very similar to that obtained in Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022). As in § 3.2, we ensure that the background stratification influences the fast dynamics of 
$b'$ by requiring 
$b'/\alpha = O(w') = O(\alpha ^{1/2})$, so 
$b' = O(\alpha ^{3/2})$. This implies that 
$b_{00}', b_{01}', b_{02}' = 0$, but 
$b_{03}' \neq 0$.
Substituting the ansatz (3.14) into the mean and fluctuation buoyancy equations and using the information collected so far we therefore have
$$\begin{gather} \frac{\partial }{\partial \zeta}(\overline{w'_{01} b'_{03}}) + \bar{w}_{02} = \frac{\partial ^2 \bar{b}_{12}}{\partial \zeta^2}, \end{gather}$$
$$\begin{gather}\frac{\partial b'_{03}}{\partial t_f} + \bar{\boldsymbol{u}}_{{\perp} 00} \boldsymbol{\cdot}\boldsymbol{\nabla}_f b'_{03} + w'_{01} = \frac{\alpha}{Pe_b} \left(\nabla_f^2 b'_{03} + \frac{\partial ^2 b'_{03}}{\partial \zeta^2}\right), \end{gather}$$where the diffusion term for the fluctuations can be retained to regularise the equation, but is formally higher order.
We note that the reduced buoyancy fluctuation equation in this regime differs slightly from the one derived by Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) given in (3.12b), because it does not contain a term of the form 
$w'_{01} \partial \bar b/\partial \zeta$, which is formally higher order when 
$Pe_b \ll 1$. Also, we see that the mean buoyancy equation in that regime differs from the asymptotic low Péclet number (LPN) equation of Lignières (Reference Lignières1999), and from the one derived by Skoutnev (Reference Skoutnev2023), which do not contain a turbulent flux term. This discrepancy arises because their derivation assumes that all dynamics is diffusive, whereas in this intermediate regime the fluctuation dynamics is not and can therefore influence the mean buoyancy field at leading order.
Finally, we examine the vertical component of the momentum equation. Based on past experience in the 
$Pe_b \ge O(1)$ case (see § 3.2), one would naively expect to recover hydrostatic equilibrium at leading order in the mean vertical momentum equation (3.8a). Because 
$\bar b = O(\alpha Pe_b)$, this would imply 
$\alpha = (Fr^2 /Pe_b)^{1/2}$ as in the single-scale low-
$Pe_b$ equations of § 2.2. However, that choice leads to an irreconcilable inconsistency in the fluctuation equation (3.8b): the fluctuation pressure gradient term is 
$O(\alpha ^{-1/2})$, as is the fast inertial dynamics of 
$w'$, but the fluctuation buoyancy term is 
$O(\alpha ^{3/2}/Fr^2) = O(\alpha ^{-1/2}/Pe_b)$, which is formally much larger than any other term and is therefore unbalanced. In other words, we cannot reconcile hydrostatic equilibrium at leading order for the mean flow with a balanced equation for the fluctuation 
$w'$.
The solution to this conundrum is to insist instead that the equation for 
$w'$ be balanced, in which case 
$O(\alpha ^{-1/2} ) = O( \alpha ^{3/2}Fr^{-2})$, thereby recovering the standard scaling relationship 
$\alpha = Fr$ (Billant & Chomaz Reference Billant and Chomaz2001; Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007; Chini et al. Reference Chini, Michel, Julien, Rocha and Caulfield2022). With this choice, the leading-order vertical pressure gradient in the mean equation is asymptotically small (e.g. as for the wall-normal pressure gradient in laminar boundary-layer theory of Batchelor Reference Batchelor1967). This perhaps unexpected result is discussed in § 4 and below. Substituting (3.14) into (3.8), we obtain the reduced mean and fluctuating vertical component of the momentum equation at leading order:
$$\begin{gather} \frac{\partial \bar p_{00}}{\partial \zeta} = 0, \end{gather}$$
$$\begin{gather}\frac{\partial w'_{01}}{\partial t_f} + \bar{\boldsymbol{u}}_{{\perp} 00} \boldsymbol{\cdot}\boldsymbol{\nabla}_{f} w'_{01} ={-} \frac{\partial p'_{01}}{\partial \zeta} + b'_{03} + \frac{\alpha}{Re_b} \left(\nabla_f^2 w'_{01} + \frac{\partial ^2 w'_{01}}{\partial \zeta^2}\right), \end{gather}$$where the higher-order viscous term is added to regularise the fluctuation equation.
The set of equations formed by (3.15), (3.16), (3.18) and (3.19) are the intermediate regime analogues of the reduced model given in Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022). They are valid as long as 
$Re_b \ge O(1)$, and 
$\alpha \ll Pe_b \ll 1$. Given that 
$\alpha = Fr$ in this regime, this inequality constraint is equivalent to requiring that 
$Fr \ll 1$ (so 
$\alpha \ll 1$), 
$Re \ge Fr^{-2}$ and 
$Fr^{-1} \ll Pe \ll Fr^{-2}$. It is evident that, in the absence of fluctuations, these equations nearly recover the single-scale equations derived in § 2.2, except that 
$\partial \bar p_{00}/\partial \zeta = 0$. Considering higher-order terms in the vertical momentum equation, we successively find that 
$\partial _\zeta \bar {p}_{01} = 0$ as well; it is only at the next order that buoyancy influences the mean flow, i.e. 
$\partial _\zeta \bar {p}_{10} = \bar {b}_{12}$. By replacing 
$\bar {p}_{00}$ in the vertical mean momentum equation with a composite pressure, 
$\bar {p}_c = \bar {p}_{00} + \alpha ^{1/2} \bar {p}_{01} + Pe_b \bar {p}_{10}$, we recover (in the absence of fluctuations) the single-scale low-
$Pe_b$ equations and those in Skoutnev (Reference Skoutnev2023). For consistency, the corresponding composite mean horizontal momentum equation should be derived accurate to 
$O(\alpha ^{1/2},Pe_b)$ in terms of a composite mean horizontal velocity 
$\bar {\boldsymbol {u}}_{\perp c}$. This exercise, however, does not yield any structurally new terms that do not involve fluctuations (but only higher-order corrections to existing mean terms in (3.16a)). Consequently, in the singular limit of vanishing fluctuations, the multiscale model with a composite mean pressure reduces to the single-scale low-
$Pe_b$ equations in § 2.2. The physical implications of (3.15), (3.16), (3.18) and (3.19), and their potential caveats, are discussed in § 4.
3.3.2. Case 2: $Pe_b \ll \alpha$ (the fully diffusive regime)
We now consider the regime where 
$Pe_b \ll \alpha$, in which both mean and fluctuating buoyancy fields are dominated by diffusion. Accordingly, we refer to this part of parameter space as the fully diffusive regime. Inspection of (3.7b) shows that the diffusion term in the fluctuation equation is unbalanced unless 
$b' = O(Pe_b w')$. We previously found that the vertical velocity fluctuations are 
$O(\alpha ^{1/2})$, which implies here that 
$b'= O(\alpha ^{1/2} Pe_b)$. We conclude that 
$b'_{0i} = 0$ 
$\forall i$, and that 
$b'_{10} = 0$ as well. Combined with the results obtained from analysis of the mean buoyancy equation in the diffusive limit, we conclude that 
$b_{10} = 0$ while 
$b_{11} = b'_{11}$.
Using this information and substituting (3.14) into (3.7), we obtain at lowest order
$$\begin{gather} \bar{w}_{02} = \frac{\partial ^2 \bar{b}_{12}}{\partial \zeta^2}, \end{gather}$$
$$\begin{gather}w'_{01} = \nabla^2_f b'_{11} + \frac{\partial ^2 b'_{11} }{\partial \zeta^2}, \end{gather}$$which is as expected from the LPN dynamics central to this regime (Lignières Reference Lignières1999). The buoyancy equation is linear and does not contain any time dependence, it instead instantaneously couples the vertical velocity and buoyancy fields. The validity of Lignières’ LPN equations was verified numerically by Cope et al. (Reference Cope, Garaud and Caulfield2020), for instance.
As usual, the last step of the derivation involves analysis of the vertical component of the momentum equation. As in § 3.3.1, requiring hydrostatic equilibrium for the mean flow would imply 
$\alpha ^2 = Fr^2 / Pe_b$ (Lignières Reference Lignières2020; Skoutnev Reference Skoutnev2023), but leads to an inconsistency in the fluctuation equation, where the buoyancy term would be unbalanced. To see this, note that 
$b'/Fr^2 = O( \alpha ^{1/2} Pe_b/Fr^2) = O( \alpha ^{-3/2})$ with that choice for 
$\alpha$, while the fluctuating pressure term and all other dominant terms in the equation are only 
$O(\alpha ^{-1/2})$. As in the previous section, the resolution to this inconsistency is to insist that the buoyancy term in the fluctuation equation for 
$w'$ should be balanced instead. Here, this implies
\begin{equation} O\left(\frac{p'}{\alpha}\right) = O\left(\frac{b'}{Fr^2}\right) \rightarrow \alpha = \frac{Fr^2}{Pe_b}, \end{equation}
using the fact that 
$p' = O(\alpha ^{1/2})$ and 
$b' = O(\alpha ^{1/2} Pe_b)$. Recalling that 
$Pe_b = \alpha ^2 Pe$, we then recover the crucial scaling relationship
\begin{equation} \alpha = \left(\frac{Fr^2}{ Pe}\right)^{1/3} = Fr_M^{4/3}, \end{equation}which had originally been proposed by Cope et al. (Reference Cope, Garaud and Caulfield2020) based on their DNS data. Our multiscale analysis therefore provides a sound theoretical basis for their empirical results.
After substituting (3.14) into (3.8) and using the available information, we obtain
$$\begin{gather} \frac{\partial \bar p_{00}}{\partial \zeta} = 0, \end{gather}$$
$$\begin{gather}\frac{\partial w'_{01}}{\partial t_f} + (\bar{\boldsymbol{u}}_{{\perp} 00} \boldsymbol{\cdot}\boldsymbol{\nabla}_{f}) w'_{01} ={-} \frac{\partial p'_{01}}{\partial \zeta} + b'_{11} + \frac{\alpha}{Re_b} \left(\nabla_f^2 w'_{01} + \frac{\partial ^2 w'_{01}}{\partial \zeta^2}\right), \end{gather}$$
which is very similar to the system obtained in the intermediate regime studied in § 3.3.1, except for the appearance of the buoyancy fluctuation term 
$b'_{11}$ instead of 
$b'_{03}$. As before, formally higher-order viscous terms are retained to regularise the fluctuation equation.
The set of equations formed by (3.15), (3.16), (3.20) and (3.23) are the fully diffusive regime analogues of the reduced model derived by Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022). They are valid as long as 
$Re_b \ge O(1)$, and 
$Pe_b \ll \alpha \ll 1$. Given that 
$\alpha = (Fr^2/Pe)^{1/3}$ in this regime, this parameter constraint is equivalent to requiring that 
$Fr^2 \ll Pe$ (to ensure 
$\alpha \ll 1$), 
$Pe \ll Fr^{-1}$ (to ensure 
$Pe_b \ll \alpha$) and 
$Pe \ge Pr^3 Fr^{-4}$ (to ensure that 
$Re_b \ge O(1)$). The three conditions demarcate a triangle in logarithmic parameter space in which the equations are valid – see § 4.2. As an important self-consistency check, we see that the 
$Pe_b = \alpha$ transition between the fully diffusive and intermediate regimes is the same (
$Pe = Fr^{-1}$) whether the transition is approached from the former or latter part of parameter space. Here, too, (weak) buoyancy effects can be included in the mean dynamics by considering higher-order pressure contributions to the vertical momentum equation. We successively find that 
$\partial _\zeta \bar {p}_{01} =0$ and 
$\partial _\zeta (\overline {w'_{01} w'_{01}}) = - \partial _\zeta \bar {p}_{02} + \bar {b}_{12}$. Hence, when fluctuations are absent from the system of multiscale equations given above, the resulting equations with an appropriately defined composite mean pressure 
$\bar {p}_c = \bar {p}_{00} + \alpha ^{1/2} \bar {p}_{01} + \alpha \bar {p}_{02}$ recover hydrostatic balance and the single-scale low-
$Pe_b$ equations given in § 2.2. In the next section (§ 4), we consider the physical implications of equations (3.15), (3.16), (3.20) and (3.23) as well as potential caveats on their applicability.
4. Discussion
In this study, we have extended the results of Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) by performing a multiscale asymptotic analysis of stratified turbulence at low Prandtl number. Our work demonstrates the existence of several different regimes depending on the strength of the stratification (quantified by the inverse Froude number) and the rate of buoyancy diffusion (quantified by the inverse Péclet number). In each regime, the asymptotic analysis self-consistently yields a slow–fast system of quasilinear equations describing the concurrent evolution of a highly anisotropic, slow, large-scale mean flow and isotropic, fast, small-scale fluctuations. The small scales are self-consistently excited by an instability of the emergent, local vertical shear. When they are not, the equations for the large-scale anisotropic flow recover exactly those obtained in § 2 (with a composite mean vertical momentum equation in the low 
$Pe_b$ case).
The large-scale anisotropy is characterised by the aspect ratio 
$\alpha$ (the ratio of the vertical to horizontal scales of the large-scale flow), whose functional dependence on 
$Fr$ and 
$Pe$ naturally emerges from the analysis. The various regime boundaries are locations in parameter space where relevant Reynolds and Péclet numbers arising in the mean flow and fluctuation equations are 
$O(1)$, signifying transitions between a viscous and (formally) inviscid dynamics, and/or a diffusive and non-diffusive dynamics. We now summarise our findings in each regime and then discuss the model assumptions as well as the implications of our results for low Prandtl number fluids.
4.1. Synopsis of multiscale equations and their validity
The first regime is characterised by 
$Fr \ll 1$ and 
$Pe_b, Re_b \ge O(1)$, where 
$Re_b = \alpha ^2 Re$ and 
$Pe_b = \alpha ^2 Pe$. In this regime, we recover the reduced model of Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) and confirm that 
$\alpha = Fr$. Recalling that, at leading order, 
$w'_1 = w' / \alpha ^{1/2}$, 
$b'_3 = b' / \alpha ^{3/2}$, etc., that 
$\partial /\partial t_f = \alpha \partial /\partial t$ (and similarly 
$\boldsymbol {\nabla }_f = \alpha \boldsymbol {\nabla }_\perp$), and finally that 
$\partial /\partial \zeta = \alpha \partial /\partial z$, we can rewrite (3.10), (3.11), (3.12) and (3.13) as the following quasilinear system for mean and fluctuations, expressed in the original variables 
$(x,y,z,t)$.
Mean flow equations:
$$\begin{gather} \frac{\partial \bar{\boldsymbol{u}}_{{\perp}} }{\partial t} + \bar{\boldsymbol{u}}_{{\perp}} \boldsymbol{\cdot} \boldsymbol{\nabla}{_\perp} \bar{\boldsymbol{u}}_{{\perp}} {+ \bar{w} \frac{\partial \boldsymbol{\bar{u}}_\perp}{\partial z}} ={-} \boldsymbol{\nabla}_{{\perp}} \bar{p} - \frac{\partial }{\partial z} (\overline{w' {\boldsymbol u'}_{\!\!\perp}}) + \frac{1}{Re} \frac{\partial ^2 \bar{\boldsymbol{u}}_{{\perp}}}{\partial z^2} + {\bar{\boldsymbol{F}}_h}, \end{gather}$$
$$\begin{gather}\frac{\partial \bar{p}}{\partial z} = \frac{\bar b}{Fr^2},\quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\bar{u}} = 0, \end{gather}$$
$$\begin{gather}\frac{\partial \bar{b}}{\partial t} + \bar{\boldsymbol{u}}_{{\perp}} \boldsymbol{\cdot}\boldsymbol{\nabla}{_{{\perp}}} \bar{b} + \bar{w} \frac{\partial \bar{b}}{\partial z} + \bar{w} ={-} \frac{\partial }{\partial z} (\overline{w' b'}) + \frac{1}{Pe} \frac{\partial ^2\bar{b}}{\partial z^2}. \end{gather}$$Fluctuation equations:
$$\begin{gather} \frac{\partial \boldsymbol{u}_{{\perp}}'}{\partial t} + \bar{\boldsymbol{u}}_{{\perp}} \boldsymbol{\cdot}\boldsymbol{\nabla}_{{\perp}} \boldsymbol{u}_{{\perp}}' + w' \frac{\partial \bar{\boldsymbol{u}}_{{\perp}}}{\partial z} ={-} \boldsymbol{\nabla}_{{\perp}} p' + \frac{1}{Re} \nabla^2 \boldsymbol{u}_{{\perp}}' , \end{gather}$$
$$\begin{gather}\frac{\partial w'}{\partial t} + \bar{\boldsymbol{u}}_{{\perp}} \boldsymbol{\cdot}\boldsymbol{\nabla}_{{\perp}} w' ={-} \frac{\partial p'}{\partial z} +\frac{b'}{Fr^2} + \frac{1}{Re} \nabla^2 w', \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}'= 0, \end{gather}$$
$$\begin{gather}\frac{\partial b'}{\partial t} + \bar{\boldsymbol{u}}_{{\perp}} \boldsymbol{\cdot}\boldsymbol{\nabla}_{{\perp}} b' + w' \frac{\partial \bar{b}}{\partial z} + w' = \frac{1}{Pe} \nabla^2 b' . \end{gather}$$ This system, which is equivalent to (2.28)–(2.35) in Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022), is valid when 
$Fr \ll 1$, and 
$Re,Pe \ge O(Fr^{-2})$, corresponding to 
$Re_b, Pe_b \geq O(1)$. Noting that the overbar denotes an intermediate-averaging operation (see (3.3)), our reduced model has a generalised quasilinear (GQL) form (Marston, Chini & Tobias Reference Marston, Chini and Tobias2016; Marston & Tobias Reference Marston and Tobias2023), albeit here derived in physical space (whereas the GQL reduction usually is performed in Fourier space). When slow horizontal variation is suppressed and the overbar denotes a strict horizontal average, the above reduced equations have the same form as an ad hoc quasilinear (QL) reduction of the Boussinesq equations (e.g. Garaud Reference Garaud2001; Fitzgerald & Farrell Reference Fitzgerald and Farrell2018, Reference Fitzgerald and Farrell2019). In this work, we have demonstrated that this (G)QL form is in fact a natural outcome of the slow–fast asymptotic expansion and is therefore asymptotically exact in the given distinguished limit. As detailed in Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022), in two dimensions the fluctuation equations without diffusion can be reduced to the Taylor–Goldstein equation (see for example Craik Reference Craik1985) as the mean fields are independent of the fast variables 
$\boldsymbol {x}_{\perp f}$ and 
$t_f$ and hence are effectively frozen during the fast evolution of the fluctuations. Thus, the fluctuation equations admit as solutions all of the ‘classical’ stratified shear instabilities of parallel shear flows, including in particular Kelvin–Helmholtz instabilities. This mathematical reduction thus may explain the occurrence of Kelvin–Helmholtz billows within fully turbulent but strongly stratified flows.
Dimensionally, this regime has a characteristic vertical scale 
$H^* = \alpha L^* = U^*/N^*$. The characteristic vertical velocity is 
$W^* = \alpha ^{1/2} U^* = (U^{*3}/N^*L^*)^{1/2}$, and 
$w$ is dominated by small-scale fluctuations. The characteristic buoyancy scale is 
$N^{*2} H^*$, and 
$b$ is dominated by large scales. Note that the horizontal field 
$\boldsymbol {u}_\perp$ is predominantly large scale, by assumption, but contains small scales with (dimensional) amplitudes 
$O(Fr^{1/2}U^*)$. The scaling law for 
$H^*$ has been validated in numerical and laboratory experiments (e.g. Holford & Linden Reference Holford and Linden1999; Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007; Oglethorpe et al. Reference Oglethorpe, Caulfield and Woods2013), with at least suggestive observational evidence of this LAST regime being obtainable from seismic oceanography surveys, as already noted above (Falder et al. Reference Falder, White and Caulfield2016). The scaling law for 
$W^*$ was first tentatively identified by Maffioli & Davidson (Reference Maffioli and Davidson2016). More recently Garaud et al. (Reference Garaud, Chini, Cope, Shah and Caulfield2024a) measured the root-mean-squared vertical velocity of the flow in DNS of strongly stratified simulations where the turbulence is spatio-temporally intermittent. Using different weight functions to distinguish averages taken within and outside turbulent regions of the flow, they confirmed that 
$W^* \propto Fr^{1/2} U^*$ within turbulent patches. In quiescent regions of the flow outside these turbulent patches, where small-scale fluctuations are suppressed, Garaud et al. (Reference Garaud, Chini, Cope, Shah and Caulfield2024a) found that 
$W^* \propto FrU^*$, consistent with the predictions of Billant & Chomaz (Reference Billant and Chomaz2001) and Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007).
The scenario in which 
$Fr \ll 1$, 
$Re_b \ge O(1)$ and 
$\alpha \ll Pe_b \ll 1$ is an intermediate regime where the mean flow is diffusive while the fluctuations are not. In this regime, we also find that 
$\alpha = Fr$. The slow–fast system of (3.15), (3.16), (3.18) and (3.19), expressed in the original variables 
$(x,y,z,t)$, becomes the following QL system.
Mean flow equations:
$$\begin{gather} \frac{\partial \bar{\boldsymbol{u}}_{{\perp}} }{\partial t} + \bar{\boldsymbol{u}}_{{\perp}} \boldsymbol{\cdot} \boldsymbol{\nabla}{_\perp} \bar{\boldsymbol{u}}_{{\perp}} + \bar{w} \frac{\partial \boldsymbol{\bar{u}}_\perp}{\partial z} ={-} \boldsymbol{\nabla}_{{\perp}} \bar{p} - \frac{\partial }{\partial z} (\overline{w' {\boldsymbol u'}_\perp}) + \frac{1}{Re} \frac{\partial ^2 \bar{\boldsymbol{u}}_{{\perp}}}{\partial z^2} + {\bar{\boldsymbol{F}}_h}, \end{gather}$$
$$\begin{gather}\frac{\partial \bar{p}}{\partial z} = 0,\quad \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\bar{u}} = 0 , \end{gather}$$
$$\begin{gather}\bar{w} ={-} \frac{\partial }{\partial z} (\overline{w' b'}) + \frac{1}{Pe} \frac{\partial ^2\bar{b}}{\partial z^2}. \end{gather}$$Fluctuation equations:
$$\begin{gather} \frac{\partial \boldsymbol{u}_{{\perp}}'}{\partial t} + \bar{\boldsymbol{u}}_{{\perp}} \boldsymbol{\cdot}\boldsymbol{\nabla}_{{\perp}} \boldsymbol{u}_{{\perp}}' + w' \frac{\partial \bar{\boldsymbol{u}}_{{\perp}}}{\partial z} ={-} \boldsymbol{\nabla}_{{\perp}} p' + \frac{1}{Re} \nabla^2 \boldsymbol{u}_{{\perp}}', \end{gather}$$
$$\begin{gather}\frac{\partial w'}{\partial t} + \bar{\boldsymbol{u}}_{{\perp}} \boldsymbol{\cdot}\boldsymbol{\nabla}_{{\perp}} w' ={-} \frac{\partial p'}{\partial z} + \frac{b'}{Fr^2} + \frac{1}{Re} \nabla^2 w', \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}'= 0, \end{gather}$$
$$\begin{gather}\frac{\partial b'}{\partial t} + \bar{\boldsymbol{u}}_{{\perp}} \boldsymbol{\cdot}\boldsymbol{\nabla}_{{\perp}} b' + w' = \frac{1}{Pe} \nabla^2 b' . \end{gather}$$ These equations are valid for 
$Re \ge Fr^{-2}$ (so 
$Re_b \ge O(1)$) and 
$Fr^{-1} \ll Pe \ll Fr^{-2}$ (so 
$\alpha \ll Pe_b \ll 1$). By writing them in the original isotropic variables, we now see that the correct QL equations at this order are almost the same as in the Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) regime, except that terms in 
$\bar {b}$ are dropped (in the mean hydrostatic balance and in the buoyancy perturbation equation) because they are formally of higher order. In addition, the mean buoyancy equation takes the LPN form of Lignières (Reference Lignières1999), modified by the fluctuation-induced buoyancy flux. Given that the mean fields are independent of fast variables, the fluctuation subsystem again can be treated as a (two-dimensional) eigenvalue problem on fast scales that, in this case, admits as solutions linear instability modes of low-
$Pr$ parallel shear flows (e.g. Jones Reference Jones1977; Lignières et al. Reference Lignières, Califano and Mangeney1999).
Dimensionally, the vertical length scale and vertical velocity scale are the same as in the non-diffusive regime. The buoyancy field, however, is now dominated by large scales if 
$Pe_b \ge \alpha ^{1/2}$ and by small scales if 
$Pe_b \le \alpha ^{1/2}$. Testing these scaling laws numerically will be very difficult, unfortunately, because the range of the intermediate region is very small (holding 
$Pe$ constant while varying 
$Fr^{-1}$ or vice versa) unless 
$Pr$ is itself very small.
Finally, the case where 
$Fr\ll 1$, 
$Re_b \ge O(1)$ and 
$Pe_b \ll \alpha$ corresponds to a fully diffusive regime in which both the mean flow and fluctuations are dominated by diffusion and satisfy the LPN balance derived by Lignières (Reference Lignières1999). In this regime, we have demonstrated that 
$\alpha = (Fr^2/Pe)^{1/3}$. The slow–fast system of (3.15), (3.16), (3.20) and (3.23), written in the original variables 
$(x,y,z,t)$, is given below.
Mean flow equations:
$$\begin{gather} \frac{\partial \bar{\boldsymbol{u}}_{{\perp}} }{\partial t} + \bar{\boldsymbol{u}}_{{\perp}} \boldsymbol{\cdot} \boldsymbol{\nabla}{_\perp} \bar{\boldsymbol{u}}_{{\perp}} + \bar{w} \frac{\partial \boldsymbol{\bar{u}}_\perp}{\partial z} ={-} \boldsymbol{\nabla}_{{\perp}} \bar{p} - \frac{\partial }{\partial z} (\overline{w' {\boldsymbol u'}_{\!\!\perp}}) + \frac{1}{Re} \frac{\partial ^2 \bar{\boldsymbol{u}}_{{\perp}}}{\partial z^2} +{\bar{\boldsymbol{F}}_h}, \end{gather}$$
$$\begin{gather}\frac{\partial \bar{p}}{\partial z} = 0,\quad \boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{\bar{u}} = 0 , \end{gather}$$
$$\begin{gather}\bar{w} = \frac{1}{Pe} \frac{\partial ^2\bar{b}}{\partial z^2}. \end{gather}$$Fluctuation equations:
$$\begin{gather} \frac{\partial \boldsymbol{u}_{{\perp}}'}{\partial t} + \bar{\boldsymbol{u}}_{{\perp}} \boldsymbol{\cdot}\boldsymbol{\nabla}_{{\perp}} \boldsymbol{u}_{{\perp}}' + w' \frac{\partial \bar{\boldsymbol{u}}_{{\perp}}}{\partial z} ={-} \boldsymbol{\nabla}_{{\perp}} p' + \frac{1}{Re} \nabla^2 \boldsymbol{u}_{{\perp}}', \end{gather}$$
$$\begin{gather}\frac{\partial w'}{\partial t} + \bar{\boldsymbol{u}}f_{{\perp}} \boldsymbol{\cdot}\boldsymbol{\nabla}_{{\perp}} w' ={-} \frac{\partial p'}{\partial z} + \frac{b'}{Fr^2} + \frac{1}{Re} \nabla^2 w', \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}'= 0, \end{gather}$$
$$\begin{gather}w' = \frac{1}{Pe} \nabla^2 b' . \end{gather}$$ This system is valid for 
$Re_b \ge O(1)$, and 
$Pe_b \ll \alpha \ll 1$. As 
$\alpha = (Fr^2/Pe)^{1/3}$, these constraints are equivalent to 
$Fr^2 \ll Pe$ (so that 
$\alpha \ll 1$), 
$Pe \ll Fr^{-1}$ (so that 
$Pe_b \ll \alpha$) and 
$Pe \ge Pr^3 Fr^{-4}$ (so that 
$Re_b \ge O(1)$).
Dimensionally, the characteristic vertical length scale 
$H^* = \alpha L^* = (Fr^2/Pe)^{1/3} L^* = (U^* \kappa ^*/N^{*2})^{1/3}$. The characteristic vertical velocity 
$W^* = \alpha ^{1/2} U^* =(Fr^2/Pe)^{1/6} U^* = (U^{*7}\kappa ^*/N^{*2} L^{*3})^{1/6}$ and is dominated by small-scale fluctuations. The characteristic buoyancy scale is 
$Pe (Fr^2/Pe)^{5/6} L^* N^{*2} = (U^{*11/6} \kappa ^{*-1/6} N^{*-5/3} L^{*-3/2}) L^* N^{*2}$ and is dominated by small-scale fluctuations as well. These scalings have been validated by the DNS of Cope et al. (Reference Cope, Garaud and Caulfield2020) in fully turbulent flows. Garaud et al. (Reference Garaud, Chini, Cope, Shah and Caulfield2024a) showed that they also remain valid within turbulent patches for DNS that are in the intermittent regime of Cope et al. (Reference Cope, Garaud and Caulfield2020). In quiescent regions outside of the turbulent patches, predictions from the single-scale equations apply (cf. § 2.2 and Skoutnev Reference Skoutnev2023).
Crucially, we find that in all three regimes the characteristic vertical velocity 
$W^* = \alpha ^{1/2} U^*$ is significantly larger than that predicted from the single-scale equations given in §§ 2.1 and 2.2, where 
$W^* = \alpha U^*$. This larger scaling has implications for turbulent vertical transport of buoyancy and passive scalars in stellar interiors (see § 5).
Finally, note that, in all of these regimes, we have assumed 
$Re_b \ge O(1)$, which then implies that 
$Re_b \gg \alpha$ since 
$\alpha \ll 1$. Recalling that the fluctuation equations have an effective Reynolds number 
$Re_b/\alpha$, this condition is necessary to ensure that the small-scale fluctuations can develop without being suppressed by viscosity, and is therefore key to the multiscale expansions derived here. When 
$Re_b < 1$, viscous effects become important and could strongly affect our conclusions. We do not pursue this issue further here, deferring discussion of viscous regimes to future work. However, we note that Garaud et al. (Reference Garaud, Chini, Cope, Shah and Caulfield2024a) confirmed that 
$Re_b = O(1)$ correctly predicts the boundary of the region of parameter space where the multiscale equations and their predicted scalings apply.
While this study focuses on the scaling relationships that emerge from the derivation of these multiscale equations and their implications for turbulence regimes (see § 4.2 below), solving these equations is, of course, also an area of interest not just from a scientific but also an algorithmic standpoint. Accordingly, using existing algorithms for slow–fast systems (Michel & Chini Reference Michel and Chini2019; Chini et al. Reference Chini, Michel, Julien, Rocha and Caulfield2022; Ferraro Reference Ferraro2022), our ongoing work involves numerically solving the multiscale partial differential equations for a range of values of 
$Re_b$, which will allow numerical identification of the parameter values at which these behaviours are realised (Shah Reference Shah2022).
4.2. Regimes of stratified stellar turbulence
Our findings partition parameter space into various regimes of stratified turbulence, which are illustrated in figure 1 for (a) a 
$Pr=1$ fluid , (b) a 
$Pr=10^{-2}$ fluid and (c) a 
$Pr=10^{-6}$ fluid. Panel (b) is typical of liquid metals and the bottom of some stellar interiors (see Garaud Reference Garaud2021). In all panels, the horizontal axis shows the inverse Froude number 
$Fr^{-1}$, so that stratification increases to the right. The vertical axis shows the outer-scale Péclet number 
$Pe$, so that diffusive effects decrease upward. In both panels, the grey region shows where 
$\alpha = O(1)$, so the large-scale flow is isotropic and effectively unstratified; this regime is not the focus of our study. For stronger stratification, the large-scale flow becomes anisotropic, and the possible regimes of stratified turbulence depend on 
$Pr$.
Figure 1. Regime diagram for a fluid with (a) 
$Pr = 1$, (b) 
$Pr = 10^{-2}$ and (c) 
$Pr=10^{-6}$. The Péclet number is on the vertical axis and inverse Froude number is on the horizontal axis (such that stratification increases to the right). The unstratified regime is shown in grey. Regions in white are viscously controlled. The 
$Pe_b \ge O(1)$ regime is marked in green. In each panel, the blue dashed line marks the transition 
$Pe_b = O(1)$ where the mean flow becomes diffusive, and the red solid lines mark the viscous transition where 
$Re_b = O(1)$. With 
$Pr \ll 1$ (b,c), the solid blue line marks 
$Pe_b = O(\alpha )$ where the small-scale fluctuations become diffusive. A region of parameter space opens up between 
$Pe_b =O(1)$ and 
$Re_b =O(1)$, where the two new regimes identified in this work exist: the intermediate regime (yellow), and the fully diffusive regime (purple). (c) A parameter set indicative of the solar tachocline is marked with a black circle, using typical parameters given by Garaud (Reference Garaud2020) (
$Pe=10^8$, 
$Fr^{-1}=10^{3}$).
For 
$Pr = O(1)$, shown in the upper plot, the partitioning of parameter space for 
$Fr^{-1} \gg 1$ is straightforward. When 
$Pe_b, Re_b \ge O(1)$ (green region), viscosity and diffusion play a secondary role in both the mean and fluctuation dynamics. Billant & Chomaz (Reference Billant and Chomaz2001) and Lindborg (Reference Lindborg2006) put forward the basis for this regime, for which asymptotic theory was developed by Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022). By contrast, if 
$Pe_b, Re_b \ll 1$ (white region) then both effects strongly influence the flow dynamics. The transition takes place when 
$Re,Pe \approx Fr^{-2}$.
At low 
$Pr$, diffusion becomes important long before viscosity does, so the 
$Pe_b = O(1)$ transition is distinct from the 
$Re_b = O(1)$ transition. This distinction opens up parameter space to the two new regimes discussed in this work: the intermediate regime (yellow region, the dynamics of which is described in § 3.3.1) and the fully diffusive regime (purple region, the dynamics of which is described in § 3.3.2). We now clearly see that the fully diffusive regime is confined to a triangle in log–log parameter space for a given Prandtl number. It is delimited from above by the intermediate regime (yellow region), from below by the isotropic regime (grey region), and from the right by the viscous regime (white region). More specifically, it is bounded by the points A 
$(Fr^{-1} = 1, Pe = 1)$, B 
$(Fr^{-1} = Pr^{-1}, Pe = Pr^{-1})$ and C 
$(Fr^{-1} = Pr^{-1/2}, Pe = Pr)$ and thus becomes increasingly wide as 
$Pr$ decreases, but shrinks towards the point A as 
$Pr \rightarrow 1$ (e.g. compare (b,c)). Note that the isotropic (grey) region in (b,c) contains an additional triangular area which arises from the scaling relationship (3.22): 
$\alpha = 1$ corresponds to 
$Pe = Fr^{2}$ there, in contrast to the 
$Pr \simeq \mathit {O}(1)$ case in (a) where 
$\alpha = 1$ implies 
$Fr = 1$. In all panels, the non-turbulent regime where 
$Re < 1$ is marked by 
$Pe < Pr$, which has been left white.
A major implication of our results is that, depending on the choice of 
$Pe$, different regimes are encountered as the stratification increases. We now discuss horizontal transects through figure 1(c), for different values of 
$Pe$. We first consider the case 
$Pr < Pe < 1$, which is the regime discussed in Cope et al. (Reference Cope, Garaud and Caulfield2020). As 
$Fr^{-1}$ increases (while holding 
$Pe$ and 
$Re$ constant), our model predicts that the turbulence ought to be isotropic until 
$Fr^{-1} = Pe^{-1/2}$ (interestingly, because diffusion partially relaxes the effects of stratification when 
$Pe \ll 1$). As 
$Fr^{-1}$ continues to increase, the turbulence enters the fully diffusive anisotropic regime, and remains in that regime until 
$Fr^{-1} = (Pe/Pr^3)^{1/4}$, at which point viscosity begins to affect the mean flow. This series of regime transitions is qualitatively consistent with that observed in the low Prandtl number DNS of Cope et al. (Reference Cope, Garaud and Caulfield2020).
At the other extreme, let us consider a transect for 
$Pe > Pr^{-1}$, which is the case considered by Garaud (Reference Garaud2020), who primarily analysed simulations for which 
$Pe = 60$ and 
$Pr = 0.1$. This regime is relevant for strongly sheared layers in stellar interiors, such as the solar tachocline. For moderate stratification, namely 
$1 \le Fr^{-1} \le Pe$, our analysis shows that the turbulence is expected to be both anisotropic and non-diffusive, and its properties should be captured by the model of Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022). As stratification increases past 
$Fr^{-1} = Pe$, the turbulence is predicted to enter the intermediate regime, where the mean flow is dominated by diffusion but the fluctuations are not. Beyond 
$Fr^{-1} = \sqrt {Pe /Pr}$, viscous effects should become important. Note how, at these large values of the Péclet number, the fully diffusive inviscid regime discussed in § 3.3.2 is not accessible. Instead, viscosity begins to influence the mean flow before diffusion influences the fluctuations.
This series of regime transitions is qualitatively consistent with the simulations reported in Garaud (Reference Garaud2020). Specifically, she found that the turbulence is mostly isotropic for 
$Fr^{-1} < 1$, then becomes anisotropic with little effect of diffusion for intermediate values of 
$Fr^{-1}$. However, her empirically derived scaling laws (
$H^* \propto Fr^{2/3} L^*$, 
$W^* \propto Fr^{2/3} U^*$) do not match those predicted by the Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) theory. Garaud et al. (Reference Garaud, Chini, Cope, Shah and Caulfield2024a) have now resolved this discrepancy by showing that the apparent 
$W^* \propto Fr^{2/3} U^*$ scaling law is an artefact of the coexistence of turbulent patches where 
$W^* \sim Fr^{1/2} U^*$ and laminar regions where 
$W^* \sim Fr U^*$, whose volume fraction respectively decreases and increases as stratification increases.
For even larger values of 
$Fr^{-1}$, Garaud (Reference Garaud2020) found that diffusion becomes important, and that the dynamics is governed by the LPN equations. However, at this point the turbulence is also viscously affected, a situation qualitatively consistent with the right-most regime transition in figure 1.
The case where 
$1 < Pe < Pr^{-1}$ is relevant for weaker shear layers in stellar interiors. As 
$Fr^{-1}$ increases, our multiscale analysis predicts that the turbulence ought to be isotropic until 
$Fr^{-1} = 1$, at which point the regime analysed by Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) sets in and both the mean and fluctuation dynamics are non-diffusive. As the stratification continues to increase, the vertical eddy scale decreases gradually until the mean flow becomes dominated by diffusion at 
$Fr^{-1} = \sqrt {Pe}$ and the intermediate regime is manifest. The flow dynamics remains essentially unchanged in that regime until diffusion also begins to affect the fluctuations as well, at 
$Fr^{-1} = Pe$, at which point the turbulence becomes fully diffusive. Finally, viscous effects begin to be important when 
$Fr^{-1} = (Pe/Pr^3)^{1/4}$. To date, no DNS have been published probing this range of 
$Pe$. Verifying the existence of these successive transitions with DNS is one of our areas of active work.
4.3. Mathematical considerations
Having constructed a map of parameter space based on the results of our multiscale analysis, we now discuss important consequences and caveats of the resulting reduced models. As in Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022), we find that the reduced equations in each regime form a closed set describing the concurrent evolution of a highly anisotropic large-scale mean flow together with isotropic small-scale fluctuations. Crucially, the fluctuation equations derived in all regimes identified are linear in the fluctuation fields, but feed back nonlinearly on the mean flow evolution. The reduced models thus all have a QL form, which emerges self-consistently from the asymptotic analysis rather than being imposed a priori. This emergent quasilinearity in the slow–fast limit therefore appears to be an inherent property of stratified turbulence regardless of the Prandtl number. Our findings can be used as a theoretical basis not only for using the QL approximation to create reduced models of stratified turbulence, but also to know precisely which terms to keep in each region of parameter space (see e.g. Marston & Tobias Reference Marston and Tobias2023). Of course, the fundamental premise upon which the analysis is predicated is that the stratification drives a scale-separated dynamics characterised by anisotropic large-scale and roughly isotropic small-scale flow structures. To the extent that spectrally non-local interactions between these disparate scales of motion play a crucial role in stratified turbulence, the conclusions and predictions of our multiscale analysis are likely to be valid.
In QL systems subject to fast instabilities, nonlinear saturation requires the feedback from the finite-amplitude fluctuations to render the mean fields marginally stable to disturbances of any horizontal wavenumber. It is straightforward to verify that this condition of approximate marginal stability indeed characterises each regime. This deduction is consistent with the empirical observation of ‘self-organised criticality’ noted in § 1, giving further credence to the appropriateness of the multi-scale approach used here. In the 
$Pe_b \ge O(1)$ and intermediate regimes, the fluctuation equations are non-diffusive at leading order, and therefore describe the growth (or decay) of perturbations due to a standard vertical shear instability. The dimensional vertical shear 
$S^*$ of the mean flow (which emerges once velocity layers develop) can be estimated to be roughly 
$U^*/H^*$ since 
$|\boldsymbol {\bar u}| = O(1)$. The gradient Richardson number
\begin{equation} J = \frac{N^{*2}}{S^{*2}} = O\left(\frac{N^{*2}H^{*2}}{U^{*2}}\right) = O\left(\frac{\alpha^2}{Fr^2}\right) = O(1) \end{equation}
in these regimes, which is indicative of marginal stability according to the Richardson and Miles–Howard criteria (Richardson Reference Richardson1920; Howard Reference Howard1961; Miles Reference Miles1961). Note that Garaud et al. (Reference Garaud, Khan and Brown2024b) recently verified that 
$J$ is indeed 
$O(1)$ in DNS of stratified turbulence at 
$Pe_b = O(1)$.
In the 
$Pe_b \ll \alpha$ regime, by contrast, the fluctuation equations are inherently diffusive, and one therefore expects the vertical shear instability to be of the diffusive kind (Zahn Reference Zahn1974; Jones Reference Jones1977; Lignières et al. Reference Lignières, Califano and Mangeney1999). It has been shown, at least for sinusoidal shear flows (Garaud, Gallet & Bischoff Reference Garaud, Gallet and Bischoff2015a), that the condition for marginal linear stability is 
$J Pe_S = O(1)$ where 
$Pe_S = U^* H^*/ \kappa ^*$ is the Péclet number based on the vertical shear profile. We can easily check that the scalings found in § 3.3.2 satisfy this criterion. Indeed
\begin{equation} J Pe_S = O\left(\frac{N^{*2}H^{*2}}{U^{*2}}\frac{U^* H^*}{\kappa^*} \right) = O\left(\alpha^3 \frac{Pe}{Fr^2} \right) = O(1) , \end{equation}
as required. Taken together with the empirical scaling laws obtained from the DNS data of Cope et al. (Reference Cope, Garaud and Caulfield2020), this evidence confirms that (3.22) is the correct expression for 
$\alpha$ in stratified turbulence at very low 
$Pe_b$.
A somewhat less intuitive result of our analysis is the requirement that 
$\partial \bar p_{00}/\partial \zeta = 0$, while 
$\bar p_{00} = O(1)$, in both the intermediate and fully diffusive regimes. Physically, this condition can be understood by noting that in the 
$Pe_b \ll 1$ scenario, departures from the mean background stratification formally are exceedingly small (owing to the strong thermal diffusion), which explains why they do not affect the assumed background hydrostatic balance. Furthermore, recalling that 
$\zeta$ has been rescaled by the vertical length scale 
$H^* = \alpha L^*$ with 
$\alpha \ll 1$, we note that 
$\partial _\zeta \bar {p}_{00} = 0$ only applies on that scale, and does not preclude 
$\bar p_{00}$ from potentially varying on larger scales. Indeed, the possibility exists that two vertical scales may be incorporated into the dynamics, with the leading-order mean pressure being hydrostatically coupled to the leading-order mean buoyancy on the larger vertical scale, although we leave exploration of that possibility for future work. Alternatively, as discussed in § 3.3, higher-order mean pressure terms, which couple to the mean buoyancy, can be incorporated using a composite mean vertical momentum equation. Then, in the singular limit of no fluctuations, the sets of multiscale equations for the intermediate (4.2) and fully diffusive (4.3) regimes recover the single-scale equations in § 2.2. Regardless, we emphasise that buoyancy anomalies do affect the fluctuation dynamics, which in turn modifies the mean flow.
5. Conclusions
In this study, inspired by numerical evidence of scale separation and flow anisotropy, we have conducted a formal multiscale asymptotic analysis of the Boussinesq equations governing the dynamics of SST at low Prandtl number. A key outcome of our work is a new map of parameter space, shown in figure 1, that demarcates different regimes of stratified turbulence. Crucially, we find that new regions of parameter space open up at low Prandtl number in which diffusive turbulent flows exist, scenarios that are not possible at 
$Pr = \mathit {O}(1)$. For each of these new regimes, scaling laws for the vertical velocity and vertical length scale of turbulent eddies naturally emerge from the analysis. These scaling laws are summarised in § 4.2 and recover previous findings by Chini et al. (Reference Chini, Michel, Julien, Rocha and Caulfield2022) and Cope et al. (Reference Cope, Garaud and Caulfield2020) in appropriate distinguished limits. Finally, recent work by Garaud et al. (Reference Garaud, Khan and Brown2024b) has demonstrated numerically that the presence of a mean vertical shear (ignored in this work) has little impact on these scalings as long as its amplitude is smaller than the small-scale emergent vertical shear 
$U^*/H^*$. This finding can be proved more formally using the asymptotic tools developed here, suggesting that the new regime diagram is robust (at least if other effects such as rotation and magnetic fields are absent; see below).
These results have important implications for low 
$Pr$ fluids, such as liquid metals (where 
$Pr \sim 0.01\unicode{x2013}0.1$), planetary interiors (where 
$Pr \sim 0.001\unicode{x2013}0.1$) and stellar interiors (where 
$Pr \sim 10^{-9}\unicode{x2013}10^{-2}$). In particular, the derived scaling laws enable us to propose simple parameterisations for the turbulent diffusion coefficient (
$D^*$) of a passive scalar in each identified regime by multiplying the characteristic vertical length and vertical velocity scales.
(i) In the
$Pe_b \geq \mathit {O}(1)$ and the intermediate regimes
(5.1)\begin{equation} D^* \propto H^* W^* \propto Fr^{3/2} L^* U^* \propto U^{* \, 5/2}/(L^{* \, 1/2}N^{*3/2}). \end{equation}(ii) In the fully diffusive regime
(5.2)\begin{equation} D^* \propto H^* W^* \propto (Fr^2/Pe)^{1/2} L^* U^* \propto U^{* \, 3/2} \kappa^{*1/2} / (L^{* \, 1/2} N^*). \end{equation}
As such, our work challenges current understanding of stratified turbulence in stars. Indeed, the most commonly used model for stratified turbulence in stellar evolution calculations is the model of Zahn (Reference Zahn1992), which proposes a vertical turbulent diffusivity equivalent to (5.2). However, we have demonstrated that this expression is only valid in a relatively small region of the parameter space (see figure 1). In particular, Zahn's model assumes that the turbulence is always fully diffusive, but this is not the case in the intermediate and 
$Pe_b \geq \mathit {O}(1)$ regimes, which are appropriate for more strongly sheared fluid layers such as the solar tachocline. Zahn's model also assumes that viscosity is negligible, which is not the case for sufficiently stratified flows.
Our multiscale analysis can also be used to estimate the magnitude of the dimensional turbulent buoyancy flux 
$\mathcal {B}^* = |\overline {w'b'}| N^{*2} U^* L^*$, using the characteristic sizes of 
$w'$ and 
$b'$ found in each regime.
(i) In the
$Pe_b \ge O(1)$, $Re_b \gg 1$ regime, $w' = O(\alpha ^{1/2})$ and $b' = O(\alpha ^{3/2})$ with $\alpha = Fr$, so
(5.3)\begin{equation} \mathcal{B}^* \propto \alpha^{2} N^{*2} U^* L^* \propto U^{*3}/L^*. \end{equation}(ii) In the
$Pe_b \ll \alpha, Re_b \gg 1$ regime, $w' = O(\alpha ^{1/2})$ and $b' = O(Pe_b \alpha ^{1/2})$ with $\alpha = (Fr^2/Pe)^{1/3}$ so we also find in this case that
(5.4)\begin{equation} \mathcal{B}^* \propto Pe_b \alpha N^{*2} U^* L^* \propto Pe \alpha^3 N^{*2} U^* L^* \propto U^{*3}/L^*. \end{equation}
This demonstrates that both regimes have a flux coefficient 
$\varGamma = \mathcal {B}^* / \varepsilon ^*$ (where 
$\varepsilon ^* \propto U^{*3}/L^*$ estimates the kinetic energy dissipation rate) and a mixing efficiency 
$\eta = \mathcal {B}^*/(\mathcal {B}^* + \varepsilon ^*)$ that are 
$\mathit {O}(1)$ and independent of stratification (provided the stratification is sufficiently strong that 
$\alpha \ll 1$). This conclusion is consistent with the analysis of Maffioli, Brethouwer & Lindborg (Reference Maffioli, Brethouwer and Lindborg2016) and with the DNS of Cope et al. (Reference Cope, Garaud and Caulfield2020) at 
$Pr \ll 1$.
Further work is needed to investigate the effects of rotation and magnetisation, which are also important in stellar and planetary interiors. Both effects are likely to stabilise the horizontal turbulence to some extent, which will therefore affect the emergent vertical shear instability as well.
Nonetheless, this work has already demonstrated the power of formal multiscale asymptotic analysis for discovering and validating the existence of new regimes of stratified turbulence in stellar and planetary interiors.
Acknowledgements
The authors gratefully acknowledge the Geophysical Fluid Dynamics Summer School (NSF 1829864), particularly the 2022 program.
Funding
P.G. acknowledges funding from NSF AST-1814327. K.S. acknowledges funding from the James S. McDonnell Foundation.
Declaration of interests
The authors report no conflict of interest.
Open access
