The dynamics of stratified horizontal shear flows at low P\'eclet number

We consider the dynamics of a vertically stratified, horizontally-forced Kolmogorov flow. Motivated by astrophysical systems where the Prandtl number is often asymptotically small, our focus is the little-studied limit of high Reynolds number but low P\'eclet number (which is defined to be the product of the Reynolds number and the Prandtl number). Through a linear stability analysis, we demonstrate that the stability of two-dimensional modes to infinitesimal perturbations is independent of the stratification, whilst three-dimensional modes are always unstable in the limit of strong stratification and strong thermal diffusion. The subsequent nonlinear evolution and transition to turbulence is studied numerically using direct numerical simulations. For sufficiently large Reynolds numbers, four distinct dynamical regimes naturally emerge, depending upon the strength of the background stratification: the unstratified turbulent regime; the stratified turbulent regime; the intermittent regime and the viscous regime. By considering dominant balances in the governing equations, we are able to derive scaling laws for each regime which explain the numerical data.


Introduction
Statically stable stratified flows, where the background equilibrium fluid density decreases (at least on average) upwards in a gravitational field, are ubiquitous. Examples in geophysics include atmospheres, oceans and lakes, while they also occur on astrophysical scales in planetary and stellar interiors. A key physical process in such flows is that fluid parcels perturbed vertically from their equilibrium position experience a restoring 'buoyancy force'. Furthermore, it is generic that the fluid will also have a spatially and temporally varying background velocity distribution that is expected to interact with the background 'stable' stratification.
In many cases, this interaction develops into so-called stratified turbulence, key to transport processes in geophysical flows (Ivey et al. 2008;Ferrari & Wunsch 2009) and also thought to play a crucial role in stellar interiors (Zahn 1974(Zahn , 1992Spiegel & Zahn 1992). Stratified turbulence exhibits a wide range of dynamical, and often counter-intuitive behaviours, not least because it leads to complex, and still controversial, irreversible energetic exchange pathways between the kinetic, potential and internal energy reservoirs. Understanding and modelling those pathways, in particular the 'efficiency' of mixing (associated with the irreversible conversion of kinetic energy into potential energy) is of great importance for larger-scale descriptions of geophysical flows, such as weather forecasts, ocean circulation simulation or indeed climate models, and astrophysical flows that regulate planetary and stellar evolution. In what follows, we first describe the current understanding of stratified turbulence in geophysical flows, and explain why these results need to be revisited in the astrophysical context, which is the purpose of this work.

Stratified turbulence in geophysical flows
While the subject of astrophysical stratified turbulence remains in its infancy (see more on this below), there has been a great amount of research into transition, turbulence and mixing in stratified flows with focus on the relevance to atmospheric and oceanic flows (e.g. Ivey et al. 2008). Within this context, the simplest idealised (yet commonly considered) situation has three fundamental modelling assumptions: that the fluid velocity is solenoidal, i.e. ∇ · u = 0; that the density differences within the flow are sufficiently small for the 'Boussinesq approximation' with a linear equation of state to be an appropriate model; and that the density variations in the fluid are associated with a single stratifying agent, avoiding the occurrence of 'double diffusive' phenomena (which may still be very important in a variety of different circumstances, see for example the reviews of Schmitt 1994;Radko 2013;Garaud 2018). Without loss of generality, the density field ρ may be assumed to be a function of temperature T alone, such that where ρ 0 and T 0 are reference densities and temperatures, and α is the thermal expansion coefficient. Since temperature satisfies an advection-diffusion equation where κ is the thermal diffusivity, the density fluctuations also satisfy the same advectiondiffusion equation. Both irreversible mixing and turbulent viscous dissipation, leading respectively to irreversible changes in the potential energy and internal energy of the flow, rely inherently on the action of diffusive processes. Under the three simplifying assumptions above, the stratified fluid under consideration has only two relevant diffusivities: the kinematic viscosity ν quantifying the diffusivity of momentum; and κ, quantifying the diffusivity of density. Together with these diffusivities, there are at least three additional dimensional parameters required to describe a stratified flow: a characteristic velocity scale U c , a characteristic length scale L c , and a characteristic buoyancy frequency N c associated with the background buoyancy frequency profile N b (z), defined as where g is the acceleration due to gravity, and ρ b and T b are background profiles of The dynamics of stratified horizontal shear flows at low Péclet number 3 density and temperature, respectively. Note that the Boussinesq approximation requires that the total variation of a background scalar quantity q b (z) satisfies L c |dq b /dz| q 0 . A natural set of non-dimensional parameters can be constructed as: a Reynolds number Re quantifying the relative magnitude of inertia to momentum diffusion by viscosity; a Péclet number P e quantifying the relative magnitude of inertia to the diffusion of the density; and a Froude number F r quantifying the relative magnitude of the inertia to the stratification, defined as where P r = ν/κ is the Prandtl number. Note that for vertically sheared flows, the Froude number is related to a bulk Richardson number as (1.5) Also note that we have implicitly restricted our focus on a regime where the scales of motion are sufficiently small and fast so that the effects of rotation can be ignored, otherwise an additional parameter is necessary. As discussed in detail in Riley & Lelong (2000) and Brethouwer et al. (2007), oceanic and atmospheric flows are often very strongly stratified, in the specific sense that if both L c and U c are identified with typical scales of horizontal motions, then F r 1 (Ri 1). Nevertheless, turbulence still occurs, at least in spatio-temporally varying patches (Portwood et al. 2016). This has profound implications for understanding the dynamics of such flows. Brethouwer et al. (2007), following Billant & Chomaz (2001), demonstrated that when both Re 1 and F r 1, several different flow regimes are possible. Each regime can be understood as a distinct dominant balance between various terms in the Navier-Stokes equations, dependent on their relative sizes. Of central significance to these balances, however, are two additional geophysically-motivated parameter choices, both of which we wish to revisit in this manuscript which aims to extend this work to astrophysicallyrelevant flows. The first of these parameter choices is motivated by the expectation (and empirical observation) that 'strong' stratification leads to anisotropy in the flow, so the characteristic vertical length scales L v are expected to be very different from characteristic horizontal length scales L h ≡ L c . The second relies on the fact that the Prandtl number is of order unity or larger in geophysical flows. Typically, P r ∼ O(1) for gases (e.g. P r 0.7 for air) while for fresh water P r ∼ O(10) (with some variability with temperature and pressure, although the canonical value is chosen to be 7). If the density variations are due to salinity with diffusivity D rather than temperature differences, the analogous ratio of diffusivities, known as the Schmidt number Sc = ν/D ∼ O(1000), is even higher.
With these two further choices, Brethouwer et al. (2007) discussed three particular regimes which are worthy of comment. The first, originally considered by Lilly (1983) (also see Riley & Lelong (2000) for further discussion) has Re 1 and F r 1, yet L v /L c F r and also L v /L c 1/ √ Re. With these scalings, all terms involving vertical derivatives (specifically diffusive terms and advective terms involving vertical velocity) are insignificant in the Navier-Stokes equations, and so the governing equations reduce to the evolution equations for an incompressible and inviscid 'two-dimensional' horizontal velocity u h (x, y, t). Furthermore, since P r O(1), diffusive terms in the density equation can also be ignored, and quasi-two dimensional (though possibly layerwise) flow evolution is expected.
The other two regimes discussed in detail by Brethouwer et al. (2007) still rely essentially on the fact that P r O(1). They also exploit the insight of Billant & Chomaz (2001) that the vertical length scale should not be externally imposed, but should emerge as a property of the flow dynamics. In that respect, as presented in detail by Brethouwer et al. (2007), a key parameter is the quantity commonly referred to as the 'buoyancy Reynolds number' Re b , defined as (1.6) When Re b O(1), (but still with P r O(1)), a viscously affected regime is expected, where horizontal advection is balanced by viscous diffusion, specifically associated with vertical shearing. This regime, much more likely to be relevant in experiments (or simulations) rather than in geophysical applications, has L v /L h ∼ 1/ √ Re, and does not exhibit a conventional turbulent cascade, but rather exhibits the effects of viscosity (and density diffusion) even at large scales.
Conversely, Brethouwer et al. (2007) showed that when Re b 1, viscous effects are insignificant (as is density diffusion since P r O(1)) and the remaining terms (including the advection by the vertical velocity) become self-similar with respect to zN c /U c , with z being the vertical coordinate aligned with gravity. This suggests strongly that L v ∼ U c /N c , or equivalently that the Froude number based on the vertical scale L v , defined as should be of order one, so L v L c . Such a vertical layer scale has been commonly observed in a wide variety of sufficiently high Reynolds number stratified flows (e.g Park et al. 1994;Holford & Linden 1999;Billant & Chomaz 2000;Godeferd & Staquet 2003;Brethouwer et al. 2007;Oglethorpe et al. 2013;Lucas et al. 2017;Zhou & Diamessis 2019) and appears to be a generic property of high Re b and high P r stratified turbulence. This regime is characterised not only by anisotropic length scales but also by anisotropy in the velocity field, and hence the associated turbulence, leading Falder et al. (2016) to refer to this flow regime as the 'layered anisotropic stratified turbulence' (LAST) regime.
The vertical layering on the scale L v is key to understanding how turbulence can be maintained in the LAST regime despite the fact that F r 1. Indeed, these 'layers' in the density distribution consist of relatively weakly stratified wider regions separated by relatively thinner 'interfaces' with substantially enhanced density gradient. As such, local values of the buoyancy frequency can vary widely from the characteristic value N c . When the local vertical shear is sufficiently strong compared to the local density gradient, then the gradient Richardson number Ri g , defined as can drop to values low enough for shear instabilities to be able to develop. If in addition the Reynolds number is sufficiently large for inertial effects to be dominant, this allows for the possibility of turbulence through the breakdown of shear instabilities, albeit with both spatial and temporal intermittency. It is crucial to appreciate that this LAST regime relies inherently on the assumption that P r O(1), as high Reynolds number thus implies high Péclet number, so localized turbulent events can erode the stratification and in turn participate in the formation or maintenance of the layers. Although appropriate for the atmosphere and the ocean, this fundamental assumption most definitely does not apply in the astrophysical context, The dynamics of stratified horizontal shear flows at low Péclet number 5 where P r 1 (see below). As we now demonstrate, density layering is prohibited in that case, suggesting that LAST dynamics cannot occur. This raises the interesting question of whether analogous or fundamentally different regimes exist when P r 1.

Stratified turbulence in astrophysical flows
The fluid from which stars and gaseous planets are made is a plasma comprised of photons, ions, and free electrons. As a result, one of the main differences between astrophysical and geophysical flows is the value of the Prandtl number, which is much smaller than one as mentioned above. In typical stellar radiative zones, for instance, Pr usually ranges between 10 −9 and 10 −5 (e.g. Garaud et al. 2015b). The microphysical explanation for this difference is that heat can be transported by photons efficiently while momentum transport usually requires collisions between ions (which comprise most of the mass), so ν κ and Pr 1. This crucially introduces the possibility of a new regime of flow dynamics where Re 1 while P e = Pr Re 1, which is never realized in geophysics.
Astrophysical fluids are also not incompressible. However, under a set of assumptions that are almost always satisfied sufficiently far below the surface of stars and gaseous planets, the Spiegel-Veronis-Boussinesq approximation (Spiegel & Veronis 1960) can be used to reduce the governing equations to a form that is almost equivalent to that used for geophysical flows. In particular, where c p is the specific heat at constant pressure. In comparison with (1.2), the new term wg/c p accounts for compressional heating (or cooling) as the parcel of fluid contracts (or expands) to adjust to the ambient pressure as it moves downwards (or upwards) in a gravitational field. As a result, the background buoyancy frequency profile N b (z) is modified from (1.3) to (1.10) from which a new characteristic buoyancy frequency N c can be defined. Interest in stratified shear instabilities at low Prandtl number and/or low Péclet number in stars dates back to Zahn (1974). In this regime, thermal dissipation greatly mitigates and modifies the effect of stratification in comparison to flows with P r O(1). In particular, as demonstrated by Lignières (1999) (see also Spiegel 1962;Thual 1992), a dominant balance emerges in the temperature equation whereby (at least to leading order in P e −1 ), showing that temperature fluctuations and vertical velocity fluctuations are slaved to one another (see more on this in section 2). Mass conservation, combined with appropriate boundary conditions, then generally implies that the horizontal average of T should be zero. Physically, this simply states that due to the very rapid diffusion of the temperature fluctuations (and hence density), perturbations cannot modify the background. Density layering is therefore prohibited, as stated above, so the local buoyancy frequency remains close to the characteristic value N c everywhere. Another important consequence of this highly diffusive limit (Lignières 1999) is that the Péclet and Froude (or Richardson) numbers are no longer independent control parameters for the system dynamics, but always appear together as P e/F r 2 or RiP e. Zahn (1974) argued that, as a result, the threshold for vertical shear instability should be RiP e Re/Re c where Re c is the critical Reynolds number for instability in unstratified, unbounded shear flows (which he estimated would be O(1000)). Zahn's criterion for instability is now commonly written as RiPr K Z , where K Z ∼ O(10 −3 ). This was recently independently confirmed using direct numerical simulations (DNSs) by Prat et al. (2016) (see also Prat & Lignières 2013, 2014 and Garaud et al. (2017), who found that K Z 0.007. With the aforementioned estimates for Pr , we see that shear-induced turbulence in low P e vertical shear flows is therefore likely for Ri up to ∼ 10 2 or higher. On the other hand, for astrophysical flows with RiPr K Z , or for horizontally-sheared flows (see below), one may naturally ask whether any pathway to turbulence exists, since the density layering that is central to the LAST regime is not possible here. This paper aims to answer this question for the case of horizontally-sheared flows.
Before proceeding, however, it is useful to briefly review the most commonly used model of stratified turbulence in stars (see Lignières 2018, for a more comprehensive review of the topic). Zahn (1992) considered successively both vertically-sheared flows and horizontally-sheared flows. For a vertically-sheared flow with characteristic shearing rate S c , he argued based on work by Townsend (1958) and Dudis (1974) that the largest unstable vertical scale in the flow would satisfy RiP e l ∼ O(1), where here Ri = N 2 c /S 2 c and where P e l ≡ S c l 2 /κ is an eddy-scale Péclet number. This defines the characteristic Zahn scale l Z as (1.12) Using dimensional analysis, Zahn (1992) then proposed a simple expression for a turbulent diffusion coefficient, namely (1.13) The relevance of the Zahn scale to the dynamics of low Péclet number stratified turbulence in vertically-sheared flows was confirmed by Garaud et al. (2017) using direct numerical simulations (DNSs). They also verified that (1.13) applies for flows that have both low Péclet number and sufficiently high Reynolds number, as long as l Z is much smaller than the domain scale, and RiPr K Z (see also Prat & Lignières 2013, 2014Prat et al. 2016).
In the horizontally-sheared case, Zahn (1992) postulated (following an argument attributed to Schatzman & Baglin 1991), that while the turbulence would be mostly twodimensional on the large scales owing to the strong stratification, it could become threedimensional below a scale L c where thermal dissipation becomes important. This scale is by definition the Zahn scale, and is therefore given by (1.12) where here S c = U c /L c (and U c is the characteristic velocity of eddies on scale L c ). Since Pr 1, this scale is also unaffected by viscosity, so one would expect a turbulent cascade with well-defined kinetic energy transfer rate of order U 3 c /L c . If, in addition, dissipative irreversible conversions into the potential energy reservoir are negligible, then U 3 c /L c = ε where ε is the viscous energy dissipation rate. Solving for L c and U c yields (see Lignières 2018, for an alternative derivation of these scalings): (1.14) The dynamics of stratified horizontal shear flows at low Péclet number 7 from which a turbulent diffusion coefficient can then be constructed as (1.15) The Zahn (1992) model for stratified turbulence driven by horizontal shear at low Péclet number and/or low Prandtl number has, to our knowledge, never been tested. In addition to verifying (1.14) and (1.15), we are also interested in testing the assumption that all energy dissipation is exclusively viscous. Although this assumption is superficially plausible, there is growing evidence (Maffioli et al. 2016;Garanaik & Venayagamoorthy 2019) for flows with P r O(1) that non-trivial irreversible mixing converting kinetic energy into potential energy continues to occur even in the limit F r → 0 of extremely strong stratification.
In what follows we therefore study the simplest possible model of a stratified horizontal shear flow, and focus in this paper on the limit where thermal diffusion is important, or equivalently, the low Péclet number limit. As we shall demonstrate, distinct parameter regimes naturally emerge, each with its own characteristic scalings, determined not through heuristic physical reasoning but through the identification of different dominant balances in the relevant governing equations. Section 2 presents the model setup, and section 3 summarizes the results of a linear stability analysis of the problem. Section 4 describes the results of a few characteristic simulations and identifies four separate regimes. These are then systematically reviewed in section 5, where we study the dominant balances for each regime and derive pertinent scaling laws that are then compared with the numerical data. We discuss these results and draw our conclusions in section 6.

Mathematical model
We consider an incompressible, body-forced, stably stratified flow with streamwise velocity field aligned with the x-axis. In accordance with the Spiegel-Veronis-Boussinesq approximation (Spiegel & Veronis 1960), we assume that the basic state comprises a linearized temperature distribution T b (z) given by T b (z) = T 0 + z(dT b /dz), where T 0 is a reference temperature, along with a body-forced laminar velocity field u L (y). The total temperature field, T , includes perturbations T (x, y, z, t) away from the basic state such that T = T b (z) + T (x, y, z, t). As discussed in section 1.2, the density fluctuations ρ and temperature fluctuations T are related according to the linearized equation of state where ρ 0 is a reference density and α = −ρ −1 0 (∂ρ/∂T ) is the coefficient of thermal expansion. The three-dimensional velocity field is given by u(x, y, z, t) = ue x +ve y +we z . For numerical efficiency, we impose triply-periodic boundary conditions on the body force F and the variables T and u such that (x, y, z) ∈ [0, L x ) × [0, L y ) × [0, L z ). A suitable candidate for the applied force is a monochromatic sinusoidal forcing driving a horizontal Kolmogorov flow: This choice of forcing is computationally straightforward to implement and also reveals interesting dynamics (see Lucas et al. 2017). Figure 1 illustrates the basic laminar state. The governing Spiegel-Veronis-Boussinesq equations (Spiegel & Veronis 1960) for this model setup are: where ν is the kinematic viscosity, κ is the thermal diffusivity, χ is the forcing amplitude, p is the pressure, c p is the specific heat at constant pressure and gravity g acts in the negative z-direction. In this study, we specify that L y = L z while L x may vary continuously such that the aspect ratio of the domain is given by λ = L x /L y . λ > 1 corresponds to domains which are longer in the streamwise direction.

Non-dimensionalization and model parameters
In equilibrium, we anticipate a balance between the body force and fluid inertia such that u · ∇u ∼ χ sin (2πy/L y ) e x in the streamwise direction. For a characteristic length scale L y /2π, this gives a characteristic velocity scale χL y /2π and a characteristic time scale L y /2πχ. Combined with the vertical temperature gradient scale (dT b /dz + g/c p ), we use the equivalent non-dimensionalization as in Lucas et al. (2017) to give the following system of equations, in which all quantities are non-dimensional: The dynamics of stratified horizontal shear flows at low Péclet number 9 We thus have three non-dimensional numbers: the Reynolds number Re; the buoyancy parameter B; and the Prandtl number P r, which determine the dynamics of the system: where N b is the dimensional buoyancy frequency defined in (1.10), which is now constant by construction. Note that B is related to the Froude number as (2.10) It is also convenient to introduce the Péclet number P e, defined as (2.11) Both sets of parameters, (Re, B, P r) or (Re, B, P e), uniquely define the system and will be used interchangeably throughout this study. In all that follows, the domain is a cuboid such that ( 2π), and variables p, T and u have triply-periodic boundary conditions. This system, defined by (2.6), (2.7) and (2.8), will henceforth be referred to as the standard system of equations.

Low-Péclet number approximation
As discussed in section 1.2, when the thermal diffusion timescale is much shorter than the advective timescale, a quasi-static regime is established where temperature fluctuations are slaved to the vertical velocity field. Motivated by astrophysical applications such as stellar interiors, we consider the standard set of equations (2.6)-(2.8) in the asymptotic limit of low Péclet number (LPN). This limit was studied by Spiegel (1962) and Thual (1992) in the context of thermal convection, and more recently by Lignières (1999) in the context of stably stratified flows. Lignières proposed that the standard equations can be approximated by a reduced set of equations called the "low-Péclet number" equations (LPN equations hereafter), in which the density fluctuations are slaved to the vertical velocity field: (2.14) These can be derived by assuming a regular asymptotic expansion of T in powers of P e, i.e. T = T 0 + T 1 P e + O(P e 2 ), and by assuming that the velocity field is of order unity. At lowest order (P e −1 ), we get ∇ 2 T 0 = 0 implying that T 0 = 0 is required to satisfy the boundary conditions, while at the next order (P e 0 ), the equations yield w = ∇ 2 T 1 ≈ P e −1 ∇ 2 T as required.
Noting that (2.13) can be re-written formally as T = P e∇ −2 w, we derive the reduced set of LPN equations:

15)
∇ · u = 0. (2.16) These equations explicitly demonstrate that under the LPN approximation (and in contrast to the standard equations), there are only two non-dimensional parameters governing the flow dynamics, notably the Reynolds number Re and the product of the buoyancy parameter and the Péclet number, BP e = P eF r −2 . This combined parameter, which we consider to be a measure of the stratification, can take any value (even for small Péclet numbers) because B can be arbitrarily large, or equivalently F r can be arbitrarily small, as the stratification becomes strong.
There are advantages of studying the LPN equations rather than the standard equations. For example, this reduced set of equations allows for the derivation of mathematical results such as an energy stability threshold that explicitly depends on BP e (see Garaud et al. 2015a). Throughout this study, we will discuss both systems of equations, verifying the validity of the LPN equations where possible.

