# Turbulent Transport in a Strongly Stratified Forced Shear Layer with Thermal Diffusion

Pascale Garaud, Logithan Kulenthirarajah

► To cite this version:

Pascale Garaud, Logithan Kulenthirarajah. Turbulent Transport in a Strongly Stratified Forced Shear Layer with Thermal Diffusion. The Astrophysical Journal, 2016, 821, 10.3847/0004-637X/821/1/49. insu-03675422

HAL Id: insu-03675422

https://insu.hal.science/insu-03675422v1

Submitted on 21 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# TURBULENT TRANSPORT IN A STRONGLY STRATIFIED FORCED SHEAR LAYER WITH THERMAL DIFFUSION

PASCALE GARAUD<fnref n="1" /> AND LOGITHAN KULENTHIRARAJAH<fnref n="2" />

## ABSTRACT

This work presents numerical results on the transport of heat and chemical species by shear-induced turbulence in strongly stratified, thermally diffusive environments. The shear instabilities driven in this regime are sometimes called “secular” shear instabilities, and can take place when the Richardson number of the flow is large, provided the Péclet number is small. We have identified a set of simple criteria to determine whether these instabilities can take place or not. Generally speaking, we find that they may be relevant whenever the thermal diffusivity of the fluid is very large (typically larger than $10^{14} \text{ cm}^2 \text{ s}^{-1}$), which is the case in the outer layers of high-mass stars ($M \ge 10 M_\odot$), for instance. Using a simple model setup in which the shear is forced by a spatially sinusoidal, constant-amplitude body-force, we have identified several regimes ranging from effectively unstratified to very strongly stratified, each with its own set of dynamical properties. Unless the system is in one of the two extreme regimes (effectively unstratified or completely stable), however, we find that (1) only about 10% of the input power is used toward heat transport, while the remaining 90% is viscously dissipated; (2) that the effective compositional mixing coefficient is well-approximated by the model of Zahn, with $D \approx 0.02 \kappa_T/J$ where $\kappa_T$ is the thermal diffusivity and $J$ is the Richardson number. These results need to be confirmed, however, with simulations in different model setups and at higher effective Reynolds number.

**Key words:** hydrodynamics – instabilities – stars: massive – turbulence

## 1. INTRODUCTION

The continued progress in stellar spectroscopy, combined with the advent of asteroseismology, have opened new prospects for stellar astrophysics, providing more reliable ways of testing the accuracy of stellar evolution models. These challenge our understanding of the micro- and macro-physical processes that participate in all aspects of a star’s life, from the deep interior to the surface. Today, the greatest source of uncertainty in stellar modeling is arguably turbulent mixing in radiative regions. The need for non-canonical mixing (i.e., turbulent mixing of chemical species or angular momentum that is not related to convection or overshoot) in stellar evolution has long been recognized: despite the general agreement between theory and observations of stars, some discrepancies remain that can only be reconciled by invoking additional mixing in stellar radiation zones (see the review by Pinsonneault 1997, for instance). Examples of observational evidence for missing compositional mixing comes from surface abundances (see, for instance, Skumanich 1972; Fliegner et al. 1996; Gratton et al. 2000; Sestito & Randich 2005; Charbonnel & Zahn 2007) and helioseismology (e.g., Elliott & Gough 1999) while that related to missing angular-momentum mixing comes for instance from asteroseismology of red giant branch (RGB) stars (Eggenberger et al. 2012; Mosser et al. 2012; Deheuvels et al. 2014).

Solving the “missing mixing” problem begins with identifying possible sources of hydrodynamical or magneto-hydrodynamical instabilities, determining under which conditions they operate, and then quantifying the amount of heat, angular momentum, and compositional transport they induce. A commonly discussed mixing mechanism is shear-induced turbulence (e.g., Zahn 1974, 1992; Schatzman 1977; Endal & Sofia 1978; Barnes et al. 1999). Shear—in particular rotational

shear—is, to various levels, always present in a star. Whether shear-induced turbulence develops or not, and how much mixing it causes, depends on the amplitude of the viscous force and of the buoyancy force acting on the fluid, which both have a tendency to quench the instability. In what follows, we now discuss in turn what is known about the conditions for which shear instabilities can develop in stars, which will lead us to discuss the effects of thermal diffusion, and then summarize some of the main shear-induced mixing models that have been proposed.

### 1.1. *Under Which Conditions do Shear Instabilities Occur?*

The nature of the criterion for the onset of shear instabilities has long been the subject of debate, starting from the work of Richardson (1920) who argued on energetic principles that shear-induced turbulence can only be sustained provided $J < J_c$, where $J = N^2/S^2$ is the local gradient Richardson number, $N$ is the local Brunt-Väisälä frequency (often called buoyancy frequency), $S = |du/dr|$ is the local shearing rate of the flow field $u$ and $J_c$ is a universal constant which he argues must be of order one. Later, Miles (1961) and Howard (1961) formally proved that a necessary<fnref n="3" /> condition for linear instability to occur in an inviscid and non-diffusive shear flow is that $J$ must drop below 1/4 somewhere in the fluid. Unfortunately, $J$ rarely ever drops below one in stellar interiors—typical values are far greater, suggesting that the shear should always be stabilized by the stratification (both from linear theory, and on energetic principles).

The effect of thermal diffusion, however, can relax this so-called Richardson criterion by damping the buoyancy force (Townsend 1958; Dudis 1974; Gage & Miller 1974; Zahn

<footnote n="1">¹ Department of Applied Mathematics and Statistics, Baskin School of Engineering, University of California at Santa Cruz, 1156 High Street, Santa Cruz CA 95064, USA</footnote>

<footnote n="2">² Institut de Recherche en Astrophysique et Planétologie (IRAP), 14, avenue Edouard Belin, F-31400 Toulouse, France
*Received 2015 December 29; accepted 2016 February 25; published 2016 April 7*</footnote>

<footnote n="3">³ But not always sufficient.</footnote>
1974; Jones 1977). Indeed, in cases where the background temperature gradient within the star is the main source of stratification (by contrast with the possibility of chemical stratification which we ignore in this work), then the potential energy cost of overturning motions is reduced if the moving fluid has time to adjust thermally to its surroundings. Whether it does or not is typically measured by the eddy Péclet number

$$Pe_l = \frac{Sl^2}{\kappa_T}, \qquad (1)$$

where *l* is the vertical lengthscale of the eddy considered, and $\kappa_T$ is the thermal diffusivity. If $Pe_l \ll 1$, then the motion is strongly diffusive, the eddy is always in thermal equilibrium with its surrounding, and the influence of thermal stratification effectively disappears—at least on the scale *l*.

While a strict yet simple and general linear stability criterion akin to the Miles-Howard theorem does not exist in the case of diffusive shear flows, a combination of mathematical and physical arguments backs the notion that thermal diffusion vastly increases the critical gradient Richardson number $J_c$ above which a flow is stable. Indeed, Townsend (1958) and Dudis (1974) first performed a linear stability analysis of stratified atmospheric shear flows in the optically thin limit, where radiative losses play a similar role as thermal diffusion in damping the buoyancy force. In both cases, they found that thermal radiation begins to influence the stability of the fluid when $St_{cool}$ roughly drops below one (where $t_{cool}$ is the thermal cooling timescale associated with the optically thin radiation law). The inviscid linear instability criterion is no longer $J < J_c = 1/4$, but becomes $J < J_c \sim O(St_{cool})^{-1}$ instead, which implies that strongly stratified shear flows with $J \gg 1$ can still be destabilized provided $St_{cool} \ll 1$.

Zahn (1974) used these results to conjecture on the stability of optically thick shear flows, which are much more relevant for stellar interiors. The equivalent cooling timescale $t_{cool}$ for an eddy of size *l* can be approximated by $t_{cool} = l^2/\kappa_T$, so the condition $St_{cool} < O(1)$ is equivalent to the low Péclet number condition discussed earlier,

$$Pe_l < O(1), \qquad (2)$$

where $Pe_l$ was defined in (1). Proceeding with the analogy with the optically thin case, Zahn argued that the linear inviscid criterion for instability $JSt_{cool} < O(1)$ should therefore become

$$JPe_l < (JPe)_c, \qquad (3)$$

where $(JPe)_c$ is a constant of order unity. This new criterion clearly allows for instability at large $J$ provided $Pe_l$ is small enough.

By contrast with the optically thin limit, however, this criterion depends on the so-far-unspecified lengthscale *l* of the turbulent eddies, and naively suggests that one could always pick *l* small enough for (3) to apply. In other words, shear instabilities should always be possible as long as there is an infinitesimal amount of thermal diffusion. In practice, however, this is of course not the case, since viscosity cannot be ignored on very small scales. Zahn (1974) therefore argued that *l* should be set equal to the smallest scale that is not stabilized by viscosity, $l_c$, which satisfies $Sl_c^2/\nu = Re_c$ where $\nu$ is the kinematic viscosity and the critical Reynolds number $Re_c$ is a universal constant which Zahn suggests is of the order of a thousand. On that scale, $Pe_{l_c} = Pr Re_c$, where $Pr = \nu/\kappa_T$ is the

Prandtl number and the criterion for instability (3) becomes

$$JPr Re_c < (JPe)_c \Leftrightarrow \frac{JPe}{(JPe)_c} < \frac{Re}{Re_c}, \qquad (4)$$

where Pe and Re can now be defined according to any lengthscale *L*, instead of $l_c$ specifically. Note how while (3) only addressed the question of the stability of a shear flow to perturbations of a given size *l*, (4) more generally addresses the question of the stability of shear flows to perturbations at arbitrary lengthscales. This new criterion clearly allows for the development of instabilities for values of *J* much larger than one provided the Prandtl number Pr is sufficiently low, and is sometimes referred to as a criterion for “secular shear instabilities” (Endal & Sofia 1978). We note, however, that Zahn’s criterion is not a formal mathematical result on either linear or nonlinear stability of stratified diffusive shear flows.<fnref n="4" /> Rather, it should be viewed as a plausible heuristic criterion whose validity remains to be verified. Furthermore, given that this criterion is in principle only valid for low Péclet number flows, its domain of applicability also needs to be clarified (see below for more on this topic).

Significant progress has been made in this direction in the past 20 years. Lignières (1999) and Lignières et al. (1999) showed that the dynamics of low Péclet number flows in general can be studied with a reduced set of equations in which the temperature field is slaved to the velocity field (see Section 2, and see also Thual 1992). In these equations, the Richardson number and the Péclet number always appear together as a product called the “Richardson-Péclet” number hereafter. This formally shows that any criterion for stability or instability to shear at low Péclet number should involve the Richardson-Péclet number rather than the Richardson number alone, which is indeed the case in (3) and (4) for instance. Using these so-called “low Péclet number” (LPN) equations, Lignières et al. (1999) then showed that the linear stability properties of a viscous stratified hyperbolic tangent shear layer do indeed satisfy a criterion of the kind written in (4) (see also Jones 1977). First results in the nonlinear regime were obtained by Prat & Lignières (2013), who showed numerically that a viscous stratified linear shear layer is nonlinearly unstable to finite amplitude perturbations that have a mean lengthscale *l* provided $JPe_l < (JPe)_c \approx 0.426$, hence verifying (3). Their experimental protocol (see below for more detail) however was unable to address the more general and more crucial question of the validity of (4).

More recently, Garaud et al. (2015a) showed that any periodic stratified shear flow in the limit of asymptotically low Péclet number is formally stable to any initial perturbations (of arbitrary amplitude and shape) provided the Richardson-Péclet number is greater than a factor of order unity times the Reynolds number, which is equivalent to Zahn’s criterion (4) except with $Re_c \sim O(1)$ instead of $O(10^3)$. This implies that flows with $J$ as large as $Pr^{-1}$ could in principle be subject to finite amplitude instabilities, potentially expanding the region of parameter space prone to shear instabilities even further compared with Zahn’s predictions, albeit in the limit of small Péclet numbers. In practice, Garaud et al. (2015a) found that

<footnote n="4">⁴ If anything, Zahn’s argument is necessarily one that pertains only to the nonlinear stability of the flow, since for small enough scales, the background shear will always be very nearly linear, and linear shear flows are inherently linearly stable.</footnote>
only a subset of that region is nonlinearly unstable, but their results remain to be confirmed with more extensive simulations as well as with non-periodic flows. In summary, all of these studies point in the same direction, namely that Zahn's criterion given in (4) appears to hold for low Péclet number flows both in the linear sense and in the nonlinear sense, although the actual values of $(JPe)_c$ and $Re_c$ remain to be determined in each case, and depend on the global shape of the shear flow.

## 1.2. What is a Low Péclet number Flow?

The domain of validity of the LPN equations, and by proxy, the domain of applicability of criterion (4), is clearly limited to low Péclet number flows. However, since one can define a number of different Péclet numbers depending on the length- and velocity-scales of interest, it is worth asking which Péclet number has to be small for these equations to be valid. The original derivation of Lignières (1999) argues that one should use the rms velocity of the flow, $u_{rms}$, together with the typical vertical scale of energy-bearing eddies, $l_v$, to do so. In other words, the LPN equations and criterion (4) should apply as long as the turbulent Péclet number

$$Pe_t \equiv \frac{u_{rms} l_v}{\kappa_T} \ll 1. \qquad (5)$$

This seems to be confirmed by the numerical simulations of Prat & Lignières (2013; see also Prat & Lignières 2014), and of Garaud et al. (2015a). As mentioned above, Prat & Lignières (2013) considered strictly linear stratified shear flows, which are known to be linearly stable (Knobloch 1984). Each one of their simulations has a well-defined background shearing rate $S$, and is initiated with finite amplitude perturbations that have a given velocity spectrum defining a dominant eddy size $l_v$, and a typically velocity scale $u_{rms}$ (which together define their initial turbulent Péclet number $Pe_t = u_{rms}l_v/\kappa_T$). By comparing the temporal evolution of the total kinetic energy in the various simulations, they were able to identify nonlinearly stable cases (i.e., cases where the kinetic energy gradually decays to zero) from nonlinearly unstable ones (i.e., where the total kinetic energy increases), and thus measured numerically the critical threshold for nonlinear instability. They found that the critical Richardson-Péclet number for the onset of instability as measured in simulations that use the LPN equations was an excellent match to those obtained using the full Navier-Stokes equations in the limit where $Pe_t$ was smaller than one, and that it is given by $(JPe)_c \approx 0.426$ in that limit. In a related study, Garaud et al. (2015a) studied the development of shear instabilities in sinusoidally forced stratified shear flows. They defined their Péclet number based on the macroscopic lengthscale $L$ and velocity $U_L$ of the laminar solution of the forced problem, $Pe = U_L L/\kappa_T$. As for Prat & Lignières (2013), they found that the critical Richardson-Péclet number for the onset of instability as measured in simulations that use the LPN equations was an excellent match to the one obtained with the full equations provided $Pe < 1$. However, they also found that the domain of validity of the LPN equations extended to larger $Pe$ as well, a result they attributed to the fact that the turbulent Péclet number $Pe_t$ of the simulation could be smaller than one even if $Pe$ was larger than one.

The results of Garaud et al. (2015a) illustrate that, while correct, Lignieres' criterion for the applicability of the LPN equations and of criterion (4) (namely $Pe_t < 1$) is not always practical for stellar evolution calculations since it requires

a priori knowledge of the typical eddy scale $l_v$ and velocity $u_{rms}$ of the turbulent flow, but these can only be obtained a posteriori from full hydrodynamic simulations. Instead, it would be more desirable to have a criterion that predicts under which conditions the LPN equations will apply, based only on known global parameters. With this knowledge, we would then be able to determine fairly easily whether there are stars for which the low Péclet number limit is relevant, and then apply Zahn's criterion more specifically to determine which stellar regions may be subject to secular shear instabilities at large Richardson number.

## 1.3. How Much Mixing do Shear Instabilities Cause?

Knowing when secular shear instabilities are likely to occur is only half of the problem—one also needs to quantify the transport of momentum, heat and composition induced, in order to make progress toward solving the missing mixing problem. Zahn (1974) did not seem to address the question, and the first reference to the mixing rate of shear instabilities in the presence of strong diffusion can be found in the work of Endal & Sofia (1978) instead. They argued, as is commonly done, that the turbulent mixing coefficient $D$ can be estimated by multiplying a typical lengthscale $l$ to a typical velocity $u \sim Sl$, and proposed that $l$ be the smallest possible scale $l_c$ for which $Sl_c^2/\nu = Re_c$ (i.e., the same scale discussed earlier in the context of the secular shear instability criterion). This then yields

