The Combined Effects of Vertical and Horizontal Shear Instabilities in Stellar Radiative Zones

Shear instabilities can be the source of significant amounts of turbulent mixing in stellar radiative zones. Past attempts at modeling their effects (either theoretically or using numerical simulations) have focused on idealized geometries, where the shear is either purely vertical or purely horizontal. In stars, however, the shear can have arbitrary directions with respect to gravity. In this work, we use direct numerical simulations to investigate the nonlinear saturation of shear instabilities in a stably stratified fluid, where the shear is sinusoidal in the horizontal direction and either constant or sinusoidal in the vertical direction. We find that in the parameter regime studied here (nondiffusive, fully turbulent flow), the mean vertical shear does not play any role in controlling the dynamics of the resulting turbulence, unless its Richardson number is smaller than 1 (approximately). As most stellar radiative regions have a Richardson number much greater than 1, our result implies that the vertical shear can essentially be ignored in the computation of the vertical mixing coefficient associated with shear instabilities for the purpose of stellar evolution calculations, even when it is much larger than the horizontal shear (as in the solar tachocline, for instance).


INTRODUCTION
Quantifying vertical mixing by shear instabilities in stellar radiative zones is a long-standing question that dates back to the 1970s and the works of Zahn (1974) (see also Schatzman 1969;Spiegel & Zahn 1970).Indeed, shear is ubiquitous in stars.For instance, the gradual spin-down of the stellar surface by magnetized winds naturally generates some level of radial shear on a global scale via angular momentum transport (Spada et al. 2016).Similarly, large-scale meridional flows advecting the star's mean angular momentum poleward or equatorward naturally drive the emergence of a global-scale latitudinal shear (Zahn 1992;Spiegel & Zahn 1992).More localized shear layers can be maintained in the vicinity of a differentially-rotating convection zone, as evidenced by the helioseismic discovery of the solar tachocline, which has both strong radial and latitudinal shear (Christensen-Dalsgaard & Schou 1988;Thompson et al. 1996).Internal gravity waves can in some cases amplify pre-existing shear and create strong localized radial shear layers in a star (Kumar & Quataert 1997;Ringot 1998;Kumar et al. 1999;Charbonnel & Talon 2005).
Under the right conditions, the kinetic energy of the shear can drive and maintain turbulence via shear instabilities (Richardson 1920).It is therefore important to understand what these conditions are and to quantify the resulting vertical mixing.This problem turns out to be particularly complex, however, because shear instabilities can take many different forms depending on the orientation of the shear with respect to the density stratification, as well as the strength and even the shape of the shear and density profiles (cf.Drazin & Reid 2004;Caulfield 2021).Shear instabilities can also be affected by diffusion (viscous, thermal, or compositional), rotation and magnetic fields (Chandrasekhar 1961).For these reasons, one must take a very gradual approach in studying the problem, adding only one new physical effect at a time to fully grasp its impact on the properties of the developing turbulence.Furthermore, while other instabilities can often be studied quite successfully using linear stability theory (i.e.stability to infinitesimal perturbations), this is rarely the case with shear instabilities.As long as viscosity is small, they can almost always be destabilized by suitably chosen finite-amplitude perturbations even when the background flow is linearly stable (see, e.g.Grossmann 2000;Garaud et al. 2015;Avila et al. 2023), and thus, the outcome sensitively depends on the initial conditions applied.As a result, studying the conditions for instability and resulting turbulent mixing often requires a combination of heuristic arguments and heavy-duty computational tools.
For all of these reasons, progress in modeling mixing by shear instabilities in stars has been relatively slow, focusing for now on understanding simple geometries (vertical shear only, or horizontal shear only), ignoring the effects of magnetic fields (except in linear stability analyses) and mostly including the effect of rotation only insofar as it is the source of the shear in stars with differential rotation, but it is ignored thereafter.Within that limiting set of assumptions, however, starting with the work of Zahn (1974) and Zahn (1992), combined with the more recent theoretical work and numerical experiments performed over the past 10 years (see the reviews of Lignières (2020) and Garaud (2021), as well as Sections 1.1 and 1.2 below), we now have a growing understanding of the criterion for instability.We also have models for the turbulent mixing associated with (non-rotating, non-magnetic) vertical and horizontal shear instabilities when these processes are considered independently.
The main goal of this paper is to expand on this previous work to determine what happens to mixing by shear instabilities when vertical and horizontal shear are both present at the same time in a (non-rotating, non-magnetic) stably stratified region.With this in mind, we begin by briefly reviewing what is known about mixing by horizontal and vertical shear instabilities separately, which will also serve as a pedagogical introduction to the concepts needed for the interpretation of the red numerical simulations presented in this paper.
1.1.What is known about (non-rotating, non-magnetic) vertical shear instabilities Loosely defined, vertical shear instabilities are processes that extract energy from the local vertical shear to generate and maintain some amount of vertical mixing.To do so in a stellar radiative zone, the shear must be strong enough to overcome the stabilizing effects of both viscosity and stratification.Because the viscosity is very small in stars, one could naively assume that it is never relevant -and that is often, but not always the case (see below).Stratification, on the other hand, plays a dominant role in quenching vertical shear instabilities.Indeed, in order to tap into the mean shear kinetic energy reservoir, the turbulent eddies must move vertically but this incurs a corresponding energy transfer to the potential energy reservoir of the stratification when doing so adiabatically.The competition between the kinetic energy gain and potential energy cost is loosely captured by the so-called Richardson number J, defined here as the square of the ratio of the mean buoyancy frequency N to the mean vertical shearing rate S m , i.e.
Assuming viscous and thermal diffusion are both negligible, the rate at which eddies can gain kinetic energy from the mean flow is proportional to S 2 m , while the rate at which they must provide potential energy to the stratification (by irreversibly lifting dense fluid parcels and lowering buoyant fluid parcels) is proportional to N 2 .Vertical shear instabilities can therefore develop (linearly or nonlinearly) provided J is (roughly) smaller than one (Richardson 1920).A similar criterion that describes a necessary condition for the onset of linear instability can be derived more formally (cf. Howard 1961); in that case, the existence of an extremum in the shearing rate, or in other words, an inflection point in the flow profile, is generally also required (Fjørtoft 1950;Drazin & Reid 2004).
It is quite unusual, however, for stars to have extended regions that satisfy J ≤ 1 for long periods of time, except very close to the edge of a convective zone (where N → 0), or when an external mechanism is present to sustain an intense shear (such as accretion of material from a disk, see MacDonald 1983;Avila et al. 2023).As such, the 'dynamical' shear instabilities listed above do not appear to play a significant role in stellar evolution.
All of the aforementioned arguments have assumed adiabatic perturbations, which is not necessarily the case if the stratification is primarily due to temperature and thermal diffusion is fast compared with the local dynamics.In the early 1970s, progress in understanding shear instabilities in the Earth's atmosphere (Townsend 1958), where perturbations are not adiabatic due to radiative losses, led Spiegel & Zahn (1970) to the realization that rapid thermal diffusion would have a similar effect in stellar interiors.Zahn (1974Zahn ( , 1992) ) then proposed the first model of thermally diffusive vertical shear instabilities (now often referred to as 'secular' shear instabilities), where he argued that they mix momentum and chemical species vertically with a diffusivity and where P r = ν/κ T is the Prandtl number, ν is the viscosity and κ T is the thermal diffusivity.A series of direct numerical simulations (DNS hereafter) by Prat & Lignières (2013, 2014); Prat et al. (2016), as well as Garaud & Kulenthirarajah (2016) and Garaud et al. (2017) confirmed the validity of (2) as long as J ≫ 1 and thermal diffusion is important on the scale of the shear layer.
In practice, however, D mix is always fairly small in secular shear instabilities (at most, a couple of orders of magnitude larger than ν, see Garaud 2021).For this reason, they are not usually a significant source of mixing unless no other instabilities are present.

