V Witzke

L J Silvers

B Favier

Evolution and characteristics of forced shear flows in polytropic atmospheres: Large and small Péclet number regimes

Keywords: methods, numerical -stars, interiors -hydrodynamics -instabilities -turbulence

Complex mixing and magnetic field generation occurs within stellar interiors particularly where there is a strong shear flow. To obtain a comprehensive understanding of these processes, it is necessary to study the complex dynamics of shear regions. Due to current observational limitations, it is necessary to investigate the inevitable small-scale dynamics via numerical calculations. Here, we examine direct numerical calculations of a local model of unstable shear flows in a compressible polytropic fluid primarily in a two-dimensional domain, where we focus on determining how key parameters affect the global properties and characteristics of the resulting saturated turbulent phase. We consider the effect of varying both the viscosity and the thermal diffusivity on the non-linear evolution. Moreover, our main focus is to understand the global properties of the saturated phase, in particular estimating for the first time the spread of the shear region from an initially hyperbolic tangent velocity profile. We find that the vertical extent of the mixing region in the saturated regime is generally determined by the initial Richardson number of the system. Further, the characteristic quantities of the turbulence, i.e. typical length-scale and the root-mean-square velocity are found to depend on both the Richardson number, and the thermal diffusivity. Finally, we present our findings of our investigation into saturated flows of a 'secular' shear instability in the low Péclet number regime with large Richardson numbers.

INTRODUCTION

Developing a complete model of stellar dynamics requires a comprehensive understanding of the microphysical and macrophysical processes present in all stellar regions. Current stellar evolution models are challenged by some discrepancies between theory and observations, which can only be resolved by introducing additional mixing (see Pinsonneault 1997, and references therein). In stars, a possible source for such mixing processes is shear-induced turbulence [START_REF] Zahn | Union I. A[END_REF]Schatzman 1977;Endal & Sofia 1978). In addition, complex gas dynamics in stellar interiors is not only important for mixing processes but also plays a crucial role in magnetic field generation (see Miesch & Toomre 2009;Jones, Thompson & Tobias 2010). Therefore, ongoing research focuses on understanding possible hydrodynamical and magnetohydrodynamical instabilities leading to turbulence in various stellar regions (Rüdiger, Kitchatinov & Elstner 2012;Witzke, Silvers & Favier 2015;Garaud & Kulenthirarajah 2016).

E-mail:Veronika.Witzke.1@city.ac.uk (VW); Lara.Silvers.1@city.ac.uk (LJS); Favier@irphe.univ-mrs.fr (BF)

Main-sequence stars have common dynamical elements, such as large-scale shear flows resulting from differential rotation. However, these shear flows can occur in different regions with different characteristics i.e. different transport coefficients, thermal stratification, etc. One important body is the Sun, where considerable effort has recently been directed towards extending our understanding of the tachocline and its role in the solar dynamo (see Silvers 2008, and references therein). In the Sun there is also the near-surface shear layer that is located at the upper boundary of the convection zone. The near-surface shear layer is believed to be important for the magnetic field generation (Brandenburg 2005;Cameron, R. H. & Schüssler, M. 2017), but it has significantly different transport coefficients (Thompson et al. 1996;Miesch & Hindman 2011;[START_REF] Barekat | [END_REF]. Thus non-dimensional numbers such as the Reynolds number or the Péclet number can vary in astrophysical objects. The Péclet number, which is the ratio of advection to temperature diffusion, plays an important role in the linear and non-linear stability of shear flows. In particular, a low Péclet number, i.e. less than unity, affects the stability threshold of shear flows [START_REF] Zahn | Union I. A[END_REF]Lignières et al. 1999;Garaud et al. 2015), as it destabilises the system and facilitates an instability. Such a low Péclet number regime can be reached in stellar regions where there is high thermal diffusivity, which can occur deep in stellar interiors or in envelopes of massive stars (Garaud & Kulenthirarajah 2016). In the classical stability analysis of a vertical shear flow in an ideal fluid, the Richardson number, which corresponds to the ratio of the Brunt-Väisälä frequency and the turnover rate of the shear, gives the measure for stability. A low Richardson number, less than 1/4, is required for the system to become unstable. However, in diffusive systems, where the Péclet number is low, a possible 'secular' shear instability can be present. This type of shear flow instability develops even for high Richardson numbers but only if the thermal diffusivity is large enough to weaken the stable stratification [START_REF] Zahn | Union I. A[END_REF]Lignières, Califano & Mangeney 1999;Garaud, Gallet & Bischoff 2015;Witzke, Silvers & Favier 2015). So far investigations of low Péclet number shear flows have made use of the Boussinesq approximation (Prat et al. 2016), which does not permit the study of a system that is larger than a pressure scale height as is present in stellar interiors. Numerical investigations of the dynamics in stellar interiors are challenging because the viscosity as well as the Prandtl number present in such regions are very low, which leads to high Reynolds numbers and low Péclet numbers on small scales. The considerable computational cost of such numerical calculations is not accessible with current resources. Using considerably larger values of the Prandtl number, or the viscosity, in order to obtain a tractable system is a way forward but this will necessarily lead to some differences in the dynamics that would be found in real astrophysical shear flows. Moreover, the Richardson number in the tachocline is approximated to be O(10). While this approximation is obtained from spatially and time-averaged measurements, turbulent motions can be present on smaller length-scales and time-scales. Recent investigations have focused on the conditions under which low Péclet number flows become linearly or non-linearly unstable (see Prat & Lignières 2013;Garaud et al. 2015;Garaud & Kulenthirarajah 2016) and mixing can occur. However, the relevant parameters that determine whether thermal diffusivity has significant impact on the non-linear dynamics of the system are dictated by the typical length-scales and velocities in the turbulent regime and so they are not known a priori, but can only be determined from fully non-linear hydrodynamical calculations. In stellar interiors these parameters evolve due to complex mechanisms, which are not entirely understood. Moreover, previous numerical studies used a periodic domain in the vertical direction (Garaud & Kulenthirarajah 2016) or a linear shear profile (Prat & Lignières 2014). Both approaches correspond to modelling a very localised part of a larger shearing region. While they have extended our understanding, it remains unclear if the findings of these approaches persist when investigating shear transition regions in their entirety. In this paper we examine a local region of a fully compressible, stratified fluid to improve our understanding of shear regions in stellar interiors. Using a hyperbolic tangent velocity profile permits us to model a larger region of differential rotation, where a sharp shear flow is localised at the middle. Such a shear flow allows for different dynamics to occur compared to the aforementioned studies. Moreover, our setup is not periodic in the vertical direction and highly stratified, which differs significantly from a standard approach where all dimensions are periodic. Our investigations address the following questions: Does the viscosity have an impact on the spread of a shear flow instability and how are the turbulent length-scales affected? What effect does the thermal diffusivity have on the resulting shear region and turbulence characteristics present there? To what extent is the spread of the shear instability controlled by the Richardson number? While we expect viscosity to have an effect on the smallest possible length-scales it is not obvious if the typical length-scale will be significantly affected. We will primarily focus on classical shear instabilities in a two-dimensional setup with different Péclet numbers in order to investigate the resulting saturated regime and what affects the extent of the shear region after saturation. The paper will proceed as follows. In Section 2 the governing equations are given together with the numerical methods used. In Section 3 we examine how unstable shear flows saturate and evolve into quasi-static states when transport coefficients are changed. Focusing on a low Richardson number instability, an overview of the effects of varying the Reynolds number and initial Péclet number is given. Finally, the possibility of a 'secular' instability is investigated, where large Richardson numbers are considered. Here, the non-linear regime is compared to a low Péclet number turbulence induced by a classical low Richardson number instability.

