© ESO, 2016

1. Introduction

Stellar evolution plays a central role in astrophysics by linking observable quantities to fundamental parameters of stars, such as mass, radius, and age. These parameters also provide constraints on galaxies, which are made of stars, and on planetary systems, which orbit around them. In addition, stellar models are used to calibrate many measurements for other domains of astrophysics, such as cosmological distance measurements based on models of Cepheids and type-Ia supernovae.

Stellar evolution theory involves many physical processes with very different lengths and occuring at very different time scales. In particular, rotationally induced transport processes, as well as overshoot and various magneto-hydrodynamical (MHD) instabilities, often refered to as “extra mixing”, have a significant impact on the distribution of angular momentum and chemical elements, and thus on the structure and the evolution of stars. The predictive power of stellar evolution codes is limited by large uncertainties related to our poor knowledge of these (magneto-)hydrodynamical processes and by the fact that, being fundamentally 3D, they cannot be included self-consistently in 1D computations. Some evolution codes include rotation (Meynet et al. 2013) or other 3D processes, such as overshoot (Moravveji et al. 2015), semi-convection (Ding & Li 2014), thermohaline mixing (Wachlin et al. 2014), or the magneto-rotational instability (Wheeler et al. 2015). However, these processes, including convection itself, are described by phenomenological models which involve free parameters.

Another problem of stellar evolution is the fact that all the information we have about stars comes from their surfaces (except for the Sun where we also have information from neutrinos) and is generally not spatially resolved. This implies that there is no straightforward way to constrain the transport generated by hydrodynamical processes inside stars. However, we do have some constraints, mainly from spectroscopy and asteroseismology.

First, the analysis of absorption lines in stellar spectra gives access to surface chemical abundances of stars. Some of these abundances are strongly affected by mixing processes occurring in certain regions of stars. For instance, elements produced (helium and nitrogen) or depleted (carbon) in the core during the CNO cycle are good indicators of mixing between the core and the surface (see for example Hunter et al. 2009; Martins et al. 2015). One major issue is that measured abundances depend on the initial chemical composition of observed stars, which is usually not known.

In contrast, asteroseismology allows us to probe the interiors of stars by analysing their oscillation spectra. Indeed, regular spectral patterns can be related to the stellar structure, including the rotation profile. These techniques have been used so far mainly for the Sun and a few red giants and subgiants (Beck et al. 2012; Mosser et al. 2012; Deheuvels et al. 2012, 2014, 2015), but a great quantity of new data from other stars is now available thanks to the CoRoT and Kepler satellites. This is clearly a very promising source of additional constraints on transport processes (see for example Moravveji et al. 2015).

The so-called shear mixing is one of these processes. It occurs when differential rotation is strong enough to overcome the effect of stable stratification in a radiative zone. Stellar radiative zones are characterised by a very efficient radiative thermal diffusion, which tends to weaken the stabilising effect of thermal stratification. This has been formalised by Zahn (1992) using phenomenological arguments. Since then, several related models including additional physics such as chemical stratification or horizontal turbulence have been proposed and implemented in stellar evolution codes (Maeder & Meynet 1996; Talon & Zahn 1997). More rigourous models based on mean-field equations exist. However, these models usually express the mixing rate as a function of quantities which are not available in stellar evolution codes (e.g. Lindborg & Brethouwer 2008), and additional assumptions such as closure models (e.g. Canuto 1998, 2011) are thus required, introducing uncertainties. Although it would be interesting to test such more complex models with numerical simulations, we focus here on Zahn’s model, which is commonly used in stellar evolution calculations.

An increasing number of studies suggest that the transport predicted by Zahn’s model (or related ones) is insufficient to fit some obervational constraints. For example, synthetic stellar populations computed by Brott et al. (2011) and Potter et al. (2012) cannot reproduce several groups of stars observed by Hunter et al. (2009), notably nitrogen-rich slow rotators. Note that the existence of these stars might be explained by rotational mixing if some braking is active at the surface, for instance due to a surface magnetic field (Meynet et al. 2011). Other studies show that additional mixing is needed to explain observed rotation profiles of red giants or subgiants (Eggenberger et al. 2012; Ceillier et al. 2012; Marques et al. 2013). Moreover, even adding other transport processes is still not sufficient to explain the observations (Cantiello et al. 2014, for the Tayler-Spruit dynamo; Fuller et al. 2014; Belkacem et al. 2015, for internal gravity waves). Constraints from the spin of neutron stars and white dwarfs at birth also indicate that the transport as predicted by Zahn’s model is not sufficient (Heger et al. 2005, for neutron stars; Suijs et al. 2008, for white dwarfs). Additional constraints are needed to determine whether these discrepancies point to an incorrect description of shear mixing or to the presence of another transport process.