1.2.
What is known about (non-rotating, non-magnetic) horizontal shear instabilities While vertical shear instabilities are either rare (in the dynamical regime) or not very efficient at mixing (in the secular regime), Zahn (1992) realized that horizontal shear instabilities could be a much more reliable source of mixing.Horizontal shear instabilities extract energy from the horizontal shear without any potential energy cost as long as the motions remain horizontal.As such, as long as magnetic fields and rotation are ignored, a horizontal shear is almost always unstable to purely horizontal perturbations (Balmforth & Young 2002;Arobone & Sarkar 2012;Cope et al. 2020;Park et al. 2020), that drive a meandering of the mean flow, and/or the formation of large-scale 'pancake'like vortices.The low viscosity then implies that these meanders and vortices decouple in the vertical direction, thus leading to the emergence of strong vertical shear on small vertical scales, even in the absence of any mean vertical shear.Small-scale vertical shear instabilities, and associated vertical mixing, can thus be triggered by the presence of an unstable large-scale horizontal shear.This general picture has been validated in DNS of Cope et al. (2020) and Garaud (2020a) (see more on this below).
By contrast with the case of vertical shear instabilities where models generally agree with each other, and with the data from numerical simulations, quantifying mixing by horizontal shear instabilities is a topic of ongoing investigations where no consensus has yet been reached.Several theoretical models have been put forward, starting with Zahn (1992), followed by the more recent works of Cope et al. (2020), Garaud (2020a), Lignières (2020), Chini et al. (2022), Skoutnev (2023) and Shah et al. (2023).These models vary significantly in their detailed predictions because of differences in their respective founding assumptions and in the mathematical approach selected, and because their domains of validity are sometimes limited.It is not the purpose of this section to discuss each model in turn, as it would be difficult to do so without any bias (see Shah et al. 2023, for a rapid review of the subject).Instead, we summarize here what general agreement the models have, with the hope that these are concepts that can be relied upon with reasonable confidence.
Generally speaking, all models attempt to relate the vertical mixing coefficient D mix to properties of the horizontal turbulence.Zahn (1992), for instance, directly related D mix to the horizontal turbulent diffusivity ν h , while the other models relate D mix to the assumed properties of large-scale horizontal eddies (which have characteristic large-scale horizontal velocity U h and lengthscale L h ).It is also generally agreed that as long as the flow is fully turbulent, D mix ∝ l z w rms (where l z is a characteristic vertical eddy lengthscale and w rms is a characteristic vertical eddy velocity), and that both of these quantities must decrease with the increasing stratification.Finally, we know that there ought to be different regimes depending on whether the turbulent Péclet number based on w rms and l z , namely P e t = w rms l z /κ T , is larger than 1 (adiabatic regime) or smaller than 1 (diffusive regime), and that the transition between the two regimes is a function of the stratification and of the thermal diffusivity.It is also worth noting that viscosity is expected to become relevant once the stratification exceeds some threshold, although not all models take this into account.Because of this, the partitioning of parameter space into various regimes is relatively complex (see, for instance, Shah et al. 2023).
The main disagreements between the models emerge in the derivations of w rms and l z as functions of the properties of the horizontal turbulence.This also affects where the boundaries of the various regimes (adiabatic, diffusive, viscous) lie in parameter space.While the theoretical controversy is ongoing, there is some hope of gaining insight from DNS data.Cope et al. (2020) and Garaud (2020a) recently ran a series of DNS of horizontally sheared, stratified Kolmogorov flows at low Prandtl number.In these simulations, a constant-in-time body force F = F 0 sin(k h y)e x is applied to drive a vertically invariant, sinusoidal flow.The anticipated horizontal eddy scale is therefore an input parameter, L h = O(k −1 h ), while the anticipated horizontal velocity scale is also an input parameter, obtained from balancing the nonlinear terms and the forcing in the streamwise component of the momentum equation: U h = O( F 0 /ρ m k h ) where ρ m is the mean density of the fluid.Cope et al. (2020) explored cases where thermal diffusion is important at all scales, so the Péclet number based on the large scales, P e = U h L h /κ T , is smaller than 1.Meanwhile Garaud (2020a) explored cases where thermal diffusion is negligible on the large scales (P e ≫ 1).
In the thermally diffusive regime (P e < 1), Cope et al. (2020) found, rather robustly, that there exists an intermediate regime of stratified turbulence where This regime is valid as long as the vertical scale of the turbulent eddies is much smaller than the domain size, but yet large enough to not be influenced by viscosity.The scaling law (3) is consistent1 with the model predictions from Zahn (1992), Lignières (2020), Skoutnev (2023) and Shah et al. (2023) in their respective thermally-diffusive regimes, which is rather surprising as the scalings for w rms and l z individually differ in each model.However, the constraint P e < 1 is rarely satisfied in stars.In the case where P e ≫ 1, Garaud (2020a) found that (3) does not always apply and that two regimes exist depending on whether the turbulent Péclet number P e t defined above is larger or smaller than one.More specifically, D mix is independent of the thermal diffusivity when P e t ≫ 1, but recovers the diffusive results of Cope et al. (2020) when P e t ≪ 1.Unfortunately, as P e t is an emergent quantity (rather than an input parameter), a model is needed to predict when the transition from non-diffusive to thermally diffusive turbulence occurs as a function of the input parameters.Here, the disagreement between the models themselves is significant; for instance, Skoutnev (2023) finds that the solar tachocline would lie close to the regime transition while Shah et al. (2023) argue instead that it lies squarely in the non-diffusive regime.Finally, we note that running DNS in the regime P e ≫ 1 and P r ≪ 1 is numerically very challenging (as these require extremely large Reynolds numbers).As a result, numerical experiments have yet to reach sufficiently high values of P e to be able to test model predictions in that regime, and establish which one (if any) is correct.

Outline of the proposed work
While we will need to wait for future DNS to validate and/or invalidate existing models, we now return to the original question posed earlier, namely what happens when mean vertical and horizontal shear are combined, which is the more relevant scenario in stellar interiors.Indeed, the interplay between rotation, meridional circulations, spin-down by stellar winds, and shear instabilities naturally contributes to the maintenance of a large scale rotational shear in both radius and latitude, clarified in the seminal paper by Zahn (1992).The solar tachocline, for example, supports a large-scale horizontal differential rotation (with the equator rotating approximately 30 percent faster than the poles), as well as the well-known sharp vertical shear that is at the origin of its name (Spiegel & Zahn 1992).Red Giant Branch (RGB) stars are known to have substantial radial shear between the core and envelope (see Aerts et al. 2019, and references therein) which must necessarily be associated with horizontal shear by Zahn's argument.Conversedly, stars with outer convective zones can exhibit substantial surface latitudinal differential rotation (Reinhold et al. 2013) which, by analogy with the solar tachocline, suggests the likely existence of strong radial shear below.
In this study, we continue past investigations by Garaud & Kulenthirarajah (2016); Garaud et al. (2017), Cope et al. (2020) and Garaud (2020a), and run a series of DNS experiments which combine vertical and horizontal shear to address the effects of the former on the latter.Our goal is to establish a set of guiding principles to determine if and when one type of instability clearly dominates over the other, or when both interact in such a way that they cannot be modeled separately.We focus on the high Péclet number regime (P e ≫ 1) defined earlier, which is arguably the more relevant situation for the majority of stars (Garaud 2020a).Section 2 describes the model setup, in which two different vertical shear profiles (constant or sinusoidal) are tested.Section 3 presents our results, both qualitatively and then more quantitatively.Finally, Section 4 summarizes our results while discussing important caveats of our model assumptions.