MODEL

We consider an ideal monatomic gas with constant dynamic viscosity, µ, constant thermal conductivity, κ, constant heat capacities c p at constant pressure and c v at constant volume, and with an adiabatic index γ = c p /c v = 5/3. In this study we chose to fix the dynamic viscosity rather than the kinematic viscosity. The domain throughout most of this paper is a x-z plane, where the horizontal xdirection is periodic and the depth in the z-direction is d. The depth of the domain in dimensional units is given by z, where the two boundaries are located at z = 0 and z = d. Most of our calculations were performed in this two-dimensional setup in order to reduce the computational cost. Thus we consider a two-dimensional set of differential equations, i.e. we use three-dimensional arrays, but we neglect the spatial variations in a possible third y-direction, which is perpendicular to the x-z plane, and assume that all quantities are zero in this possible third direction. The full set of dimensionless, differential equations is:

∂ρ ∂t = -∇• (ρu) (1) ∂(ρu) ∂t = σC k ∇ 2 u + 1 3 ∇(∇•u) -∇• (ρuu) -∇p + θ(m + 1)ρ ẑ + F (2) ∂T ∂t = C k σ(γ -1) 2ρ |τ| 2 + γC k ρ ∇ 2 T -∇• (T u) -(γ -2)T ∇• u (3)
where ρ is the density, u the velocity field, T the temperature, p is the pressure, θ denotes the temperature difference across the layer and ẑ is the unit vector in the z-direction. In the dimensionless equations above, all lengths are given in units of the domain depth d, where z = z/d is the non-dimensional vertical length. The temperature and density are recast in units of T t and ρ t , the temperature and density at the top of the layer, and we take the sound-crossing time, which is given by t = d/[(c pc v )T t ] 1/2 , as the reference time.

There are two additional dimensionless numbers in the set of equations above: the Prandtl number, σ = µc p /κ, which is the ratio of viscosity to thermal diffusivity and the thermal diffusivity parameter C k = κt/(ρ t c p d 2 ). The strain rate tensor in equation (3) has the form

τ i j = ∂u j ∂x i + ∂u i ∂x j -δ i j 2 3 ∂u k ∂x k , (4)
where its tensor norm is |τ| 2 = τ ji τ i j . For the basic state a polytropic relation between pressure and density is taken. In this paper the polytropic index m always satisfies the inequality m > 1/(γ -1) = 3/2, such that the atmosphere is stably stratified. Note, that the gravitational acceleration, g, is derived from the hydrostatic equilibrium. It is assumed constant throughout our domain due to the local assumption of a polytropic atmosphere, which leads to g = θ(m + 1), and the Cowling approximation (Cowling 1941). Our parameter choices for all calculations presented in this paper are summarised in Table 1 andTable 2.

The boundary conditions at the top and the bottom of the domain are impermeable and stress-free velocity and fixed temperature:

u z = ∂u x ∂z = 0 at z = 0 and z = 1, (5)
T = 1 at z = 0 and

T = 1 + θ at z = 1. ( 6
)
The dimensionless initial temperature and density profiles are of the form:

T (z) = (1 + θz) (7) ρ(z) = (1 + θz) m . ( 8
)
This basic state corresponds to an equilibrium state if the fluid is at rest. However, we assume that an external force, denoted F in equation (2), sustains the following initial background velocity profile

U 0 = (u 0 (z), 0, 0) T = U 0 2 tanh 2 L u (z -0.5) êx (9)
where U 0 is the shear amplitude and L u is the width of the shear profile. A hyperbolic tangent shear profile was chosen to minimize the boundary effects. Although the relevant parameter values for the investigations were chosen to keep the spread of the instability confined in the middle domain, it is impossible to avoid completely any boundary effects. Furthermore, the boundary conditions introduced in equation ( 5) restrict the shear profile to values of L u that will result in a low enough value of the z-derivative of u 0 at the boundaries. A visualisation of the general form of shear, density and temperature profiles used in this paper can be found in Witzke, Silvers & Favier (2015).

The force term, F, in equation ( 2) aims to model external forces resulting from large-scale global effects (such as Reynolds stresses associated with thermal convection in global-scale calculations for example) that are not included in our local approach. In order to balance the viscous dissipation associated with the initial shear flow profile given by equation ( 9), the force

F = -σC k ∇ 2 U 0 (10)
was included in equation (2). For the set of equations ( 1)-(3) with the viscous forcing given by equation (10) the basic state described above with the initial velocity profile is only an equilibrium state if viscous heating, the first term on the right-hand side of equation (3), is neglected. For all the cases that we consider, the time-scale of the shear instability is at least two orders of magnitude smaller than the viscous heating time-scale of the system. Thus there is a time-scale separation between the growth rate of the instability and the viscous evolution of the shear flow, so that assuming we have an equilibrium is reasonable. This method has been broadly applied to model forced shear flows and the dynamics of the solar tachocline (e.g. Miesch 2003;Silvers et al. 2009). Note that this method only balances the viscous diffusion of momentum associated with the target profile and does not depend on the actual non-linear solution. It is always the case that any kind of forcing method has an effect on the saturated regime. Moreover, most forcing methods reach a state where the injected energy and the dissipated energy are in balance. We chose a method that ensures a balance and provides a local forcing that has minimal effect on the changes of the background profile induced by the instability. Thus it is suitable to study the characteristics of turbulent motion that are triggered by an instability in a viscous fluid. For a more detailed discussion of the appropriate forcing method see Witzke, Silvers & Favier (2016).

Our calculations were initialised by adding a small random temperature perturbation to the equilibrium state including the additional shear flow in equation ( 9). In order to evolve the system in time, equations ( 1)-(3) were solved using a hybrid finitedifference/pseudo-spectral code (see Matthews, Proctor & Weiss 1995;Silvers, Bushby & Proctor 2009;Favier & Bushby 2012, 2013). In addition to conducting fully non-linear direct numerical calculations, we also considered the linear stability analysis of this system, as detailed in Witzke, Silvers & Favier (2015). For this the eigenvalue-problem was numerically solved on a one-dimensional grid in the z-direction that is discretised uniformly, and this method is adapted from Favier et al. (2012).

For the characterisation of the initial state it is convenient to introduce several dimensionless numbers. Applying our nondimensionalisation on the Brunt-Väisälä frequency, as derived in Andrews (2000, p. 33), it becomes

N 2 (z) = θ (m + 1) T pot ∂T pot ∂z , (11)
where T pot = T P 1/γ-1 is the potential temperature. Then, the minimum value of the Richardson number, Ri, across the layer is defined as