$$D \approx Sl_c^2 \approx \nu Re_c, \qquad (6)$$

which implies that $D$ should always be larger than $\nu$ by a factor of $Re_c$ regardless of the level of stratification and regardless of the thermal diffusivity, as long as the criterion for secular shear instabilities is met.

Later, Zahn (1992) similarly argued that $D$ should be the product of a characteristic lengthscale $l$ and a characteristic velocity $u$, but that the relevant lengthscale $l$ should instead be the largest possible scale for which $JPe_l$ is still in the nonlinearly unstable regime, in other words, one for which

$$JPe_l = J \frac{ul}{\kappa_T} = (JPe)_c, \qquad (7)$$

where $(JPe)_c$ is the same universal constant of order one discussed in the context of Equation (3) (see Prat & Lignières 2013, for a discussion of this formula). This then implies that the diffusion coefficient $D$ must take the form

$$D = \beta ul = \beta (JPe)_c \frac{\kappa_T}{J} = \nu \beta \frac{(JPe)_c}{JPr}, \qquad (8)$$

where $\beta$ is another proportionality constant of order one. This shows that $D$ should simply be proportional to $\kappa_T/J$. If the flow is marginally stable to secular shear instabilities—or in other words, if the flow just satisfies the stability criterion (4) so that $JPr \approx (JPe)_c/Re_c$, then the estimates from Endal & Sofia (1978) and from Zahn (1992) are roughly equivalent. However, allowing for the possibility of flows that are significantly beyond marginal stability, that is, flows for which $JPr \ll Re_c^{-1}$, then the mixing rate proposed by Zahn (1992) can be orders of magnitudes larger than the one proposed by Endal & Sofia (1978), and does depend both on the local shearing rate and stratification (through $J$) as well as on the local thermal
diffusivity $\kappa_T$. Note that Zahn's prediction for $D$ must clearly fail in the strict limit where $J \to 0$. However, as discussed earlier, stellar shear layers are rarely, perhaps never, that weakly stratified, so the proposed coefficient is probably reasonable for all practical purposes.

A first attempt to test numerically which of these two models for $D$ is more accurate was recently reported by Prat & Lignières (2013, 2014). Their experimental protocol, which was discussed above, led them to focus only on flows that were marginally stable to finite amplitude perturbations with a prescribed characteristic lengthscale $l$. For these “fixed lengthscale” marginal solutions, they found that $D \simeq 0.058 \kappa_T/J$ in the limit where $Pe_l \le 1$, as predicted by Zahn (1992). At a first glance, this appears to be a remarkable validation and calibration of Zahn's theory which also rules out the Endal & Sofia (1978) prescription. It is crucial to note, however, that Prat & Lignières (2013) effectively test Zahn's model precisely under the conditions for which it was derived (i.e., $J$ is selected to be at the marginal nonlinear stability threshold for the selected lengthscale $l$), so it is not entirely surprising to see that the model works well in that case. Whether these particular conditions would naturally arise in stars where the shear and the scale of the turbulent eddies are free to evolve remains to be determined, and therefore so is the validity of Zahn's model.

In summarizing past work, we have therefore raised three questions. (1) What is the correct criterion to apply when trying to establish whether a particular stellar region is undergoing shear-induced mixing or not, taking into account the effect of thermal diffusion (but ignoring compositional stratification, at least for now) (2) Can one determine from first principles and without resorting to numerical calculations whether thermal diffusion will have a strong effect on the dynamics of stellar shear layers or not, and (3) What controls the turbulent mixing rates in stars where the shear layer and the resulting shear-induced turbulence are both free to evolve naturally? While the answer to question (1) has already been discussed fairly comprehensively in Section 1.1 (and references therein), this paper focusses on questions (2) and (3). In Section 2 we present our model setup, and briefly discuss the LPN equations. Section 3 presents typical outcomes of the numerical experiments, both in the low Péclet number limit and in the high Péclet number limit, for weakly and strongly stratified shear flows. Section 4 addresses and answers question (2), namely under which simple conditions the LPN equations are valid approximations of the full Navier-Stokes equations. Section 5 then discusses which stars are most likely to harbor regions where the LPN equations are valid, or equivalently, where secular shear instabilities might be relevant. Section 6 looks more quantitatively at the available numerical data, and draws preliminary conclusions on the transport properties of secular shear instabilities thus partially addressing question (3). Section 7 finally discusses the implication of the results for astrophysical modeling, and suggests future avenues of investigation.

## 2. THE MODEL

We consider a region of a stellar radiation zone located around the radius $r_0$, with a local background density profile $\bar{\rho}(r) \simeq \rho_0 + (r - r_0)d\bar{\rho}/dr + \dots$, and a local background temperature profile $\bar{T}(r) = T_0 + (r - r_0)d\bar{T}/dr + \dots$. The region also has a background adiabatic temperature gradient $dT_{ad}/dr \simeq (T_0/P_0)\nabla_{ad}$ where $P_0$ is the background pressure at $r = r_0$.

We then consider a small domain around $r_0$, and model it using a Cartesian coordinate system with gravity defining the $z$ direction (so $z = r - r_0$, for instance). We assume that the domain is much thinner than a pressure scaleheight, and that any fluid flow is slow enough to use the Boussinesq approximation (Spiegel & Veronis 1960). As in the work of Garaud et al. (2015a), we assume that the shear is created by a body-force applied in the $x$-direction, whose amplitude varies in the $z$-direction but is otherwise constant with time:

$$F = F(z)e_x. \qquad (9)$$

This “constant forcing” setup could model, for instance, the role of Euler’s force (e.g., the fictitious acceleration that arises when a rotating frame of reference is spinning up or down) in the generation of rotational shear by the very slow differential expansion and contraction of stellar layers. It could also model the effects of quasi-static tides in binary systems. For simplicity, we assume that $F(z) = F_0 \sin(kz)$ is spatially sinusoidal. While this may not be particularly realistic, the periodic nature of $F$ allows us to use a fully spectral code, which assumes that all the perturbations (to the temperature, velocity field, pressure) are triply periodic. Section 7 discusses the potential pros and cons of using a sinusoidal force versus other possible ways of forcing the shear.

This system is described by the following dimensional equations (which neglect the role of compositional stratification):

$$\begin{aligned} \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} = & -\frac{1}{\rho_0} \nabla p + \alpha g T e_z + \nu \nabla^2 \mathbf{u} \\ & + \frac{F_0}{\rho_0} \sin(kz) e_x, \\ \nabla \cdot \mathbf{u} = & 0, \\ \frac{\partial T}{\partial t} + \mathbf{u} \cdot \nabla T + \left( \frac{dT}{dr} - \frac{dT_{ad}}{dr} \right) w = & \kappa_T \nabla^2 T, \end{aligned} \qquad (10)$$

where $p$ and $T$ are the pressure and temperature perturbations away from the hydrostatic, linearly stratified background quantities $\bar{\rho}(r)$ and $\bar{T}(r)$, and $\mathbf{u} = (u, v, w)$ is the velocity field. The diffusivities $\kappa_T$ and $\nu$, as well as the thermal expansion coefficient $\alpha = -\rho_0^{-1}(\partial\rho/\partial T)$ and the local gravity $g$ are all assumed to be constant.

As discussed by Garaud et al. (2015a), these equations have a laminar solution with no temperature fluctuations, where the velocity field is given by:

$$\mathbf{u}_L = \frac{F_0}{k^2 \nu \rho_0} \sin(kz) e_x \equiv U_L \sin(kz) e_x. \qquad (11)$$

The amplitude and lengthscale of this velocity field can be used to define a new unit system where

$$[u] = U_L = \frac{F_0}{k^2 \nu \rho_0} \quad \text{is the unit velocity,}$$

$$[l] = k^{-1} \quad \text{is the unit length,}$$

$$[t] = \frac{[l]}{[u]} = \frac{k \nu \rho_0}{F_0} \quad \text{is the unit time,}$$

$$[T] = k^{-1} \left| \frac{dT}{dr} - \frac{dT_{ad}}{dr} \right| \quad \text{is the unit temperature,}$$

$$[p] = \rho_0 [u]^2 = \frac{F_0^2}{\rho_0 k^4 \nu^2} \quad \text{is the unit pressure.} \qquad (12)$$
This system of units is traditionally used in studies of stratified shear instabilities, as in Lignières et al. (1999) and Garaud et al. (2015a) for instance. The non-dimensionalized equations become:

$$ \frac{\partial \tilde{u}}{\partial t} + \tilde{u} \cdot \nabla \tilde{u} = - \nabla \tilde{p} + \text{Ri} \tilde{T} e_z + \frac{1}{\text{Re}} \nabla^2 (\tilde{u} - \sin(z) e_x), \quad (13) $$

$$ \nabla \cdot \tilde{u} = 0, \quad (14) $$

$$ \frac{\partial \tilde{T}}{\partial t} + \tilde{u} \cdot \nabla \tilde{T} + \tilde{w} = \frac{1}{\text{Pe}} \nabla^2 \tilde{T}, \quad (15) $$

where $\tilde{u}$, $\tilde{p}$ and $\tilde{T}$ are the non-dimensional velocity, pressure and temperature fields, where the differential operators as well as the independent variables x, y, z and t have been implicitly non-dimensionalized as well, and where

$$ \begin{aligned} \text{Re} &= \frac{U_L}{k\nu} = \frac{F_0}{\rho_0 \nu^2 k^3}, \\ \text{Ri} &= \frac{N^2}{U_L^2 k^2} = \frac{\alpha g \left| \frac{d\tilde{T}}{dr} - \frac{dT_{\text{ad}}}{dr} \right| \rho_0^2 \nu^2 k^2}{F_0^2}, \\ \text{Pe} &= \frac{U_L}{k\kappa_T} = \frac{F_0}{\rho_0 \nu k^3 \kappa_T}. \end{aligned} \quad (16) $$

The laminar solution (11) now takes the dimensionless form $\tilde{u}_L = \sin(z) e_x$. It is worth bearing in mind that, being based on the laminar solution of the equations, the Reynolds, Péclet and Richardson numbers just defined are not necessarily relevant to the dynamics of the turbulent solution.

As discussed in Section 1, we will specifically be interested in low Péclet number flows, for which thermal diffusion is significant. According to Lignières (1999), in this limit the temperature fluctuations are slaved to the vertical velocity fluctuations as

$$ \tilde{w} = \frac{1}{\text{Pe}} \nabla^2 \tilde{T}, \quad (17) $$

and the flow dynamics can be modeled using the reduced low Péclet number momentum equation:

$$ \begin{aligned} \frac{\partial \tilde{u}}{\partial t} + \tilde{u} \cdot \nabla \tilde{u} = & - \nabla \tilde{p} + \text{RiPe} \nabla^{-2} \tilde{w} e_z \\ & + \frac{1}{\text{Re}} \nabla^2 (\tilde{u} - \sin(z) e_x). \end{aligned} \quad (18) $$

There are two clear advantages in using this system compared with the standard equations at low Péclet number (Lignières 1999): on the one hand, there is one less variable, and therefore one less equation to solve for, and on the other hand this reduced asymptotic system bypasses the complications that may arise from having to follow the evolution of several fields that evolve on vastly different timescales (i.e., the stiffness problem). Meanwhile, the inverse Laplacian is dealt with trivially in our triply periodic spectral code, noting that $\tilde{w} = \frac{1}{\text{Pe}} \nabla^2 \tilde{T}$ becomes an algebraic equation in spectral space and that we are only seeking periodic solutions for $\tilde{T}$ that have zero horizontal and vertical means (which eliminates the additive linear solutions). As discussed by Garaud et al. (2015a), the LPN equations implicitly assume that the

horizontally averaged temperature profile must always be the same as the background profile, which is not the case for the full equations. We shall compare in this work the predictions of the standard equations and of the LPN equations in more detail, using Direct Numerical Simulations.

Note that there is an alternative non-dimensionalization of the system that will turn out to be very useful (see Section 4.2 and beyond). Indeed, dimensionally speaking one can also create a typical velocity based on the imposed force as:

$$ U_F = \left( \frac{F_0}{k \rho_0} \right)^{1/2}. \quad (19) $$

This new velocity scale $U_F$ no longer depends on any diffusivities, but instead only depends on the characteristics of the forcing (and of the background density). As such, it is a quantity that is more likely to be relevant at high Reynolds numbers than the laminar flow amplitude $U_L$. A more detailed physical interpretation of this velocity will be given in Section 4.2. If we use $U_F$ to create a new system of units, exactly as in (12) but with $[u] = U_F$, we can define new Péclet, Reynolds and Richardson numbers as:

$$ \begin{aligned} \text{Re}_F &= \frac{U_F}{k\nu} = \left( \frac{F_0}{k^3 \rho_0 \nu^2} \right)^{1/2} = \text{Re}^{1/2}, \\ \text{Pe}_F &= \frac{U_F}{k\kappa_T} = \left( \frac{F_0}{k^3 \rho_0 \kappa_T^2} \right)^{1/2} = \text{Re}^{-1/2} \text{Pe}, \\ \text{Ri}_F &= \frac{N^2}{k^2 U_F^2} = \frac{N^2 \rho_0}{k F_0} = \text{ReRi}. \end{aligned} \quad (20) $$

In this new system of units, the velocity, pressure and temperature fields are denoted as $\hat{u}$, $\hat{p}$ and $\hat{T}$, and the new system of governing equations is:

$$ \frac{\partial \hat{u}}{\partial t} + \hat{u} \cdot \nabla \hat{u} = - \nabla \hat{p} + \text{Ri}_F \hat{T} e_z + \frac{1}{\text{Re}_F} \nabla^2 \hat{u} + \sin(z) e_x, \quad (21) $$

$$ \nabla \cdot \hat{u} = 0, \quad (22) $$

$$ \frac{\partial \hat{T}}{\partial t} + \hat{u} \cdot \nabla \hat{T} + \hat{w} = \frac{1}{\text{Pe}_F} \nabla^2 \hat{T}. \quad (23) $$

Note that the unit length and the unit temperature have not changed, so $\hat{T} = \hat{T}$. Meanwhile, $\hat{u} = \text{Re}_F \hat{u} = \text{Re}^{1/2} \hat{u}$.

### 3. NUMERICAL SIMULATIONS

#### 3.1. The Numerical Model