Standard equations
We begin by considering the stability of a laminar flow to infinitesimal perturbations, with initial focus on the standard set of equations (2.6)-(2.8). The background flow u L (y), which satisfies Re −1 ∇ 2 u L + sin(y)e x = 0, is given by (3.1) Note that if one wishes to consider a basic state with generic amplitude aRe instead of amplitude Re, it is straightforward to apply a rescaling using the method described in appendix A. For small perturbations u (x, y, z, t) away from this laminar flow, i.e. letting u = u L (y) + u (x, y, z, t), the linearised perturbation equations are: In this set of partial differential equations (PDEs), the coefficients are periodic in y but independent of x, z and t. Consequently, and in the conventional fashion, we consider normal mode disturbances of the form: where q ∈ (ũ ,ṽ ,w , T , p) and k x and k z are the perturbation wavenumbers in the x and z-directions respectively. The geometry of the model set-up requires that k x ∈ R and k z ∈ Z. We seek periodic solutions forq(y) given bŷ Substituting this ansatz into equations (3.2)-(3.4) and using the orthogonality property of complex exponentials, we obtain a 5 × (2L + 1) = (10L + 5) algebraic system of equations for the u l , v l , w l , T l , p l for l ∈ (−L, L): This system can be re-formulated as a generalised eigenvalue problem for the complex growth rates σ, .., p L ), A and B are (10L + 5) × (10L + 5) square matrices and B i,j = {δ ij , i, j (8L + 4); 0, otherwise}. Equation (3.12) has (10L+5) eigenvalues σ. For perturbation wavenumbers k x and k z and system parameters Re, B and P r, the eigenvalue with the largest real part determines the growth rate of the linear instability. The eigenvalue problem can be solved numerically, with L chosen such that convergence is achieved. For example, for a 2D perturbation such that k x = 1 and k z = 0 and system parameters Re = 100, B = 100 and P r = 1, we find that max(Re{σ}) = σ = −0.01, demonstrating that this particular mode is stable to infinitesimal perturbations.
3.1.1. Comparison with previous results at P r = 1 We first consider the case of P r = 1 for ease of comparison with previous work. Both Deloncle et al. (2007) and Arobone & Sarkar (2012) considered the linear stability of horizontal shear layers with somewhat different base flows, and Lucas et al. (2017) considered the linear stability of the specific horizontally-sheared Kolmogorov flow considered here, exclusively for P r = 1. Letting B = 100, we consider the linear stability of the basic state flow u L (see (3.1)) across a range of Reynolds numbers for both 2D and 3D perturbation modes.
Figure 2(a) shows the neutral stability curves (σ = 0) for varying vertical wavenumbers k z ∈ (0, ..., 6) in the (Re, k x ) space. Our results are in good agreement with those of Lucas et al. (2017). Stability (σ < 0) is found to the left and above the curves while instability (σ > 0) occurs to the right and below. The black curve illustrates the 2D (k z = 0) mode. This neutral stability curve intercepts the x-axis when Re = 2 1/4 1.19, implying that the system is linearly stable when Re < 2 1/4 (in agreement with Beaumont 1981;Balmforth & Young 2002, once the correct rescaling is applied (see section 3.3)). For large Re, it asymptotes to k x = 1 but, in agreement with Lucas et al. (2017), always lies below this line, leading to the conclusion that domains such that λ = L x /L y 1 are linearly stable to the 2D mode.
The coloured curves show the neutral stability curves for the first six 3D modes (k z ∈ (1, ..., 6)). The onset of instability in the 3D modes is found to occur for higher Reynolds numbers than the 2D mode, with the critical Reynolds number for instability of these 3D modes increasing monotonically with increasing k z . For a range of Re ∼ O(100) (corresponding to P e ∼ O(100)), the 3D curves actually cross the line k x = 1 implying that these modes are unstable for domains where λ = 1, i.e. cubic domains.
Figures 2(b)-(c) further analyse the information in figure 2(a) by computing, for each Reynolds number and k z , the largest (positive) growth rate across all values of k x , σ max , and the value of k x for which that maximum is achieved, k x,max . We see that, as well as being the mode that becomes unstable first, the 2D mode is always the fastest growing one. In addition, the ratio of the growth rate of the 2D mode to that of the 3D modes increases with Re. We therefore predict that the 2D mode would strongly influence the dynamics when it is unstable (i.e. for domain sizes such that λ > 1). Finally, we note Variation with Reynolds number for a collection of kz wavenumbers of: (b) the largest growth rate σmax maximised across all horizontal wavenumbers kx; (c) the associated horizontal wavenumber kx,max. The curves plotted include kz = 0 (black) and kz = 1, 2, 3, 4, 5, 6 (coloured) and the standard equations were used with B = 100 and P r = 1 fixed.
that the corresponding streamwise wavenumbers of the fastest growing 3D modes satisfy k x,max → 0 in the limit Re → ∞, while those of the fastest growing 2D mode remain constant.