Ri min = min 0 z 1       N(z) 2 ∂u 0 (z) ∂z 2       = min 0 z 1          θ 2 L 2 u (m + 1) m+1 γ -m (1 + θz) U 0 -4u 0 (z) 2 /U 0 2          , (12)
where the derivative of the background velocity profile, defined in equation ( 9), with respect to z corresponds to a local turnover rate of the shear. In most cases the minimum Ri value is at z = 0.5, but for some parameter choices where there is a large temperature gradient, θ, and a broad shear width, the minimum is shifted towards greater z. For selected cases the Ri(z) was plotted further below in Fig. 5 and Fig. 6. The 1/4 criterion is a necessary, but not sufficient, requirement for instability in an incompressible fluid. Thus, in order to verify that the cases we considered are unstable, the linear stability problem was solved. We considered unstable shear flows with a Richardson number, Ri, less than 1/4 at a point in the domain. However, for some investigations, systems with a minimum Richardson numbers greater than 1/4 were considered in order to study the 'secular' instability. Furthermore, to characterise the system we used the initial Péclet number at the top of the domain, i.e. z = 0, which we define as

Pe = U 0 L u C k ρ(0) , (13)
where U 0 and L u are as defined in equation ( 9). The Péclet number is also useful for the so-called 'secular' instabilities. In Subsection 3.4 we considered shear instabilities induced by the destabilizing effect of thermal diffusion for which larger values of Ri can be used (Dudis 1974;[START_REF] Zahn | Union I. A[END_REF]Lignières et al. 1999). Finally, the initial Reynolds number, defined at the top of the domain, is

Re = U 0 L u σC k ρ(0). ( 14
)
Our choice of fundamental units is useful when a fully compressible, polytropic atmosphere is studied, and naturally differs from previous Boussinesq studies (see for example Jones 1977;Lignières et al. 1999;Peltier & Caulfield 2003). In the following we discuss the link between the free parameters that appear in the equations ( 1)-( 3) and the Reynolds number, the Richardson number and the Péclet number, which are typically used when incompressible shear flows are studied. Varying the Prandtl number, σ, corresponds to a change in the Reynolds number only. The Péclet number can be associated with the thermal diffusivity parameter, C k . However, in order to vary only the Péclet number, it is necessary to keep the viscosity fixed, i.e. σ needs to be adjusted accordingly. Finally, the Richardson number can be varied independently from Re and Pe by changing the temperature gradient, θ, or the polytropic index, m. To characterise the initial setup the initial Pe number, Re number and minimum Ri number were calculated, where the Pe and Re numbers were obtained at the top but using the typical length-scale from the shear profile. Since the domain is strongly stratified, and the initial configuration is a polytropic state, the pressure scale-height changes across the domain. The maximum and minimum pressure scale-heights at the bottom and top layer respectively were calculated and a rounded value is given for all the considered cases in Table 1 and2.

RESULTS

When focusing on small-scale dynamics such as shear-induced turbulence, we would like to understand whether there is a relation between the transport coefficients, such as viscosity and thermal diffusivity, and the characteristic length-scales and the velocities of the resulting turbulent saturated state. Furthermore, observations of relevant shear flows, such as the tachocline for example, only provide spatial and time averaged measurements (see for example Kosovichev 1996;Charbonneau et al. 1999). Therefore, understanding what controls the global properties of the resulting mean flow after saturation, can help to draw a connection between observations of shear flows in astrophysical objects and numerical calculations. For this we investigated the effect of varying the values of the transport coefficients on non-linear dynamics. Note, that the viscosity value is a limiting factor for the smallest scales that can be achieved numerically in the turbulent regime. We focus on varying Re in Subsection 3.1, where we also introduce all calculated quantities that are investigated for all studied cases. Then, the Pe number was varied in Subsection 3.2 and we investigate how varying Ri affects the saturated phase in Subsection 3.3. These investigations were conducted in two dimensions, but some representative three-dimensional cases were performed and compared to two-dimensional cases in Subsection 3.5. All cases that were considered are summarised in Table 1. Finally, in Subsection 3.4 we discuss the so-called 'secular' shear instability (Endal & Sofia 1978) in a two-dimensional setup. Here, we investigated if the resulting state in the non-linear regime shows significantly different behaviour compared to that of a low Péclet number regime triggered by a 'classical' shear flow instability. Parameters for the 'secular' instability cases are summarised in Table 2.

The effect of varying the Reynolds number

In order to understand the effect of viscosity on shear induced turbulence, investigations in a two-dimensional domain with spatial resolution 480 × 512 were performed. This study is summarised in Table 1. Here, a highly supercritical system, i.e. Ri 1/4, was chosen, where Ri = 0.006. The Mach number remains less than 0.1 throughout the domain, which was achieved by taking U 0 /2 = 0.095, 1/L u = 60 and considering a temperature difference of θ = 1.9. The viscosity was varied by several orders of magnitude while the thermal diffusivity was fixed in such a way that the Péclet number is greater than unity (cases A to C in Table 1).

In general, before the system reaches a quasi-static state the shear flow instability grows exponentially and evolves throughout the saturation phase. After the system saturates, it enters a regime where statistical quantities fluctuate around a mean value. These regimes can be identified from the time evolution of the volume averaged vertical root-mean-square (rms) velocity in two dimensions

w = 1 N x N z Nx i=1 Nz j=1 w(i, j) 2 . ( 15
)
Since the overall evolution of the volume averaged vertical velocity is similar in all unstable systems, the time evolution of w for cases A to C is shown in Figure 1. We chose these three cases because they have significantly different viscosities and illustrate the effect of viscosity on the saturation. In Fig. 1 (a) the exponential growth of the instability is displayed, which is almost identical for the three cases because only the viscosity was changed. During the exponential growth phase, the instability growth rate does not vary for the cases in Fig. 1. After approximately nine sound crossing times the system starts to saturate. This phase persists longer for lower viscosities, as can be seen in Fig. 1 (b). Finally, when the volume averaged vertical velocity fluctuates around a mean value the system reaches a statistically steady state. Case A enters the statistically steady state after 150 sound crossing times, but for case C the statistically steady state is approximately after 450 sound crossing times. All statistics presented in this paper were time-averaged over a sufficiently long time interval during the statistical steady state. In order to ensure that the system is evolved for a sufficiently long time we considered the largest diffusive time-scale t σC k ,0 = 1/(σC k ), where 1 is the non-dimensionalised length of the domain. Our calculations were evolved for at least a significant fraction of this timescale, > 0.2 times the largest diffusive time-scale, which is long enough to give meaningful statistics. Note, the statistically steady state in this investigation is, to some extent, affected by the forcing method we used. The forcing method we chose here reaches a state where the work done by the forcing on the system and the viscous dissipation rate of momentum in the system are in balance (Witzke, Silvers & Favier 2016).

We begin with a qualitative observation of the flow, by considering the vorticity component perpendicular to the x-z-plane just after the exponential growth phase when the billows start to overturn, and during the statistical steady state. The dynamics alter significantly with decreasing viscosity. In Fig. 2 snapshots for cases A and C are displayed: As the viscosity decreases, vorticity structures are generated on much smaller spatial scales, as expected. However, it also becomes evident that the height of the horizontal layer in which mixing occurs changes only slightly as the viscosity is decreased. In order to estimate the extent of the effective shear region the horizontally averaged velocity in x-direction was calculated as