In this section, we present and discuss results using the traditional non-dimensionalization based on the laminar flow velocity unless specifically mentioned. In the majority of the simulations discussed in this paper, the set of Equation (10) is solved in a triply periodic domain of non-dimensional size $L_x = 10\pi$, $L_y = 4\pi$ and $L_z = 2\pi$, using the pseudo-spectral code originally developed by S. Stellmach (Stellmach et al. 2011; Traxler et al. 2011) and modified for the purposes of this work to include the body forcing. Simulations are either initiated with small random fluctuations in the velocity field, or started from the results of a previous run at different parameter values. Table 1 shows a record of all simulations run in this
<table><caption>**Table 1**
Summary of the Main Results for All the Runs with the Standard Equations</caption><thead><tr><th>0.01Pe<br>Pe<sub>F</sub></th><th>10<sup>4</sup>Ri<br>Ri<sub>F</sub></th><th>100&#x0302;<sub>rms</sub><br>&#x0302;<sub>rms</sub></th><th>I&#x0302;<sub>v</sub><br>I&#x0302;<sub>v</sub></th><th>-100(&#x0302;&#x0307;T&#x0307;)<br>-(&#x0302;&#x0307;T&#x0307;)</th><th>-100(&#x0302;&#x0307;C&#x0307;)<br>-(&#x0302;&#x0307;C&#x0307;)</th><th>100(&#x0302;&#x0307; &middot; &#x0307;F&#x0307;)<br>&#x0302;&#x0307; &middot; &#x0307;F&#x0307;</th></tr></thead><tbody><tr><td>0.1</td><td>10<sup>-3</sup></td><td>2.21 &plusmn; 0.06</td><td>0.791 &plusmn; 0.01</td><td>1.48 &plusmn; 0.2</td><td>12.7 &plusmn; 2</td><td>0.369 &plusmn; 0.04</td></tr><tr><td>0.1</td><td>10<sup>-2</sup></td><td>2.2 &plusmn; 0.05</td><td>0.802 &plusmn; 0.01</td><td>1.21 &plusmn; 0.2</td><td>10.1 &plusmn; 1</td><td>0.370 &plusmn; 0.03</td></tr><tr><td>0.1</td><td>10<sup>-1</sup></td><td>2.4 &plusmn; 0.05</td><td>0.78 &plusmn; 0.017</td><td>0.48 &plusmn; 0.1</td><td>3.34 &plusmn; 0.6</td><td>0.624 &plusmn; 0.04</td></tr><tr><td>0.1</td><td>1</td><td>2.46 &plusmn; 0.09</td><td>0.7 &plusmn; 0.016</td><td>0.1 &plusmn; 0.03</td><td>0.94 &plusmn; 0.1</td><td>0.908 &plusmn; 0.07</td></tr><tr><td>0.1</td><td>10</td><td>3.39 &plusmn; 0.04</td><td>0.6 &plusmn; 0.006</td><td>0.01 &plusmn; 0.001</td><td>0.28 &plusmn; 0.01</td><td>1.86 &plusmn; 0.04</td></tr><tr><td>0.1</td><td>10<sup>2</sup></td><td>5.12 &plusmn; 0.02</td><td>0.494 &plusmn; 0.003</td><td>0.004 &plusmn; 0.002</td><td>0.12 &plusmn; 0.003</td><td>3.31 &plusmn; 0.02</td></tr><tr><td>0.1</td><td>10<sup>3</sup></td><td>11.3 &plusmn; 0.04</td><td>0.355 &plusmn; 0.002</td><td>9.6 &times; 10<sup>-4</sup> &times; 4 &times; 10<sup>-5</sup></td><td>0.03 &plusmn; 0.001</td><td>7.86 &plusmn; 0.03</td></tr><tr><td>1</td><td>10<sup>-4</sup></td><td>2.21 &plusmn; 0.06</td><td>0.796 &plusmn; 0.01</td><td>8.5 &plusmn; 2</td><td>11.8 &plusmn; 3</td><td>0.366 &plusmn; 0.06</td></tr><tr><td>1</td><td>10<sup>-3</sup></td><td>2.20 &plusmn; 0.06</td><td>0.797 &plusmn; 0.01</td><td>7.81 &plusmn; 1</td><td>10.7 &plusmn; 2</td><td>0.368 &plusmn; 0.04</td></tr><tr><td>1</td><td>10<sup>-2</sup></td><td>2.40 &plusmn; 0.06</td><td>0.775 &plusmn; 0.02</td><td>2.87 &plusmn; 0.5</td><td>3.67 &plusmn; 0.6</td><td>0.613 &plusmn; 0.04</td></tr><tr><td>1</td><td>10<sup>-1</sup></td><td>2.30 &plusmn; 0.05</td><td>0.775 &plusmn; 0.02</td><td>1.59 &plusmn; 0.2</td><td>1.71 &plusmn; 0.2</td><td>0.7 &plusmn; 0.05</td></tr><tr><td>1</td><td>1</td><td>2.92 &plusmn; 0.1</td><td>0.674 &plusmn; 0.02</td><td>0.25 &plusmn; 0.06</td><td>0.39 &plusmn; 0.05</td><td>1.42 &plusmn; 0.1</td></tr><tr><td>1</td><td>10</td><td>4.90 &plusmn; 0.03</td><td>0.518 &plusmn; 0.005</td><td>0.04 &plusmn; 0.002</td><td>0.125 &plusmn; 0.004</td><td>3.13 &plusmn; 0.03</td></tr><tr><td>1</td><td>10<sup>2</sup></td><td>11.0 &plusmn; 0.04</td><td>0.384 &plusmn; 0.003</td><td>0.01 &plusmn; 3 &times; 10<sup>-4</sup></td><td>0.03 &plusmn; 0.001</td><td>7.61 &plusmn; 0.03</td></tr><tr><td>10</td><td>10<sup>-4</sup></td><td>2.23 &plusmn; 0.06</td><td>0.8 &plusmn; 0.01</td><td>14.3 &plusmn; 2</td><td>13.8 &plusmn; 2</td><td>0.367 &plusmn; 0.03</td></tr><tr><td>10</td><td>10<sup>-3</sup></td><td>2.23 &plusmn; 0.09</td><td>0.8 &plusmn; 0.01</td><td>10.3 &plusmn; 1</td><td>10.0 &plusmn; 1</td><td>0.376 &plusmn; 0.04</td></tr><tr><td>10</td><td>10<sup>-2</sup></td><td>2.26 &plusmn; 0.1</td><td>0.819 &plusmn; 0.02</td><td>7.83 &plusmn; 1</td><td>7.23 &plusmn; 1</td><td>0.446 &plusmn; 0.04</td></tr><tr><td>10</td><td>10<sup>-1</sup></td><td>2.35 &plusmn; 0.09</td><td>0.762 &plusmn; 0.02</td><td>1.35 &plusmn; 0.3</td><td>1.32 &plusmn; 0.2</td><td>0.741 &plusmn; 0.07</td></tr><tr><td>10</td><td>1</td><td>2.74 &plusmn; 0.2</td><td>0.738 &plusmn; 0.03</td><td>0.22 &plusmn; 0.8</td><td>0.21 &plusmn; 0.9</td><td>0.86 &plusmn; 0.3</td></tr><tr><td>10</td><td>10</td><td>6.05 &plusmn; 1</td><td>0.750 &plusmn; 0.1</td><td>0.18 &plusmn; 0.4</td><td>0.14 &plusmn; 0.3</td><td>3.14 &plusmn; 1</td></tr><tr><td>10</td><td>10<sup>2</sup></td><td>13.9 &plusmn; 3</td><td>0.893 &plusmn; 0.1</td><td>0.07 &plusmn; 0.3</td><td>0.05 &plusmn; 0.3</td><td>9.36 &plusmn; 3</td></tr></tbody></table>

Note. All runs have Re = 10⁴ (equivalently, Re<sub>F</sub> = 100), and a resolution of 960 × 384 × 192 effective mesh points (equivalently, 320 × 128 × 64 Fourier modes). The first column reports Re⁻¹/²Pe = 0.01Pe = Pe<sub>F</sub>; the second column reports ReRi = 10⁴Ri = Ri<sub>F</sub>. The third column is the rms velocity of the flow defined in (25), the fourth column is the vertical scaleheight defined in (27), the fifth column is minus the time- and volume-averaged turbulent heat flux, the sixth column is minus the time- and volume-averaged turbulent compositional flux, and the seventh column is the time- and volume-averaged input power (see Section 6.4).

format, together with some of the mean properties of the turbulent solutions. We have also modified the code to solve instead the LPN momentum Equation (18) together with the continuity equation, and have run a number of simulations with this new setup in a slightly smaller domain of size L<sub>x</sub> = 10π, L<sub>y</sub> = 2π and L<sub>z</sub> = 2π. The latter are summarized in Table 2. Note that the onset of instability in these simulations was already discussed by Garaud et al. (2015a), although we have also added new ones at larger values of the Richardson-Péclet number.

## 3.2. Typical Results at Low Péclet number

We first look at typical results in the limit of low Péclet numbers, for Re = 10⁴ or equivalently Re<sub>F</sub> = 100 (our largest Reynolds number value). As we shall demonstrate in Section 4 the behavior of all runs with the standard equations at Re = 10⁴ and Pe ≤ 10 is actually statistically indistinguishable from that of the LPN equations at the same Re. In all unstable cases, after a short transient period, the flow settles into a well-defined statistically stationary state. The behavior of the solutions in that state can loosely be classified into three categories: those with Ri → 0, which are all very similar to the unstratified limit Ri = 0, those where Ri is somewhat larger, for which the effects of stratification becomes important, and those where Ri exceeds unity, which are well into the linearly stable region of parameter space, and for which new dynamics emerge.

### 3.2.1. The Nearly Unstratified Limit (Ri → 0)

Figure 1 shows snapshots of the vertical velocity field $\tilde{u}$, of the horizontal velocity field $\tilde{u}$ and of the temperature perturbations away from the linearly stratified background, $\tilde{T}$, once the forced

sheared system has reached a statistically stationary state. In this simulation, Pe = 10 and Ri = 10⁻⁶. The very weak stratification, compounded by the relatively important effect of thermal diffusion, implies that the vertical motions proceed un-impeded, and have more-or-less the same scaleheight as that of the imposed forcing (see Section 4 for a more quantitative estimate of that scaleheight). The scale of the temperature perturbations is commensurate with that of the domain. The horizontal velocity perturbations are significantly stronger than the mean flow, which is hard to recognize in the snapshot. Instead, we see a shear layer that meanders spatially and temporally. The horizontally averaged and time-averaged mean velocity, defined as

$$ \tilde{u}(z) = \frac{1}{t_2 - t_1} \int_{t_1}^{t_2} \langle \tilde{u} \rangle_h(z, t) dt \quad \text{where} \\ \langle \tilde{u} \rangle_h(z, t) = \frac{1}{L_x L_y} \iint \tilde{u}(x, y, z, t) dxdy \quad (24) $$

(where the time interval [t₁, t₂] over which the data is averaged is taken once the simulation has reached a statistically stationary state) is however well-defined. It is very close to being perfectly sinusoidal, is in phase with the forcing, and its amplitude is of order Re⁻¹/².

### 3.2.2. The Stratified Regime

Figure 2 shows snapshots of the same fields as in Figure 1 for a simulation with Pe = 10 and Ri = 0.01 so RiPe = 0.1. The results are significantly different from those obtained in the nearly unstratified limit. We see for instance that the scale of the vertical velocity fluctuations is much smaller than before, and the same is true for the temperature fluctuations. The horizontal velocity perturbations are, this time, significantly
<table><caption>**Table 2**
Summary of the Main Results for All the Runs with the LPN Equations</caption><thead><tr><th>Re<sup>1/2</sup><br>Re<sub>F</sub></th><th>Re<sup>1/2</sup>RiPe<br>Ri<sub>F</sub>Pe<sub>F</sub></th><th>N<sub>x</sub>, N<sub>y</sub>, N<sub>z</sub></th><th>Re<sup>1/2</sup> ˙u<sub>rms</sub><br>  ˙u<sub>rms</sub></th><th>l̇<sub>v</sub><br>l̇<sub>v</sub></th><th>-Re⟨ŵT⟩/Pe<br>-⟨ŵT̂⟩/Pe<sub>F</sub></th><th>Re<sup>1/2</sup> ˙S<br>  ˙S</th></tr></thead><tbody><tr><td>10</td><td>10<sup>-3</sup></td><td>480, 96, 96</td><td>2.40 ± 0.008</td><td>0.991 ± 0.002</td><td>23.7 ± 0.351</td><td>0.756 ± 0.1</td></tr><tr><td>10</td><td>10<sup>-2</sup></td><td>480, 96, 96</td><td>2.06 ± 0.002</td><td>0.996 ± 3 × 10<sup>-5</sup></td><td>16.5 ± 1.20 × 10<sup>-2</sup></td><td>0.92 ± 1 × 10<sup>-5</sup></td></tr><tr><td>10</td><td>10<sup>-1</sup></td><td>480, 96, 96</td><td>2.12 ± 0.2</td><td>0.944 ± 0.03</td><td>1.84 ± 0.9</td><td>1.81 ± 1</td></tr><tr><td>10</td><td>1</td><td>480, 96, 96</td><td>3.74 ± 0.4</td><td>0.982 ± 0.009</td><td>0.943 ± 0.504</td><td>4 ± 2</td></tr><tr><td>10</td><td>1.5</td><td>480, 96, 96</td><td>5.16 ± 0.0001</td><td>0.991 ± 3 × 10<sup>-6</sup></td><td>0.568 ± 9.00 × 10<sup>-5</sup></td><td>6.67 ± 1 × 10<sup>-5</sup></td></tr><tr><td>10</td><td>2</td><td>480, 96, 96</td><td>6.40 ± 0.005</td><td>0.993 ± 3 × 10<sup>-5</sup></td><td>0.193 ± 4.00 × 10<sup>-3</sup></td><td>8.7 ± 0.1</td></tr><tr><td>10</td><td>2.2</td><td>480, 96, 96</td><td>6.85 ± 0.001</td><td>0.994 ± 1 × 10<sup>-6</sup></td><td>6.30 × 10<sup>-3</sup> ± 0.0003</td><td>9.6 ± 0.1</td></tr><tr><td>33.3</td><td>10<sup>-3</sup></td><td>480, 96, 96</td><td>2.18 ± 0.1</td><td>0.901 ± 0.02</td><td>13.1 ± 4</td><td>0.86 ± 0.3</td></tr><tr><td>33.3</td><td>10<sup>-2</sup></td><td>480, 96, 96</td><td>2.18 ± 0.1</td><td>0.910 ± 0.02</td><td>9.53 ± 3</td><td>1.05 ± 0.5</td></tr><tr><td>33.3</td><td>10<sup>-1</sup></td><td>480, 96, 96</td><td>2.34 ± 0.1</td><td>0.850 ± 0.02</td><td>1.64 ± 0.6</td><td>1.83 ± 0.9</td></tr><tr><td>33.3</td><td>1</td><td>720, 144, 144</td><td>3.42 ± 0.1</td><td>0.731 ± 0.01</td><td>0.202 ± 0.04</td><td>3.5 ± 0.8</td></tr><tr><td>33.3</td><td>3</td><td>720, 144, 144</td><td>4.25 ± 0.1</td><td>0.665 ± 0.01</td><td>9.10 × 10<sup>-2</sup> ± 0.01</td><td>4.67 ± 0.5</td></tr><tr><td>33.3</td><td>6</td><td>720, 144, 144</td><td>5.42 ± 0.3</td><td>0.617 ± 0.01</td><td>5.80 × 10<sup>-2</sup> ± 0.02</td><td>6.2 ± 1</td></tr><tr><td>33.3</td><td>7</td><td>720, 144, 144</td><td>16.1 ± 7</td><td>0.745 ± 0.2</td><td>6.20 × 10<sup>-2</sup> ± 0.3</td><td>18 ± 10</td></tr><tr><td>33.3</td><td>8</td><td>720, 144, 144</td><td>23.5 ± 0.01</td><td>0.992 ± 1 × 10<sup>-7</sup></td><td>7.31 × 10<sup>-3</sup> ± 1 × 10<sup>-5</sup></td><td>33 ± 0.001</td></tr><tr><td>50</td><td>10<sup>-3</sup></td><td>480, 96, 96</td><td>2.19 ± 0.1</td><td>0.857 ± 0.02</td><td>11.7 ± 3</td><td>0.841 ± 0.3</td></tr><tr><td>50</td><td>10<sup>-2</sup></td><td>480, 96, 96</td><td>2.36 ± 0.1</td><td>0.845 ± 0.02</td><td>5.68 ± 2</td><td>1.28 ± 0.4</td></tr><tr><td>50</td><td>10<sup>-1</sup></td><td>480, 96, 96</td><td>2.42 ± 0.1</td><td>0.798 ± 0.02</td><td>1.57 ± 0.5</td><td>1.86 ± 0.5</td></tr><tr><td>50</td><td>1</td><td>480, 96, 96</td><td>3.30 ± 0.1</td><td>0.678 ± 0.01</td><td>0.17 ± 0.03</td><td>3.33 ± 1</td></tr><tr><td>50</td><td>3</td><td>720, 144, 144</td><td>3.99 ± 0.1</td><td>0.627 ± 0.009</td><td>8.89 × 10<sup>-2</sup> ± 0.008</td><td>4.25 ± 0.7</td></tr><tr><td>50</td><td>6</td><td>720, 144, 144</td><td>4.73 ± 0.1</td><td>0.586 ± 0.008</td><td>6.00 × 10<sup>-2</sup> ± 0.005</td><td>5.32 ± 1</td></tr><tr><td>50</td><td>10</td><td>720, 144, 144</td><td>5.58 ± 0.1</td><td>0.552 ± 0.007</td><td>4.40 × 10<sup>-2</sup> ± 0.004</td><td>6.4 ± 1</td></tr><tr><td>50</td><td>15</td><td>720, 144, 144</td><td>6.51 ± 0.1</td><td>0.522 ± 0.006</td><td>3.30 × 10<sup>-2</sup> ± 0.004</td><td>7.6 ± 1</td></tr><tr><td>50</td><td>20</td><td>720, 144, 144</td><td>7.34 ± 0.2</td><td>0.500 ± 0.006</td><td>2.70 × 10<sup>-2</sup> ± 0.005</td><td>8.5 ± 1</td></tr><tr><td>50</td><td>25</td><td>720, 144, 144</td><td>10.9 ± 0.1</td><td>0.487 ± 0.007</td><td>2.70 × 10<sup>-2</sup> ± 0.003</td><td>9.47 ± 1</td></tr><tr><td>50</td><td>30</td><td>720, 144, 144</td><td>11.8 ± 0.1</td><td>0.473 ± 0.006</td><td>2.26 × 10<sup>-2</sup> ± 0.02</td><td>10.2 ± 1</td></tr><tr><td>50</td><td>35</td><td>720, 144, 144</td><td>12.5 ± 0.1</td><td>0.46 ± 0.006</td><td>1.97 × 10<sup>-2</sup> ± 0.002</td><td>11 ± 1</td></tr><tr><td>50</td><td>40</td><td>720, 144, 144</td><td>13.3 ± 0.1</td><td>0.452 ± 0.006</td><td>1.73 × 10<sup>-2</sup> ± 0.002</td><td>11.5 ± 1</td></tr><tr><td>50</td><td>45</td><td>720, 144, 144</td><td>14 ± 0.2</td><td>0.442 ± 0.006</td><td>1.51 × 10<sup>-2</sup> ± 0.004</td><td>12.4 ± 1</td></tr><tr><td>100</td><td>10<sup>-4</sup></td><td>960, 162, 162</td><td>2.27 ± 0.07</td><td>0.794 ± 0.02</td><td>14.8 ± 4</td><td>0.83 ± 2</td></tr><tr><td>100</td><td>10<sup>-3</sup></td><td>960, 162, 162</td><td>2.27 ± 0.1</td><td>0.796 ± 0.02</td><td>11.9 ± 3</td><td>0.86 ± 0.4</td></tr><tr><td>100</td><td>10<sup>-2</sup></td><td>960, 162, 162</td><td>2.39 ± 0.1</td><td>0.839 ± 0.02</td><td>13 ± 2</td><td>1.1 ± 0.4</td></tr><tr><td>100</td><td>10<sup>-1</sup></td><td>960, 162, 162</td><td>2.42 ± 0.1</td><td>0.740 ± 0.02</td><td>1.73 ± 0.5</td><td>1.84 ± 0.6</td></tr><tr><td>100</td><td>1</td><td>960, 162, 162</td><td>3.41 ± 0.05</td><td>0.601 ± 0.01</td><td>0.148 ± 0.02</td><td>3.5 ± 0.6</td></tr><tr><td>100</td><td>10</td><td>960, 162, 162</td><td>5.11 ± 0.04</td><td>0.494 ± 0.01</td><td>4.00 × 10<sup>-2</sup> ± 0.03</td><td>5.78 ± 0.7</td></tr><tr><td>100</td><td>50</td><td>960, 162, 162</td><td>8.51 ± 0.03</td><td>0.395 ± 0.003</td><td>1.50 × 10<sup>-2</sup> ± 0.0005</td><td>10 ± 1</td></tr><tr><td>100</td><td>100</td><td>960, 162, 162</td><td>11.3 ± 0.05</td><td>0.354 ± 0.03</td><td>9.50 × 10<sup>-3</sup> ± 0.0004</td><td>13.1 ± 1</td></tr><tr><td>100</td><td>120</td><td>960, 162, 162</td><td>12.3 ± 0.07</td><td>0.343 ± 0.003</td><td>8.30 × 10<sup>-3</sup> ± 0.0005</td><td>14.1 ± 1</td></tr><tr><td>100</td><td>130</td><td>960, 162, 162</td><td>12.8 ± 0.06</td><td>0.338 ± 0.002</td><td>7.80 × 10<sup>-3</sup> ± 0.0004</td><td>14.6 ± 1</td></tr><tr><td>100</td><td>140</td><td>960, 162, 162</td><td>13.3 ± 0.08</td><td>0.334 ± 0.003</td><td>7.40 × 10<sup>-3</sup> ± 0.0005</td><td>15 ± 1</td></tr><tr><td>100</td><td>150</td><td>960, 162, 162</td><td>13.7 ± 0.1</td><td>0.330 ± 0.002</td><td>7.00 × 10<sup>-3</sup> ± 0.0006</td><td>15.7 ± 0.9</td></tr><tr><td>100</td><td>160</td><td>960, 162, 162</td><td>14.5 ± 0.7</td><td>0.326 ± 0.003</td><td>6.70 × 10<sup>-3</sup> ± 0.003</td><td>16.1 ± 1</td></tr><tr><td>100</td><td>170</td><td>960, 162, 162</td><td>20.8 ± 0.04</td><td>0.325 ± 0.002</td><td>7.40 × 10<sup>-3</sup> ± 0.0003</td><td>17.5 ± 1</td></tr><tr><td>100</td><td>180</td><td>960, 162, 162</td><td>21.3 ± 0.03</td><td>0.320 ± 0.002</td><td>7.00 × 10<sup>-3</sup> ± 0.0003</td><td>17.6 ± 1</td></tr></tbody></table>