Stability at low P r
Astrophysical applications motivate an understanding of the effects of the stratification parameter B and Prandtl number P r on the linear stability of the basic state. Consequently, in figures 3(a)-(c) we plot the neutral stability curves in exactly the same fashion as in figure 2(a), for three different Prandtl numbers: P r = 0.1 (first column), P r = 0.01 (second column) and P r = 0.001 (third column), keeping B = 100 constant. Whilst the neutral stability curves for the 2D mode are identical, clear trends exist for the 3D modes. A reduction in the value of P r shifts the critical Reynolds numbers for the onset of instability of the 3D modes towards higher values, thereby making these modes less unstable. This result is consistent with Arobone & Sarkar (2012), who investigated the stability of a stratified, horizontally-sheared hyperbolic flow. We also note that the same trend is found by letting B → 0 and keeping P r constant (not plotted). Thus, B → 0 (at fixed P r) and P r → 0 (at fixed B) have the same effect: the 3D modes of instability are suppressed while the 2D mode remains unstable. The explanation for this emerges from consideration of (2.7). As the Prandtl number tends to zero (keeping the Reynolds number finite), the Péclet number becomes small and so the buoyancy diffusion becomes important. In this case, a parcel of fluid that is advected into surrounding fluid of a different density adjusts very rapidly to its surroundings, thereby reducing the buoyancy force and so approximating an unstratified system. However, it is important to note that another distinguished limit exists in which B → ∞ and P r → 0, while the product BP r remains finite. This limit is relevant to stellar interiors, and behaves quite differently from the case where B is fixed while P r → 0, as we now demonstrate.
The dynamics of stratified horizontal shear flows at low Péclet number