(z) = 1 N x Nx i=1 u x (i, z), (16)
where the overbar denotes that the quantity u x is horizontally averaged, and N x is the resolution in x-direction. Then, the effective shear width, L e f f , was obtained by fitting the function

f (z) = U e f f 2 tanh 2 L e f f (z -0.5) ( 17
)
to the resulting time averaged ūx (z). Since our aim is to approximate the spread, for some cases the middle of the domain was excluded from the fit. To obtain the fit we applied a non-linear least squares method using a trust region algorithm (Moré, J. J. & Sorensen, D. C. 1983) to find U e f f and L e f f . As the goodness of the fit achieves a very small root-mean-squared error of order 10 -3 for all cases, we used the 95% confidence boundaries for the obtained coefficients to estimate the error in L e f f . Note that, we also tested a cut-off method as an alternative approach, where we determined the region in which the averaged velocity drops below 95% of the maximum value of the velocity amplitude, L cut . The results are summarised in Table 1 and are qualitatively the same.

The averaged velocity profiles, as calculated in equation ( 16), are displayed in Fig. 3, where asymmetries can occur due to the fact that we considered a strongly stratified system. For all of the cases considered the confinement of the effective shear region is due to the dynamics, where the boundaries have only a negligible effect on the form of the horizontally averaged profiles. In Fig. 3 (a) ūx is shown for cases A to C and shows that all of the cases are very similar to each other despite the difference of three orders of magnitude in viscosity. Interestingly, for cases with smaller diffusivities, the form of the averaged velocity profile developes a 'staircase like' profile. However, upon checking the corresponding averaged density profiles, we determined that no 'staircase like' behaviour is present.

The effective shear width, L e f f , for these cases, summarised in Table 1, confirms the previous observation that the vertical extent of the region, where mixing occurs, is barely affected by viscosity:

The effective width decreases only slightly, as viscosity is changed over two orders of magnitude.

During the evolution of an unstable shear flow the horizontally averaged profiles for density, temperature and velocity can be modified. Therefore, the effective minimal Ri number, which we define as

min Ri e f f = min

          -θ (m + 1) γ ∂ū(z) ∂z 2 γ -1 ρ(z) ∂ ρ(z) ∂z - 1 T (z) ∂ T (z) ∂z           , (18)
where the overbar denotes horizontally averaged quantities, Understanding the relevant parameters affecting the turbulent characteristics, can provide a comprehensive picture of the possible dynamics in stellar interiors. For the characteristic properties of the turbulent regime the root-mean-square velocity of the perturbations and the typical turbulent length-scales were calculated. The systems that we were considering are stratified such that most quantities will change with depth, z, throughout the domain. Therefore, investigating horizontally-averaged profiles varying with depth, be-fore averaging over depth, provides further insight in the dynamics during the saturated regime. The horizontal turbulent length-scale of the overall velocity can be defined as

l t (z) = 2π E(k x , z)/k x dk x E(k x , z) dk x , (19)
where k x is the horizontal wave number. The corresponding energy spectrum E(k x , z) takes the form

E(k x , z) = 1 4 u(k x , z) • ρu * (k x , z) + ρu(k x , z) • u * (k x , z) , (20)
where denotes the Fourier transform and the * is used to indicate the complex conjugate. In previous studies the vertical scale of the vertical motion was calculated (Garaud & Kulenthirarajah 2016).

Due to inherently inhomogeneous nature of the system considered, it is impossible to calculate exactly the same quantity. However, to obtain a comparable quantity we calculated the typical horizontal scale of the vertical motion by taking only the vertical velocity into account in equation ( 20), such that the energy spectrum E w (k x , z) was obtained. The corresponding turbulent length-scale is then given by

l w (z) = 2π E w (k x , z)/k x dk x E w (k x , z) dk x . (21)
However, the resulting typical length-scales, l w , are always smaller than the overall turbulent length-scale, l t , but show the same trends with varied Re, Pe and Ri in our investigations. The root-meansquare of the fluctuating velocity u rms (z) we calculated as

u rms (z) = Nx x=1 (u(x, z) -ūx (z)) 2 /N x , (22)
where ūx (z) is the horizontally averaged velocity in x-direction as defined in equation ( 9). Here, we averaged over the horizontal layers after the root-mean-square velocity was obtained at each position. The u rms reveals that the actual turbulent region after saturation is more confined as indicated by effective shear width. When the instability is spread the effective shear region is enlarged, but during the saturated regime the outer layers of this region become no longer turbulent. Thus the effective shear region is an upper bound for the turbulent region. In addition, we calculated a local turbulent Reynolds number

Re t (z) = ρ(z)l t (z)u rms (z)/(σC k ), (23)
where ρ(z) is the horizontally averaged density. Since our domain it inhomogeneous in the vertical direction the turbulent Reynolds number varies across the domain.

It can be seen in Table 1 that Re t increases with decreasing σ, as expected. However, the minimum l w decreases from case A to case B but is found to increase slightly for case C. This indicates that if viscosity is further decreased then the smallest typical lengthscales present at the middle of the domain might converge towards a certain value. The u rms increases with increasing viscosity and the maximum u rms at the middle of the domain as well. Therefore, we conclude that varying Re by changing the Prandtl number does not lead to significant changes of the global characteristics, but affects the typical length-scales of the turbulence as expected.

Varying Péclet numbers