Note. The first column reports $Re_F = Re^{1/2}$; the second column reports $Ri_F Pe_F = Re^{1/2}RiPe$. The third column is the number of equivalent mesh points in each direction. The fourth column is the rms velocity of the flow defined in (25). The fifth column is the vertical scaleheight defined in (27). The sixth column is minus the time- and volume-averaged turbulent heat flux divided by $Pe_F$ (see Section 6.4). The seventh column is the mean shear at $z = \pi$ (see Section 6.2).

weaker than the mean flow, and the latter is clearly recognizable in the snapshot. The amplitude of the mean flow is now much larger and its shape is no longer sinusoidal, but instead appears to tend to a piecewise linear profile. It is still essentially in phase with the forcing, however.

### 3.2.3. The Very Strongly Stratified Limit

As discussed by Garaud et al. (2015a), turbulent solutions exist for values of Ri beyond the threshold for linear stability (which is such that $RiPe \simeq 0.25$ at very large Reynolds numbers). These solutions have been obtained by careful

continuation of previous solutions, progressively increasing Ri (equivalently, RiPe in the LPN equations) across the threshold and beyond. Using the LPN equations, we have in fact been able to push into the linearly stable regime somewhat further than Garaud et al. (2015a) did, and found that, for very strongly stratified flows ($RiPe > 1.6$ at $Re = 10^4$), the system dynamics change once again quite dramatically. This is illustrated in Figure 3 (see also Figure 9), which shows that the mean flow no longer has the same symmetries as the imposed force, and adopts instead a new skewed state where the minima and maxima are shifted away from their original positions. This shift effectively enlarges one of the regions of near-constant
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Simulations snapshots for a run using the standard equations with Re = 10⁴, Pe = 10, and Ri = 10⁻⁶ (equivalently, Re<sub>F</sub> = 100, Pe<sub>F</sub> = 0.1, and Ri<sub>F</sub> = 0.01). From top left to bottom right, we see the vertical velocity fluctuations, the temperature fluctuations, the horizontal velocity fluctuations in the x direction, and the horizontally averaged mean flow in the x direction. The dashed green line in the fourth quadrant shows the function sin(z) times the amplitude u<sub>max</sub> = max(û) = Re¹/² max(û̃) of the mean flow for comparison with the form of the imposed force.</figcaption></figure>

shear, and creates two thinner (and therefore stronger) shear layers on either side. Surprisingly, we see from the flow snapshots that the stronger/thinner shear layers become laminar, while the turbulence is confined to the weaker/wider one. Why this symmetry breaking occurs, and what stabilizes the region of strongest shear, remain to be determined.

### 3.3. Typical Results at High Péclet Number

We now look at typical results in the limit of very large Péclet numbers, here for Pe = 1000. The nearly unstratified limit in that case is still achieved whenever RiPe < 10⁻³ (e.g., for Ri < 10⁻⁶ if Pe = 1000), and is the same as the one for low Péclet numbers described in Section 3.2.1. For larger Ri, however, we observe that instead of settling into a relatively regular statistically stationary state with weak fluctuations around the mean, the system adopts a quasi-periodic behavior that cycles between intense mixing events that destroy the existing shear and render it too weak to maintain turbulence, and fairly quiescent periods during which the system is close to laminar and where the forcing gradually amplifies the shear (see Figure 4). During these laminar periods, the instantaneous rms flow velocity grows linearly with time. This regime will be discussed in more detail in a forthcoming publication. Indeed, while interesting for many reasons, it has only been observed to exist so far in the case where Pe is extremely large, and is therefore somewhat outside of the scope of this paper. It is worth noting, however, that similar quasi-periodically bursting

stratified shear flows have recently been observed in laboratory experiments by Meyer & Linden (2014).

### 4. THE VALIDITY OF THE LPN EQUATIONS

We now turn to a more quantitative analysis of our simulations. In this section, we focus on comparing the results obtained using the LPN equations with those obtained using the standard equations.

#### 4.1. Comparison between the LPN Equations and the Full Equations

In order to compare simulations using the LPN equations and those using the standard equations more quantitatively, we now look at typical global properties of the turbulent flow, such as the rms velocity, and the typical vertical lengthscale of the vertical fluid motions in each of the available simulations.

From here on, we define û<sub>rms</sub> as the non-dimensional rms velocity of the flow,

$$ \tilde{u}_{\text{rms}} = \langle |\tilde{u}(\mathbf{x}, t)|^2 \rangle^{1/2}, \qquad (25) $$

where the ⟨·⟩ notation will be used hereafter to denote any volume- and time-average, as in

$$ \langle q \rangle = \frac{1}{t_2 - t_1} \int_{t_1}^{t_2} \left[ \frac{1}{L_x L_y L_z} \iiint q(\mathbf{x}, t) d^3\mathbf{x} \right] dt \qquad (26) $$
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Simulation snapshot for a run at Re = 10⁴, Pe = 10, and Ri = 10⁻² (equivalently, Re<sub>F</sub> = 100, Pe<sub>F</sub> = 0.1, and Ri<sub>F</sub> = 100). From top left to bottom right, we see the vertical velocity fluctuations, the temperature fluctuations, the horizontal velocity fluctuations in the x direction, and the horizontally averaged mean flow in the x direction. The dashed green line in the fourth quadrant shows the function sin(z) times the amplitude u<sub>max</sub> = max(û) = Re¹/² max(û) of the mean flow for comparison with the form of the imposed force.</figcaption></figure>

where the time interval [t₁, t₂] is chosen as in Equation (24). We also define ˜l<sub>v</sub> as the typical non-dimensional vertical scale of the energy-bearing vertical fluid motions (see, for instance, Batchelor 1953), with

$$ \tilde{l}_v = \frac{1}{t_2 - t_1} \int_{t_1}^{t_2} \left[ \frac{\sum \tilde{k}_x \sum \tilde{k}_y \sum \tilde{k}_z \neq 0 E_w(\tilde{k}, t) \tilde{k}_z^{-1}}{\sum \tilde{k}_x \sum \tilde{k}_y \sum \tilde{k}_z \neq 0 E_w(\tilde{k}, t)} \right] dt, \quad (27) $$

where E<sub>w</sub>(k̃, t) is the instantaneous kinetic energy of the vertical velocity field associated with wavenumber k̃ = (k̃<sub>x</sub>, k̃<sub>y</sub>, k̃<sub>z</sub>). Note that there are other ways of defining ˜l<sub>v</sub> (see Section 6 for detail). In particular, we could have defined it using a weighted average of k̃<sub>z</sub><sup>-1</sup> over the total kinetic energy spectrum instead, rather than the one based on the vertical motions only. However, that definition has the disadvantage of being contaminated by the vertical scale of the mean flow, or, should they be present, by vortices whose structure is invariant in z and whose vorticity is aligned with the z axis. Neither of these flows participate in vertical transport processes, however. By focussing on vertical fluid motions in Equation (27), by contrast, we only take into account flows that affect vertical transport. In what follows, we shall refer to ˜l<sub>v</sub> for simplicity as the “vertical eddy lengthscale,” bearing in

mind that it really is the vertical lengthscale of vertical fluid motions.

The results are presented in Tables 1 and 2, and shown in Figure 5, for the standard equations at various Péclet numbers, and for the LPN equations. For ease of comparison between the datasets, we plot ˜l<sub>v</sub> (equivalently, ˆl<sub>v</sub>) and 100 ˜u<sub>rms</sub> (equivalently, ˆu<sub>rms</sub>) as functions of the Richardson–Péclet number. Figure 5 reveals a number of interesting facts.

First, we see that for each curve (corresponding to each value of Peinvestigated), both ˜l<sub>v</sub> and ˆu<sub>rms</sub> asymptote to the respective constants ˜l<sub>0</sub> ≃ 0.8 and ˆu<sub>0</sub> ≃ 0.022 in the limit of very low RiPe, and these constants appear to be independent of the Péclet number. They simply represent the turbulent properties of an unstratified sinusoidally forced shear flow at Re = 10⁴. Second, we see that the effect of the stratification becomes relevant even for the LPN equations when RiPe is greater than about 10⁻³. As RiPe increases, the typical vertical eddy scale decreases, and the rms velocity increases. The system appears to be stable to finite-amplitude instabilities for RiPe greater than a few, as discussed by Garaud et al. (2015a). We shall study the properties of the turbulent solutions in Section 6 in more detail.

In both cases (for ˜l<sub>v</sub> and for ˆu<sub>rms</sub>), we see quite clearly that the LPN equations are a good-to-excellent approximation of the standard equations up to Pe = 100, but no longer for Pe = 1000. This shows that the domain of validity of the LPN equations extends well beyond Pe = 1. As discussed by
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Left: snapshots of $\tilde{u}$ and $\tilde{w}$ (as a function of $(x, z)$ in the $y=0$ plane), obtained with LPN equations at $Re = 10^4$ and $RiPe = 1.8$ (equivalently, $Re_F = 100$ and $Ri_F Pe_F = 180$). Right: the mean flow in the $x$ direction (red line), and for comparison, the corresponding laminar solution (green dashed line). The mean flow is clearly no longer symmetric but becomes skewed. The regions of strongest mean shear are essentially laminar, while the turbulence subsists in the region of weaker shear (near the middle of the box).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Left: snapshots of the horizontal velocity field $\tilde{u}$ at selected times, in a simulation with $Re = 10^4$, $Pe = 10^3$, and $Ri = 0.01$ (equivalently, $Re_F = 100$, $Pe_F = 10$, and $Ri_F = 100$). Right: evolution of the instantaneous rms flow velocity $\tilde{u}_{\text{rms}}(t) = [(L_x L_y L_z)^{-1} \iiint \|\tilde{u}(x, t)\|^2 d^3x]^{1/2}$, as well as the instantaneous vertical turbulent temperature flux $\tilde{F}_T(t) = (L_x L_y L_z)^{-1} \iiint \tilde{w} \tilde{T} d^3x$. This is a typical example of high Péclet number bursting behavior. The time intervals where $\tilde{u}_{\text{rms}}(t)$ grows linearly correspond to the times where the fluid is essentially laminar, while the sharp drops correspond to brief turbulent mixing events. The snapshot illustrates the dynamics around the second mixing event.</figcaption></figure>

Garaud et al. (2015a), this is not entirely surprising. Indeed, $Pe$ is not a particularly good estimate of the actual turbulent Péclet number $Pe_t = u_{\text{rms}}l_v/\kappa_T = \tilde{u}_{\text{rms}}\tilde{l}_v$. $Pe$ of the unstable flow once a statistically stationary state has been reached. This is demonstrated in Figure 6. We see that $Pe_t$ is significantly smaller than one for $Pe = 10$, and remains of order one even

for $Pe = 100$. This naturally explains the good correspondence of the LPN simulations and of the standard simulations up to $Pe = 100$ seen in Figure 5, and confirms the argument of Lignières (1999), that the LPN solutions are a good approximation to the true solution provided the turbulent Péclet number $Pe_t$ is of order one or smaller.
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Mean values of the vertical eddy lengthscale $\tilde{l}_v$ (equivalently, $\hat{l}_v$) and of the rms velocity $100\hat{u}_{rms}$ (equivalently, $\hat{u}_{rms}$), as a function of RiPe and RiFPeF, for Re = 10⁴ (equivalently, ReF = 100). The error bars represent the rms fluctuations around the mean value. The LPN equations are a good approximation to the full equations for up to Pe = 100 (equivalently, PeF = 1). In the limit where RiFPeF ≪ 1, we find that $\tilde{l}_v \simeq 0.8$ and $\hat{u}_{rms} \simeq 0.022$ (equivalently, $\hat{u}_{rms} = 2.2$).</figcaption></figure>