MODEL SETUP
In what follows, we consider a stably-stratified (i.e., radiative) region of a star.Gravity defines the vertical direction (g = −ge z ), and the vertical extent of the domain is assumed to be much smaller than a density, pressure or temperature scaleheight.With this assumption, the local Brunt-Väisälä frequency N can be assumed to be constant.We use the Boussinesq approximation for gases (Spiegel & Veronis 1960), noting that fluid flows in the deep interiors of stars are much slower than the local sound speed.A mean horizontal flow is driven by application of a body force F = F e x , and/or the use of a shearing box (see below for more detail).We then study the development of shear instabilities on this mean flow, ignoring in this work the effects of rotation and magnetic fields.This assumption is made for simplicity and is briefly discussed in Section 4. Future work will investigate their effects separately.
The governing equations describing this model are where u = (u, v, w) is the velocity field, T is the potential temperature fluctuation away from the stably-stratified background state, and p is the pressure fluctuation away from hydrostatic equilibrium.In addition, ρ m is the mean density of the fluid, α is the coefficient of thermal expansion, ν is the viscosity, and κ T is the thermal diffusivity.All of these fluid properties are assumed to be constant in the small region considered.The body force F = F e x varies sinusoidally in the y direction, with a horizontal wavenumber k h = 2π/L y where L y is the domain width.
We then non-dimensionalize the equations as in Cope et al. (2020) and Garaud (2020a).We use k h as the inverse unit length.Anticipating a balance between the forcing and the inertial terms in the horizontal component of the momentum equation, we define where F 0 is the amplitude of the sinusoidal body force, and use U F as the unit velocity.The unit time is then defined as (k h U F ) −1 and the unit temperature as N 2 /αgk h .Note that since the lengthscale used in ( 7) is k −1 h , this non-dimensionalization more appropriately captures the scalings of anticipated horizontal shear instabilities than of vertical shear instabilities (which would likely care about the vertical lengthscale of the shear or the domain height L z ).
The non-dimensional governing equations corresponding to (4)-( 6) are then: In that follows, hatted quantities are non-dimensional, while those without hats are dimensional, except for the independent variables x, y, z, t and associated operators which are assumed to be nondimensional (hats in this case are removed to avoid crowding the notation).Several parameters are introduced in these equations: the domain-scale Reynolds number Re, which is the ratio of the viscous diffusion timescale to the turbulent advection timescale; the domain-scale Péclet number P e, which is the ratio of the thermal diffusion timescale to the turbulent advection timescale; and the buoyancy parameter B, which is the square of the ratio of the Brunt-Väisälä frequency N to the horizontal shearing rate k h U F of the anticipated mean flow: As we are interested in studying the combination of vertical and horizontal shear on the development of shear instabilities, we note that there are different possible ways of driving the shear.In this work, we always drive the horizontal shear using a body force as in Cope et al. (2020) and Garaud (2020a), so the results can be compared with their work in the limit where the vertical shear is zero.To apply vertical shear, however, there are several options.One of them is to use the standard shearing box formalism (as in, e.g., Prat et al. 2016) and drive a constant vertical shear.Another option is to use a body force that is sinusoidal in both horizontal and vertical directions.Both types of shear may be expected in stars.For instance, the constant shear case would be appropriate to model stars where the vertical shear varies slowly with depth.On the other hand, the sinusoidal shear case could be appropriate to model stars with complex vertical shear profiles, such as the ones that can be driven by the nonlinear interaction between gravity waves and differential rotation.
In what follows, we use and compare both options, which serves two purposes.First, it might help establish the role of the vertical shear lengthscale on the emergent turbulence.In the sinusoidal model proposed below, the vertical shear lengthscale is finite while it is infinite in the constant shear case.Second, and likely more importantly, a sinusoidal vertical shear has an inflection point, and can excite a linear vertical shear instability provided the Richardson number is small enough.In the constant shear case, by contrast, the vertical shear is not (on its own) linearly unstable even at low Richardson number.One might therefore expect to obtain rather different outcomes with and without the inflection point.
These two different options require two different codes, but are both based on the PADDI code that was already used in the works of Cope et al. (2020) and Garaud (2020a).PADDI is a highperformance, triply-periodic pseudo-spectral DNS code, described in Traxler et al. (2011) (see their Appendix A).

Sinusoidal shear in the horizontal direction and uniform shear in the vertical direction
In a first model setup, we apply the same non-dimensional body force in equation ( 8) as in Cope et al. (2020) and Garaud (2020a), namely, In addition, we drive a constant vertical shear using the shearing-box formalism.The numerical methodology is, overall, quite similar to the PADDI code, and is described in detail in Brown & Radko (2021).In this formalism, the simulation grid deforms over time, following the flow of the (prescribed) uniform non-dimensional vertical shear, Ŝ, which is described by the following coordinate transformation: The Fourier modes of the simulation are expanded as complex exponentials in the primed coordinate system, and thus, in addition to horizontal periodicity in x and y, the fields are assumed to be periodic in z ′ instead of z.When the inclination of the grid, ∂z ′ /∂x ′ exceeds Lz /2 Lx , the simulation is remapped to a new coordinate system with an inclination of ∂z ′′ /∂x ′′ = ∂x ′ /∂z ′ − L z /L x .This transformation is performed with the appropriate dealiasing (Delorme 1985).Transforming the system in this way ensures minimal cell deformation and retains periodicity in the new coordinate system.Note that in these simulations, the velocity field û represents the flow perturbations around a state of constant shear, rather than the full flow field.In all simulations, the non-dimensional size of the domain is fixed to Lx = 4π and Ly = Lz = 2π.Selecting a domain that is at least twice a long in the x direction as in the y direction ensures that it contains at least one wavelength of the fastest-growing mode of the primary horizontal shear instability (Cope et al. 2020).The resolution is fixed to have 384 equivalent grid points per interval of length 2π in each direction, which ensures that they are all well-resolved at the parameters selected (see below).
The initial conditions for the simulation are either û(x, y, z, 0) = sin(y) + μ(x, y, z) where μ is some small amplitude white noise, or, the simulation is continued from the end of another simulation at different values of B (this greatly helps reduce compute time).The final statistically stationary state achieved by the simulations seems to be independent of the initial conditions, at least for the cases run here2 .

Sinusoidal shear in the horizontal and vertical directions
In a second model setup, we force the mean shear with a body force that varies sinusoidally in both horizontal and vertical directions, such that As in the case with constant vertical shear described above, the non-dimensional horizontal dimensions of the domain are Lx = 4π and Ly = 2π.To ensure periodicity in the vertical direction, we choose Lz = 2π/ Ŝ unless otherwise specified.The resolution is 384 equivalent grid points per interval of length 2π in each direction, as above, so the results of the constant vertical shear case and the sinusoidal vertical shear case can easily be compared to one another.For the largest shearing   rate considered in our simulations, namely Ŝ = 4, taking Lz = π/2 results in a rather short vertical domain.Worrying that this could artificially impact the results, we also ran simulations with Lz = 2π, but found that the time-averaged mean properties of the flow extracted in the shorter and taller domain sizes are very similar (see Table 3.1).We therefore kept the larger Lz = 2π runs for the Ŝ = 4 case (since they were available) but were satisfied with a shorter domain height Lz = π for the Ŝ = 2 simulations.
The initial conditions for the simulations are either û(x, y, z, 0) = sin(y + Ŝz) + μ(x, y, z) or the simulation is restarted from the end of another simulation at different values of B. Again, the final statistically stationary state achieved appears to be independent of the initial conditions with the same caveat as above.