Low-Péclet number equations
We now examine the linear stability of the LPN equations, given by equations (2.15) and (2.16). We follow the same steps as in the previous section, however we find ourselves working this time with a reduced set of four equations rather than five. We obtain a 4 × (2L + 1) = (8L + 4) algebraic system of equations for the u l , v l , w l , p l for l ∈ (−L, L): As before, this can be re-formulated as a generalised eigenvalue problem for the complex growth rates σ,  In order to test the validity of the LPN equations, we first compare the results of the linear stability analysis in the LPN limit to that obtained using the standard equations. Figure 3 illustrates this comparison. The top row (as already discussed) shows the neutral stability curves from the standard equations and the bottom row shows the equivalent results from the LPN equations. The value of BP r = BP e/Re in the bottom row decreases from left to right by two orders of magnitude, in line with reductions in the value of BP e at fixed Re in the standard equations in the top row. As demonstrated in section 2.3, the LPN equations are asymptotically correct in the limit where P e → 0. Figure 3 shows that they remain valid up to P e 0.1 (i.e. within the regions shown in grey). Outside of these regions, increasingly large differences emerge, especially as P e increases above one. In particular, the neutral stability curves for the 3D modes never cross the line k x = 1 in the LPN equations, suggesting that horizontal shear instabilities do not arise for cubic domains when the LPN approximation is used.
The LPN system of equations depend on the combined parameter BP r = BP e/Re. In the limit of strong stratification (B → ∞) and strong thermal diffusion (P r → 0), this parameter remains finite and is not necessarily small. As can be seen in figure 3, the 3D modes remain unstable in this limit, in agreement with the results from the standard system of equations.
We now focus on the case when BP r = 1. By way of comparison with the standard equations at P r = 1, figures 4(b)-(c) show, for each Reynolds number and k z wavenumber, the largest (positive) growth rate across all values of k x , σ max , and the value of k x for which that maximum is achieved, k x,max . As before, we observe that the 2D mode is both the first mode to become unstable, and is always the fastest growing mode. There are, however, two significant differences between high and low Prandtl number dynamics. Firstly, in the LPN limit, figure 4(b) shows that the growth rates of the fastest growing 3D modes increase in line with those of the fastest growing 2D mode. Secondly, the corresponding values of k x,max remain constant as Re → ∞. Consequently, the 3D modes remain important relative to the 2D mode and we therefore predict that, in contrast to the case when P r = 1, both the 2D and 3D modes would strongly influence the dynamics in this limit. These results, combined with the fact that the 3D modes remain unstable in the limit of strong stratification and strong thermal diffusion, have important consequences, as we shall see section 4.

Critical Reynolds number for instability
An important finding in this study, determined numerically across a broad spectrum of parameters, is the fact that the critical Reynolds number, Re c , for the onset of linear instability, as given by the 2D mode, is independent of both the stratification and Prandtl numbers, being fixed at Re c = 2 1/4 . This result holds for both the standard equations and the LPN equations and differs quite substantially from that obtained in Garaud et al. (2015a) for the case of a vertically orientated shear, where stratification was found to be able to stabilize a system.
To see why this is the case, we observe in equations (3.7)-(3.11) (or the equivalent LPN equations (3.13)-(3.16)) that setting k z = 0 reduces the problem to the study of equations (3.7), (3.8) and (3.11) (or equivalently (3.13), (3.14) and (3.16)). This reduced problem is well studied (Beaumont 1981;Balmforth & Young 2002), being the linear stability of an unstratified (B = 0) flow. The critical Reynolds number for instability around a basic state of u L (y) = sin(y)e x has been shown to be √ 2 (Beaumont 1981). As detailed in appendix A, a simple transformation given by relations (A 15) (for a = Re −1 ) gives the corresponding critical Reynolds number for 2D modes in this study to be 2 1/4 , which corresponds to the result obtained numerically. Consequently, we note that horizontally-sheared Kolmogorov flows with Re > 2 1/4 are always unstable, irrespective of the stratification, and thus form a convenient basis from which to study the subsequent nonlinear evolution of stratified, low-Péclet number flows.

Direct numerical simulations
We now present results from a series of DNSs of horizontal shear flows at low Péclet number following the model setup and equations described in section 3. As we shall demonstrate, the system presents a rich ecosystem of instabilities that feed on each other, leading to a number of distinct dynamical regimes that will be further characterized in section 5.