## 4.2. Estimating the Turbulent Péclet Number

As discussed in Section 1, we would like to estimate without having to run numerical simulations, for a given dimensional forcing amplitude $F_0$, shear lengthscale $k^{-1}$, viscosity $\nu$ and thermal diffusivity $\kappa_T$, what PeT may be. As it turns out, this is actually fairly easy. Indeed, assuming that in the statistically stationary state characteristic of low Péclet number simulations there is a balance in the momentum equation between the inertial term and the forcing, such that

$$ \rho_0(\mathbf{u} \cdot \nabla \mathbf{u}) \cdot \mathbf{e}_x \sim F_0 \quad (28) $$

in dimensional terms, then we can define a new characteristic dimensional flow amplitude $U_F$ as

$$ U_F = \left( \frac{F_0}{k\rho_0} \right)^{1/2} . \quad (29) $$

This is precisely the velocity introduced in Section 2, which is at the heart of the proposed alternative non-dimensionalization of the system. We now revisit the simulation results, in the units based on $U_F$ (where the non-dimensional quantities are hatted).

Figure 5 shows $l_v$ and $u_{rms}$ plotted in this new set of units and against RiFPeF (see the top axis). As we can see, $\hat{u}_{rms}$ is now of order unity for most simulations at low enough RiFPeF, and only grows slowly with RiFPeF for RiFPeF ≫ 1. In other words, $U_F$ seems to be a good predictor for the dimensional rms velocity of the turbulent shear flow. Furthermore, the apparent transition from weakly stratified to strongly stratified regimes now occurs for RiFPeF ~ 0.1–1, which shows that RiFPeF is a more meaningful bifurcation parameter than RiPe (for which the same transition occurs around 0.01). Since $U_F$ is a good predictor for $u_{rms}$, and since $l_v$ is not too different from $k^{-1}$, we then predict that $Pe_F = U_F/k\kappa_T$ should be a good predictor for PeT = $u_{rms}l_v/\kappa_T$. We see in Figure 6 that this is indeed the case. This is especially true in the weakly stratified limit, where $Pe_T \sim 1.75Pe_F$ but remains also true within an order of magnitude for larger RiFPeF. This is a very convenient result since $Pe_F$ is calculated purely from quantities that are

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Turbulent Péclet number $Pe_T = u_{rms}l_v/\kappa_T = \hat{u}_{rms}\tilde{l}_vPe = \hat{u}_{rms}\hat{l}_vPe_F$ as a function of RiPe. The horizontal lines show the corresponding value $Pe_F$ for each set of simulations, showing that $Pe_T$ is a reasonably good estimate of $Pe_F$.</figcaption></figure>

readily known, rather than quantities that have to be determined a posteriori from the numerical simulations. We can now use it to determine which stars are likely to harbor low Péclet number shear layers.

## 5. RELEVANCE OF DIFFUSIVE SHEAR INSTABILITIES IN STARS

As found in Section 4, in order to determine whether diffusive—or secular—shear instabilities could be relevant in stellar interiors, one merely has to calculate $Pe_F$ using the available information (strength of perturbing force, local diffusivities and local Brunt–Väisälä frequency) and see when the latter is smaller than one. The diffusivities and the local Brunt–Väisälä frequency are known from the background properties of the stellar model. The perturbing force, however, depends on the situation considered—whether the shear is induced by the differential contraction and expansion of the
star, by tidal forces or other—and can vary by orders of magnitude accordingly. Bearing this in mind, we can get very rough estimates of the typical values of Pe<sub>F</sub>, Re<sub>F</sub> and Ri<sub>F</sub> in the interiors of solar-type stars as follows:

$$
\begin{align}
\mathrm{Pe}_F &= 10^9 \left( \frac{F_0/\rho_0}{10^5 \mathrm{cm} \, \mathrm{s}^{-2}} \right)^{1/2} \left( \frac{k^{-1}}{10^9 \mathrm{cm}} \right)^{3/2} \left( \frac{\kappa_T}{10^7 \mathrm{cm}^2 \, \mathrm{s}^{-1}} \right)^{-1}, \nonumber \\
\mathrm{Re}_F &= 10^{15} \left( \frac{F_0/\rho_0}{10^5 \mathrm{cm} \, \mathrm{s}^{-2}} \right)^{1/2} \left( \frac{k^{-1}}{10^9 \mathrm{cm}} \right)^{3/2} \left( \frac{\nu}{10 \mathrm{cm}^2 \, \mathrm{s}^{-1}} \right)^{-1}, \nonumber \\
\mathrm{Ri}_F &= 0.01 \left( \frac{N^2}{10^{-6} \, \mathrm{s}^{-2}} \right) \left( \frac{k^{-1}}{10^9 \mathrm{cm}} \right) \left( \frac{F_0/\rho_0}{10^5 \mathrm{cm} \, \mathrm{s}^{-2}} \right)^{-1}, \tag{30}
\end{align}
$$

where the numerical values chosen for comparison for $k^{-1}$, $\kappa_T$ and $\nu$ are order-of-magnitude estimates of the solar tachocline properties (Gough 2007), a well-known shear layer at the interface between the convection zone and the radiative zone of the Sun (Hughes et al. 2007). Meanwhile $F_0/\rho_0$ is compared with its local gravity $g \simeq 10^5 \mathrm{cm} \, \mathrm{s}^{-2}$.

This clearly shows that, unless $F_0/\rho_0$ is many orders of magnitude smaller than $g$, a tachocline-like shear cannot be viewed as a low Péclet number flow. And, supposing that $F_0/\rho_0$ were to be of the order of $10^{-13} \mathrm{cm} \, \mathrm{s}^{-2}$, then $Ri_F$ would be of the order of $10^{20}$. This shows that diffusive shear instabilities are unlikely to play a role in the dynamics of the solar tachocline, and that the latter is most probably stable to its vertical shear.<fnref n="5" /> Another way to have a low $Pe_F$ would be to consider the possibility of very thin shear layers. Using again the tachocline value of $\kappa_T \sim 10^7 \mathrm{cm}^2 \, \mathrm{s}^{-1}$, and taking $F_0/\rho_0 \sim 0.1 \mathrm{cm} \, \mathrm{s}^{-2}$, would require the shear layer to be no thicker than $k^{-1} \sim 10^6 \mathrm{cm}$ (10 km) to be diffusive. Either way, it appears that diffusive shear instabilities are unlikely in the interior of solar-type stars.

In much more massive stars, however, the thermal diffusivity can be much larger. This is illustrated in Figure 7, which shows $\kappa_T$ as a function of radius for a range of stars on the Main Sequence. Here, $\kappa_T$ is calculated as in Garaud et al. (2015b). We see that for stars larger than about $10 M_\odot$, an increasingly large fraction of their outer envelope has very high thermal diffusivities, with values of $10^{14} \mathrm{cm}^2 \, \mathrm{s}^{-1}$ or larger. The reason why $\kappa_T$ increases with stellar mass is mostly due to the fact that higher-mass stars typically have lower densities and higher temperatures than lower-mass stars at the same fractional mass coordinate.

With the same scalings as before, but with $\kappa_T$ compared with $10^{14} \mathrm{cm}^2 \, \mathrm{s}^{-1}$ instead, we get

$$
\mathrm{Pe}_F = 100 \left( \frac{F_0/\rho_0}{10^5 \mathrm{cm} \, \mathrm{s}^{-2}} \right)^{1/2} \left( \frac{k^{-1}}{10^9 \mathrm{cm}} \right)^{3/2} \left( \frac{\kappa_T}{10^{14} \mathrm{cm}^2 \, \mathrm{s}^{-1}} \right)^{-1}. \quad (31)
$$

This time, $Pe_F$ values smaller than one can be achieved with a reasonably thin shear layer ($k^{-1} \sim 10^7 \mathrm{cm}$), or with a reasonable force of $F_0/\rho_0 \sim 10 \mathrm{cm} \, \mathrm{s}^{-2}$. This implies that diffusive shear instabilities could be relevant in the outer layers of high-mass stars. At the same time, the values of $N^2$ in the radiative layers of high-mass stars range between $10^{-8}$ and $10^{-7} \mathrm{cm}^2$, so

<figure><img src="image_10.png" /><figcaption>Figure 7. Thermal diffusivity $\kappa_T$ (in cm$^2$ s$^{-1}$) as a function of radius for stars of various masses $M_\star$ ranging from 2 $M_\odot$ to 30 $M_\odot$.</figcaption></figure>

the bifurcation parameter

$$
\mathrm{Ri}_F \mathrm{Pe}_F = 0.1 \left( \frac{N^2}{10^{-7} \mathrm{s}^{-2}} \right) \left( \frac{k^{-1}}{10^9 \mathrm{cm}} \right)^{5/2} \times \left( \frac{F_0/\rho_0}{10^5 \mathrm{cm} \, \mathrm{s}^{-2}} \right)^{-1/2} \left( \frac{\kappa_T}{10^{14} \mathrm{cm}^2 \, \mathrm{s}^{-1}} \right)^{-1}, \quad (32)
$$

will likely be in the range $10^{-3}–10^3$ depending on the exact value of $\kappa_T$, $N^2$, $F_0$ and $k$ used. In other words, $Ri_F Pe_F$ could in principle span most of parameter space between the nearly unstratified limit and the strongly stratified limit. With this in mind, we now go back to the simulations and look at their basic transport properties in more detail.

## 6. GLOBAL PROPERTIES OF THE TURBULENT LPN SOLUTIONS AT LOW AND HIGH $Ri_F Pe_F$

We now restrict our analysis to flows that have a low predicted Péclet number (with $Pe_F \le 1$), focussing on the results obtained using the LPN equations unless otherwise specified. In all that follows, we adopt the non-dimensionalization based on the forcing defined in Section 2, where the velocity is expressed in units of $U_F$, and use $Ri_F Pe_F$ and $Re_F$ as our basic input parameters.

### 6.1. Properties of the Turbulent Flow

Figure 8(a) shows the rms velocity of the flow as a function of $Ri_F Pe_F$ for various values of $Re_F$. As already discussed in Section 3, we find that the rms velocity becomes independent of $Ri_F Pe_F$ in the limit of very low stratification ($Ri_F Pe_F \ll 1$). This is expected since the buoyancy force plays an insignificant role in the momentum equation in this limit, so the characteristics of the system should only depend on the remaining parameter $Re_F$. Interestingly, we find that viscosity is already essentially negligible for $Re_F \ge 33.3$, with $\hat{u}_{\mathrm{rms}}(Ri_F Pe_F \to 0, Re_F) \simeq 2.2$. Dimensionally, this implies that $u_{\mathrm{rms}} \simeq 2.2(F_0/k\rho_0)^{1/2}$ in high $Re_F$, low $Ri_F Pe_F$ shear flows.

Figure 8(a) also shows that the rms velocity increases monotonically (albeit not always smoothly, see below) with $Ri_F Pe_F$ for all $Re_F \ge 33.3$. This trend can easily be explained:

<footnote n="5">⁵ On the other hand it may still be unstable to its horizontal shear (Watson 1981; Garaud 2001), which could in turn affect its stability to vertical shear, as discussed by Talon & Zahn (1997).</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Left: rms velocity of the flow for the LPN equations for various Re<sub>F</sub>. For each value of Re<sub>F</sub>, the vertical arrow marks the limit for linear stability of the flow, while the horizontal arrow marks the value of ̂u<sub>rms</sub> corresponding to the laminar solution. Note how the lower Re<sub>F</sub> simulations smoothly tend to the laminar solution while the higher Re<sub>F</sub> simulations follow a nonlinearly unstable turbulent branch that ignores the laminar solution entirely. Right: partition between the different contributions to the total kinetic energy of the flow, for Re<sub>F</sub> = 100 and various Ri<sub>F</sub>Pe<sub>F</sub>, in the LPN equations. In both figures, the quantities plotted are in the units based on U<sub>F</sub>.</figcaption></figure>

as the stratification increases, it becomes gradually more
difficult to destabilize the forced horizontal flow, so the shear
must grow to larger amplitudes before turbulence can set in. As
a result, the total kinetic energy in the fluid in the turbulent state
also increases. However, it is interesting to note this increase is
solely due to the increase in the amplitude of the mean flow.
Figure 8(b) demonstrates this by comparing the mean flow
kinetic energy, given by

$$
\hat{E}_{\text{mean}} = \frac{1}{L_z} \int \frac{1}{2} \hat{u}^2(z) dz, \quad (33)
$$

(where ̂u(z) is horizontal mean flow defined as ̃u(z) in
Equation (24) but expressed in the non-dimensionalization
based on U<sub>F</sub>) with that of entire flow field ̂u, given by

$$
\hat{E}_{\text{tot}} = \frac{1}{2} \langle |\hat{u}(\mathbf{x}, t)|^2 \rangle, \quad (34)
$$

and that of the perturbations, given by

$$
\hat{E}_{\text{pert}} = \hat{E}_{\text{tot}} - \hat{E}_{\text{mean}}. \qquad (35)
$$

In the limit Ri<sub>F</sub>Pe<sub>F</sub> ≪ 1, we find that the kinetic energy in the mean flow is only about 1/16 of the total kinetic energy. In this limit, most of the energy is in the perturbations. By contrast, when Ri<sub>F</sub>Pe<sub>F</sub> ≫ 1 the mean flow is the major contributor to the total kinetic energy in the system. The properties of the mean flow are discussed in more detail in Section 6.2.

Going back to Figure 8(a), we see that the manner in which
û<sub>rms</sub> increases with Ri<sub>F</sub>Pe<sub>F</sub> is notably different between the low
Reynolds number runs (Re<sub>F</sub> = 10 and Re<sub>F</sub> = 33.3) and the
high Reynolds number runs (Re<sub>F</sub> = 50 and Re<sub>F</sub> = 100). For
Re<sub>F</sub> = 10 and Re<sub>F</sub> = 33.3, û<sub>rms</sub> is a smooth function of
Ri<sub>F</sub>Pe<sub>F</sub>, and the respective curves gradually approach the point
of marginal linear stability (which is marked by the intersection
of the horizontal and vertical arrows). In other words, the
turbulence gradually dies away as the system approaches that
marginal state.

For Re<sub>F</sub> = 50 and Re<sub>F</sub> = 100, however, the linear stability
threshold is no longer relevant, and turbulent solutions extend
significantly into the linearly stable region of parameter space

(Garaud et al. 2015a). For both sets of simulations, we see that
after crossing the threshold for marginal linear stability, the
û<sub>rms</sub> (Ri<sub>F</sub>Pe<sub>F</sub>) curve has a very sharp step, followed by another
gradual incline.<fnref n="6" /> This step corresponds to the transition
discussed in Section 3.2.3 from a symmetric to a skewed flow
profile, the latter having a significantly higher total kinetic
energy. The shape of the step being somewhat reminiscent of
half a hysteresis curve, we have tried to find evidence of
multiple equilibria in this system. However, we have not found
any: the functional dependence of û<sub>rms</sub> on Ri<sub>F</sub>Pe<sub>F</sub> is the same
whether Ri<sub>F</sub>Pe<sub>F</sub> is gradually increased or gradually decreased
across the step.<fnref n="7" /> Finally, Figure 8(a) reveals that the position of
the step depends on Re<sub>F</sub>, showing that viscosity plays a role in
this transition. This is not a surprising result, given that the
transport of momentum across the strongly sheared, laminar
portion of the flow (see Section 3.2.3) has to be diffusive.