Because of the degeneracy of observational constraints and the impossibility of performing laboratory experiments in the stellar regime, numerical simulations are a second way to obtain new constraints on specific transport processes. The purpose of this work is to test Zahn’s family of models for shear mixing with local direct numerical simulations, and potentially to propose new prescriptions including new physical ingredients. Whereas in the previous papers (Prat & Lignières 2013, 2014, thereafter refered to as PL13 and Paper I, respectively) we focused on thermal diffusion and the dynamical effect of chemical stratification, the present paper is intended to test the robustness of the constraints on Zahn’s model that one can obtain from shear mixing simulations. We achieve this by comparing the results of simulations performed with two different numerical codes and different radial boundary conditions, and by studying the effect of viscosity and the size of the numerical domain.

The paper is organised as follows. We summarise the model of shear mixing proposed by Zahn (1992) in Sect. 2, and we describe the two numerical codes used in this paper in Sect. 3. The results of our simulations are presented in Sect. 4. Finally, we discuss the robustness of numerical constraints on Zahn’s model in Sect. 5.

2. Model by Zahn (1992)

In the laminar, adiabatic, non-viscous case, the classical Richardson criterion for linear shear instability reads (Miles 1961) $\mathit{Ri}\mathit{<}\frac{\mathrm{1}}{\mathrm{4}}\mathit{,}$(1)where Ri = N²/S² is the Richardson number comparing the stabilising effect of thermal stratification (given by the Brunt-Väisälä frequency N) and the destabilising effect of differential rotation (given by the shear rate S = dU/ dz of a velocity field U(z) varying along the spatial coordinate z). Typical Richardson number values in stellar radiative zones are between 10² and 10⁴. Following Townsend (1958), Zahn (1974) introduced the effect of thermal diffusion on thermal stratification. The main assumptions underlying Zahn’s model are: (i) in the regime of very high thermal diffusivities, the Richardson criterion is replaced by $\mathit{RiP}{\mathit{e}}_{\mathit{\ell }}\mathit{<}{\mathit{R}}_{\mathrm{crit}}\mathit{,}$(2)where Pe_ℓ = wℓ/κ is the turbulent Péclet number (which can be seen as the ratio between the diffusive cooling time and the dynamical time of turbulence), w and ℓ turbulent velocity and length scales, κ the thermal diffusivity of the fluid, and R_crit the critical value of RiPe_ℓ, of the order of one; (ii) turbulent flows tend to reach a statistical steady state which is marginally stable (i.e. RiPe_ℓ = R_crit); (iii) the turbulent diffusion coefficient D_t is proportional to wℓ; and (iv) the turbulent Reynolds number (i.e. the ratio between the viscous time and the dynamical time of turbulence, which typically is between 10² and 10⁵ in stellar radiative zones) Re_ℓ = wℓ/ν, where ν is the viscosity of the fluid, is larger than a critical value Re_crit of the order of 10³. These assumptions, especially the first one, which arises from the hypothesis that a linear treatment can provide some approximate treatment of turbulent flows, are not rigourously justified, and hence need to be tested with numerical simulations such as those presented in this work.

From the first two assumptions, we can deduce $\mathit{P}{\mathit{e}}_{\mathit{\ell }}\mathrm{=}\frac{{\mathit{R}}_{\mathrm{crit}}}{\mathit{Ri}}\mathit{,}$(3)and by definition of the turbulent Péclet number, we have $\mathit{w\ell }\mathrm{=}\frac{\mathit{\kappa }{\mathit{R}}_{\mathrm{crit}}}{\mathit{Ri}}\mathrm{·}$(4)Finally, the third assumption gives us ${\mathit{D}}_{\mathrm{t}}\mathrm{\propto }\frac{\mathit{\kappa }}{\mathit{Ri}}\mathrm{·}$(5)It is important to note that in the asymptotic regime described by Zahn (1992), the turbulent diffusion coefficient does not depend on viscosity. However, the existence of the fourth assumption suggests that when it is not verified, the turbulent diffusion coefficient may depend on the value of the turbulent Reynolds number. In the asymptotic regime, from Eq. (4), we get $\mathit{R}{\mathit{e}}_{\mathit{\ell }}\mathrm{\propto }\mathrm{(}\mathit{RiPr}{\mathrm{)}}^{-1}\mathit{,}$(6)where Pr = ν/κ is the Prandtl number (comparing the relative effects of viscosity and thermal diffusion), which is typically of the order of 10^-6 or smaller in stellar interiors. This suggests that RiPr is a proxy for the turbulent Reynolds number. We can thus describe the dependence of the turbulent diffusion coefficient on the turbulent Reynolds number as a dependence on RiPr. In particular, one can check how the quantity D_t/ (κRi^-1), which is constant in Zahn’s model according to Eq. (5), varies with RiPr. This will be done in Sect. 4.3. Such a dependence is easier to implement in stellar evolution codes, where RiPr can be computed, but not Re_ℓ.