Numerical algorithm
The DNSs are performed using the PADDI code first introduced by Traxler et al.
and Stellmach et al. (2011) to study double-diffusive fingering. The code has since then been modified to study many different kinds of instabilities, including bodyforced vertical shear instabilities, using both the standard equations and the LPN approximation (Garaud et al. 2015a;Garaud & Kulenthirarajah 2016;Gagnier & Garaud 2018;Kulenthirarajah & Garaud 2018). PADDI is a triply-periodic pseudo-spectral algorithm that uses pencil-based Fast Fourier Transforms, and third order backwarddifferencing Adams-Bashforth adaptive timestepping (Peyret 2002) in which diffusive terms are treated implicitly while all other terms are treated explicitly. The velocity field is made divergence-free at every timestep by solving the relevant Poisson equation for the pressure. Two versions of the code exist, one that solves the standard equations (2.6)-(2.8), and one that solves the LPN equations (2.15)-(2.16).
Based on the linear stability analysis performed in section 3, we have selected a domain size such that L y = L z = 2π and L x = 4π. This allows for the natural development of a single 2D mode of instability (for which k x = 0.5), without being computationally prohibitive at high Reynolds number (see below). A comparison of simulation outcomes for different domain lengths is presented in Cope (2019, section 4.2) but only for Re = 50, for which only two dynamical regimes exist. A systematic exploration of the effect of domain aspect ratio at high Reynolds number will be the subject of future work.
Tables 1 and 2 present all the runs that have been performed with this model setup, using equations (2.6)-(2.8) and (2.15)-(2.16), respectively. The number of Fourier modes used in each direction (after dealiasing) depends on the Reynolds number Re selected; these are summarized in tables 1 and 2. The same resolution is used regardless of the values of B and P e. To save on computational time, only one of the simulations at each Reynolds number is initiated from the original initial conditions (i.e. u = sin(y)e x plus some small amplitude white noise). All the others are restarted from the end point of a simulation at nearby values of P e or B. In all cases, we have run the simulations until they reach a statistically stationary state, except where explicitly mentioned. Note that for very large values of B or very small values of P r, we have found it necessary to decrease the value of the maximum allowable timestep substantially. This is because the system of equations becomes increasingly stiff and is otherwise susceptible to the development of spurious elevator modes (i.e. modes that are invariant in the vertical direction). To save on computational time, we only ran simulations using the standard equations in that limit.

Typical simulations: early phase
We begin by presenting the early phases of development of the horizontal shear instability, in two typical simulations at moderately large Reynolds number (Re = 300), high stratification (B = 30 000 and 300 000, respectively) and relatively low Péclet number (P e = 0.1). Both simulations were initialized with u = sin(y)e x plus small amplitude white noise. Figures 5(a) and 5(d) show the root mean square values of the streamwise (u rms ), spanwise (v rms ) and vertical (w rms ) velocities for each simulation, computed at each instant in time as q rms (t) = q 2 1/2 , (4.1) where the angular brackets denote a volume average such that For both values of B, we clearly see the growth of the streamwise flow due to the forcing. Spanwise and vertical fluid motions first decay, until the onset of the 2D mode of instability (i.e. whereby v rms begins to grow while w rms continues to decay), rapidly followed by the 3D mode, for which w rms finally also begins to grow. Snapshots of the streamwise velocity fields near the saturation of these instabilities are presented in figures 5(b)-(c) and 5(e)-(f). In both cases, the snapshot at t 1 illustrates the early development of the 2D and 3D modes of instability. The 2D mode causes a meandering of the background flow, and the 3D mode causes a vertical modulation of the position of the meanders. We also see that the 3D mode has a substantially smaller vertical scale for larger B. The snapshot at t 2 shows how the instability further evolves with time: the meanders and their vertical shifts both grow in amplitude, leading to the development of substantial vertical shear of the streamwise flow.
While similar early-time dynamics are observed at all parameter values (assuming the 2D mode is unstable), what happens beyond that depends on Re, B and P e. We now describe in turn the various regimes that can be found.  Table 2. Summary of all the runs obtained using the low Péclet number equations, with parameters Re, and BP e. Quantities in columns 3-5 are computed in the manner described in section 4.4. In terms of equivalent grid points, the resolution used is 576×288×288 (Re = 100 and Re = 300 runs), and 768×384×384 (Re = 600 runs).

Typical simulations: nonlinear saturation
The nonlinear saturation of this body-forced horizontal shear flow depends crucially on the selected value of the stratification parameter B. In what follows, we investigate the effect of varying B. Snapshots of the streamwise velocity, vertical velocity and local viscous dissipation rate, taken during the statistically stationary state, are presented in figure 6. In all but the last row, Re = 300, and P e = 0.1. For the last row, Re = 50.
For very large values of B (bottom row in figure 6), the vertical scale of the 3D mode of instability is relatively small. Even though substantial shear develops between successive meanders of the streamwise jets, this shear is too small to overcome the stabilizing effect of viscosity, and remains stable. The resulting flow takes the form of thin layers, crucially in the velocity field, each of which presents a meandering jet with its own distinct phase. These jets are weakly coupled in the vertical direction through viscosity. The vertical velocity field is non-zero, however, and is presumably generated by the weak horizontal divergence of the flow within each jet.
As B decreases (i.e. moving up in figure 6), the reduced stratification now allows for the development of secondary vertical shear instabilities between the meanders, albeit only intermittently. Spatially localized overturns can be seen in figures 6(g)-(i). These become more numerous and more frequent as B continues to decrease. The viscous dissipation is clearly enhanced in the turbulent regions compared with the laminar regions.
For intermediate values of B (see figures 6(d)-(f)), the flow becomes fully turbulent. The vertical scale of the eddies remains relatively small, however, consistent with stratification playing a role in shaping the dynamics of the turbulence. The meandering streamwise jets are still clearly visible. The dissipation rate snapshot shows that the scale of the turbulent eddies is small in both horizontal and vertical directions.
Finally, for low values of B, the scale of the eddies is now the domain scale, and the turbulence is unaffected by stratification. In fact, this system is very similar to the one obtained in weakly stratified vertically sheared flows (see Garaud & Kulenthirarajah 2016), except for the horizontally-averaged mean flow (which varies with y instead of z).
These observations therefore suggest the existence of at least four distinct dynamical regimes: unstratified turbulence for very low B; stratified turbulence for intermediate values of B; intermittent turbulence for higher values of B; and finally, viscouslydominated stratified laminar flow for the highest values of B. We will now proceed to characterize these different regimes more quantitatively.

Data extraction
Each of the simulations we have performed was integrated until the system had reached a statistically stationary state. This can take a long time, especially for the very strongly stratified systems, so data in that limit is scarce except for the lowest values of Re. Once in that statistically stationary state, we compute the time average, and deviations around that average, of w rms (t) and T rms (t), where q rms (t) for any quantity q was defined in (4.1). These are reported as w rms ± δw rms and T rms ± δT rms , respectively, in tables 1 and 2.
We also compute the temperature flux as for DNSs that use the standard equations and P e −1 F T (t) = w∇ −2 w (4.4) for DNSs that use the LPN equations, where the angular bracket was defined in (4.2). We finally compute the viscous energy dissipation rate as (t) = Re −1 |∇u| 2 . (4.5) Note that is the non-dimensional version of ε introduced in section 1.2. This shows that the rate at which the body force does work on the flow, u sin(y) , is The dynamics of stratified horizontal shear flows at low Péclet number 21 Figure 7. Autocorrelation function Aw(l, t) as defined in (4.10) computed at six randomly selected times during the statistically stationary state of a simulation with parameters Re = 300, B = 10 000 and P e = 0.1. Note how Aw(l, t) has a well-defined first zero, whose time-average defines the vertical eddy scale lz.
partitioned between energy that is dissipated viscously (through ), and energy that is converted into potential energy (at a rate BF T ). The fate of the latter can be established by multiplying the temperature equation by T and integrating over the domain, which reveals that ∂ ∂t for the full equations (while in the LPN limit, the time derivative simply disappears). This shows that BF T is ultimately dissipated thermally at a rate P e −1 |∇T | 2 . From these considerations, it is common to define a so-called instantaneous mixing efficiency (see e.g. Maffioli et al. 2016) at a given point in time, which measures the efficiency with which kinetic energy, produced by the applied forcing, is converted into potential energy as opposed to being dissipated viscously. We have computed η(t) for all simulations produced, and report its time average and deviation from that average, while in a statistically stationary state, as η ± δη in tables 1 and 2. Finally, another useful diagnostic of the flow is the typical vertical scale of the turbulent eddies. As discussed in Garaud & Kulenthirarajah (2016) and Garaud et al. (2017), there are many different ways of extracting such a length scale, either from weighted averages over the turbulent energy spectrum, or from spatial autocorrelation functions of the velocity field. Garaud et al. (2017) compared these different methods and concluded that the spatial autocorrelation function was a more physical and reliable way of extracting the vertical length scale. In what follows, we therefore compute the function A w (l, t) = 1 L x L y L z w(x, y, z, t)w(x, y, z + l, t)dxdydz (4.10) at each timestep for which the full fields are available, using periodicity of w to deal with points near the domain boundaries. Sample functions for six randomly selected times during the statistically stationary state are shown in figure 7 for a simulation with parameters Re = 300, B = 10 000 and P e = 0.1 (a simulation from the stratified intermittent regime, snapshots from which are shown in figures 6(g)-(i)). We clearly see that A w (l, t) has a well-defined first zero at each timestep, which we call l z (t). The vertical eddy scale thus obtained is then averaged over all available timesteps during the statistically stationary state to obtain the mean vertical eddy scale l z and its standard deviation δl z .