Figure 8(b) studies the various contributions to the total kinetic energy in the system, for the runs with Re<sub>F</sub> = 100. We find that even though the total kinetic energy increases with Ri<sub>F</sub>Pe<sub>F</sub>, the kinetic energy of the perturbations E<sub>pert</sub> remains more-or-less constant for all simulations. This is a rather remarkable result since Ri<sub>F</sub>Pe<sub>F</sub> varies by more than six orders of magnitude, which appears to be one of the defining properties of forced shear flows at low Péclet numbers. A closer inspection of the contributions to the turbulent kinetic energy from both streamwise horizontal motions, cross-stream horizontal motions and vertical motions, shows that even though their total remains constant, the kinetic energy in the vertical motions drops slowly with increasing Ri<sub>F</sub>Pe<sub>F</sub>, while that in the horizontal motions increases slightly to compensate. What controls the partitioning of the turbulent kinetic energy, which is always of order one in this system, into vertical and horizontal perturbations respectively, remains to be determined.

<footnote n="6">⁶ The incline is difficult to see in the case where Re<sub>F</sub> = 100, because there are only two points close to one another before the turbulent solution is lost, but it is indeed present.</footnote>

<footnote n="7">⁷ It remains possible that there are indeed multiple equilibria, but if they exist, they only do so over a range of Ri<sub>F</sub>Pe<sub>F</sub> too narrow for us to find.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Variation of the mean flow profile and amplitude as a function of Ri<sub>F</sub>Pe<sub>F</sub> for Re<sub>F</sub> = 100, using the LPN equations.</figcaption></figure>

## 6.2. Mean Flow Properties

As discussed in Section 3, for fixed Reynolds number, both the amplitude and shape of the mean flow change significantly as Ri<sub>F</sub>Pe<sub>F</sub> increases from the unstratified limit up to the critical value above which turbulent solutions cease to exist. This is illustrated in Figure 9. For low Ri<sub>F</sub>Pe<sub>F</sub>, the profiles are very close to being perfectly sinusoidal and are in phase with the forcing. As Ri<sub>F</sub>Pe<sub>F</sub> increases beyond one, the profiles become more and more triangular, but the maxima and minima of $\hat{u}(z)$ remain aligned with the maxima and minima of the imposed forcing (i.e., at $z = \pi/2$ and $z = 3\pi/2$). Finally, for Ri<sub>F</sub>Pe<sub>F</sub> = 170 (not shown) and Ri<sub>F</sub>Pe<sub>F</sub> = 180 (shown), the profiles become asymmetric, as discussed in Section 3.2.3.

Interestingly, the shear at $z = \pi$ remains a smooth function of Ri<sub>F</sub>Pe<sub>F</sub> despite the dramatic change of behavior associated with the loss of symmetry. This can be seen more clearly in Figure 10, which shows the product of Pe<sub>F</sub> with the gradient Richardson number J as a function of Ri<sub>F</sub>Pe<sub>F</sub>. J is calculated from the amplitude of the mean shear at $z = \pi$, as

$$J = \frac{Ri_F}{\hat{S}^2} \quad \text{where, here, } \hat{S} = \left| \frac{d\hat{u}}{dz} \right|_{z=\pi} \quad (36)$$

We see that for the larger Re<sub>F</sub> runs (which are the only ones we believe to be relevant for stellar interiors), JPe<sub>F</sub> grows smoothly as Ri<sub>F</sub>Pe<sub>F</sub> increases. For very low Ri<sub>F</sub>Pe<sub>F</sub>, JPe<sub>F</sub> is proportional to Ri<sub>F</sub>Pe<sub>F</sub>, which can easily be explained from the fact that the mean flow amplitude, and therefore the mean shear, are independent of Ri<sub>F</sub>Pe<sub>F</sub> in that limit (see Figure 5(b)). For intermediate values of Ri<sub>F</sub>Pe<sub>F</sub>, JPe<sub>F</sub> seems to scale like $(Ri_F Pe_F)^{1/2}$, a result which remains to be explained.

For large Ri<sub>F</sub>Pe<sub>F</sub>, JPe<sub>F</sub> appears to saturate at a value that is proportional to Re<sub>F</sub>. This result could be consistent with the notion that a strongly stratified shear flow adapts itself to satisfy a theoretical marginal stability<fnref n="8" /> criterion of the kind JPr = (JPr)<sub>c</sub> where (JPr)<sub>c</sub> is independent of Re<sub>F</sub>, as in (4) for instance. The two data points available (for Re<sub>F</sub> = 50 and 100) suggest that for large enough stratification, the flow satisfies

$$JPr \approx (JPr)_c \approx 0.006 \quad (37)$$

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Variation of JPe<sub>F</sub> as a function of Ri<sub>F</sub>Pe<sub>F</sub>. The horizontal lines mark the apparent saturation values for each curve in the nonlinearly unstable regime (Re<sub>F</sub> ≥ 50 and large Ri<sub>F</sub>Pe<sub>F</sub>), and are consistent with JPe<sub>F</sub> ≃ 0.006Re<sub>F</sub>. The line fitting the data at very low Ri<sub>F</sub>Pe<sub>F</sub> has JPe<sub>F</sub> = 1.4Ri<sub>F</sub>Pe<sub>F</sub>; the line fitting the data at intermediate values of Ri<sub>F</sub>Pe<sub>F</sub> has JPe<sub>F</sub> = 0.1(Ri<sub>F</sub>Pe<sub>F</sub>)<sup>1/2</sup>. The errorbars represent rms fluctuations around the mean JPe<sub>F</sub>, and are much larger for the low Ri<sub>F</sub>Pe<sub>F</sub> runs where the energy in the fluctuations is much larger than the energy in the mean flow (see Section 6.1).</figcaption></figure>

This is in fact remarkably similar to the criterion proposed by Zahn (1974) from heuristic energy stability arguments (see Equation (4)). It is important to note however that this result only applies to the strongly stratified, nonlinearly unstable regime: for low and intermediate values of Ri<sub>F</sub>Pe<sub>F</sub>, the observed shear layer appears to be far from marginal stability as defined by Zahn's criterion.

It is unfortunate that we have so few data points in the strongly stratified limit, which prevents us from firmly establishing the validity of (37). The forced shear layer rapidly becomes fully stable for larger values of Ri<sub>F</sub>Pe<sub>F</sub> (at least in the case of this sinusoidal forcing, see Section 7 for detail), so we are unable to determine whether (37) continues to hold for even stronger levels of stratification. Numerical simulations at larger Re<sub>F</sub> may help resolve this problem, since they are likely to exhibit instability for larger values of Ri<sub>F</sub>Pe<sub>F</sub>. However, they are computationally prohibitive to date.

Finally, we note that while the transition from JPe<sub>F</sub> ~ Ri<sub>F</sub>Pe<sub>F</sub> to JPe<sub>F</sub> ~ (Ri<sub>F</sub>Pe<sub>F</sub>)<sup>1/2</sup> is independent of viscosity and always occurs around Ri<sub>F</sub>Pe<sub>F</sub> ~ 0.01, the transition from JPe<sub>F</sub> ~ (Ri<sub>F</sub>Pe<sub>F</sub>)<sup>1/2</sup> to JPe<sub>F</sub> ~ (JPe<sub>F</sub>)<sub>c</sub> occurs at progressively larger values of Ri<sub>F</sub>Pe<sub>F</sub> as Re<sub>F</sub> increases. For this particular problem, we find that, very roughly,

$$ \begin{align} JPe_F &\approx 1.4Ri_F Pe_F \quad \text{for } Ri_F Pe_F \le 0.005, \nonumber \\ JPe_F &\approx 0.1(Ri_F Pe_F)^{1/2} \quad \text{for } 0.005 \le Ri_F Pe_F \le (0.06 Re_F)^2, \nonumber \\ JPe_F &\approx 0.006 Re_F \quad \text{for } (0.06 Re_F)^2 \le Ri_F Pe_F \nonumber \\ &\le (Ri_F Pe_F)_c \sim Re_F, \tag{38} \end{align} $$

where the numerical constants were determined by fitting the data, and the ranges of validity were determined by requiring approximate continuity of JPe<sub>F</sub> with Ri<sub>F</sub>Pe<sub>F</sub>. The upper limit for the existence of turbulent solutions, (Ri<sub>F</sub>Pe<sub>F</sub>)<sub>c</sub>, was determined numerically to be a constant of order unity times Re<sub>F</sub>, for the sinusoidal forcing selected (Garaud et al. 2015a). Taken at face value, Equation (38) would imply that the third regime, where JPe<sub>F</sub> ≃ 0.006Re<sub>F</sub>, may disappear altogether for

<footnote n="8">⁸ It is crucial to remember that this has nothing to do with marginal stability to linear perturbations, since the system continues to be turbulent well-into the region of linear stability.</footnote>
large enough Re<sub>F</sub>. However, given our current lack of explanation for the scalings observed in the intermediate regime, and the uncertainties in the determination of (Ri<sub>F</sub>Pe<sub>F</sub>)<sub>c</sub> (see Garaud et al. 2015a for detail), we caution the reader against using the formulas given in Equation (38) too far outside of the range of the available data. Within that range, however, (38) can be used to predict *J* (and hence the mean shear) resulting from a given forcing, and a given set of fluid parameters (viscosity, thermal diffusivity and stratification).

### 6.3. Typical Scale and Shape of the Turbulent Eddies

We now study in more detail the scale and shape of the turbulent eddies and how they vary with input parameters. In Section 4, we defined the vertical scale $\hat{l}_v$ (or equivalently, $\hat{l}_v$ since the two are the same) of the turbulent eddies, and showed their variation with Ri<sub>F</sub>Pe<sub>F</sub> in the LPN simulations in Figure 5. By analogy, we can also define the typical horizontal scales $\hat{l}_{h,x}$ and $\hat{l}_{h,y}$

$$ 
\begin{aligned}
\hat{l}_{h,x} &= \frac{1}{t_2 - t_1} \int_{t_1}^{t_2} \left[ \frac{\sum_{\hat{k}_x \neq 0} \sum_{\hat{k}_y} \sum_{\hat{k}_z} \hat{E}_u(\hat{k}, t) \hat{k}_x^{-1}}{\sum_{\hat{k}_x \neq 0} \sum_{\hat{k}_y} \sum_{\hat{k}_z} \hat{E}_u(\hat{k}, t)} \right] dt, \\
\hat{l}_{h,y} &= \frac{1}{t_2 - t_1} \int_{t_1}^{t_2} \left[ \frac{\sum_{\hat{k}_x} \sum_{\hat{k}_y \neq 0} \sum_{\hat{k}_z} \hat{E}_v(\hat{k}, t) \hat{k}_y^{-1}}{\sum_{\hat{k}_x} \sum_{\hat{k}_y \neq 0} \sum_{\hat{k}_z} \hat{E}_v(\hat{k}, t)} \right] dt, \quad (39)
\end{aligned}
$$

where $\hat{E}_u(\hat{k}, t)$ and $\hat{E}_v(\hat{k}, t)$ are the kinetic energies associated with fluid motion in the x- and y-directions respectively, with wavenumber $\hat{k}$. From their definitions, we see that $l_{h,x}$ is really the x-lengthscale of fluid motions in the x-direction, while $l_{h,y}$ is the y-lengthscale of fluid motion in the y-direction.

Figure 11 shows the variation of $\hat{l}_{h,x}$, $\hat{l}_{h,y}$ and $\hat{l}_v$ with Ri<sub>F</sub>Pe<sub>F</sub> for Re<sub>F</sub> = 100. As discussed in Section 3.2, $\hat{l}_v$, as well as $\hat{l}_{h,x}$ and $\hat{l}_{h,y}$, become independent of Ri<sub>F</sub>Pe<sub>F</sub> for sufficiently low Ri<sub>F</sub>Pe<sub>F</sub>. The fact that $\hat{l}_v \to 0.8$ at low Ri<sub>F</sub>Pe<sub>F</sub> suggests that the vertical eddy scale is indeed very similar to the imposed forcing lengthscale in this limit, which is not surprising since this is the only available non-diffusive scale in the system when buoyancy is negligible. In the opposite limit, we clearly see that the vertical eddy scale gradually decreases as stratification becomes more important. This trend agrees with the common notion that vertical overturning becomes more and more difficult, so only eddies whose vertical scale is small enough for diffusion to mitigate the effect of stratification are allowed. It is interesting, and somewhat surprising, to see that both $\hat{l}_{h,x}$ and $\hat{l}_{h,y}$ also decrease with increasing Ri<sub>F</sub>Pe<sub>F</sub>, so the anisotropy of the eddies does not vary as strongly with Ri<sub>F</sub>Pe<sub>F</sub> as one may naively expect: across most of the range, $\hat{l}_{h,y} \simeq \hat{l}_v$, while $\hat{l}_{h,x}/\hat{l}_v \sim 3.0-3.4$. In other words, the eddies take the form of “cigares” of more-or-less circular cross sections in the (y, z) plane, and three times longer in the x-direction. Whether this result still holds at much larger values of Ri<sub>F</sub>Pe<sub>F</sub> remains to be determined.

Since the vertical scale of the eddies is sometimes used (as in Zahn 1992, for instance) to make predictions on the effective mixing rate induced by the turbulence, one would ideally like to create a quantitative model for the variation of $\hat{l}_v$ with Ri<sub>F</sub>Pe<sub>F</sub> and Re<sub>F</sub>. However, we first note that the estimated vertical

<figure><img src="image_14.png" /><figcaption>Figure 11. Vertical lengthscales $\hat{l}_v$ and $\hat{l}_v'$, as well as the horizontal lengthscales $\hat{l}_{h,x}$ and $\hat{l}_{h,y}$ of the eddies as a function of Ri<sub>F</sub>Pe<sub>F</sub>, for Re<sub>F</sub> = 100, in the LPN equations (see the text for details).</figcaption></figure>

eddy scale depends somewhat on the definition of $\hat{l}_v$ adopted. Indeed, we could equivalently have chosen to use the definition

$$ \hat{l}'_v = \frac{1}{t_2 - t_1} \int_{t_1}^{t_2} \left[ \frac{\sum_{\hat{k}_x} \sum_{\hat{k}_y} \sum_{\hat{k}_z} \hat{E}_w(\hat{k}, t) \hat{k}_z}{\sum_{\hat{k}_x} \sum_{\hat{k}_y} \sum_{\hat{k}_z} \hat{E}_w(\hat{k}, t)} \right]^{-1} dt, \quad (40) $$

which is just as plausible as the one advocated in Equation (27). The two definitions would in fact be equivalent (within a constant factor) if the kinetic energy spectrum were a strict power law. However, this is not the case in our simulations, and as a result $\hat{l}_v$ and $\hat{l}_v'$ are not proportional to one another. We find that, for large Ri<sub>F</sub>Pe<sub>F</sub>,

$$ 
\begin{aligned}
\hat{l}_v &\sim (\mathrm{Ri}_F \mathrm{Pe}_F)^{-1/6}, \\
\hat{l}'_v &\sim (\mathrm{Ri}_F \mathrm{Pe}_F)^{-2/9}.
\end{aligned}
\quad (41) $$

The fact that the two power laws are different suggests that there are significant changes in the shape of the power spectrum as Ri<sub>F</sub>Pe<sub>F</sub> increases, which in turn implies that it will be difficult to create a simple a priori model for the variation of the eddy scale with stratification. At best, one can look at the empirical data and infer that $\hat{l}_v$ varies approximately as

$$ \hat{l}_v \sim (\mathrm{Ri}_F \mathrm{Pe}_F)^{-\alpha_l} \quad \text{for } \mathrm{Ri}_F \mathrm{Pe}_F \gg 1, \quad (42) $$

for large Ri<sub>F</sub>Pe<sub>F</sub>, where the constant of proportionality is of order one, and where $\alpha_l \simeq 0.19 \pm 0.03$ (depending on the definition adopted). Again, this is a tentative measurement, which would benefit from being confirmed with simulations at higher Re<sub>F</sub> and higher Ri<sub>F</sub>Pe<sub>F</sub>, so we caution the reader against using it too far outside of the range of parameters for which the formula has been established.

### 6.4. Energy Budget and Heat Transport

In order to study the kinetic energy budget, we start from the original set of Equations (10), and express them in the non-dimensionalization based on the forcing. We then dot the momentum equation with $\hat{u}$ and multiplying the thermal equation
with $\hat{T}$, and take a spatial average over the domain to find that

$$ \frac{1}{2} \frac{\partial}{\partial t} \langle |\hat{u}|^2 \rangle = \mathrm{Ri}_F \langle \hat{w} \hat{T} \rangle - \frac{1}{\mathrm{Re}_F} \langle |\nabla \hat{u}|^2 \rangle + \langle \hat{u} \cdot \hat{F} \rangle, \quad (43) $$