Here the main focus is to investigate different turbulent systems that have different Péclet numbers but where both Re and Ri are fixed. It is important to distinguish between two Péclet number limits: In the large Péclet number limit Pe 1, which means that the typical time-scale on which advection occurs is shorter than the time-scale on which thermal diffusion acts. The small Péclet number limit starts around Péclet number of order unity, where the time-scales for advection and diffusion are of the same order, and it continues for all Pe < 1. Moreover, thermal diffusion becomes important in systems where the thermal diffusion time-scale is shorter than the buoyancy time-scale, which is the case in a system with Pe < 1. Thus in the small Péclet number limit thermal diffusion weakens the stable stratification, i.e. the system becomes less constrained by buoyancy such that vertical transport is enhanced. We will focus our discussion here on five cases with different C k and σ, where C k was chosen so that the initial Péclet number increases by four orders of magnitude over the cases we considered. In order to examine how the non-linear dynamics change with varying the Péclet number we started with a qualitative comparison of the smallest and largest Péclet numbers considered for cases D to H. A visualisation of the vorticity after the system has saturated is shown in Fig. 4. For case D, displayed in Fig. 4 (a), with a Péclet number of order 10 -2 , very strong positive vorticity in a narrow region around the middle plane is present. Here, patches are stretched along the x-axis with a few small interruptions of negative vorticity. Positive vorticity regions are stretched outwards from the middle plane at z = 0.5 and are overturning. A significantly different pattern is present in Fig. 4 (b), case H, where the Péclet number is of order 10 2 . Here the vorticity amplitude is notably less than in the other case and a vertically extended turbulent region is present. The turbulent region is more isotropic with positive and negative vorticity patches. Furthermore, smaller scale vortices are present in case H compared to case D. From this we conclude that case D, which is in the small Péclet number regime, leads to a different turbulent state, where larger fluid parcels are present compared to case H, which is in the large Péclet number regime. This indicates a complex effect of Pe number on the length-scales present in the turbulent regime. Since the Péclet number depends on the typical length-scale, for perturbations on sufficiently short lenght-scales the small Péclet number regime is reached. Thus for such perturbations the stabilising effect of the stratification is weaker, which was first noted by [START_REF] Zahn | Union I. A[END_REF], and they can develop. So when decreasing C k the length-scales that are destabilised become smaller and thus smaller fluid parcels are present. Since the Pe number is varied significantly in cases D to H, another quantity, the temperature fluctuations around the initial background temperature, δT , is of interest. A visualisation of δT for the cases D and H is shown in Fig. 4 (c) and Fig. 4 (d) respectively. The absence of small scale fluctuations in Fig. 4 (c) is a natural consequence of a greater thermal diffusivity. Note, when increasing the Pe number, from D to H, the extent of the turbulent layer after saturation looks visually larger. This is confirmed by the effective shear width, L e f f in Table 1. Since the dynamical viscosity is fixed, these effects result solely from different Pe. Therefore, the Péclet number, which is associated with the thermal diffusion, plays an important role in the non-linear dynamics and in particular on the vertical spread of the shear induced turbulence. We find that in the limit of large Péclet numbers the effective spread increases with increasing Péclet numbers as the effective shear width, L e f f , becomes larger for cases F to H as shown in Table 1. This shows that a decreased Pe number significantly damps the spread of perturbations. In the large Pe regime the temperature of a fluid parcel will adjust faster to the surrounding temperature when Pe is decreased. Therefore, part of the kinetic energy is irreversibly converted into internal energy quicker and so the further propagation of the fluid parcels is hindered. However, in the limit of small Péclet numbers the opposite trend is observed, where the effective shear width decreases with increasing Pe. This suggests a complex interplay between the effect of thermal diffusivity and energy contained in the system, which we investigate in the next subsection. We compared the typical turbulent length-scale, l w , and the rootmean-square of the velocity perturbations, u rms , for different initial Péclet numbers (cases D to H in Table 1). These cases reveal that the Pe number significantly affects the turbulence scale l w , where the smallest l w is found for the case with Pe = 60. Furthermore, the maximum root-mean-square of the velocity perturbations (see Table 1) reduces with increasing thermal diffusion in the limit of large Péclet numbers. This trend is consistent with the recent Boussinesq investigations of the low Péclet number regime (Garaud & Kulenthirarajah 2016) where the Richardson number was varied. However, in our setup the observed turbulent length-scales always exceed the minimal pressure scale-height, H min , by at least a factor of two, whereas for some cases the maximum pressure scale-height, H max is slightly larger than the turbulent length-scale. So that the Boussinesq approximation is not valid for our system.

Varying the Richardson number

To understand how the energy provided by the background flow will affect the saturated regime, the Richardson number was varied independently of all other parameters. This was achieved by changing θ separately for a low Pe number (cases D1, D2, D3 ), an intermediate Pe number (cases FF1, FF2, FF3), and for a large Pe number (cases G1, G2, G3), which are summarised in Table 1. For Pe = 0.04 and Pe = 6 we changed the Richardson number by one order of magnitude via varying the thermal stratification between θ = 0.8 and θ = 6. Fig. 5 shows the Ri across the middle of the domain for the cases D1 to D3. The region where this number remains less than 1/4 decreases slightly with increasing Ri. As a result the effective shear width significantly decreases with increasing Richardson number for the cases D1 to D3 and FF1 to FF3. The Richardson number is proportional to the ratio of the potential energy that is needed to overcome the stabilising stratification and the available kinetic energy of the background flow. Since a perturbation loses its initial energy when moving vertically, in the ideal case it will continue to spread as long as it has more energy than needed to overcome the stratification. Therefore, the vertical extent of the turbulent region increases with decreasing Richardson numbers. However, in the limit of large Pe the opposite is observed, and so the effective shear widths increases with increasing Ri, which indicates a more complex effect. Therefore, we conclude that there are two parameters affecting the extent of the effective shear region generated by an unstable shear flow in a stratified system. The first is the Richardson number, since it provides information on the ratio between available kinetic energy and the potential energy. The second parameter that controls the extent of the effective shear region is the Péclet number, which affects the vertical motion of fluid parcels.

The product of the Richardson number and the Péclet number, RiPe, has been used as an input parameter in Boussinesq calculations (Prat & Lignières 2014;Garaud & Kulenthirarajah 2016) and so it is natural to consider if this product can be used to quantify the system in these compressible calculations. Therefore, we now focus on comparing the typical length-scale and root-mean-square of the velocities obtained in the cases D1 to G3 with results obtained in Garaud & Kulenthirarajah (2016). For the small Pe regime the decrease of the typical length-scale with increasing RiPe, is recovered. However, the root-mean-square velocities in our calculations decrease with increasing RiPe, whereas in Garaud & Kulenthirarajah (2016) they increase with increasing RiPe. Interestingly, in the large Pe regime the typical length-scale increases with RiPe, if the Péclet number is increased, but decreases if the Ri number is increased. This result suggests that in stratified systems, and also in the large Pe regime the product of Ri and Pe can not be used to characterise the saturated dynamics of the system. It is necessary to consider each dimensionless number separately.

Secular Instability

In the previous section we focused only on classical shear flow instabilities that can be present within the stellar interiors where Table 2. For the 'secular' instability cases O to R the temperature gradient, θ = 2.0, and the polytropic index, m = 3.0, are fixed, but the initial Ri number changes from 0.4 for cases O to Q to 1.0 for case U. For all cases the pressure scale-height is the same, where H max = 0.4 and H min = 0.15. The initial Péclet number,Pe, Reynolds number, Re, and minimum Richardson number, Ri, are listed. The effective shear width is calculated after saturation and Re t , and l w are averaged over the whole domain. Additionally, we provide a shear width, L cut obtained by a cut-off method at 95%.

Secular Instability

Resolution N x = 512, N z = 480

Case: σ C k U 0 L u Pe Re Ri L e f f L cut Re t ūrms max u rms l w min(l w )
O 0.001 0.05 0.1 0.029 0.06 5.7 × 10 1 0.4 0.076 ± 0.002 0.20 4.8 × 10 1 8.5 × 10 -04 6.7 × 10 -3 0.58 0.20 P 0.0006 0.05 0.1 0.029 0.06 9.5 × 10 1 0.4 0.094 ± 0.002 0.24 1.6 × 10 2 1.1 × 10 -03 7.4 × 10 -3 1.16 0.26 Q 0.0002 0.05 0.1 0.029 0.06 2.9 × 10 2 0.4 0.097 ± 0.001 0.25 8.1 × 10 2 1.2 × 10 -03 6.5 × 10 -3 1.12 0.36 R 0.0002 0.05 0.07 0.033 0.05 2.3 × 10 2 1.0 0.063 ± 0.001 0.17 7.8 × 10 2 9.8 × 10 -04 6.3 × 10 -3 1.4 0.82 Figure 6. Richardson number with depth around the middle of the domain for the secular instability cases and the most of the cases from Table 1.