Nonlinear saturation: scaling regimes
In our quest for a quantitative description of the four dynamical regimes described in section 4.3, we endeavour to derive scaling laws that explain our data, in an analogous fashion to the approach of Brethouwer et al. (2007) in which the focus was on geophysically-relevant parameters (P r O(1)). Consistent with our goal to study systems in which the Péclet number P e is small, we have run a range of simulations using the standard equations (2.6)-(2.8) for three different Péclet numbers (0.01, 0.1 and 1), which we compare alongside simulations using the LPN equations (2.15) and (2.16), noting excellent agreement. Using both sets of equations, we consider four different Reynolds numbers (50, 100, 300, 600) and investigate a wide range of background stratifications.

Effects of stratification on mixing and the vertical scale of eddies
The first flow diagnostic that we discuss is the vertical eddy scale l z , computed using the method described in section 4.4. Figure 8(a) shows l z as a function of BP e, consistent with our expectations on the potential relevance of this parameter for low Péclet number flows (as discussed in section 2.3). For all but the largest values of BP e (which corresponds to the viscous regime discussed in section 4.3), we confirm that BP e is indeed the relevant parameter, and that l z is independent of Re. As a result, all the data collapses on a single universal curve. For weak stratification, which we refer to as the unstratified regime, the vertical eddy scale is invariant with respect to both stratification and Reynolds number. We find that l z 2, which is of the order of the size of the periodic domain. For intermediate values of BP e, corresponding to the stratified turbulent regime described in section 4.3, we find that l z 2(BP e) −1/3 . (5.1) Finally, for very strong stratification in the stratified viscous regime, the vertical eddy scale becomes independent of BP e and now depends solely on Reynolds number, with the empirical relationship given by: analogously to the viscously-affected stratified regime considered by Brethouwer et al. (2007) and discussed in the introduction (since l z in (5.2) is non-dimensional and scaled by a characteristic horizontal length scale). While only three clear regimes are evident in this plot, data from the stratified intermittent regime discussed in section 4.3 lies in the region of parameter space between the l z ∼ (BP e) −1/3 and l z ∼ Re −1/2 regimes, as the DNSs begin to feel the effects of viscosity, and hence Re.
It is also of interest to observe how the mixing efficiency η, discussed in section 4.4, depends on the stratification BP e and Reynolds number Re. Figure 8(b) shows η as a function of BP e for each of our simulations. This time, the four regimes can be clearly identified. For the unstratified regime, the mixing efficiency depends only on BP e, and is given by: As the stratification increases, the mixing efficiency increases until it reaches a plateau at η 0.4 which, as we argue below, is a defining property of the stratified turbulent regime. The range of values of BP e for which η is approximately constant is very small for Re = 50, but clearly increases with Re, and is quite substantial for Re = 600. However, in all cases, a threshold is reached where η begins to decrease again. To understand why this is the case, note that the vertical eddy scale decreases rapidly (as discussed above) as the stratification increases and inevitably reaches a point where the effects of viscosity begin to play a role. This is manifest in the fact that η begins to depends on Re. The system enters the intermittently turbulent regime, where we observe the new empirical scaling: Finally, for even larger values of BP e, our DNSs suggest a fourth regime for very large stratification, where we observe the scaling η 0.25Re 2 (BP e) −1 (5.5) which we described earlier (see section 4.3) as characteristic of the stratified viscous regime. These empirical observations inspire us to attempt to derive scaling laws using ideas of dominant balance in the governing equations.

Derivation of scaling regimes
In the following analysis, and consistent with our study of low Péclet number systems, we always assume a LPN balance in equation (2.7) such that Our approach, therefore, is to consider the dominant balance between terms in the momentum equation (2.6), specifically the relative importance of stratification, inertia and viscosity.

Unstratified regime
We begin by considering the unstratified regime, described in section 4.3 and illustrated in figures 6(a)-(c). Motivated by the qualitative observation of the domain-filling eddies in figures 6(a) and 6(b), we make the assumptions that each of the three velocity components and eddy length scales are approximately isotropic with u rms , v rms , w rms ∼ O(1); l x , l y , l z ∼ O(1). (5.7) These assumptions for w rms and l z are confirmed in figures 8(a) and 8(c), indicated by the red lines. By combining the LPN approximation (5.6) with assumptions (5.7), we find a scaling for the typical temperature perturbations: In terms of the mixing efficiency η, (5.7) implies u sin(y) ∼ O(1). Thus η ∼ B wT u sin(y) ∼ Bw rms T rms ∼ BP e. (5.9) The theoretically derived scalings (5.8) and (5.9) are consistent with the empirical scalings determined using our DNSs, shown using the red lines in figures 8(d) and 8(b) respectively. The lack of Re-dependence affirms the irrelevance of viscosity in this regime. Finally, it is of interest to compute the condition of validity for this unstratified regime. In the vertical momentum equation (2.6), we have assumed that stratification is weak relative to fluid inertia, such that BT u · ∇w. Using the scalings derived above, we find that this is true when BP e O(1).
(5.10) Condition (5.10) combined with the condition for linear instability (Re > 2 1/4 ) defines the The dynamics of stratified horizontal shear flows at low Péclet number 25 region of parameter space in which we would expect to observe this regime of unstratified turbulence.

Stratified turbulent regime
As the stratification increases, the system transitions into the stratified turbulent regime, first presented in section 4.3. This regime is defined by a constant mixing efficiency. We anticipate that the velocity components will no longer be isotropic, but have no reason to assume the same about the length scales of the eddies, which are now small-scale and are typically generated through localized shear-driven Kelvin-Helmholtz type instabilities. Consequently, we assume that u rms , v rms ∼ O(1); l x ∼ l y ∼ l z . (5.11) With this assumption, the LPN approximation becomes w rms ∼ P e −1 T rms l 2 z . (5.12) In the vertical momentum equation, we now expect a dominant balance between inertia and stratification such that u · ∇w ∼ BT , implying that u rms w rms l −1 z ∼ BT rms .
( 5.13) This, combined with (5.12) and u rms ∼ O(1), leads to the vertical eddy length scale l z ∼ (BP e) −1/3 , (5.14) which is confirmed by the yellow line in figure 8(a). Empirically, we find that the prefactor is close to 2, and confirm that this scaling is independent of the Reynolds number. As mentioned earlier, η ∼ O(1) is a defining property of the stratified turbulent regime, with a roughly equal partitioning between viscous dissipation and thermal dissipation. The yellow line in figure 8 There is strong evidence for both of these scalings as illustrated by the yellow lines in figures 8(d) and 8(c) respectively, and the empirical data are consistent with the associated prefactors being close to one in each case. Once again, we highlight the lack of dependence on Re in (5.17).
In this regime, we can finally estimate a generic non-dimensional turbulent diffusivity for vertical transport of a passive scalar as D turb ∼ w rms l z ∼ (BP e) −1/2 , (5.18) with a prefactor that is expected to be of order unity. This result can be compared with the mixing coefficient expected in low Péclet number stratified turbulence caused by vertical shear, which scales as (RiP e) −1 instead (see (1.13), when cast in non-dimensional form). We see that D turb decreases much less rapidly with increasing stratification in horizontally-sheared flows than in vertically sheared flows, at least while the system is in this stratified turbulent regime. The assumptions that we made in the vertical momentum equation balance, i.e. that the viscous terms are negligible (Re −1 ∇ 2 w BT ), along with scalings for l z , u rms and T rms , lead to the condition BP e Re 2 . This suggests that the stratified turbulent regime scalings should apply when: 1 BP e Re 2 .
(5.19) Condition (5.19), computed more precisely in section 5.2.4, uniquely defines the region of parameter space in which we would expect to observe this particular type of stratified turbulence in flows at low Péclet number.