3. Physical and numerical configurations

The simulations presented in this paper have been performed with two different codes. The first one is the Balaïtous code, that we used in PL13 and Paper I. This corresponds to what we call the forcing configuration. The second one is the Snoopy code (Lesur & Longaretti 2005), associated with the shearing-box configuration.

The physical configurations used in this paper are essentially the same, as in PL13 and Paper I, and the two codes also have some numerical aspects in common, which are listed in Sect. 3.1. The forcing configuration is described in Sect. 3.2 and we present the specificities of the shearing-box configuration in Sect. 3.3.

3.1. Common aspects

All our simulations are performed using the Boussinesq approximation, which neglects density fluctuations except in the buoyancy term, and the so-called small-Péclet-number approximation (SPNA, see Lignières 1999), which allows us to explore the regime of very high thermal diffusivities at a reasonable computational cost. Local mean gradients are characterised by a uniform velocity shear, and uniform temperature and chemical composition gradients in the vertical/radial direction, as described in Fig. 1.

Fig. 1 Sketch of the flow configuration.

We choose L, S^-1, ΔU = SL, ΔT = LdT/ dz, Δc = LdC/ dz, and ρS²L² as length, time, velocity, temperature, concentration, and pressure units (where L is the size of the numerical domain, T(z) and C(z) are the background temperature and concentration profiles, and ρ is the background density). In this framework, the dimensionless Navier-Stokes equations read $\begin{array}{ccc}& & {\nabla }\mathrm{·}{u}\mathrm{=}\mathrm{0}\mathit{,}\\ & & \\ & & {\mathit{u}}_{\mathit{z}}\mathrm{=}\mathrm{\Delta }\mathit{\psi ,}\end{array}$with ψ = θ/Pe, where u, P, and θ are the dimensionless velocity, pressure, and temperature fluctuations around the background profile, and Pe = SL²/κ is the Péclet number, assumed to be very small. These equations highlight the two main control parameters: (i) the product of the Richardson and Péclet numbers RiPe = N²L²/ (Sκ) (where N² = αgΔT/L, α being the thermal expansion coefficient of the fluid and g the local gravity), also called the Richardson-Péclet number, which characterises the effect of stratification affected by thermal diffusion; and (ii) the Reynolds number Re = SL²/ν, which characterises the effect of viscosity. These control parameters are clearly dependent on the size of the numerical domain L. One can define their turbulent counterparts RiPe_ℓ = N²wℓ/ (S²κ) and Re_ℓ = wℓ/κ, based on the turbulent velocity and length scales, which are diagnostic variables. The only relevant control parameter that does not depend on the size of the numerical domain is the product of the Richardson and Prandtl numbers RiPr = RiPe/Re = N²ν/ (S²κ).

To these equations, we add the equation for passive advection/diffusion of the chemical concentration: $\frac{\mathit{\partial c}}{\mathit{\partial t}}\mathrm{+}{u}\mathrm{·}{\nabla }\mathit{c}\mathrm{=}\frac{\mathrm{1}}{\mathit{P}{\mathit{e}}_{\mathrm{c}}}\mathrm{\Delta }\mathit{c,}$(10)which introduces a new control parameter, the chemical Péclet number Pe_c = SL²/D_m, which characterises the effect of the chemical molecular diffusity D_m.