Richardson numbers become smaller than 1/4. However, another well-known type of shear instability can occur even when the Richardson number is significantly greater than 1/4 (Dudis 1974;[START_REF] Zahn | Union I. A[END_REF]Garaud et al. 2015), but only in the low Péclet number regime.

Low Péclet numbers are most likely to be present in upper regions of massive stars, where differential rotation is present and 'secular' shear instabilities can develop. Previous investigations have considered diffusive systems either by using the Boussinesq approximation (Lignières et al. 1999;Prat & Lignières 2013;Garaud et al. 2015) or restrict the study to linear stability analysis when employing the fully compressible equations (Witzke et al. 2015). In order to investigate the evolution that occurs during the saturated phase, non-linear studies are required. However, it is extremely difficult to consider small Péclet number regimes using the full set of equations as the time-stepping is restricted by stability constraints related to the diffusion time. Therefore, it has previously been impossible to numerically investigate small Péclet number shear flows without using approximations (as in Prat & Lignières 2013). However, we show that it is possible to calculate low Péclet and large Richardson number cases in a two-dimensional domain. Here we present some calculations, with a spatial resolution of 480×512 to investigate the differences between the classical and 'secular' instabilities during the saturated phase.

To ensure that the classical KH instability is not triggered, we set Ri = 0.4, which was achieved by taking θ = 2, m = 3.0, U 0 = 0.1 and 1/L u = 70. Taking C k = 0.05 leads to an initial Péclet number of order 10 -2 , such that a secular instability can develop due to the destabilising effect of low Péclet numbers. Fig. 6 shows how the Ri number changes in the middle of the domain comparing the profile to the small Ri = 0.006 cases investigated above. For cases O, P and Q the Reynolds number was varied via the Prandtl number in order to investigate its effect on the turbulent length-scale. In order to study how different Ri affect the dynamics of a secular instability, for case R the Richardson number was increased to Ri = 1.0. This was achieved by changing U 0 = 0.07 and 1/L u = 60, and the dynamical viscosity was fixed to be 10 -5 . We compared the growth rate and the most unstable mode predicted by means of a linear stability analysis as used in Witzke et al. (2015) to that found in all cases of the secular instability. Furthermore, we checked that the instability is a consequence of the destabilising mechanism at low Péclet numbers, by conducting test cases with the same dynamical viscosity as used in cases O to R but C k = 0.0002. This result in an initial Pe ≈ 12 and we find that for both Ri numbers the system remains stable and the initial perturbations decay.

When the instability starts to saturate in any of the cases O to R we observe very little overturning, which is unlikely to develop into a turbulent regime at least for the Reynolds numbers considered here (see Fig. 7). For all cases the vertical spread of the flow is significantly smaller than observed in the previous cases (see Fig. 3 (c) and Table 2). This is because the Richardson numbers are very large and the available kinetic energy for the perturbations is decreased.

We find for all cases here that the vertical spread of the temperature perturbations is slightly greater than it is for the vertical velocity during the saturation phase, which can be seen from Fig. , where the Péclet number is of the same order, we see that for the classical instability no layering occurs. Turning to the characteristic length-scale during the saturated regime we find that, from the data in Table 2, decreasing the viscosity results in greater typical turbulent length-scales. The opposite trend was observed in the previous study of an unstable system at large Péclet number. The trend for the 'secular' instability cases can be explained by the peculiar pattern observed for the vertical velocity and temperature fluctuations as shown in Fig. 7. Since the fluctuations are sheared out the typical length-scale between the up and down motions is increased. When the viscosity is decreased it becomes easier for the horizontal movement of the background flow to elongate the fluctuation pattern even more, such that the typical length increases. However, such a pattern does not transit to developed turbulence for the cases considered. In a test case with Ri = 0.1, a similar layering is observed when the Péclet number dropped below unity. However, case D, which has a Pe = 0.04 and a Ri = 0.006, does not show any similar layering. Therefore, we conclude that such layering can develop in systems that are not very far from the stability threshold. Such a behaviour is not only a feature of the secular instability, but can be present in a system with Richardson number close to the stability threshold, but with a sufficiently small Péclet number.

To summarize, we find that the 'secular' shear flow instability in a fully compressible, stratified fluid, shows the expected trend for the spread of the instability. Although the cases studied do not become turbulent, a turbulent regime will be eventually reached when con-sidering larger Reynolds numbers. We observe a difference in the trend for typical length-scales during the saturated regime, when compared to cases where a classical KH instability was triggered. Moreover, a peculiar pattern for the temperature fluctuations and vertical velocity fluctuations are found to be present.

3.5 Saturated regime using three-dimensional calculations Two-dimensional calculations may alter the non-linear dynamics from what could occur in three dimensions, since vortex stretching is suppressed. Therefore, it is crucial to investigate two representative cases in three dimensions to confirm the non-linear dynamics.

Therefore, here we will discuss cases where the full set of dimensionless, three-dimensional, differential equations is used. An additional horizontal y-direction prependicular to the z-x-plane is considered, which is periodic and its length is normalised by the depth d.

We chose two different Péclet number regimes, because both small and large Péclet number regimes can occur in stellar interiors depending on the region considered and the type of star. Here, case I represents the large Péclet number regime, Pe > 1 and case J the low Péclet number regime with Pe < 1. The calculations were evolved over a sufficient fraction of the largest thermal diffusion time-scale in the system, which is given by the domain depth squared divided by the thermal diffusivity, t thermal = 1/C k , where 1 is the non-dimensional depth of the domain. The low Pe case was even evolved for several thermal times, t thermal . The spatial resolution for these calculations is N x = 256, N y = 256 and N z = 320. The spatial extent of the horizontal dimensions needs to be taken larger for case J, because in a smaller box the secondary instability that propagates in the y-direction is suppressed. The exact parameters for these cases are summarized in Table 1, where the corresponding two-dimensional cases K and L have the same parameters and resolution as their corresponding three-dimensional case. Note, that this comparison study was performed at a greater viscosity, σC k , than all other cases due to computational cost.

Here, we compared the characteristics of the three-dimensional calculations to two-dimensional calculations. Furthermore, direct comparisons of global properties obtained in the saturated regime for both the small and large Péclet number regime calculations were performed. Fig. 8 shows contour plots of the vertical velocity, w, for case I and case J. These plots were taken a long time after the system has saturated for both cases. Visually the three-dimensional calculations show similar dynamics: The turbulent regions have approximately the same vertical extent, as well as the overturning regions, for both considered cases. This is explained by the fact that the initial Péclet numbers of the two cases are close to the threshold between small and large Péclet number regime, such that the dynamics in both regimes are similar. Similar dynamics are also observed in two-dimensional calculations for the small Péclet number regime. However, here in the three-dimensional calculations, secondary instabilities lead to turbulent motions in the y-direction. Fig. 8 clearly shows a well developed turbulence in both horizontal directions. When comparing the three-dimensional cases to their corresponding two-dimensional cases K and L, the effective shear width L e f f is similar for the large Pe calculations. However, in the small Pe regime the three-dimensional case has a reduced shear width L e f f = 0.40 compared to the corresponding two-dimensional case, where it is L e f f = 0.48. Although the values differ between the twoand three-dimensional cases, it is important to note that a similar increase in the spread for larger Pe persists. The discrepancies in the confinement in three-dimensional calculations can be explained by the different non-linear dynamics where, for example, a secondary instability can evolve perpendicular to the x-direction. In Fig. 9 (a) the horizontally averaged velocity in the x-direction is shown for all four cases. While for the low Pe number cases J and L a 'staircase like' profile is found, the velocity profile in the large Pe regime is best described as a hyperbolic tangent profile.