Stratified viscous regime
As discussed in section 4.3, for very strong stratification we observe the formation of thin and viscously coupled layers, each containing almost two-dimensional flow. Consequently, we expect that horizontal and vertical velocity components and length scales will both be strongly anisotropic. Denoting horizontal length scales as l h , we make the following assumptions: In what follows, we split the velocity field into a horizontal and vertical component, u = u h + we z , with a corresponding decomposition of the gradient operator ∇ = (∇ h , ∂/∂z). The momentum equation can be split into its horizontal and vertical components as: If we assume a dominant balance between viscosity and the forcing in the horizontal momentum equation (5.22), then Re −1 ∂ 2 z u h ∼ sin(y)e x ∼ O(1). This balance, combined with u h ∼ O(1), leads to the classical viscous scaling for the vertical length scales (cf. Brethouwer et al. 2007): Substantial evidence for this scaling is visible in figure 8(a), where the series of blue lines correspond to l z 2Re −1/2 for each individual Reynolds number. Note that these strongly stratified DNSs exhibit large amplitude quasi-time-periodic behaviour, a feature that we believe to be an intrinsic property of such flows. We consequently attribute the large error bars associated with some DNSs to this observation. In the vertical momentum equation, we assume that the dynamics are hydrostatic, therefore ∂ z p ∼ BT implies pl −1 z ∼ BT rms . This approximation, combined with the requirement from the balance in the horizontal momentum equation that p ∼ O(1), and with the scaling (5.24) for l z , gives us a scaling for temperature perturbations: T rms P e ∼ Re 1/2 (BP e) −1 . (5.25) This stratified viscous regime is considerably more challenging to simulate than the other three regimes, a consequence of the very small time steps required and long integration times. However, we see in figure 8(d) that the blue lines, which represent the scalings in (5.25), fit the few available data points well, once again with a prefactor close to one. The LPN approximation (5.12), combined with (5.24) and (5.25), leads to a scaling for the vertical velocity field: w rms ∼ Re 3/2 (BP e) −1 . (5.26) Again we see a good correspondence between the blue curves in figure 8(c), which represent this scaling, and the data, with a prefactor of 0.25. Using these results, we finally find that with a prefactor of 0.25 for consistency with the data obtained for w rms and T rms . This is consistent with observations for Re = 50 and Re = 100 in figure 8(b). We can also estimate a generic non-dimensional turbulent diffusivity for vertical transport of a passive scalar as with a prefactor that is again expected to be of order unity. The viscous regime is achieved in the opposite limit to the one derived in (5.19) for the stratified turbulent regime, namely when BP e Re 2 . Thus we find that the system parameters must satisfy 2 1/2 < Re 2 BP e (5.29) when combined with the condition for linear instability. Condition (5.29), computed more precisely in section 5.3, defines the region of parameter space in which we would expect to observe this stratified viscous regime. We note for consistency that each of the scalings obtained here do depend on the value of the Reynolds number, as one would expect.

Stratified intermittent regime
There exists a fourth regime, visible both in the DNSs and the results presented in figures 8(a)-(d). This final regime is a transitional regime that occurs between the stratified turbulent regime and the stratified viscous regime. As discussed in section 4.3, it is inherently intermittent in the sense that we observe spatially and temporally localised patches of small-scale turbulence generated via vertical shear instabilities, surrounded by more laminar, viscously dominated flow. Whilst we have been unable to derive satisfactory scalings for this regime, we can nevertheless deduce some of them empirically from figures 8(b) and 8(c).
For instance, we see in figure 8(b) that the onset of this stratified intermittent regime (indicated by the green lines) is characterised by a sudden change in the dependence of the mixing efficiency η on BP e, from the constant value of 0.4 observed in the stratified turbulent regime to a regime where η is given by (5.4). It is interesting and perhaps reassuring to note that the parameter group BP e/Re 2 , which controls η in this regime, is the same parameter group that appears in the viscous regime. Note that for η 0.1, we observe a temporary flattening of this scaling just before the onset of the viscous regime. It is certainly possible that this feature is an artefact of inherent variability in the simulations (and therefore the measurement of η has larger associated error bars). It is interesting to note, however, that this "knee" in the curve does occur for flows with at least three different Reynolds numbers and at the same value of the key parameter BP e in each case. Figure 9. Regime diagram, applicable in the LPN limit, illustrating five dynamical regimes across system parameters BP e (horizontal axis) and Re (vertical axis). Each regime is associated with a colour: linearly stable (purple); unstratified (red); stratified turbulent (yellow); stratified intermittent (green); stratified viscous regime (blue). The four example DNSs presented in figure  6 are associated with parameters corresponding to the red, yellow, green and blue squares.
In addition, figure 8(c) suggests that w rms scales as w rms 0.05Re 3/4 (BP e) −1/2 . (5.30) No clear scalings for l z or T rms /P e appear to be deducible from the numerical results. Figure 9 summarises the four regimes of nonlinear saturation that were described in section 5.2 along with the inclusion of the linearly stable regime that was discussed in section 3. The unstratified regime, indicated in red in figure 9, occurs when BP e O(1); Re > 2 1/4 . (5.31)

Regime diagram
The stratified turbulent regime, indicated by the yellow region in figure 9, and the stratified viscous regime, indicated by the blue region, exist when BP e Re 2 and BP e Re 2 respectively, with the stratified intermittent regime (green) lying at the transition. Greater precision on these regime boundaries, permitting the identification of the domain of validity of this stratified intermittent regime, can be determined from figure 8(a).
If we assume that, for each Reynolds number, the transition between the stratified turbulent and stratified intermittent regimes occurs when η 0.4, then the boundary is given by BP e 0.0016Re 2 . This provides the more precise condition for the stratified turbulent regime: 1 BP e 0.0016Re 2 , (5.32) labelled in figure 9. We note that this regime does not intersect the region of linear stability, indicating that for certain Reynolds numbers for which instability occurs (2 1/4 < Re < 25) this particular type of stratified turbulence does not exist. From figure 8(b) we can also estimate that the transition between the stratified intermittent regime and the stratified viscous regime approximately occurs when η 0.05 irrespective of the Reynolds number, which would imply that the boundary is given by BP e 4.6Re 2 . Thus a more precise condition for the stratified viscous regime is given by  Figure 9 shows that the stratified intermittent regime can exist for any value of Re, provided that the system is linearly unstable.

Discussion
As summarized in section 5.3, our numerical experiments have revealed that stratified horizontal Kolmogorov flows at high Reynolds number but low Péclet number exhibit (at least) four different non-trivial dynamical regimes depending on the respective values of the parameters BP e and Re (where B, P e and Re were defined in (2.9) and (2.11)). In all but one of these regimes, well-defined dominant balances in the momentum equation lead to simple scaling laws for the turbulent properties of the flow. We now first compare our results with prior studies of stratified mixing in the geophysical context, and then discuss the implications of our findings for stratified mixing in stars, whose understanding motivated this study.