Simulations parameters
In what follows, we present a large suite of DNS of horizontally and vertically sheared stratified flows.In almost all cases, we have run simulations at the same parameters for the case of constant and sinusoidal vertical shear, which enables us to distinguish between their effects in controlling the properties of the turbulence.The results presented have a Reynolds number Re = 600, which is the largest value we can realistically achieve while at the same time running a large number of simulations.As shown by Garaud (2020a) in the case without vertical shear, Re = 600 is large enough to ensure that the effects of viscosity on the flow dynamics are essentially negligible, except for the largest stratification selected below (B = 100), where they are partially suppressing the turbulence, leading to spatial intermittency in the turbulence intensity.The Prandtl number is fixed to be P r = 0.1, which is substantially smaller than one, while at the same time ensuring that the Péclet number based on the large-scale horizontal flow remains very large (P e = 60).This parameter regime, where both Re and P e are large, but P r is small, is consistent with what one would expect in stellar interiors, but the numerical values achieved are much less extreme than in reality (where Re and P e are generally many orders of magnitude larger, and P r is many orders of magniture smaller, see Garaud 2020a).
Having fixed Re = 600, P e = 60, we then ran a grid of simulations for B ∈ {1, 10, 30, 100} and Ŝ ∈ {0, 0.577, 1, 2, 4}.The range of B was selected to span regimes where the turbulence is effectively unstratified (B = 1), where stratification controls but does not inhibit the turbulence (B = 10 and B = 30) and where stratification is sufficiently strong to intermittently suppress the turbulence (B = 100); see Cope et al. (2020) and Garaud (2020a) for a description and characterization of these various regimes.Note that the threshold values of B that respectively delimit the three regimes (when Ŝ = 0) depend on the Reynolds and Péclet number, so the choices made here are only valid for Re = 600, P e = 60 (see Shah et al. 2023, for a more complete map of parameter space).The range of Ŝ was selected to span regimes of zero and weak vertical shear ( Ŝ = 0.577), moderate vertical shear ( Ŝ = 1, for which the typical horizontal and vertical shear are the same) and large vertical shear (up to Ŝ = 4).We did not run simulations for Ŝ larger than 4, as these require much smaller timesteps and become too computationally expensive.The model parameters, and salient results, are presented in Table 1.
This choice of parameters implies that the vertical shear instabilities expected of the mean shear alone would be of the dynamical kind rather than the diffusive kind, should they take place.Indeed, we can define a Péclet number based on the vertical shear to be where S is the dimensional shearing rate, and L S is the dimensional vertical shear lengthscale.In the models described above, L S = ( Ŝk h ) −1 for the sinusoidally-sheared case (so P e S = P e/ Ŝ), and L S is technically infinite in the constant shear case, but in practice is related to size of the domain L z (so P e S = ŜP e).With the choice of parameters selected here, P e S is always much greater than one.Based on the work of Garaud & Kulenthirarajah (2016), we then anticipate that instabilities of the mean applied vertical shear would drive adiabatic motions rather than thermally diffusive ones in the absence of horizontal shear.We now define the quantity J in as the 'input' Richardson number based on the imposed vertical shear: This quantity is listed in Table 1 for each simulation.As explained in Section 1.1, when the mean Richardson number is significantly smaller than 1, adiabatic perturbations gain more energy from the mean vertical shear than they lose from mixing the stratification (Richardson 1920).Turbulence can be sustained in this manner, and we therefore expect the vertical shear to play an important role in the flow dynamics.By contrast, when the Richardson number (and the Péclet number) are both significantly greater than 1, turbulence cannot be sustained by the mean vertical shear only, and one must rely on the horizontal shear instabilities instead to drive it.Note that most stars generally have J ≫ 1 in their stably stratified radiative zones, even in regions of strong vertical shear (in the bulk of the solar tachocline, for instance, J ∼ 10 3 − 10 4 and only drops to zero very close to the base of the convection zone).Nevertheless, for the sake of completeness, we have run simulations that span a wide range of values of J in including cases where J in < 1. 6.25 1.31 ± 0.08 0.19 ± 0.02 0.14 ± 0.04 Table 1.Input parameters and salient results for the constant vertical shear case (columns 4-6) and the sinusoidal shear case (columns 7-9).In all of the simulations, Re = 600, P e = 60.The rms velocities ûrms and ŵrms , and the vertical eddy scale lz , are extracted from the DNS as described in Section 3.3.(a) denotes simulations from Garaud (2020a).(b) denotes sinusoidal shear simulations in a domain of height Lz = 2π rather than Lz = 2π/ Ŝ.