Numerically, both codes use a Fourier colocation method associated with periodic boundary conditions in the horizontal (x- and y-) directions. They also both allow us to perform direct numerical simulations (DNS), in which all the physically relevant scales of the problem are resolved, typically down to the dissipation scale. The main differences between the codes lie in the numerical method and the boundary conditions used in the vertical (z-)direction, and in the way background profiles (of shear, temperature, and composition) are imposed.

3.2. Forcing configuration

The forcing configuration uses compact finite differences in the vertical direction, associated with impenetrable boundary conditions and imposed shear, temperature and chemical composition at the upper and lower boundaries. Due to these boundary conditions, the flow is not rigourously statistically homogeneous in the vertical direction close to the upper and lower boundaries, as mentioned in PL13 and Paper I.

The linear background profiles are imposed by forcing terms f_v and f_c added to Eqs. (8)and (10), respectively. These terms are defined as ${{f}}_{\mathrm{v}}\mathrm{=}\frac{{U}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{-}\overline{{u}}}{\mathit{\tau }} \mathrm{and} {\mathit{f}}_{\mathrm{c}}\mathrm{=}\frac{\mathit{C}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{-}\overline{\mathit{c}}}{\mathit{\tau }}\mathit{,}$(11)where U(z) and C(z) are the background velocity and composition profiles, the overline denoting the horizontal average, and τ is the forcing time scale. Normally a similar term should be added to the temperature equation, but in the SPNA, only infinitesimally small temperature fluctuations are considered.

Tests have shown that, provided that the forcing time scale τ is much smaller than a few hundred shear times, its exact value does not affect the results. As a consequence of the forcing method, the horizontally-averaged profiles (here $\overline{{u}}$ and $\overline{\mathit{c}}$) remain very close to the background profiles (here U and C).

3.3. Shearing-box configuration

The shearing-box configuration uses a Fourier colocation method also in the vertical direction, associated with shear-periodic boundary conditions at the upper and lower boundaries. This is made possible by the fact that the equations are solved for the fluctuations of all quantities around their background profiles and not for the quantities themselves. The resulting flow is rigourously statistically homogeneous in all directions, in constrast with the forcing configuration.

There is no explicit forcing in the shearing-box configuration, but the fact that fluctuations are shear-periodic in the vertical direction ensures that the spatial averages of the shear and of the temperature and composition gradients are equal to their background values. Equivalently, it means that the horizontal averages of velocity, temperature and composition differences between the upper and lower boundaries are fixed. Consequently, horizontally averaged gradients are not necessarily uniform at a given time and mean profiles are able to evolve under the action of the instability.

4. Results of the simulations

We performed five series of simulations at constant chemical Péclet number Pe_c = 4 × 10⁴, with different values of RiPe: one with Balaïtous at Re = 4 × 10⁴, and four with Snoopy with Re ranging from 2 × 10⁴ to 1.6 × 10⁵. In each simulation, we estimate the turbulent diffusion coefficient as ${\mathit{D}}_{\mathrm{t}}\mathrm{=}\mathrm{-}\frac{\mathrm{⟨}{\mathit{c}}^{\mathrm{\prime }}{\mathit{u}}_{\mathit{z}}\mathrm{⟩}}{\mathrm{d}\mathit{C}\mathit{/}\mathrm{d}\mathit{z}}\mathit{,}$(12)where ⟨⟩ denotes the temporal and spatial average, assumed to be equivalent to the statistical average, c′ is the fluctuation of chemical composition about its background profile, and u_z is the vertical component of the velocity. One can also estimate the turbulent viscosity coefficient, also known as eddy viscosity, as ${\mathit{\nu }}_{\mathrm{t}}\mathrm{=}\mathrm{-}\frac{\mathrm{⟨}{\mathit{u}}_{\mathit{y}}^{\mathrm{\prime }}{\mathit{u}}_{\mathit{z}}\mathrm{⟩}}{\mathrm{d}\mathit{U}\mathit{/}\mathrm{d}\mathit{z}}\mathit{,}$(13)where ${\mathit{u}}_{\mathit{y}}^{\mathrm{\prime }}$ is the fluctuation of the y-component of the velocity about the mean shear profile. The SPNA simulation described in PL13 is also considered here for comparison. Table 1 summarises the relevant quantities of our simulations.

In this section, we investigate four different aspects of our simulations. First, we compare the results given by the two codes in Sect. 4.1. Second, we analyse the effect of the limited size of the numerical domain in Sect. 4.2. Third, we study the dependence of the turbulent transport on the turbulent Reynolds number in Sect. 4.3. Finally, we constrain the relation between the coefficients of turbulent diffusion and viscosity in Sect. 4.4.