$$ \frac{1}{2} \frac{\partial}{\partial t} \langle \hat{T}^2 \rangle = - \langle \hat{w} \hat{T} \rangle - \frac{1}{\mathrm{Pe}_F} \langle |\nabla \hat{T}|^2 \rangle. \quad (44) $$

Note that the pressure term, $\langle \hat{u} \cdot \nabla \hat{p} \rangle = \langle \nabla \cdot (\hat{u} \hat{p}) \rangle = 0$ since the domain is periodic. Using similar steps, one can easily show that the same is true of the nonlinear terms $\langle \hat{u} \cdot (\hat{u} \cdot \nabla \hat{u}) \rangle$ and $\langle \hat{T} \hat{u} \cdot \nabla \hat{T} \rangle$. If, in addition, the system achieves a statistically stationary state, the thermal energy balance implies:

$$ \langle \hat{w} \hat{T} \rangle = - \frac{1}{\mathrm{Pe}_F} \langle |\nabla \hat{T}|^2 \rangle, \quad (45) $$

which shows that the turbulent heat flux $\langle \hat{w} \hat{T} \rangle$ caused by shear instabilities is negative (i.e., downward) in a stably stratified region. Substituting $\langle \hat{w} \hat{T} \rangle$ into the kinetic energy equation, we then have:

$$ \langle \hat{u} \cdot \hat{F} \rangle = \frac{\mathrm{Ri}_F}{\mathrm{Pe}_F} \langle |\nabla \hat{T}|^2 \rangle + \frac{1}{\mathrm{Re}_F} \langle |\nabla \hat{u}|^2 \rangle, \quad (46) $$

which shows how the power input into system by the force $\mathbf{F}$ is first converted into velocity and temperature fluctuations, which are then both dissipated microscopically. Note that we obtain exactly the same energy balance using the LPN equations, the only difference being that these equations assume that (45) is true at all times instead of in a quasi-stationary, domain-averaged sense (see Equation (17)). In that case, Equation (46) can also be rewritten as

$$ \langle \hat{u} \cdot \hat{F} \rangle = -\mathrm{Ri}_F \mathrm{Pe}_F \langle \hat{w} \nabla^{-2} \hat{w} \rangle + \frac{1}{\mathrm{Re}_F} \langle |\nabla \hat{u}|^2 \rangle. \quad (47) $$

It is common in the geophysical literature to measure the efficiency of stratified turbulence in mixing buoyancy through the ratio (Peltier & Caulfield 2003)

$$ \eta = - \frac{\mathrm{Ri}_F \langle \hat{w} \hat{T} \rangle}{\langle \hat{u} \cdot \hat{F} \rangle} = - \frac{\mathrm{Ri}_F \mathrm{Pe}_F \langle \hat{w} \nabla^{-2} \hat{w} \rangle}{\langle \hat{u} \cdot \hat{F} \rangle} \quad \text{in the LPN limit.} \quad (48) $$

This quantity measures the fraction of the input power that is effectively used to transport buoyancy irreversibly while $1 - \eta$ measures the fraction of the input power that is dissipated viscously (either through the dissipation of the turbulent fluctuations, or through the dissipation of the mean flow). The ratio $\eta$ is shown in Figure 12 as a function of $\mathrm{Ri}_F \mathrm{Pe}_F$ for $\mathrm{Re}_F = 100$, and for different values of $\mathrm{Pe}_F$ as well as for the LPN equations. For very low $\mathrm{Ri}_F \mathrm{Pe}_F$, $\eta$ increases linearly with $\mathrm{Ri}_F \mathrm{Pe}_F$ as $\eta \simeq 40$ $\mathrm{Ri}_F \mathrm{Pe}_F$. This is consistent with the notion that temperature is a passive scalar in that limit, and can be explained mathematically by noting that both $\langle \hat{u} \cdot \hat{F} \rangle$ and $\langle \hat{w} \nabla^{-2} \hat{w} \rangle$ (as with all other quantities in the system) are independent of $\mathrm{Ri}_F \mathrm{Pe}_F$ in effectively unstratified flows. For larger stratifications, we see that in the LPN equations $\eta$ has a first local maximum between $\mathrm{Ri}_F \mathrm{Pe}_F = 0.01$ and 0.1, then a dip around $\mathrm{Ri}_F \mathrm{Pe}_F = 1$, then another local maximum later on followed by another dip. While the dip at $\mathrm{Ri}_F \mathrm{Pe}_F = 1$ remains to be explained, the second dip corresponds to the transition to

<figure><img src="image_15.png" /><figcaption>Figure 12. Buoyancy transport efficiency factor $\eta$ at $\mathrm{Re}_F = 100$, in the LPN equations and in the standard equations at various $\mathrm{Pe}_F$. The errorbars are the rms temporal fluctuations around the mean, and are very large for the $\mathrm{Pe}_F = 10$ case at large $\mathrm{Ri}_F \mathrm{Pe}_F$ because of the bursts described in Section 3.3.</figcaption></figure>

a skewed state, and can be understood by noting that the partitioning of the flow between a turbulent zone and a laminar one effectively increases the viscous dissipation (in particular, that associated with the mean flow) and reduces the turbulent heat transport.

Generally speaking, however, we see that $\eta$ oscillates around about 0.1 (for $\mathrm{Ri}_F \mathrm{Pe}_F$ greater than about 0.01), a value which is notably smaller than the typical transport efficiency discussed in the high-Prandtl-number geophysical literature where $\eta$ is typically closer to 0.2 (see the review by Peltier & Caulfield 2003, for instance). The quantity $\eta$ being smaller at low Péclet number presumably stems from the fact that thermal diffusion plays a significant role in dissipating buoyancy fluctuations before they can cause mixing. Viewed in this light, it is in fact somewhat surprising that $\eta$ is not actually much smaller than 0.1 in the LPN equations. We also see in Figure 12 that the dependence of $\eta$ on $\mathrm{Ri}_F \mathrm{Pe}_F$ is roughly the same for the standard equations at low $\mathrm{Pe}_F$ and for the LPN equations, as expected. For higher $\mathrm{Pe}_F$, however, the results are very different and $\eta$ can be much larger than 0.1 (it can in fact be much larger than 0.2 in the limit of very strongly stratified shear flows).

These results have interesting consequences: they imply that the dimensional heat flux carried by shear-induced turbulence in the low Péclet number limit, and in the case where $\mathrm{Ri}_F \mathrm{Pe}_F \ge 0.01$, can be predicted to be

$$ \begin{align} \mathcal{F} &= \rho_0 c_p \left( -\kappa_T \frac{d\bar{T}}{dr} + \langle wT \rangle \right) \nonumber \\ &= \rho_0 c_p \left[ -\kappa_T \frac{d\bar{T}}{dr} + U_F k^{-1} \left( \frac{d\bar{T}}{dr} - \frac{dT_{\text{ad}}}{dr} \right) \langle \hat{w}\hat{T} \rangle \right] \nonumber \\ &= \rho_0 c_p \left[ -\kappa_T \frac{d\bar{T}}{dr} - \frac{\eta \langle \hat{u} \cdot \hat{F} \rangle}{\mathrm{Ri}_F} U_F k^{-1} \left( \frac{d\bar{T}}{dr} - \frac{dT_{\text{ad}}}{dr} \right) \right] \nonumber \\ &= - \rho_0 c_p \left[ \kappa_T \frac{d\bar{T}}{dr} + \frac{\eta \mathcal{P}}{\alpha \rho_0 g} \right] \tag{49} \end{align} $$

where $\rho_0$ is the local background density, $c_p$ is the specific heat at constant pressure, $\mathcal{P} = \langle \hat{u} \cdot \hat{F} \rangle$ is the dimensional power input into the system, which is always positive at steady-state,
$\alpha$ is the thermal expansion coefficient defined in Section 2, $g$ is the local gravity and $\eta$ is roughly equal to 0.1. The remainder of the power injected into the star, that is, about 90% of it, is dissipated viscously and ultimately heats the fluid locally (an effect which is however neglected in the Boussinesq equations).

## 6.5. Momentum Transport

The simulations discussed in Sections 3 and 6 have revealed a number of non-trivial results, in particular when it comes to the shape of the mean flow. Why is the sinusoidal shape preferred for low degrees of stratification, while more strongly stratified shear flows tend to have nearly piecewise linear profiles? And more curiously, why is there a sudden loss of symmetry in the flow for even more strongly stratified systems?

In order to make progress toward answering these questions, recall that the evolution of the mean flow as a function of time is given by the horizontal average of the horizontal component of the momentum equation, namely

$$ \frac{\partial \hat{u}}{\partial t} = - \frac{d}{dz} \langle \hat{u} \hat{w} \rangle_h + \frac{1}{\text{Re}_F} \frac{d^2 \hat{u}}{dz^2} + \sin(z). \quad (50) $$

When the system reaches a statistically stationary state, this equation implies a balance between the divergence of the momentum flux, the forcing, and the mean shear dissipation by viscosity. If the Reynolds stress $\langle \hat{u} \hat{w} \rangle_h$ were a known function of the local shear $d\hat{u}/dz$ (and of the input parameters $\text{Ri}_F \text{Pe}_F$ and $\text{Re}_F$) then (50) could be used in closed form to infer $\hat{u}(z)$. Ideally, one would like to use simulations such as the ones presented here to measure the dependence of $\langle \hat{u} \hat{w} \rangle_h$ on the mean shear and the input parameters, and then see if this could illuminate the flow shape problem. However, this type of constant forcing simulation does no lend itself well to this exercise. Indeed, integrating (50) in the statistically stationary state yields (within an additive constant)

$$ \langle \hat{u} \hat{w} \rangle_h = \frac{1}{\text{Re}_F} \frac{d\hat{u}}{dz} - \cos(z), \quad (51) $$

which implies that the Reynolds stress

$$ \langle \hat{u} \hat{w} \rangle_h \approx -\cos(z) + O(\text{Re}_F^{-1}) \quad (52) $$

becomes nearly independent both $\text{Ri}_F \text{Pe}_F$ and $\text{Re}_F$ for large enough $\text{Re}_F$. This contrived balance implies that momentum transport cannot easily be studied using this approach. We therefore defer the problem of understanding the momentum balance and mean flow shape to a subsequent paper, in which we shall revisit our numerical results in the light of a turbulence closure model (L. Kulenthirarajah & P. Garaud 2016, in preparation).

## 6.6. Compositional Mixing

We now finally turn to the question of mixing of a passive tracer, and whether the latter is adequately captured by the models of Endal & Sofia (1978) or Zahn (1992). For this purpose, we have measured, in a limited number of simulations using the standard equations, the vertical flux of a passive

tracer by simultaneously solving for the equation

$$ \frac{\partial \hat{C}}{\partial t} + \hat{u} \cdot \nabla \hat{C} + \hat{w} = \frac{1}{\text{Pe}_C} \nabla^2 \hat{C}, \quad (53) $$

where $C$ represents for instance the perturbation of the concentration of a particular chemical species away from a linearly stratified background. This equation is already non-dimensionalized using $[C] = k^{-1} \frac{dC}{dr}$ where $\frac{dC}{dr}$ is the dimensional constant background compositional gradient assumed to exist. The momentum equation remains unchanged since $C$ is assumed to be a passive tracer. The compositional Péclet number $\text{Pe}_C$ is defined as

$$ \text{Pe}_C = \frac{\kappa_T}{\kappa_C} \text{Pe}_F. \quad (54) $$

In all cases we have used $\text{Pe}_C = 100$, so the dependence of the results on $\text{Pe}_C$ have not been tested yet. Note that none of the mixing models currently discussed in the literature have any dependence on $\text{Pe}_C$.

Using these simulations, we extract the rate of compositional mixing by measuring the flux $\langle \hat{w} \hat{C} \rangle$. An effective turbulent diffusion coefficient can then be computed via

$$ \langle wC \rangle = -D \frac{d\bar{C}}{dr} \leftrightarrow \hat{D} \equiv \frac{kD}{U_F} = -\langle \hat{w} \hat{C} \rangle. \quad (55) $$

Our results are shown in Figures 13(a) and (b). Figure 13(a) shows $\hat{D}$ as a function of the input parameter $\text{Ri}_F \text{Pe}_F$ for simulations with $\text{Re}_F = 100$ and two different low Péclet numbers (note that we do not have compositional data for the LPN runs). As expected, the two datasets are reasonably consistent with one another. In the nearly unstratified limit ($\text{Ri}_F \text{Pe}_F \le 0.001$), we see that $\hat{D}$ tends to a constant that is roughly equal to 10. For larger values of $\text{Ri}_F \text{Pe}_F$, $\hat{D}$ appears to follow a power law with $\hat{D} \approx 0.3 (\text{Ri}_F \text{Pe}_F)^{-1/2}$. This scaling, as we shall see below, is consistent with the model of Zahn (1992), and appears to be valid from $\text{Ri}_F \text{Pe}_F \approx 0.01$ to $\text{Ri}_F \text{Pe}_F \approx 100$. Unfortunately, we do not have any data for larger $\text{Ri}_F \text{Pe}_F$, which would correspond to the very strongly stratified limit.

It is interesting to note, however, that the variation of $\hat{D}$ with $\text{Ri}_F \text{Pe}_F$ cannot be explained either by assuming that $\hat{D} \propto \hat{u}_{\text{rms}} \hat{l}_v$ or $\hat{D} \propto \hat{w}_{\text{rms}} \hat{l}_v$ (i.e., the product of a typical vertical velocity with a typical turbulent lengthscale), nor by assuming that the velocity is $\hat{S} \hat{l}_v$ so $\hat{D} \propto \hat{S} \hat{l}_v^2$ (where $\hat{S}$ is the non-dimensional mean shear at mid-layer, see Equation (36)). By contrast with $\hat{D}$, $\hat{w}_{\text{rms}} \hat{l}_v$ decreases much more weakly with $\text{Ri}_F \text{Pe}_F$, while $\hat{u}_{\text{rms}} \hat{l}_v$ and $\hat{S} \hat{l}_v^2$ actually increase slightly with $\text{Ri}_F \text{Pe}_F$. This shows that using these standard estimates for $\hat{D}$ can be very misleading.

Figure 13(b) shows the same data for $\hat{D}$, but this time plotted against $J \text{Pe}_F$ where $J = N^2/S^2 = \text{Ri}_F/\hat{S}^2$. We find in this case that $\hat{D} \approx 10$ for $J \text{Pe}_F \le 0.002$, while $\hat{D} \approx 0.02 (J \text{Pe}_F)^{-1}$ for $J \text{Pe}_F \ge 0.002$. The scaling for $\hat{D}$ in that limit is therefore consistent with the model of Zahn (1992), and suggests that, dimensionally speaking,

$$ \begin{aligned} D &\approx 10k^{-1}U_F \quad \text{for } J \text{Pe}_F < 0.002, \\ D &\approx 0.02 \frac{\kappa_T}{J} \quad \text{otherwise.} \end{aligned} \quad (56) $$
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Left: turbulent diffusion coefficient $\hat{D}$ as a function of $Ri_F Pe_F$, for $Pe_F = 0.1$ and $Pe_F = 1$. In low Péclet number flows, $\hat{D}$ appears to be well-approximated by $0.3(Ri_F Pe_F)^{-1/2}$ for $Ri_F Pe_F > 0.001$. Also shown are standard estimates of $\hat{D}$ as $\hat{w}_{rms} \hat{l}_v$ and $\hat{S} \hat{l}_v^2$, neither of which are good. Right: validation of Zahn’s model, which predicts that $\hat{D} \sim (JPe_F)^{-1}$. We find that this is the case for $JPe_F \ge 0.002$.</figcaption></figure>

The reason why Zahn’s model appears to explain our data is, however, puzzling. As discussed in Section 1, Zahn (1992) assumes that $D \propto ul$, where $l$ is the largest possible eddy size for which the shear is still nonlinearly unstable according to the criterion $JPe_l \le (JPe)_c$. As seen in Figure 13(a), however, that lengthscale cannot be $l_v$ since none of the estimates involving $l_v$ appear to work. This leaves us in the rather uncomfortable position of either trying to explain why the lengthscale that dominates the transport of passive scalars should be so very different from the vertical lengthscale of the energy-bearing eddies —something that goes against what is commonly assumed, or, to accept that the good match between our data and Zahn’s model is somewhat of a coincidence. Neither of these options are particularly satisfactory but both are equally plausible (or implausible). Future simulations using a different model setup will be needed to resolve this frustrating conundrum.