Temporal evolution of selected simulations
In this section we describe the typical temporal evolution and route to nonlinear saturation of the turbulence for selected simulations to illustrate the range of possible behaviors of the two model systems.We begin by looking in Figure 2 at a simulation which has a relatively large input Richardson number, anticipating that it should be dominated by horizontal shear instabilities.It has a mean sinusoidal shear in both y and z directions, with B = 10, Ŝ = 1 (so J in = 10), and is initialized from noise at t = 0.The top panel shows as solid lines the quantities ûrms (t), vrms (t) and ŵrms (t), defined as ûrms (t) = ⟨û 2 ⟩ 1/2 , vrms (t) = ⟨v 2 ⟩ 1/2 , and ŵrms (t where the brackets denote a volume average.Snapshots of the streamwise velocity field at selected times are shown in the bottom panels.At early times, the mean flow driven by the body force F = sin(y + Ŝz)e x simply grows linearly with time (see first snapshot).This flow is linearly unstable to shear instabilities, and perturbations begin to grow.They reach a sufficiently large amplitude to disrupt that mean flow around t = 10.The perturbations are spatially complex (see second snapshot), and their nature depends quite strongly on the dominant instability at any given point in time.While they are interesting in their own right, they are not the subject of this investigation as they rapidly evolve into fully-developed turbulence that saturates into a statistically stationary state (see third snapshot).
For this particular choice of input parameters, the statistically stationary state achieved looks quite similar to the one described in Garaud (2020a) for simulations with no vertical shear.A horizontal shear instability of the mean flow causes and maintains large-scale horizontal meanders.After saturation, these are characterized by typical horizontal velocity components that are of the same order of magnitude, i.e., ûrms ≃ vrms .The meanders decouple vertically on some short vertical lengthscale, and the resulting streamwise flow looks quite different from the one that is directly forced by F (this can be seen by comparing the first and third snapshots).This decoupling generates substantial emergent, small-scale vertical shear that seems unstable to vertical shear instabilities even though the mean shear is not and therefore drives vertical mixing on small scales.The strong stratification implies that the typical vertical velocities are much smaller than the horizontal ones, so ŵrms ≪ ûrms , vrms .
Figure 3, by contrast, shows two simulations that have a relatively low input Richardson number, namely J in = 0.625 (B = 10, Ŝ = 4).Both were restarted from other simulations at the same shearing rate but a different stratification.The solid lines show the case where the background shear is constant, while the dashed lines show the case where the background shear is sinusoidal.This time, we see that vrms is substantially smaller than ûrms in both simulations once in the statistically stationary state.We find that ŵrms ∼ vrms ≪ ûrms , consistent with the findings of Garaud & Kulenthirarajah (2016) and Garaud et al. (2017) for cases without horizontal shear.The two simulations are otherwise quite different, however.In the sinusoidal shear case, we see in the corresponding snapshot that the streamwise velocity field û is dominated by a coherent mean flow ∝ sin(y + Ŝz).This implies that it contains the necessary inflection points to directly trigger a vertical shear instability since J in is small.The constant shear case, as discussed earlier, is linearly stable to the vertical shear instability.However, the emergence of meanders from the horizontal shear instability rapidly creates the vertical structure needed to trigger the vertical shear instability (see more on this in the next section).Because the meanders are necessary for saturation, the velocity components satisfy ŵrms ≤ vrms ≪ ûrms in this case.
Finally, we note that while most simulations achieve a statistically stationary state relatively rapidly, in a few instances we have observed a long-lived transient flow dominated by vertical shear instabilities, which eventually transitions into a statistically stationary state that is dominated by horizontal shear instabilities.This can be seen in the case presented in Figure 4, which has an applied mean sinusoidal shear in both the y and z directions, with B = 30, Ŝ = 2, and was initialized at t ≃ 20 from a corresponding simulation with a lower stratification (B = 10) and identical shearing rate.At early times, the turbulence has properties that are characteristic of simulations dominated by vertical shear instabilities (with ŵrms ≃ vrms ≪ ûrms ).This state lasts for about 40 time units, and  appears to be relatively stationary in time.However, around t = 60 a complete reorganization of the flow takes place.We see that vrms increases rapidly while ûrms decreases, owing to the growth and saturation of a large-scale horizontal shear instability.A new (and this time true) statistically stationary state is established, whose properties are qualitatively similar to the ones described earlier for the B = 10, Ŝ = 1 case.Notably, we see that once again ŵrms ≪ ûrms ≃ vrms .
Based on these results, we have run all simulations for at least 50 time units once it appears that a statistically steady state has been reached (and often for much longer).

Qualitative properties of the statistically stationary state
We now investigate visually the properties of the statistically stationary state ultimately achieved by each of the simulations.Figure 5 shows typical snapshots of the streamwise component of the flow û in the plane x = 0 for many different simulations.The left-most column shows simulations with B = 1 (effectively unstratified), the middle columns show B = 10 and B = 30 (intermediate stratification) and the rightmost column shows B = 100 (strong stratification).The top row has Ŝ = 0, for reference, and is taken from the simulations of Garaud (2020a).The remainder of the figure shows cases with Ŝ = 1, 2 and 4, and two rows are shown each time: for each value of Ŝ the vertical shear is constant in the upper row, and sinusoidal in the lower row.Figure 5 reveals several important trends, with increasing stratification, and with increasing shear.In the absence of mean vertical shear (top row, Ŝ = 0, J in = ∞), we see the results of Garaud (2020a).For B = 1, the fluid is effectively unstratified.The applied force drives a mean flow (visible in the snapshot) subject to horizontal shear instabilities that rapidly evolve into 3D turbulence.The latter is almost isotropic except near the forcing scale.As the stratification increases (B = 10, 30), meanders in the streamwise flow appear (visible as a zig-zag pattern in û).These, as discussed earlier, are driven by the horizontal shear instability and become shallower as B increases.The emergent vertical shear ∂ û/∂z correspondingly increases and the flow adapts to be marginally unstable to a small-scale vertical shear instability, even though there is no mean vertical shear (see Section 3.3).For very large stratification (B = 100), the vertical scale of the flow meanders becomes so small that they begin to feel the effects of viscosity; the turbulence becomes intermittent.
When Ŝ ̸ = 0 but J in ≫ 1 (i.e. when Ŝ = 1 for any B > 1, Ŝ = 2 for B = 30 or B = 100, and Ŝ = 4 for B = 100), we see from Figure 5 that the snapshots are similar to those of the corresponding Ŝ = 0 simulation, regardless of the method by which the vertical shear is forced (constant or sinusoidal).Based on the discussion of Section 3.1, this is not surprising, as we would not expect the mean vertical shear to play a significant dynamical role, even though it is greater than the horizontal shear.
By contrast, when J in ≪ 1 (e.g. when Ŝ = 4 for B ≤ 10), the flow depends sensitively on the manner in which the vertical shear is forced, confirming what was found in the previous section.In the sinusoidal shear case, û is clearly in phase with the forcing, and is linearly unstable to the vertical shear instability, whose growth rate is O( Ŝ) since the stratification is weak.This is much larger than   the horizontal shear instability growth rate (which is of order unity in the non-dimensionalization selected), which explains why we do not see any meanders (no zigzag patterns).In the constant shear case, by contrast, the vertical shear alone is not linearly unstable, and it is the development of flow meanders that triggers the vertical shear instability.We indeed see that û has strong meanders, and moreover, we now also see that these develop on roughly the same vertical scale as those of corresponding Ŝ = 0, 1 and 2 cases with the same value of B, suggesting that they arise from the same horizontal shear instability.
Finally, when J in is close to one, the outcome is not always as easily predictable.For instance, with Ŝ = 4 and B = 30, J in = 2 and one might have expected the flow to be dominated by horizontal shear instabilities.However, the corresponding snapshot for the sinusoidal forcing case in Figure 5 looks much more similar to those of the simulations where J in ≪ 1, suggesting that it is probably dominated by vertical shear instabilities instead.Cases with J in ∼ O(1) therefore require a more quantitative inspection, which is presented in the next section.

Quantitative analysis of the results
We now take a closer look at the quantitative properties of each simulation.We first characterize the mean flow that is driven in each case and demonstrate that by constructing a Richardson number based on the actual mean vertical shear (rather than the input parameter Ŝ), it is possible to explain the qualitative behavior of the resulting turbulent flow more accurately.Garaud (2020a).For Ŝ ∈ {1, 2, 4}, the top row shows snapshots for the constant vertical shear case (Sect.2.1), and the bottom one for the sinusoidal shear case (Sect.2.2).The color bar is in units of ûrms given for each simulation in Table 3.1.All simulations were run in a domain with height Lz = 2π except for the sinusoidally-forced Ŝ = 2 case, which was run in a domain with Lz = π (two copies of the snapshot are thus shown).

Mean flow properties and associated regime
In what follows, we define the mean flow to be the component of the streamwise flow û that projects on the particular Fourier mode whose spatial structure is the same as that of the forcing.In the case of the constant shear profile, the mean flow is defined as ŪC (y, z) = Ûm sin(y) + Ŝz, where Ûm = ⟨û(x, y, z, t) sin(y)⟩ t ⟨sin 2 (y)⟩ = 2⟨û(x, y, z, t) sin(y)⟩ t , (recalling that û in this case is defined as the perturbation away from the constant shear flow) so the mean vertical shear is simply Ŝm = Ŝ.In this expression, ⟨•⟩ t denotes a volume and time average over the statistically stationary phase.
In the sinusoidal shear case, the mean flow is defined as ŪS (y, z) = Ûm sin(y + Ŝz), where Ûm = ⟨û(x, y, z, t) sin(y + Ŝz)⟩ t so we define the quantity Ŝm = Ŝ Ûm and also call it the mean vertical shear3 .Using this definition, we can construct a 'typical' Richardson number of the mean flow, J m , based on the root-mean-square (rms) shear of the mean flow (namely [ L−1 z Lz 0 (d Ū /dz) 2 dz] 1/2 , which is Ŝm for the constant shear case, but is Ŝm / √ 2 for the sinusoidal shear case): Note that choosing to define J m using the rms shear of the mean flow, rather than its amplitude Ŝm , eases the comparison of J m with the Richardson number of the full flow field J rms defined in the next section.
The left panel of Figure 6 shows J C (solid symbols) and J S (open symbols) as a function of Ŝ for different values of B. These Richardson numbers measure the strength of the typical mean vertical shear relative to the imposed stratification.Using the same symbols, the right panel of Figure 6 shows Ûm .For the constant shear case, J C = J in = B/ Ŝ2 , which is easily predictable from the input parameters (and is shown in the straight lines / full symbols).As discussed earlier, the simulations span a wide range of Richardson numbers in the interval [0.05, 500].In the sinusoidal shear case, however, J S is allowed to deviate from J in .We see that J S is close to J in for strongly stratified simulations, but can be substantially smaller than J in for more weakly stratified systems (see, in particular, B = 10 with Ŝ = 2, and B = 30 with Ŝ = 4, where J in is somewhat greater than 1, but J S < 1).This confirms our earlier conclusion based on the snapshots only (see Figure 5) that these simulations are indeed dominated by vertical shear instabilities.In Figure 6 and in all that follows, all of the simulations that have J m ≤ 1 are marked with solid and open circles, respectively, while those with J m ≫ 1 are marked with squares.We have also marked those with J m ∈ [1, 4] with triangles, to indicate that these are borderline cases.
To understand why J S can be much smaller than J in in the sinusoidal shear case, we look at the mean flow amplitude Ûm in the right panel of Figure 6.In all strongly stratified cases (square symbols), Ûm ≃ 1.This is consistent with the dynamics being dominated by horizontal shear instabilities, where in the non-dimensionalization selected here, large-scale meanders and 'pancake' eddies with horizontal lengthscale ly ∼ O(1) create O(1) horizontal Reynolds stresses to balance the O(1) forcing.For simulations dominated by the vertical shear instability (circular symbols), however, we see that Ûm can be substantially larger.This is because the turbulent eddies have ly ∼ lz ≪ 1 in that regime (Garaud & Kulenthirarajah 2016;Garaud et al. 2017).As a result, Ûm must be correspondingly larger to ensure that the Reynolds stresses match the O(1) forcing.With a large Ûm , Ŝm is also necessarily large, and J S can be substantially smaller than J in .Finally, we note that this discussion is only relevant for the interpretation of the simulation results.Observations of the mean shear in real stars (or the output of stellar evolution codes that include such information) provide Ŝm directly, so there is no ambiguity on the Richardson number.In what follows, we therefore report all the remaining results in terms of Ŝm or in terms of the corresponding Richardson number J m (J C or J S ).

Properties of the emergent shear due to the flow meanders
Having identified that two regimes are expected depending on the size of the typical Richardson number of the mean flow (J m = J C , J S ), we now investigate in turn the properties of the turbulence in each case.To do so, we begin by constructing a diagnostic of the emergent shear created by the horizontal flow meanders.We first average the streamwise flow horizontally: We then define the rms vertical shear associated with ū as ŝrms = ⟨s 2 (y, z, t)⟩ where s = Ŝ + ∂ ū/∂z in the constant shear case, and s = ∂ ū/∂z in the sinusoidal shear case.Note that we define ŝrms from the horizontally-averaged flow ū rather than local flow û, because a shear flow has to be coherent over some significant distance in the streamwise direction to drive vertical shear instabilities.With that definition, ŝrms ≃ Ŝm in the absence of meanders for both constant and sinusoidal shear cases, so the difference between ŝrms and Ŝm can be identified as the typical meander-induced shear.
The values of ŝrms extracted from each simulation are compared with Ŝm in Figure 7 (left panel).Consistent with their definition, we always have ŝrms ≥ Ŝm .Very clearly, we see that simulations that are dominated by vertical shear instabilities in the sinusoidal shear case (open circles) have ŝrms ≃ Ŝm , which is consistent with the visual impression originally gained from Figure 5 that these cases do not have substantial flow meanders.In most other cases, ŝrms ≫ Ŝm , indicating that meanders are present and play a key role in maintaining the turbulence.In cases that are dominated by horizontal shear instabilities (i.e., J m ≫ 1, open and filled square symbols) ŝrms seems to be independent of the mean vertical shear (and of the manner it is forced) but instead only depends on the stratification (via B).To better understand why this is the case, we now turn to the right panel of Figure 7, which shows the same information but in terms of the corresponding 'rms Richardson number' against J m (J S for the sinusoidal case, J C for the constant shear case).J rms can be understood as a Richardson number of emergent, meander-induced, small-scale vertical shear, by comparison of J m which is a Richardson number of the large-scale vertical shear.We see very clearly that whenever J m ≫ 1, J rms ≃ 2. In other words, the meanders adjust themselves precisely so that they are in a neutrally stable state with respect to the energetics of vertical shear instabilities.Indeed, it is easy to verify that if J rms ≃ 2 then an equivalent Richardson number based on the characteristic maximum vertical shear (instead of the rms shear) would be closer to 1.When this is the case, smallscale vertical eddies extract kinetic energy from the mean shear and provide it to the background stratification in the form of potential energy in equal amounts (Richardson 1920) Dimensionally, J rms = O(1) implies that s rms = O(N ).Since the horizontal velocity of the meanders is O(U m ), the vertical scale of the meanders l z must be O(U m /N ) to ensure that s rms = O(U m /l z ) = O(N ).In short, Figure 7 shows that the meanders of the flow adjust themselves to be close to marginal stability for the vertical shear instability, and to do so self-consistently, must have vertical structure on scales l z = O(U m /N ).This characteristic scale for stratified turbulence at high Péclet number has long been known in the geophysical literature (e.g.Holford & Linden 1999;Billant & Chomaz 2001;Brethouwer et al. 2007).It is entirely consistent (from a theoretical perspective) with the recent multiscale analysis of Chini et al. (2022) and Shah et al. (2023), who argued that this scaling is expected as long as P e ≫ B 1/2 and as long as viscous effects can be neglected.Our results suggest that it continues to hold as long as the Richardson number based on that mean vertical shear is substantially larger than one (e.g., J m ≥ 5, say).
For weaker stratification, namely J m ≤ O(1), we saw in Section 3.2.1 that the statistically stationary state achieved depends on the manner in which the flow is forced.For the sinusoidal case (open circles), we saw that meanders are essentially absent, so J rms ≃ J m .In the constant-shear case (filled circles and triangles), on the other hand, the mean flow is linearly stable to the vertical shear instability, so the latter relies on the onset of meandering via horizontal shear instabilities to set in.In that case, ŝrms must necessarily exceed Ŝm .In practice, however, we find that ŝrms is only a little larger than Ŝm or equivalently, J rms is a little smaller than J m .

Turbulent mixing
We now study vertical mixing by the turbulence in each simulation.As introduced in Section 1, an important quantity of interest for stellar evolution calculations is the so-called vertical mixing coefficient D mix , which is the turbulent diffusivity of momentum or a trace chemical element.As long as the flow is not in a viscously influenced regime where the turbulence is intermittently suppressed, this coefficient can be estimated (non-dimensionally) from the product of a characteristic vertical velocity of the turbulence ( ŵrms ) and a characteristic vertical eddy height ( lz ): where c is some constant of order unity.The quantity ŵrms is easily obtained from the simulations by taking the time average of ŵrms (t) once the system has reached a statistically stationary state.
For the one shown in Figure 2, for instance, that average is taken between t = 30 and t = 95, and is equal to ŵrms = 0.59 ± 0.03, where the ±0.03 captures the standard deviation of ŵrms (t) around the mean value 0.59.The values of ŵrms thus computed, for all available simulations, are presented in Table 1.
To calculate the eddy height lz , following Garaud (2020a), we compute the vertical auto-correlation of the vertical velocity field, defined as: The vertical lengthscale of the eddies at a given time t n is then defined as the width at half maximum of A w ( l, t n ).In other words, lz (t n ) is the solution of Finally, lz is defined as the average of all lengthscales thus computed for a given simulation once the system is in the statistically stationary state.The results are presented in Table 1.Note that the full velocity field is only saved every 10,000 timesteps at this value of Re, which, owing to the Courant condition, approximately corresponds to intervals ∆t n of 0.5 to 3 time units depending on the simulation.In general lz is then computed from an average of 20-60 snapshots.Because this number is relatively small, the statistics of lz are not as robust as those of ŵrms .Fig. 8 shows lz (top left) and ŵrms (top right) as a function of the typical vertical shearing rate of the mean flow Ŝm for various values of the stratification parameter B. Consistent with the qualitative picture revealed by Figure 5 and originally presented in Garaud (2020a) for the Ŝm = 0 dataset, we see that both the vertical eddy scale and the vertical velocity decrease with increasing stratification (increasing B).As Ŝm increases at fixed values of B = 1, 10 and 30 (but not B = 100, see below), we see that the mean shear, and the manner in which it is forced, have little influence on lz and ŵrms as long as J m ≥ 1 (open and closed square and triangular symbols), but begin to affect them when J m < 1.Overall, this is not surprising based on the discussions of the previous sections.By contrast, the effects of a strong (J m < 1) constant or sinusoidal shear on lz and ŵrms are complex (sometimes increasing them, sometimes decreasing them).Since this limit is fairly rare in stellar interiors, we do not investigate it further in this paper.
In the strongly stratified case (B = 100), we were initially surprised to find that the mean shear does affect the turbulence even though J m ≫ 1 -in this case, it does so as soon as Ŝm > 1.However, this behavior can be understood by noting that B = 100 is in a viscously influenced regime where the turbulence is spatially patchy at Ŝ = 0 (Garaud 2020a).When that is the case, the fraction of the domain covered by turbulent patches depends quite sensitively on the existence of nonlinear mechanisms that help the turbulence self-sustain and percolate through the flow (Deusebio et al. 2015;Avila et al. 2023).We hypothesize that a strong vertical shear might stretch these turbulent patches and disrupt the self-sustaining process.The effect is clearly different in the constant and sinusoidal shear cases, sometimes increasing or decreasing lz or ŵrms .We defer the analysis of the intermittent regime to future work, as we believe that it is not particularly relevant to stellar interiors.Indeed Shah et al. (2023) demonstrated that stratification is strong enough for viscosity to influence the turbulence only when B > Re (unless thermal diffusion is important, see the original paper for detail), or equivalently in dimensional terms when In stars, this corresponds to regions of very small horizontal shear, in which the corresponding vertical transport would likely be negligible anyway.
The lower two panels of Figure 8 show Dmix = ŵrms lz as functions of Ŝm (left) and J m = J C , J S (right).The values of Dmix for B = 100 are shown for reference, but should be ignored as the formula used is not reliable in that case because of the patchiness of the turbulence.Following our findings for the behavior of lz and ŵrms , we see that Dmix decreases steadily with increasing stratification and is roughly independent of the shearing rate and the manner in which the vertical shear is forced as long as J m ≥ 1 (squares and triangles).
Based on the results presented in this section we can therefore make the following statement with reasonable confidence: as long as thermal diffusion and viscosity are negligible (so the turbulence is both fully developed and adiabatic), a mean vertical shear has no effect on mixing driven by horizontal shear instabilities when its mean Richardson number is greater than one.While this is not particularly surprising in hindsight, it is worth emphasizing that this statement likely applies to the vast majority of stellar shear layers, even when the vertical shearing rate is much greater than the horizontal one, which is perhaps a little less intuitive.A good example of this is the solar tachocline, where the vertical shear is over 20 times stronger than the horizontal shear (Charbonneau et al. 1999).The fact that J = O(10 3 − 10 4 ) ≫ 1 in the tachocline (Garaud 2020a) shows that unless horizontal shear instabilities are somehow suppressed by processes not included here (see the discussion below), the turbulence they generate dominates vertical mixing, while the mean vertical shear can essentially be ignored in the computation of D mix ! 4. DISCUSSION AND CONCLUSIONS In this paper, we have used DNS to investigate the combined effects of vertical and horizontal shear on mixing by shear instabilities in the stably stratified regions of stars.We ignored rotation and magnetic fields for now.The horizontal shear was forced to be sinusoidal in the spanwise direction (here, y), and we compared two types of vertical shear profiles (constant vs. sinusoidal).The numerical simulation parameters were selected to be in a regime where thermal diffusion is negligible on the large horizontal and vertical scales, while having a low Prandtl number.This regime (rather than the numerical parameters themselves) is consistent with what is expected in a typical stellar interior (Shah et al. 2023).
Within this set of assumptions, we found that two dynamical regimes naturally emerge, depending on the Richardson number J m = N 2 /S 2 m , where S m is a typical value of the large-scale vertical shear.When J m is of order unity or smaller, vertical shear instabilities dominate the system dynamics.The form they take (and the amount of mixing they cause) depends sensitively on the geometry of the vertical shear and more specifically on whether it can trigger a linear instability or not.In addition, even though we have only considered regions with a constant stratification in this paper, there is ample evidence from the oceanographic literature (see the review by Caulfield 2021) showing that the outcome will depend sensitively on the shape of the density profile and must be studied on a case-by-case basis.Regions of stars with J m < O(1) are rare, however, so this regime is not particularly relevant to stellar evolution in general.
In the far more likely instances where J m ≫ 1, horizontal shear instabilities control the system dynamics, and the mean vertical shear has little influence on vertical mixing (unless the stratification is so strong that the turbulence itself is viscously suppressed, see Section 3.3.3).This is because the large-scale horizontal perturbations excited by the horizontal shear instability (meanders, or pancake vortices) decouple in the vertical direction on small vertical scales and drive a much stronger vertical shear.More specifically, we showed in Section 3.3 that this emergent meander-induced shear s rms satisfies the marginal stability condition J rms = N 2 /s 2 rms = O(1), so s rms = O(N ), which by definition of this regime is much larger than the mean shear S m .As a result, vertical mixing is entirely driven by the emergent rather than mean vertical shear.When that is the case, models for mixing by horizontal shear instabilities alone can be safely applied -at least once the controversy as to which model should be used is resolved, see Section 1.2 for detail.
There are, of course, several caveats to this statement.First, note that we have ignored in this work the effects of thermal diffusion by purposefully selecting parameters in which the turbulence generated is known to be relatively adiabatic.This, we believe (see Garaud & Kulenthirarajah 2016;Garaud 2020a;Shah et al. 2023), is the more common situation in stars, but that belief is somewhat controversial (see, e.g.Skoutnev 2023).Presumably, the criterion will be similar in the thermally diffusive case -i.e the mean shear can be ignored when it is much smaller than the meander-induced shear -but this hypothesis will need to be verified.More importantly, we have also ignored in this paper several effects that could suppress the horizontal shear instabilities themselves.Indeed, sufficiently strong rotation, or strong streamwise magnetic fields (both of which are present in the solar tachocline, for example) can stabilize purely two-dimensional horizontal shear instabilities (see, e.g., the discussion in Garaud 2021).When that is the case, the vertical shear may be needed to provide an alternative route to turbulence and must not be neglected.
Nevertheless, by carefully characterizing the conditions under which each instability can grow, and how much mixing they cause, we can gradually make progress in advancing our understanding of their role in stellar evolution and parametrize their effect in 1D (e.g.MESA, Paxton et al. 2013) and now 2D (cf.ESTER, Mombarg et al. 2023) stellar evolution codes.This paper has taken one more step in this direction.
y d 2 X p x 3 5 2 P e W n D y m U P 4 A + f z B 9 Y T j 1 E = < / l a t e x i t > <l a t e x i t s h a 1 _ b a s e 6 4 = " 6 B q O z 7 7 I m T M B M M 6 6 I S 3 Y k N 8 g P F A = " > A A A B 7 X i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K r 2 P Q i 8 c I 5 g H J E m Y n s 8 m Y 2 d l l p l c I S / 7 B i w d F v P o / 3 v w b J 8 k e N L G g o a j q p r s r S K Q w 6 L r f T m F l d W 1 9 o 7 h Z 2 t r e 2 d 0 r 7 x 8 0 y d 2 X p x 3 5 2 P e W n D y m U P 4 A + f z B 9 Y T j 1 E = < / l a t e x i t > <l a t e x i t s h a 1 _ b a s e 6 4 = " 6 B q O z 7 7 I m T M B M M 6 6 I S 3 Y k N 8 g P F A = " > A A A B 7 X i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K r 2 P Q i 8 c I 5 g H J E m Y n s 8 m Y 2 d l l p l c I S / 7 B i w d F v P o / 3 v w b J 8 k e N L G g o a j q p r s r S K Q w 6 L r f T m F l d W 1 9 o 7 h Z 2 t r e 2 d 0 r 7 x 8 0
t e x i t s h a 1 _ b a s e 6 4 = " 6 B q O z 7 7 I m T M B M M 6 6 I S 3 Y k N 8 g P F A = " > A A A B 7 X i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K r 2 P Q i 8 c I 5 g H J E m Y n s 8 m Y 2 d l l p l c I S / 7 B i w d F v P o / 3 v w b J 8 k e N L G g o a j q p r s r S K Q w 6 L r f T m F l d W 1 9 o 7 h Z 2 t r e 2 d 0 r 7 x 8 0T Z x q x h s s l r F u B 9 R w K R R v o E D J 2 4 n m N A o k b w W j 2 6 n f e u L a i F g 9 4 D j h f k Q H S o S C U b R S s z u k m E 1 6 5 Y p b d W c g y 8 T L S Q V y 1 H v l r 2 4 / Z m n E F T J J j e l 4 b o J + R j U K J v m k 1 E 0 N T y g b 0 Q H v W K p o x I 2 f z a 6 d k B O r 9 E k Y a 1 s K y U z 9 P Z H R y J h x F N j O i O L Q L H p T 8 T + v k 2 J 4 7 W d C J S ly x e a L w l Q S j M n 0 d d I X m j O U Y 0 s o 0 8 L e S t i Q a s r Q B l S y I X i L L y + T 5 l n V u 6 x e 3 J 9 X a j d 5 H E U 4 g m M 4 B Q + u o A Z 3 U I c G M H i E Z 3 i F N y d 2 X p x 3 5 2 P e W n D y m U P 4 A + f z B 9 Y T j 1 E = < / l a t e x i t > ˆ< l a t e x i t s h a 1 _ b a s e 6 4 = " 6 B q O z 7 7 I m T M B M M 6 6 I S 3 Y k N 8 g P F A = " > A A A B 7 X i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K r 2 P Q i 8 c I 5 g H J E m Y n s 8 m Y 2 d l l p l c I S / 7 B i w d F v P o / 3 v w b J 8 k e N L G g o a j q p r s r S K Q w 6 L r f T m F l d W 1 9 o 7 h Z 2 t r e 2 d 0 r 7 x 8 0 y x e a L w l Q S j M n 0 d d I X m j O U Y 0 s o 0 8 L e S t i Q a s r Q B l S y I X i L L y + T 5 l n V u 6 x e 3 J 9 X a j d 5 H E U 4 g m M 4 B Q + u o A Z 3 U I c G M H i E Z 3 i F N y d 2 X p x 3 5 2 P e W n D y m U P 4 A + f z B 9 Y T j 1 E = < / l a t e x i t > ˆ< l a t e x i t s h a 1 _ b a s e 6 4 = " 6 B q O z 7 7 I m T M B M M 6 6 I S 3 Y k N 8 g P F A = " > A A A B 7 X i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K r 2 P Q i 8 c I 5 g H J E m Y n s 8 m Y 2 d l l p l c I S / 7 B i w d Fv P o / 3 v w b J 8 k e N L G g o a j q p r s r S K Q w 6 L r f T m F l d W 1 9 o 7 h Z 2 t r e 2 d 0 r 7 x 8 0 T Z x q x h s s l r F u B 9 R w K R R v o E D J 2 4 n m N A o k b w W j 2 6 n f e u L a i F g 9 4 D j h f k Q H S o S C U b R S s z u k m E 1 6 5 Y p b d W c g y 8 T L S Q V y 1 H v l r 2 4 / Z m n E F T J J j e l 4 b o J + R j U K J v m k 1 E 0 N T y g b 0 Q H v W K p o x I 2 f z a 6 d k B O r 9 E k Y a 1 s K y U z 9 P Z H R y J h x F N j O i O L Q L H p T 8 T + v k 2 J 4 7 W d C J S l y x e a L w l Q S j M n 0 d d I X m j O U Y 0 s o 0 8 L e S t i Q a s r Q B l S y I X i L L y + T 5 l n V u 6 x e 3 J 9 X a j d 5 H E U 4 g m M 4 B Q + u o A Z 3 U I c G M H i E Z 3 i F Ny d 2 X p x 3 5 2 P e W n D y m U P 4 A + f z B 9 Y T j 1 E = < / l a t e x i t > Figure 2. Top: Rms velocity components ûrms , vrms , and ŵrms as a function of time, for a simulation with Re = 600, P e = 60, B = 10, Ŝ = 1 where the vertical shear is sinusoidal.Bottom: Snapshots of the streamwise velocity û at selected times (t = 9, t = 12 and t = 70, see triangular symbols for corresponding times in top panel) showing the mean flow before onset of instability, the onset of instability, and the saturated state, respectively.
t e x i t s h a 1 _ b a s e 6 4 = " 6 B q O z 7 7 I m T M B M M 6 6 I S 3 Y k N 8 g P F A = " > A A A B 7 X i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K r 2 P Q i 8 c I 5 g H J E m Y n s 8 m Y 2 d l l p l c I S / 7 B i w d F v P o / 3 v w b J 8 k e N L G g o a j q p r s r S K Q w 6 L r f T m F l d W 1 9 o 7 h Z 2 t r e 2 d 0 r 7 x 8 0T Z x q x h s s l r F u B 9 R w K R R v o E D J 2 4 n m N A o k b w W j 2 6 n f e u L a i F g 9 4 D j h f k Q H S o S C U b R S s z u k m E 1 6 5 Y p b d W c g y 8 T L S Q V y 1 H v l r 2 4 / Z m n E F T J J j e l 4 b o J + R j U K J v m k 1 E 0 N T y g b 0 Q H v W K p o x I 2 f z a 6 d k B O r 9 E k Y a 1 s K y U z 9 P Z H R y J h x F N j O i O L Q L H p T 8 T + v k 2 J 4 7 W d C J S l y x e a L w l Q S j M n 0 d d I X m j O U Y 0 s o 0 8 L e S t i Q a s r Q B l S y I X i L L y + T 5 l n V u 6 x e 3 J 9 X a j d 5 H E U 4 g m M 4 B Q + u o A Z 3 U I c G M H i E Z 3 i F Ny d 2 X p x 3 5 2 P e W n D y m U P 4 A + f z B 9 Y T j 1 E = < / l a t e x i t > ˆ< l a t e x i t s h a 1 _ b a s e 6 4 = " 6 B q O z 7 7 I m T M B M M 6 6 I S 3 Y k N 8 g P F A = " > A A A B 7 X i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K r 2 P Q i 8 c I 5 g H J E m Y n s 8 m Y 2 d l l p l c I S / 7 B i w d F v P o / 3 v w b J 8 k e N L G g o a j q p r s r S K Q w 6 L r f T m F l d W 1 9 o 7 h Z 2 t r e 2 d 0 r 7 x 8 0T Z x q x h s s l r F u B 9 R w K R R v o E D J 2 4 n m N A o k b w W j 2 6 n f e u L a i F g 9 4 D j h f k Q H S o S C U b R S s z u k m E 1 6 5 Y p b d W c g y 8 T L S Q V y 1 H v l r 2 4 / Z m n E F T J J j e l 4 b o J + R j U K J v m k 1 E 0 N T y g b 0 Q H v W K p o x I 2 f z a 6 d k B O r 9 E k Y a 1 s K y U z 9 P Z H R y J h x F N j O i O L Q L H p T 8 T + v k 2 J 4 7 W d C J S ly x e a L w l Q S j M n 0 d d I X m j O U Y 0 s o 0 8 L e S t i Q a s r Q B l S y I X i L L y + T 5 l n V u 6 x e 3 J 9 X a j d 5 H E U 4 g m M 4 B Q + u o A Z 3 U I c G M H i E Z 3 i F N y d 2 X p x 3 5 2 P e W n D y m U P 4 A + f z B 9 Y T j 1 E = < / l a t e x i t > Figure 3. Top: Rms velocity components ûrms , vrms , and ŵrms as a function of time, for simulations with Re = 600, P e = 60, B = 10, Ŝ = 4.The solid lines show the constant shear case, while the dashed lines show the sinusoidal shear case.Bottom left: Snapshot of û at t = 70 in the constant shear case.Bottom right: Snapshot of û at t = 70 in the sinusoidal shear case.
t e x i t s h a 1 _ b a s e 6 4 = " 6 B q O z 7 7 I m T M B M M 6 6 I S 3 Y k N 8 g P F A = " > A A A B 7 X i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K r 2 P Q i 8 c I 5 g H J E m Y n s 8 m Y 2 d l l p l c I S / 7 B i w d F v P o / 3 v w b J 8 k e N L G g o a j q p r s r S K Q w 6 L r f T m F l d W 1 9 o 7 h Z 2 t r e 2 d 0 r 7 x 8 0 y d 2 X p x 3 5 2 P e W n D y m U P 4 A + f z B 9 Y T j 1 E = < / l a t e x i t > ˆ< l a t e x i t s h a 1 _ b a s e 6 4 = " 6 B q O z 7 7 I m T M B M M 6 6 I S 3 Y k N 8 g P F A = " > A A A B 7 X i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K r 2 P Q i 8 c I 5 g H J E m Y n s 8 m Y 2 d l l p l c I S / 7 B i w d F v P o / 3 v w b J 8 k e N L G g o a j q p r s r S K Q w 6 L r f T m F l d W 1 9 o 7 h Z 2 t r e 2 d 0 r 7 x 8 0 y d 2 X p x 3 5 2 P e W n D y m U P 4 A + f z B 9 Y T j 1 E = < / l a t e x i t > Figure 4. Top Rms velocity components ûrms , vrms , and ŵrms as a function of time, for a simulation with Re = 600, P e = 60, B = 30, Ŝ = 2 where the vertical shear is sinusoidal.Botton: Snapshots of the streamwise velocity û at selected times (t = 40, and t = 90, see corresponding triangular symbols in the top panel.).

Figure 5 .
Figure5.Simulation snapshots of û in the x = 0 plane, for increasing shear Ŝ (from top to bottom) and stratification B (from left to right).Simulations with Ŝ = 0 are fromGaraud (2020a).For Ŝ ∈ {1, 2, 4}, the top row shows snapshots for the constant vertical shear case (Sect.2.1), and the bottom one for the sinusoidal shear case (Sect.2.2).The color bar is in units of ûrms given for each simulation in Table3.1.All simulations were run in a domain with height Lz = 2π except for the sinusoidally-forced Ŝ = 2 case, which was run in a domain with Lz = π (two copies of the snapshot are thus shown).

Figure 6 .
Figure 6.Properties of the mean flow as a function of Ŝ and B. Left: typical Richardson number J m of the mean flow, defined in equation (20) (J C in the constant shear case, solid symbols, and J S in the sinusoidal shear case, open symbols).Symbol color represents B as shown in the legend of the right panel.Simulations with J m ≤ 1 are marked with circles, and are dominated by vertical shear instabilities.Those with J m ≫ 1 are marked with squares and are dominated by horizontal shear instabilities.Intermediate cases are marked with triangles.Right: The same for the mean flow amplitude Ûm (see definition in equations 18 and 19).

Figure 7 .
Figure 7. Properties of the rms vertical shear emerging from the various flow instabilities.Left: ŝrms as a function of the typical mean shear Ŝm ( Ŝ in the constant shear case, solid symbols, and Ŝ Ûm / √ 2 for the sinusoidal shear case, open symbols).Symbol color represents B. Filled symbols are for the constant shear case, and open symbols are for the sinusoidal shear case.Right: The same data presented in terms of Richardson numbers J rms = B/ŝ 2 rms vs. J m = B/ Ŝ2 m .

Figure 8 .
Figure 8. Vertical eddy scale lz (top left), rms vertical velocity ŵrms (top right), and estimated mixing coefficient Dmix = lz ŵrms (bottom left) as a function of the typical mean vertical shearing rate Ŝm and the stratification B (indicated by the line colors, see legend).Filled symbols connected by solid lines correspond to the constant shear case, and open symbols connected by dashed lines correspond to the sinusoidal shear case.The bottom right figure shows Dmix as a function of the typical mean Richardson number J m = (J C , J S )