There is no qualitative difference in the profiles for the two-and three-dimensional calculations. We now compare the change in the typical length-scale, lw (see Table 1), from the large Pe number case to the low Pe number case obtained in three-dimensional calculations with the corresponding change in two-dimensional calculations. The typical length-scale increases in both two-and threedimensional calculations. In Fig. 9 (b) the typical length-scale, l w , is plotted with depth to illustrate the change from the middle of the domain to the boundaries. However, the values of all quantities obtained during the saturated phase, e.g. the root-mean-square velocity and typical length-scales, are different for the three-dimensional calculations compared to the corresponding two-dimensional calculations. While there are significant differences between all cases, the form of the profile around the middle of the domain is similar. A similar shift towards lower values in the large Pe regime is found for both two-and three-dimensional calculations. The differences towards the upper and lower boundaries for the two-and three-dimensional calculations indicate that there are different dynamics, which is expected due to the existence of secondary instabilities leading to a solution varying in the y-direction in threedimensional calculations. Finally, comparing the changes in the turbulent Re t , and u rms from the low Pe regime to the large Pe regime in the two-dimensional calculations with the changes obtained in the three-dimensional calculations, summarised in Table 1, we find that the trends are the same. Therefore, we conclude that while there are inevitable differences in the detailed turbulent characteristics obtained in threedimensional calculations, the overall effect obtained by varying the Pe number is qualitatively the same as in two dimensions.

CONCLUSIONS

In order to obtain a comprehensive understanding of stars in their entirety, it is crucial to understand the complex dynamics in stellar interiors by investigating shear driven turbulence. Shear driven turbulence is a promising candidate in order to explain the missing mixing problem and is important for magnetic field generation. Numerical calculations are used to obtain a comprehensive insight to the detailed small scale dynamics of shear regions. However, due to computational limitations, all calculations to date use modelling parameters that are far from the actual values in stellar interiors. Therefore, it is important that we consider how varying key properties affects our understanding of complex stellar regions.

In our study we focused on understanding the effect of different viscosities and thermal diffusivities on the saturated phase of a shear driven turbulent flow in a fully compressible polytropic atmosphere. Examining the global properties of the saturated flow of an unstable system revealed that the vertical extent of the mixing region is primarily controlled by the Richardson number, but the Péclet number also plays a key role through the time-scale on which thermal diffusivity acts on the system. For greater Richardson numbers we find that the vertical spread of the mixing decreases, which also occurs as Pe is decreased in the large Péclet number limit. In the small Péclet number limit, i.e. for Pe < 1, an opposite trend is observed, where the vertical spread increases with decreasing Pe. This increase is due to the high thermal diffusivity that weakens the effectively stratification as soon as the Pe is less than unity. This weakening effect occurs only in the small Péclet number limit. We find that viscosity does not play an important role in the formation of the global shape of the mean flow during the saturated regime. Turbulent flows can be characterised by the typical length-scale and the root-mean-square velocity of the perturbations. We showed that the typical turbulent length-scales depend on the Richardson number as well as on the Péclet number regime. Investigating strongly stratified systems, and systems in the large Péclet number regime, we find that there is a different behaviour in the turbulent characteristics depending on whether it is the Richardson number or the Péclet number that is increased. While in both cases the product RiPe was increased by the same amount, the system responds differently. Therefore, we conclude that the product of the input Péclet number and the Richardson number, RiPe, can not be used in such strongly stratified systems to extract information on the characteristics of the turbulence, and both dimensionless number should be provided. In summation, the turbulent regime of a shear flow instability depends on several parameters and these can counteract each other. The properties of the saturated regime can only be broadly predicted from the input parameters.

In the latter part of our research we focused on the low Péclet number regime. While for large Péclet numbers the initial flow requires low Ri numbers to become unstable, for low Péclet numbers it is possible to destabilize a high Richardson number shear flow (Lignières et al. 1999;Witzke et al. 2015). We examined cases where there were unstable secular shear instabilities. We found in these cases that a different dependency exists, where the typical length-scale increases with decreasing viscosity, which is not the case in large Péclet number regimes.

Having established a better understanding of what parameters significantly affect the global properties of saturated shear flow instabilities, future studies of the mixing behaviour and momentum transport of shear-driven turbulence can now be conducted. In order to gain a more comprehensive picture of the complex dynamics in stellar interiors, it is crucial to include magnetic field interactions. Future investigations with magnetic fields will help to inform how magnetic fields affect the turbulent regime. Moreover, the results obtained can be used to seek for a shear induced turbulence that is capable to drive a magnetic dynamo, which is subject to current investigations.

Figure 1 .

 1 Figure1. Volume averaged vertical velocity evolution for cases A to C. In a) the early evolution, where the exponential growth is shown, whereas in b) the long-time evolution is displayed.

Figure 2 .

 2 Figure2. Vorticity in the x-z-plane for case A and C (see Table1), where the viscosity is decreasing. At the top (a) and (b) show case A and C during saturation both at t ≈ 93 and at the bottom the same cases are displayed during the quasi-steady state (c) at t ≈ 255 and (d) at t ≈ 500 (for reference see Fig.1). For all cases the thermal diffusivity parameter is C k = 10 -4 . The dotted lines indicate the extent of the turbulent region of the saturated state as obtained from equation (17).

  ūx

Figure 3 .

 3 Figure 3. The horizontally averaged and time averaged ūx profiles are shown for cases A, B, C, D, F and H.

Figure 4 .

 4 Figure 4. The vorticity component perpendicular to the x-z-plane (top panel) and temperature fluctuations around the initial temperature profile (lower panel) for two cases shortly after the system has saturated. Case D is shown in (a) and (c) and case H is displayed in (b) and (d). The dotted lines in a) and b) indicate L e f f and the red lines in b) show the large error for L e f f , where the lower bound is indicated.

Figure 5 .

 5 Figure 5. Ri with depth for the cases D1, D2 and D3.

  7 (b) and Fig. 7 (c). A similar trend appears for a classical shear instability in the low Péclet number regime as can be seen in Fig 4. While there is some similarity in the trend in terms of the spread, there is a marked difference in the patterns observed in the saturated state between the secular and the classical instabilities. Fig. 7 (b) shows that regions of upward and downward motion are stretched along the horizontal direction. Similar pattern of the negative and positive temperature fluctuations is present in Fig. 7 (c), where layers are formed. Comparing what is seen in Fig. 7 (c) to the temperature fluctuations present in case D (see Fig. 2 (d))