Comparison with stratified mixing in geophysical flows
As we have demonstrated in this work, geophysical and astrophysical stratified turbulence is fundamentally different, because the former has a Prandtl number P r O(1) while the latter has P r 1. Therefore, crucially, in geophysically-relevant flows, a high Reynolds number flow necessarily also has a high Péclet number. Meanwhile, in astrophysics it is possible to have both Re 1 and P e 1, and the effect of thermal diffusion can become a dominant factor in the system dynamics. As demonstrated by Lignières (1999) (see also Spiegel 1962;Thual 1992), temperature and velocity fluctuations in the low Péclet number limit are slaved to one another, and density layering is prohibited (see section 1.2). This is in stark contrast with geophysical flows where density layering (or at the very least, the propensity to form alternating regions of shallower and steeper density gradients) is key to understanding the properties of stratified turbulence in the LAST regime. Indeed, the standard Miles-Howard stability criterion (Miles 1961;Howard 1961) for linear instability to vertical shear, namely Ri g < 1 4 (where Ri g is here the minimum gradient Richardson number based on the local vertical stratification and vertical shear), is at first glance incompatible with the ubiquitous presence of turbulence in most largescale stratified shear flows in geophysics (in particular in the ocean and atmosphere) where the gradient Richardson number is typically much larger than one, or indeed is irrelevant in the case of horizontally-sheared flows. However, small-scale layering releases this constraint by creating regions where the stratification is locally reduced, and the instability that is now allowed to develop continues to mix the layer, thereby allowing turbulence to sustain itself. This process, as reviewed in section 1.1, can lead to the formation of layers on the scale U c /N c , and is controlled by the buoyancy Reynolds number Re b = ReF r 2 = Re/B.
In astrophysics, typical values of the gradient Richardson number are also very large, but density layering is prohibited so this pathway to turbulence is not available. Instead, we have shown that three-dimensional perturbations of the horizontal shear (see also Deloncle et al. 2007;Arobone & Sarkar 2012;Lucas et al. 2017) cause the flow to develop layers this time in the velocity field that enhance the vertical shear (or create it when it is not initially present). For sufficiently thin velocity layers, thermal diffusion reduces the effect of stratification, allowing vertical shear instabilities to develop in between the layers. These two effects combine to drive turbulence and can cause substantial vertical mixing even when the background flow has no vertical shear. The dynamics of the system are no longer controlled by ReF r 2 , but instead, first by BP e = P e/F r 2 in the limit where BP e Re 2 , and then by the ratio BP e/Re 2 in the limit where BP e Re 2 , thus partitioning parameter space in the four different dynamical regimes discussed in section 5.
The viscous regime that we have identified (when BP e Re 2 ) is analogous to the viscously affected Re b O(1) regime discussed by Brethouwer et al. (2007), in the sense that it relies on the same dominant balances in the momentum equation. As a result, it exhibits the same scaling in terms of the vertical length scale l z ∼ Re −1/2 . It differs, however, in the treatment of the buoyancy equation, which is not surprising given the low Péclet number limit appropriate in our case. On the other hand, the stratified turbulent regime identified here bears little resemblance with the Re b 1, high P e regime of Brethouwer et al. (2007) (i.e. the LAST regime), where l z ∼ U c /L c . Indeed, for this new low Péclet number stratified turbulent regime, as found in section 5. From a dimensional analysis point of view, this new scaling can be understood as the only combination of U c and N c that can be created to form a length scale given the constraint that N 2 c and κ can only appear together as N 2 c /κ in the low Péclet number limit (as is apparent from (1.11)). But more importantly, we also saw that this scaling emerges from the assumption that the turbulent eddies are isotropic on the small scales, with l x ∼ l y ∼ l z (see section 5.2.2), which is quite different from the inherently anisotropic scalings discussed in Brethouwer et al. (2007) where l x , l y l z . In other words, the stratified turbulent regime identified here is (we believe) a genuinely new regime of turbulence, that can only exist at low Péclet number, and so we refer to it as low Péclet number stratified turbulence, LPNST.

Implications for mixing in stars
The first conclusion that can be drawn from our study is that none of the four regimes discovered support the theory proposed by Zahn (1992) for turbulence driven by horizontal shear in stellar radiation zones. Recall (see section 1.2) that the characteristic flow length scale and amplitude in his model are given by (1.14). Written in terms of the non-dimensionalization used in this work (see section 2), these are where is the non-dimensional dissipation rate (see also Lignières 2018). As he assumes that all the energy input in the system (i.e. u sin(y) , which is always of order one in the chosen units) is dissipated viscously (i.e. there is negligible irreversible conversion into the potential energy reservoir), then 1. The corresponding non-dimensional turbulent diffusivity would therefore scale as D turb ∼ (BP eRe) −1/2 . (6.3) Ignoring the unstratified turbulent regime, which is not expected ever to be relevant to stellar interiors, we see that none of the other regimes reproduce these scalings either, regardless of whether we consider the length scale, the velocity, or the turbulent mixing coefficient. A likely explanation for this discrepancy is that Zahn's theoretical scalings are based on an energetic argument, while the flow dynamics are controlled by the momentum equation (in the low Péclet number limit) itself. Since the energy equation is obtained from a projection of the momentum equation, it fails to capture information about the dominant balance of forces that must individually be satisfied in each spatial direction. In other words, (perhaps unsurprisingly) energetic arguments alone are not always sufficient to infer the correct scalings.
While we believe our results are a step forward in the study of stratified mixing in stars, they are nevertheless not yet applicable as is for a number of reasons. First and foremost is the fact that stars are actually in the high Péclet number yet low Prandtl number regime, while the simulations presented here only probe the low Péclet number regime. Indeed, a classic example of a stellar shear layer is the solar tachocline. Located just below the base of the solar convective envelope (Christensen-Dalsgaard & Schou 1988;Goode et al. 1991), this layer contains a horizontal shear flow with characteristic values of the amplitude and length scale of the base flow being U c 150 m/s and L c 5 × 10 8 m, while the buoyancy frequency is of the order of N c 10 −3 s −1 (Hughes et al. 2007). With ν 0.001m 2 /s and κ 1000m 2 /s, this implies Re ∼ O(10 14 ), P e ∼ O(10 8 ), and B ∼ O(10 7 ), with P r ∼ O(10 −6 ). Corresponding numbers for other stars are in the same parameter regime. Our low Péclet number findings are not to be casually dismissed, however: preliminary numerical results in the high Péclet and low Prandtl number regime are found to satisfy very similar scaling laws to the ones described in this paper. This is likely because the effective local Péclet number of the flow (written in terms of the actual vertical eddy scale l z instead of L c ) is low even though the Péclet number based on the global scale itself is large. A thorough study of the high Péclet and low Prandtl number regime is beyond the scope of this paper, however, and will be the subject of future work.
More crucially, however, is the fact that other effects will need to be taken into account before a comprehensive model of stratified mixing in stars can be created. The main source of shear in stars is their differential rotation, where the mean rotation rate is typically substantially larger than the shearing rate, and where the horizontal shear is usually global (i.e. with a length scale of the order of the stellar radius). This implies that the effects of curvature and angular momentum conservation must be taken into account to determine whether the horizontal shear is unstable in the first place. Two-dimensional horizontal shear flows in a rotating spherical shell were first studied by Watson (1980) (see also Garaud 2001), who found that the shearing rate must exceed a critical threshold for instability to proceed. In the context of our work, this implies that rotation could in principle inhibit the development of the primary instability. If the latter does take place, however, we anticipate that the same sequence of instabilities resulting in the development of small-scale eddies of size l z would ensue. The Rossby number based on l z is likely very large (in the tachocline, for instance, Ro ∼ U c /Ωl z ∼ O(10 4 ), where Ω ∼ 3 × 10 −6 s −1 is the mean rotation rate of the Sun), suggesting that rotation would not have a significant effect on the flow dynamics in any stratified turbulent regime. It may be relevant in the intermittent and viscous regimes on the other hand, where the horizontal eddy scale is of the order of the scale of the background flow.
In addition, stars are subject to vertical shear as well as horizontal shear, and the dynamics of shear-induced turbulence are notably different in the two cases (see section 1.2). A question of interest will therefore be to establish what controls the outcome when vertical and horizontal shear are both present. Finally, most stars are expected to be magnetized to some extent (Mestel 2012), either by the presence of a primordial magnetic field or by the action of a dynamo in a nearby convective zone. The effect of these magnetic fields will need to be taken into account to construct a truly astrophysically relevant theory of stratified turbulence.
This work was initiated as a project at the Woods Hole Geophysical Fluid Dynamics summer program in 2018. The authors thank the program for giving them the opportunity to collaborate on this topic and for their financial support. L.C. was also supported by the Natural Environment Research Council (grant number NE/L002507/1). P.G. was also funded by NSF AST 1517927. The research activity of C.P.C. was supported by the EPSRC Programme Grant EP/K034529/1 entitled Mathematical Underpinnings of Stratified Turbulence. All simulations presented here were performed on either the Hyades computer, purchased at UCSC with an NSF MRI grant, or the NSF XSEDE supercomputing facilities (Comet). The authors thank S. Stellmach for providing his code.
This system, which is now identical to the set of equations (3.7)-(3.11) except for the rescaling implicit in the hats on parameters and growth rates, can be re-formulated as a generalised eigenvalue problem for the complex growth ratesσ, A(k x , k z ,Re,B,P r)X =σBX, (A 21) and can solved using the method described in section 3. The linear stability analysis presented in section 3 considered a = 1, whereRe = Re, B = B,P r = P r andσ = σ. For a = 1, relations (A 15) provide a transformation between the original analysis and the linear stability of flows with generic amplitude au L (y).