# 7. SUMMARY AND DISCUSSION

## 7.1. Summary of Our Findings

In this paper, we have studied the dynamics of shear instabilities in stably stratified stellar regions, in the limit where the ratio of the thermal diffusion timescale to the turnover timescale of turbulent eddies is short (the so-called “low Péclet number” limit). We have shown that this limit could be relevant in the envelopes of very massive stars, where the thermal diffusivity is in excess of $10^{14} \text{ cm}^2 \text{ s}^{-1}$, but most likely does not apply for lower-mass stars, or deep within the interiors of massive stars (see Section 5). Low Péclet number shear layers can be formally studied with the LPN equations derived by Lignières (1999), which can be integrated numerically faster than the standard equations (at least in spectral codes), see Section 2.

In this first study, we have chosen to ignore the possibility of compositional stratification and horizontal shear, and focussed only on the case of thermal stratification with vertical shear. We also specifically investigated the dynamics of shear flows that arise from the application of a constant-amplitude, spatially periodic body force. In this setup, the mean shear is not prescribed, but instead is one of the quantities that we measure.

From this modeling choice, three numbers naturally emerge: $Re_F$, $Pe_F$ and $Ri_F$ defined in Equation (20), which we have shown to be good approximations to the actual turbulent Reynolds number, turbulent Péclet number and turbulent Richardson numbers of the resulting statistically stationary shear flow respectively. As such, they ideally characterize the dynamics of the shear. We first found that while the LPN equations are only formally valid in the limit $Pe_F \to 0$, they are already a good-to-excellent approximation of the full equations whenever $Pe_F \le 1$. In this limit, the only relevant non-dimensional parameters are the Richardson-Péclet product $Ri_F Pe_F$ and the Reynolds number $Re_F$. As such, instability is possible for large $Ri_F$ as long as $Pe_F$ is small enough.

As shown by Garaud et al. (2015a), for large enough Reynolds number this sinusoidal shear flow is linearly unstable provided $Ri_F Pe_F \le 0.25Re_F$, and is energy-stable (i.e., stable to any possible perturbations of any amplitude) for $Ri_F Pe_F$ greater than $\xi Re_F^3$, where $\xi$ is of order one. In practice, finite amplitude instabilities are rather difficult to trigger, and have only been found up to $Ri_F Pe_F$ of the order of a few times $Re_F$, which is not much larger than the linear instability threshold itself. Whether the turbulent solution can be continued for larger $Ri_F Pe_F$ remains to be determined. If these results can be directly applied to stars, they strongly suggest that solutions in the high $Ri_F Pe_F$ linearly stable/nonlinearly unstable region of parameter space can only be found if the shear slowly and progressively decreases over time from a point when it was linearly unstable, or if the stratification progressively increases likewise. In the opposite case, the shear flow would start as a laminar flow and remain so until the conditions are such that it becomes linearly unstable. In practice, this implies that whether a shear flow in the linearly stable/nonlinearly unstable region of parameter space is actually turbulent or not depends more on its history than on the present conditions, a classical case of hysteresis.

Using DNS, we have found that the dynamics of low Péclet number shear flows can be divided into three categories. In the limit where $Ri_F Pe_F \le 0.01$, the stratification has a negligible effect on the flow dynamics, and temperature behaves as a passive tracer. The amplitude of the resulting mean shear is of order $S \sim (kF_0/\rho_0)^{1/2}$ where $F_0$ is the amplitude of the forcing, $k^{-1}$ is the typical lengthscale of the forcing, and $\rho_0$ is the local
density. Most of the kinetic energy, however, is in the turbulent fluctuations rather than in the mean flow. Heat is transported downward, as expected from turbulence in stably stratified fluids. In stellar interiors, this means that heat is transported down the entropy gradient, but up the temperature gradient. The turbulent heat flux can be derived from (49) with $\eta \simeq 40\text{Ri}_F\text{Pe}_F$, and is dimensionally proportional to the local stratification as measured by $N^2$. The effective diffusivity of a passive tracer is given by $D \simeq 10(F_0/\rho_0 k^3)^{1/2}$.

In the opposite limit, where $1 \ll \text{Ri}_F\text{Pe}_F < (\text{Ri}_F\text{Pe}_F)_c$ (the exact lower and upper thresholds depend on the Reynolds number, as described in Equation (38)), the stratification entirely governs the dynamics of the shear flow. This region of parameter space is usually linearly stable. Most of the total kinetic energy now lies in the mean flow, and the mean shear adjusts itself to be in a state of marginal nonlinear stability, which is captured by a criterion similar to Zahn's criterion (Zahn 1974, see Equation (4)), namely $J\text{Pr} \simeq 0.006$ where $\text{Pr}$ is the Prandtl number and $J = N^2/S^2$ is the local gradient Richardson number. This implies that the mean shear is simply proportional to $N$ as $S \simeq (\text{Pr}/0.006)^{1/2}N$ in this limit. Meanwhile the turbulent heat flux is given by (49) with $\eta \simeq 0.1$, which physically means that a universal constant fraction of the total input power provided by the forcing, namely 10%, goes into buoyancy transport while the rest is dissipated viscously.

Finally, in the intermediate regime, the variation of the mean shear with stratification is more difficult to interpret (see Section 6.2), with $J\text{Pe}_F \simeq 0.1(\text{Ri}_F\text{Pe}_F)^{1/2}$. The turbulent heat flux is however still reasonably well approximated by (49) with $\eta \simeq 0.1$. In both the intermediate and in the strongly stratified regime, we find that the effective diffusivity of a passive tracer satisfies Zahn's mixing model (Zahn 1992) with $D \simeq 0.02\kappa_T/J$, even though the reason why this is the case remains elusive (see Section 6.6 for detail). In particular, we have found that contrary to what is commonly assumed in astrophysics, $D$ is not well-approximated by the product of the typical vertical velocity and typical lengthscale of the turbulent eddies.

Finally, we have found that for large enough $\text{Re}_F$ and for $\text{Ri}_F\text{Pe}_F$ approaching the nonlinear stability threshold $(\text{Ri}_F\text{Pe}_F)_c$, the system spontaneously transitions into new state which does not have the same symmetries as the imposed forcing. Instead, the mean shear becomes skewed and is partitioned between wider turbulent regions of moderate shear, and thinner laminar regions of very strong shear. The reason for this transition remains to be determined, but could be attributed to the non-monotonicity of the stress-strain relationship, in a manner that is similar to the spontaneous formation of layers in stratified systems that have non-monotonous relationships between the buoyancy flux and the buoyancy gradient (Balmforth et al. 1998). In order to study this effect in more detail, and confirm that it is a generic result rather than a peculiarity of our sinusoidal forcing assumption, we shall need to use a different model setup which is not as strongly constrained in terms of the momentum transport balance as the one we are currently using (see Section 6.5 and below for more detail). This will be the subject of a forthcoming study.

## 7.2. Discussion and Outlook

All of the aforementioned results have been obtained for simulations with moderate turbulent Reynolds numbers $\text{Re}_F$ up to 100, but need to be confirmed with future runs at higher $\text{Re}_F$ when these become feasible in a more reasonable amount of time

than what is currently possible. Specifically, we need to better determine (a) the upper limit for instability to finite amplitude perturbations $(\text{Ri}_F\text{Pe}_F)_c$ as a function of $\text{Re}_F$ and (b) the variation with $\text{Re}_F$ of the various dynamical regime thresholds discussed in Section 6.2 and the dependence of the mean shear, heat flux, and turbulent diffusivity on $\text{Ri}_F\text{Pe}_F$ within each regime.

We have, on the other hand, run a number of simulations with different domain sizes to determine the impact of the system geometry. Generally speaking, all of the simulations presented here were done in a domain that is sufficiently large to guarantee that halving it does not affect the results much, and in the few cases where we have doubled the domain width, length or height, no difference in the mean quantities larger than the error-bars quoted was noticed. The only exception was in the case of the skewed runs when the domain height was doubled. In that case, the system exhibits interesting time-dependent dynamics, whereby the mean shear oscillates between different skewed quasi-steady states. However, given our reservations about the contrived nature of this type of forcing, we defer a discussion of these results until the momentum transport through the system is better understood (see L. Kulenthirarajah & P. Garaud 2016, in preparation).

Indeed, while our choice of using a sinusoidal forcing was principally motivated by numerical convenience, this kind of model appears to suffer from two distinct issues. First, as discussed in Section 6.5, the imposed forcing strongly constrains the horizontal momentum transport equation, to the extent that it is difficult to study momentum transport with this setup. As a result, what controls the overall shape and amplitude of the mean flow still remains, to a great extent, unknown. More crucially, as noted by Garaud et al. (2015a), the dependence of the threshold for linear stability on $\text{Re}_F$ is fundamentally different in the case of a sinusoidal shear and in the case of a hyperbolic tangent shear layer (Lignières et al. 1999), the latter being much closer to the energy stability threshold than the former. This suggests that the shape of the shear itself plays a crucial role in the development of shear instabilities. Hence, how much of the results obtained here are specific to the sinusoidal shear case, and how much are generic to all types of shear layers, remains to be determined.

There are several alternatives to using a sinusoidal forcing. A commonly adopted solution is to use a uniform shear model (Rogallo 1981; Jacobitz et al. 1997; Brucker et al. 2007; Chung & Matheou 2012; Matheou & Chung 2012; Prat & Lignières 2013, 2014). In this case, the background shear is somehow imposed and maintained, and the turbulence merely adjusts itself to the given shearing rate. There are two difficulties associated with this approach, the first being the maintenance of the shear. The simplest way to create a background linear shear flow is to use, as in Prat & Lignières (2013), an adaptive force which is calculated at every timestep to drive the shear precisely back toward a linear profile should it start deviating away from it. A disadvantage of this model is that the energetics of the system are difficult to study, since there is no control over the force actually required to maintain the shear. A second possibility is to use the shearing sheet model, well-known in the context of astrophysical disks. The numerical implementation of this approach is, however, more difficult especially in a spectral code (see the discussions in Brucker et al. 2007).

The second difficulty with using a constant shear background flow is the fact that this system is well-known to be
always linearly stable (Knobloch 1984), whether unstratified or stratified, and requires finite amplitude perturbations to trigger and maintain turbulence. This, in itself, is not a technically difficult problem to overcome—one simply needs to find the unstable finite amplitude branch of solution for relatively low $Ri_F Pe_F$ and then follow it by continuation for higher and higher $Ri_F Pe_F$, as we have done already in this paper. However, it is also the case that even a small amount of curvature in the shear, together with the presence of an inflection point, is all that is needed to have instability to infinitesimal perturbations (for low enough $Ri_F Pe_F$). This raises an interesting question: all parameters being the same, is the nature of the turbulence significantly different when triggered by supercritical instabilities than when triggered by subcritical ones? In fact, there is no reason to believe that this would be the case. Even far from the onset of instability, the shape and growth timescale of the fastest-growing linearly unstable modes can still influence the dynamics of a fully turbulent flow, serving as its injection scale. Meanwhile, the same statement could perhaps apply to the finite amplitude subcritical modes, but these are likely to be very different in nature from the global linear modes. In other words, we suspect that the turbulent dynamics obtained in a simulation at constant shear are not always necessarily representative of the ones one may obtain from global simulations with spatially varying shear when linear instabilities are present. Clearly, both approaches (spatially varying body force or constant background shear) have their pros and cons, and it is our belief that much can be learned from comparing the outcomes of both types of simulations run at similar parameters. This will be the subject of future work.

This work is funded by NSF-AST1517927 and NSF-AST1412951. The authors thank Stephan Stellmach for granting us the use of his excellent code, and for many useful discussions. The authors also thank Neil Balmforth, Tobias Bischoff, Nic Brummell, Colm-cille Caulfield, Basile Gallet, and Paul Linden for their insight into this work. The authors also thank the anonymous referee for his or her outstanding suggestions, which greatly improve the clarity of the paper. The simulations presented in this paper were run on the Hyades cluster at UC Santa Cruz, purchased using an NSF-MRI grant. This paper is written in memory of Jean-Paul Zahn, whose visionary work continues to inspire us today.

## REFERENCES

* Balmforth, N. J., Smith, S. G. L., & Young, W. R. 1998, JFM, 355, 329

* Barnes, G., Charbonneau, P., & MacGregor, K. B. 1999, ApJ, 511, 466

* Batchelor, G. K. 1953, The Theory of Homogeneous Turbulence (Cambridge: Cambridge Univ. Press)

* Brucker, K. A., Isaza, J. C., Vaithianathan, T., & Collins, L. R. 2007, JCoPh, 225, 20

* Charbonnel, C., & Zahn, J.-P. 2007, A&A, 467, L15

* Chung, D., & Matheou, G. 2012, JFM, 696, 434

* Deheuvels, S., Doğan, G., Goupil, M. J., et al. 2014, A&A, 564, A27

* Dudis, J. J. 1974, JFM, 64, 65

* Eggenberger, P., Montalbán, J., & Miglio, A. 2012, A&A, 544, L4

* Elliott, J. R., & Gough, D. O. 1999, ApJ, 516, 475

* Endal, A. S., & Sofia, S. 1978, ApJ, 220, 279

* Fliegner, J., Langer, N., & Venn, K. A. 1996, A&A, 308, L13

* Gage, K. S., & Miller, J. R. 1974, PhFl, 17, 883

* Garaud, P. 2001, MNRAS, 324, 68

* Garaud, P., Gallet, B., & Bischoff, T. 2015a, PhFl, 27, 084104

* Garaud, P., Medrano, M., Brown, J. M., Mankovich, C., & Moore, K. 2015b, ApJ, 808, 89

* Gough, D. 2007, in The Solar Tachocline, ed. D. W. Hughes, R. Rosner, & N. O. Weiss (Cambridge: Cambridge Univ. Press), 3

* Gratton, R. G., Sneden, C., Carretta, E., & Bragaglia, A. 2000, A&A, 354, 169

* Howard, L. N. 1961, JFM, 10, 509

* Hughes, D. W., Rosner, R., & Weiss, N. O. 2007, The Solar Tachocline (Cambridge: Cambridge Univ. Press)

* Jacobitz, F. G., Sarkar, S., & van Atta, C. W. 1997, JFM, 342, 231

* Jones, C. A. 1977, GApFD, 8, 165

* Knobloch, E. 1984, GApFD, 29, 105

* Lignières, F. 1999, A&A, 348, 933

* Lignières, F., Califano, F., & Mangeney, A. 1999, A&A, 349, 1027

* Matheou, G., & Chung, D. 2012, PhFl, 24, 091106

* Meyer, C. R., & Linden, P. F. 2014, JFM, 753, 242

* Miles, J. W. 1961, JFM, 10, 496

* Mosser, B., Goupil, M. J., Belkacem, K., et al. 2012, A&A, 548, A10

* Peltier, W. R., & Caulfield, C. 2003, AnRFM, 35, 135

* Pinsonneault, M. 1997, ARA&A, 35, 557

* Prat, V., & Lignières, F. 2013, A&A, 551, L3

* Prat, V., & Lignières, F. 2014, A&A, 566, A110

* Richardson, L. F. 1920, RSPSA, 97, 354

* Rogallo, R. S. 1981, STIN, 81, 31508

* Schatzman, E. 1977, A&A, 56, 211

* Sestito, P., & Randich, S. 2005, A&A, 442, 615

* Skumanich, A. 1972, ApJ, 171, 565

* Spiegel, E. A., & Veronis, G. 1960, ApJ, 131, 442

* Stellmach, S., Traxler, A., Garaud, P., Brummell, N., & Radko, T. 2011, JFM, 677, 554

* Talon, S., & Zahn, J.-P. 1997, A&A, 317, 749

* Thual, O. 1992, JFM, 240, 229

* Townsend, A. A. 1958, JFM, 4, 361

* Traxler, A., Stellmach, S., Garaud, P., Radko, T., & Brummell, N. 2011, JFM, 677, 530

* Watson, M. 1981, GApFD, 16, 285

* Zahn, J.-P. 1974, in Proc. IAU Symp. 59, Stellar Instability and Evolution, ed. P. Ledoux, A. Noels, & A. W. Rodgers (Dordrecht: D. Reidel), 185

* Zahn, J.-P. 1992, A&A, 265, 115