Figure 7 .

 7 Figure 7. The vorticity component perpendicular to the x-z-plane, the vertical velocity, w, and temperature fluctuations around the initial temperature profile for case Q long after saturation. In (a) the vorticity is shown and the dotted line indicates the extent of the effective shear, L e f f , in (b) the vertical velocity is shown, and (c) shows the temperature fluctuations.

Figure 8 .

 8 Figure 8. The isosurfaces of constant vertical velocity, w, for four different velocity values in three-dimensional calculations. (a) case I at roughly t ≈ 200. (b) case J at t ≈ 110.

Figure 9 .

 9 Figure 9. Turbulence characteristics. (a) Horizontally averaged and time averaged ūx profiles. (b) Typical turbulent length-scale l w .

Table 1 .

 1 A comparison of typical length-scales, turbulent Reynolds numbers, Re t , and effective shear width, L e f f , during the saturated phase. The error for the effective shear width only accounts the fitting error. Additionally, we provide a shear width, L cut , obtained by a cut-off method at 95%. For all cases the polytropic index is m = 1.6, the shear amplitude U 0 = 0.19, and the initial shear width is L u = 0.0333. The initial Péclet number, Pe, Reynolds number, Re and minimum Richardson number, Ri, are listed. The effective shear width is calculated after saturation and Re t , and l w are averaged over the whole domain. For the cases A to H and I to L the pressure scale-heights are the same, where at the top H min = 0.2 and at the bottom of the domain H max = 0.6. For the cases D1 to G3, the maximum pressure scale-height varies from H max = 0.4 to H max = 1.2 and for the minimum pressure scale-height from H min = 0.06 to H min = 0.8.

	Case	θ	σ	C k	Pe	Re	Ri	L e f f	L cut	Re t	ūrms	max u rms	l w	min l w
			Varying Re via changing σ				Resolution N x = 512, N z = 480
	A	1.9	1	0.0001	60	6.3 × 10 1	0.006 0.61 ± 0.04 0.88 2.6 × 10 2 1.2 × 10 -2 2.4 × 10 -2 0.93	0.43
	B	1.9	0.1	0.0001	60	6.3 × 10 2	0.006 0.53 ± 0.05 0.78 5.7 × 10 2 2.5 × 10 -3 8.8 × 10 -3 0.89	0.21
	C	1.9	0.05	0.0001	60	1.3 × 10 3	0.006 0.49 ± 0.03 0.68 7.3 × 10 2 1.6 × 10 -3 6.2 × 10 -3 0.91	0.23
		Varying Pe number via changing C k				Resolution N x = 512, N z = 480
	D	1.9 0.000633	0.16	0.04 6.2 × 10 1	0.006 0.51 ± 0.02 0.89 3.2 × 10 2 8.9 × 10 -3 2.0 × 10 -2 1.04	0.60
	E	1.9	0.00633	0.016	0.4 6.2 × 10 1	0.006 0.43 ± 0.01 0.86 1.2 × 10 2 5.3 × 10 -3 2.0 × 10 -2 1.05	0.33
	F	1.9	0.0633	0.0016	4.0 6.2 × 10 1	0.006 0.43 ± 0.02 0.71 6.7 × 10 1 4.3 × 10 -3 1.9 × 10 -2 0.91	0.29
	G	1.9	0.948	0.00011	60	6.2 × 10 1	0.006 0.65 ± 0.07 0.90 3.8 × 10 1 4.0 × 10 -3 2.1 × 10 -2 0.46	0.18
	H	1.9	9.48	0.000011 600 6.2 × 10 1	0.006 0.95 ± 0.16 0.83 2.4 × 10 2 1.0 × 10 -2 2.4 × 10 -2 0.79	0.39
		Varying Ri number via changing θ				Resolution N x = 512, N z = 480
	D1	1.9 0.000633	0.16	0.04 6.2 × 10 1	0.006 0.51 ± 0.02 0.89 3.2 × 10 2 8.9 × 10 -3 2.0 × 10 -2 1.04	0.60
	D2	4.2 0.000633	0.16	0.04 6.2 × 10 1	0.018 0.47 ± 0.02 0.86 2.0 × 10 2 4.0 × 10 -3 1.6 × 10 -2 1.12	0.42
	D3	6.3 0.000633	0.16	0.04 6.2 × 10 1	0.030 0.33 ± 0.01 0.80 1.5 × 10 2 2.7 × 10 -3 1.3 × 10 -2 0.85	0.31
	FF1	0.85	0.0633	0.0016	6	6.2 × 10 1 0.0006 0.87 ± 0.09 0.89 1.7 × 10 2 1.1 × 10 -2 3.2 × 10 -2 1.23	0.53
	FF2	3.6	0.0633	0.0016	6	6.2 × 10 1	0.006 0.45 ± 0.01 0.76 1.5 × 10 2 4.7 × 10 -3 2.1 × 10 -2 1.07	0.26
	FF3	6.5	0.0633	0.0016	6	6.2 × 10 1	0.013 0.23 ± 0.01 0.67 3.1 × 10 2 5.0 × 10 -3 1.6 × 10 -2 0.68	0.32
	G1	0.49	0.948	0.00011	60	6.2 × 10 1 0.0006 0.36 ± 0.01 0.88 6.1 × 10 2 3.0 × 10 -2 4.2 × 10 -2 1.48	1.26
	G2	1.2	0.948	0.00011	60	6.2 × 10 1 0.0030 0.51 ± 0.04 0.82 2.6 × 10 2 1.3 × 10 -2 2.7 × 10 -2 1.06	0.60
	G3	1.9	0.948	0.00011	60	6.2 × 10 1 0.0060 0.65 ± 0.07 0.90 3.8 × 10 1 4.0 × 10 -3 2.1 × 10 -2 0.46	0.18
			Three-dimensional cases				Resolution N x = 256, N y = 256, N z = 320
	I	1.9	0.32	0.0016	4.0 1.2 × 10 1	0.006 0.29 ± 0.01 0.66 2.6 × 10 1 7.7 × 10 -3 3.5 × 10 -2 0.74	0.36
	J	1.9	0.016	0.032	0.2 1.2 × 10 1	0.006 0.48 ± 0.01 0.70 3.7 × 10 1 1.0 × 10 -2 4.0 × 10 -2 0.90	0.41
		Two-dimensional cases for comparison				Resolution N x = 256, N z = 320
	K (2D) 1.9	0.32	0.0016	4.0 1.2 × 10 1	0.006 0.30 ± 0.02 0.63 5.0 × 10 1 1.2 × 10 -2 3.4 × 10 -2 0.97	0.41
	L (2D)	1.9	0.016	0.032	0.2 1.2 × 10 1	0.006 0.40 ± 0.02 0.71 6.6 × 10 1 1.1 × 10 -2 3.5 × 10 -2 0.99	0.51

ACKNOWLEDGEMENTS

This research has received funding from STFC and from the School of Mathematics, Computer Science and Engineering at City, University of London. This work used the DiRAC Data Analytic system at the University of Cambridge, operated by the University of Cambridge High Performance Computing Service on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National Einfrastructure capital grant (ST/K001590/1), STFC capital grants ST/H008861/1 and ST/H00887X/1, and STFC DiRAC Operations grant ST/K00333X/1. DiRAC is part of the National E-Infrastructure.