4.1. Comparison of the two codes

To compare the behaviour of the two codes in a similar configuration, we consider here only simulations at Re = 4 × 10⁴ (indicated in Table 1 by a superscript d). Figure 2 shows the quantity D_t/ (κRi^-1) as a function of RiPr for the two codes.

Fig. 2 Dt/ (κRi-1) as a function of RiPr for simulations with Re = 4 × 104.

One can clearly see that the two codes give very similar results. For values of RiPr where both codes have been tested, the relative difference between the two estimated transport coefficients never exceeds a few percent. For values of RiPr where only one code has been tested, our results generally follow the same curve. This indicates that both codes are consistent for the study of turbulent transport in the regime of small Péclet numbers. We thus validate the shearing-box formalism for the study of the shear instability in this regime. This is an encouraging result, since the Snoopy code is more efficient and scalable than Balaïtous, mainly because of the 3D fast Fourier transform.

One may notice in Table 1 that for the simulations considered in this section, the values of ℓ/L, Re_ℓ, RiPe_ℓ, and D_t/ (wℓ) significantly differ between the codes. This results at least partly from the fact that the length ℓ is differently defined in the two codes. In Balaïtous, ℓ is the integral scale of turbulence based on horizontal 2D spectra, whereas in Snoopy it is the integral scale based on 3D spectra. As shear flows are expected to be anisotropic, it is not surprising that these two definitions give different results. To summarise, we can write $\mathit{\ell }\mathrm{=}\mathrm{2}\mathit{\pi }\frac{{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathit{E}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{/}\mathit{k}\hspace{0.17em}\mathrm{d}\mathit{k}}{{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathit{E}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{d}\mathit{k}}\mathit{,}$(14)where k is the norm of the wave vector for Snoopy, the norm of the horizontal wave vector for Balaïtous, and E(k) the kinetic energy spectrum (2D horizontal for Balaïtous and 3D for Snoopy). In both cases the turbulent velocity scale w is defined as $\sqrt{\mathrm{2}\mathit{K}}$, where K denotes the turbulent kinetic energy per unit mass.

4.2. Size of the numerical domain

According to Eq. (5), the quantity D_t/ (κRi^-1) is constant in Zahn’s model. This is clearly not what is observed in Fig. 2. Either, one of Zahn’s assumptions is not verified, or this is a numerical issue. In this section, we study the effect of the size of the box on the transport. In Sect. 4.3, we focus on the dependence of the transport on the turbulent Reynolds number.

The quantity D_t/ (κRi^-1) is plotted in Fig. 3 again as a function of RiPr, but this time for all simulations performed with Snoopy.

Fig. 3 Dt/ (κRi-1) as a function of RiPr obtained with our Snoopy simulations. The dotted line corresponds to the empirical formula given in Eq. (15). There are two overlapping symbols at RiPr = 5 × 10-4.

On the one hand, for values of RiPr down to 5 × 10^-4, there is no significant difference between the simulations at different Reynolds numbers, so that D_t/ (κRi^-1) can be considered as a function of RiPr only. On the other hand, for smaller values, a large dispersion in the results appears, larger Reynolds numbers giving larger values of D_t/ (κRi^-1). As seen previously, RiPr is a proxy for the turbulent Reynolds number, while Re depends on the size of the numerical domain. Varying Re at constant RiPr is thus equivalent to studying the effect of the numerical domain size. It means that the results of our simulations depend on the size of the box when RiPr is smaller than 5 × 10^-4, but not for larger values.

Figure 4 shows the ratio ℓ/L between the integral scale of turbulence and the size of the numerical domain as a function of RiPr.

Fig. 4 ℓ/L as a function of RiPr obtained with our Snoopy simulations. The dashed line represents the approximate limit ℓ/L = 0.35 between simulation domains that are large enough (below) and too small (above).

For simulations with RiPr ≥ 5 × 10^-4, in which the measured transport does not depend on the size of the numerical domain, we observe that the scale ratio ℓ/L is always smaller than 0.33. Simulations at low RiPr that are clearly affected by the size of the numerical domain show a scale ratio larger than 0.4. This suggests that there is a limit in ℓ/L at about 0.35 between simulations which have a sufficiently large box size and those which have a too-small numerical domain. The latter show a less vigourous turbulence, which can be related to the fact that the size of the numerical domain limits turbulence development.

We now have a reliability criterion which roughly states that there should be at least three large turbulent structures in the simulation domain to be able to give size-independent results. Assuming that this criterion is valid, two consequences arise. Firstly, the two simulations at RiPr = 2.5 × 10^-4 (runs 22 and 30) should not depend on the box size. Although we can see in Fig. 3 that they have slightly different values of D_t/ (κRi^-1), this difference is significantly smaller than the one between simulations at RiPr = 10^-4 and not much larger than in the regime of larger RiPr, for example at RiPr = 10^-3. Secondly, the simulation with Re = 1.6 × 10⁵ and RiPr = 10^-4 (run 29) is very close to the limit between size-independent simulations and size-dependent ones. Therefore, this simulation cannot be completely trusted, but the fact that it gives a value of D_t/ (κRi^-1) similar to the one in the simulation with the same Reynolds number but RiPr = 2.5 × 10^-4 (run 30) is an argument in favour of its reliability.

To conclude, this criterion should be used in future simulations to ensure that the numerical domain is large enough. Thanks to this criterion, we can now study the dependence of D_t/ (κRi^-1) on the turbulent Reynolds number, which is detailed in the next section.

4.3. Turbulent Reynolds number

Now that we have determined which simulations are physically relevant, we can focus on the dependence on the turbulent Reynolds number. If we ignore unreliable simulations (indicated in Table 1 by a superscript e), Fig. 3 shows three main features: (i) an asymptotic regime at low RiPr, which is equivalent to high turbulent Reynolds numbers; (ii) a maximum around RiPr = 3 × 10^-3; and (iii) a very sharp drop around RiPr = 7 × 10^-3 (simulations above this value showed no turbulence). Only the asymptotic regime was described by Zahn (1992), under the assumption that the turbulent Reynolds number is larger than a critical value of the order of 10³. The two other features are consequences of the transition from quasi-laminar flows dominated by viscosity to fully turbulent flows. In particular, the lower limit in Re_ℓ corresponds to the situation where the largest scale that is unstable with respect to the shear instability is too small to overcome the effect of viscosity. The fact that this limit can be associated with a critical RiPr has been predicted by Garaud et al. (2015) using a non-linear stability analysis.

At the moment, we have no physical model for the value of the turbulent diffusion coefficient outside of the asymptotic regime. However, we can propose the following empirical formula to account for the main features of Fig. 3: $\frac{{\mathit{D}}_{\mathrm{t}}}{\mathit{\kappa R}{\mathit{i}}^{-1}}\mathrm{=}\mathit{\alpha }\mathrm{+}\mathit{\beta RiPr}\mathrm{-}\mathit{\gamma }\mathrm{(}\mathit{RiPr}{\mathrm{)}}^{\mathrm{2}}\mathit{,}$(15)where α = 3.34 × 10^-2, β = 18.8, and γ = 2.86 × 10³. In the asymptotic regime (i.e. for small values of RiPr), we thus have D_t/ (κRi^-1) ≃ α. The value we find here is almost twice as low as the value 5.58 × 10^-2 found in PL13 because the asymptotic regime was not reached in that paper. The physical meaning of β/γ = 6.57 × 10^-3 is the following: for RiPr>β/γ, D_t/ (κRi^-1) drops below the asymptotic value. Since the curve is quite steep, β/γ would effectively give the order of magnitude of the critical RiPr above which there is no turbulence.

In the asymptotic regime, our new prescription predicts about 60% less transport than the model by Zahn (1992), and about 15 times less than the model by Maeder (1995). Outside of the asymptotic regime, the discrepancy between our prescription and those by Zahn and Maeder can be reduced at best by a factor two. Even then, our prescription predicts still significantly less transport than the model by Maeder (1995).

In the following figures, we kept only the relevant simulations. The relation between Re_ℓ and RiPr (see Eq. (6)) is tested in Fig. 5.

Fig. 5 Reℓ as a function of RiPr for all valid Snoopy simulations. The dashed line and the dash-dotted line represent the power laws Reℓ ∝ (RiPr)-1 and Reℓ ∝ (RiPr)− 2 / 3, respectively.

We observe that Eq. (6)is verified only when RiPr< 5 × 10^-4, that is close to the asymptotic regime (dashed line in Fig. 5). This coincides quite well with the condition Re_ℓ> 10³ given by Zahn (1992). Even at larger values of RiPr, the relation between Re_ℓ and RiPr is always a one-to-one relationship.

The fact that Eq. (6)is verified for low values of RiPr implies that RiPe_ℓ is almost constant in this regime. This is illustrated in Fig. 6.

Fig. 6 RiPeℓ as a function of RiPr for all valid Snoopy simulations.

We can see that the value of RiPe_ℓ is close to 0.65 in the considered regime. This verifies directly the assumption of Zahn’s model that in the asymptotic regime flows are characterised by a constant value of the turbulent Richardson-Péclet number. The approximate value found here is about 50% larger than 0.426, the value found in PL13. As mentioned earlier, this may result from the fact that the integral scale ℓ was somewhat differently defined in previous studies (PL13 and Paper I) compared to the present one, while the definition of w remains the same. We observe that the values of ℓ are larger in the present paper, and so are Re_ℓ and RiPe_ℓ. This may also be partly due to the fact that the SPNA simulation in PL13 was not in the asymptotic regime.

Another of Zahn’s assumptions can be tested with our simulations: the fact that the turbulent diffusion coefficient is proportional to wℓ. We observe in Fig. 7 that the ratio D_t/ (wℓ) is close to 5.5 × 10^-2 when RiPr ≤ 2.5 × 10^-4, then increases up to 6.5 × 10^-2 and finally drops to much smaller values, the value in the asymptotic regime being significantly smaller than 0.131, the one found in PL13.

Fig. 7 Dt/ (wℓ) as a function of RiPr for all valid Snoopy simulations.

Again, this is probably due to the different definitions of ℓ. We can also see that the assumption that D_t/ (wℓ) is constant is more restrictive than the previous one, since it seems to be true only for smaller values of RiPr.

4.4. Relation between turbulent diffusion and viscosity coefficients

Many stellar evolution codes consider that the diffusion coefficients for chemical elements and for angular momentum are the same, even though Zahn (1992) mentioned that they could differ by a small factor. In the weakly turbulent regime (RiPr ≳ 3 × 10^-3), our simulations show that both coefficients are small, and have nearly equal values (D_t ~ ν_t ≲ 10ν). In the fully turbulent regime, which is relevant for stellar interiors, we find that they follow the relation ${\mathit{\nu }}_{\mathrm{t}}\mathrm{\simeq }\mathrm{0.80}{\mathit{D}}_{\mathrm{t}}\mathit{,}$(16)as illustrated in Fig. 8.

Fig. 8 Scaling between νt and Dt for all Snoopy simulations. The dotted line represents the proportionality relation νt = 0.80Dt.

5. Discussion

We have shown that the two codes we have tested give similar results in the regime of small Péclet numbers. This provides a validation of the shearing-box formalism for the study of the shear instability in this regime. Our work suggests further that we do not always have a clear scale-separation between the turbulent dissipation scale, the largest shear-unstable scale, and the size of the computational domain. As a consequence, one must be cautious in the interpretation of the simulation results. In particular, we showed that it is crucial to have a numerical domain at least three times as large as the largest shear-unstable scale to ensure that the results are relevant. Our results also suggest that the asymptotic regime described by Zahn (1992) is reached when the turbulent Reynolds number is larger than 10³. This is in agreement with the fourth of Zahn’s assumptions mentioned in Sect. 2. Such high values of the turbulent Reynolds number were not reached in our previous studies (PL13 and Paper I), which were consequently not in the asymptotic regime.

The observed dependence of the turbulent diffusion on the turbulent Reynolds number (or on RiPr, which is easier to calculate in stellar evolution codes) is not trivial. The empirical prescription we propose in this paper can be implemented in stellar evolution codes and compared with other prescriptions. In addition, we have shown that it is necessary to distinguish between the turbulent viscosity coefficient used for the transport of angular momentum and the turbulent diffusion coefficient used for the mixing of chemical elements. In practice, the proportionality constant between them is close to one, which means its impact on stellar evolution calculations is expected to be minor. The fact that our new prescriptions predicts even less transport than in our previous studies does not reduce the disagreement between observations and stellar evolution models.

The validation of Snoopy in this context is an interesting result, since it is able to simulate configurations in which the velocity and entropy gradients are not aligned. Such configurations would allow us to study the horizontal transport generated by some horizontal shear and the effect of this horizontal transport on the vertical transport (see model by Talon & Zahn 1997). The horizontal transport, believed to weaken the effect of the stable stratification and thus enhance the vertical transport, is indeed poorly understood. Several prescriptions exist (see for example Zahn 1992; Maeder 2003; Mathis et al. 2004), but all are based on phenomenological arguments and have never been tested. Constraints from numerical simulations are therefore needed. One must keep in mind that the similarities in the behaviour of the two codes are not necessarily valid in the regime of large Péclet numbers. In this regime, the boundary conditions and the forcing scheme probably play an important role, since in the classical Richardson stability criterion, Ri > Ri_crit, no preferential physical scale is involved. To give reliable prescriptions of turbulent transport in this regime, which may occur in evolved stars (Hirschi et al. 2004), the impact of the differences between the two configurations must be investigated in detail.

To explore further the asymptotic regime, we must go to higher Reynolds numbers, which means better grid resolutions, but also larger box sizes, thus increasing drastically the computational cost. As mentioned in the introduction, the expected Reynolds number in stellar radiative zones is between 10² and 10⁵. Therefore, our simulations already overlap with the physically relevant regime for stellar interiors. Although it requires significantly more computational effort to reach Reynolds numbers of the order of 10⁵, such computations should be possible in the near future thanks to the continuous increase in available computational power. We note that this favourable situation contrasts with other problems in stellar hydrodynamics, for example global simulations of the Solar convective region, which are hampered by the very large Reynolds number (of the order of 10⁹) that characterises the problem, making fully resolved simulations unfeasible in the foreseeable future. For the shear problem, a cheaper alternative would be to perform large-eddy simulations (LES) in which small scales are modelled using prescriptions from DNS. But before doing so, one would first have to test this approach by comparing LES and DNS with the same control parameters.

For simplicity, the dynamical effect of chemical stratification was not taken into account in the present study, but the respective results found in Paper I are very likely to also depend on the turbulent Reynolds number. This will be the subject of a forthcoming publication. Finally, another potentially relevant ingredient is the magnetic field. In the presence of differential rotation, the magneto-rotational instability (MRI) may drive MHD turbulence and therefore provide further transport (Balbus & Hawley 1994; Spruit 1999; Menou et al. 2004; Wheeler et al. 2015; Jouve et al. 2015). The MRI in a stably stratified layer has so far mostly been described analytically, through linear analysis and phenomenological arguments. These results need to be checked with numerical simulations in a similar way to that in the present paper. Analogous simulations have been performed in the context of proto-neutron stars (Guilet & Müller 2015) and should be extended to the regime of small Péclet numbers that is relevant for stellar radiative zones.

Acknowledgments

V.P. thanks Geoffroy Lesur for his interesting ideas and his help in the use of the Snoopy code. V.P. and M.V. thank François Lignières for insightful discussions. V.P. acknowledges support from the European Research Council through ERC grant SPIRE 647383. J.G. acknowledges support from the Max-Planck-Princeton Center for Plasma Physics. This work is supported by the European Research Council through grant ERC-AdG No. 341157-COCO2CASA. The authors thank the referee Georges Meynet for his useful comments that contributed to the quality of the paper.

References

All Tables

All Figures

	Fig. 1 Sketch of the flow configuration.
In the text

	Fig. 2 Dt/ (κRi-1) as a function of RiPr for simulations with Re = 4 × 104.
In the text

	Fig. 3 Dt/ (κRi-1) as a function of RiPr obtained with our Snoopy simulations. The dotted line corresponds to the empirical formula given in Eq. (15). There are two overlapping symbols at RiPr = 5 × 10-4.
In the text

	Fig. 4 ℓ/L as a function of RiPr obtained with our Snoopy simulations. The dashed line represents the approximate limit ℓ/L = 0.35 between simulation domains that are large enough (below) and too small (above).
In the text

	Fig. 5 Reℓ as a function of RiPr for all valid Snoopy simulations. The dashed line and the dash-dotted line represent the power laws Reℓ ∝ (RiPr)-1 and Reℓ ∝ (RiPr)− 2 / 3, respectively.
In the text

	Fig. 6 RiPeℓ as a function of RiPr for all valid Snoopy simulations.
In the text

	Fig. 7 Dt/ (wℓ) as a function of RiPr for all valid Snoopy simulations.
In the text

	Fig. 8 Scaling between νt and Dt for all Snoopy simulations. The dotted line represents the proportionality relation νt = 0.80Dt.
In the text