© ESO, 2014

1. Introduction

Mixing of chemical elements in stellar interiors is a crucial feature in stellar evolution theory. Indeed, for main sequence stars, mixing can continuously draw hydrogen from outer layers to the core. For a given initial mass, stars experiencing mixing have a higher averaged mean molecular weight and thus appear brighter. The hydrogen-core-burning lifetime is increased as is the mass of the helium core. Larger helium cores modify the subsequent evolution. It is thus essential to have a good knowledge of transport processes to build reliable stellar models. Whereas microscopic processes, such as molecular diffusion, gravitational settling, and radiative acceleration, are relatively well known, the effects of macroscopic motions induced by rotation are still poorly understood. Improving such models may have a strong impact in many fields of astrophysics. This is particularly true for galaxy physics, because the properties of the stars in a galaxy provide information about the galaxy, and for planetary systems, where constraints on the stellar host help us constrain the mass and the radius of planets.

The main technique used to constrain mixing processes in stars is to determine chemical abundances on the surface of stars through measurements of the depth and the width of atomic absorption lines. These measurements may bring out some anomalies for certain stars, that is, over- or under-abundances as compared to the standard model of stellar evolution, which describes the evolution of a spherical, non-rotating star in which the only macroscopic transport process is convection. For massive stars, in which the CNO cycle is dominant, abundance anomalies of elements involved in this cycle, such as helium, carbon, or nitrogen, indicate that some deep mixing is occurring in the observed star (see e.g. Gies & Lambert 1992; Lyubimkov et al. 2000; and more recently Hunter et al. 2009 with the VLT-FLAMES survey). Martins et al. (2009, 2013) have even found evidence for the existence of stars in which mixing is so efficient that they are quasi-homogeneous.

Another sort of anomaly is linked to the fact that light elements (lithium, beryllium, and boron) are destroyed at high temperature. For massive stars, lithium and beryllium are totally depleted, whereas boron is only destroyed in the inner layers. If some mixing occurs between these layers and the outer ones, the surface abundance of boron is lower than otherwise (Venn et al. 2002; Mendel et al. 2006). For cooler stars, deep convection causes lithium and beryllium to be depleted. It follows that these elements can only be detected in a narrow range of stellar mass. Here again we find discrepancies between the standard model and the observed abundances, including for the Sun (e.g. the lithium dip, initially observed by Boesgaard & Tripicco 1986). A source of extra mixing, such as internal gravity waves (Talon & Charbonnel 2005), is thus needed.

Stellar seismology is another powerful way to constrain transport processes within stars. The analysis of stellar oscillations not only provides new global parameters that are useful for better determining fundamental parameters of stars, such as mass, radius, and age, but it also allows us to probe the internal properties of some stars, including the distribution of chemicals (for the Sun), via the speed of sound (Christensen-Dalsgaard et al. 1985; Antia & Basu 1994) and internal rotation. Beck et al. (2012) and Mosser et al. (2012) have notably determined the rotation rate in the core of red giants, and a rotation profile has even been obtained for the Sun (Thompson et al. 1996) and a subgiant (Deheuvels et al. 2012). Thanks to this promising technique, we are about to be given access both to the causes of the transport (rotation) and to its effects (chemical distribution). We still need to extend its use to a larger number of stars, which is in progress, especially for CoRoT and Kepler stars.

In many stellar evolution codes, rotationally induced mixing is taken into account thanks to a set of turbulent diffusion coefficients linked to various hydrodynamical instabilities triggered by differential rotation and meridional circulation. Among these instabilities, shear instability is thought to produce the most efficient mixing (Knobloch & Spruit 1982), sometimes called “shear mixing”. This instability normally occurs in a stably stratified sheared flow when the Richardson number Ri = (N/S)², which compares stable stratification (with the Brunt-Väisälä frequency N) and the shear (with the shear rate S = dU/ dz), is lower than a critical value Ri_c, which equals 1 / 4 in the inviscid linear stability analysis (Miles 1961). However, except in some limited layers of evolved stars (Hirschi et al. 2004), the Richardson number is generally much larger than one in stars, and this criterion would result in almost always shear-stable stellar interiors. In this high-Richardson-number or weak-shear regime, the very high thermal diffusivities found in stars are nevertheless able to reduce the amplitude of the buoyancy force through radiative losses and thus to overcome the stabilising effect of stratification (Townsend 1958; Zahn 1974).

We are interested here in one of the transport coefficients related to this instability, the radial diffusion coefficient due to radial differential rotation. It was initially derived by Zahn (1992) from the following phenomenological arguments.

• Turbulent flows generally tend to reach a statistical steady statethat is marginally stable with respect to the instability that hasgenerated them.

• The destabilising effect of the very high thermal diffusivity κ is such that the Richardson instability criterion should be replaced by RiPe<Ri_c, where Pe = UL/κ denotes the Péclet number based on velocity and length scales U and L.

• The relevant Péclet number to use in this criterion is the turbulent one Pe_ℓ = uℓ/κ, where u is the velocity scale and ℓ the length scale of turbulence.

• The turbulent diffusion coefficient is proportional to the product of these turbulent scales: D_t ≃ uℓ/ 3.

Zahn finally obtains a diffusion coefficient that reads ${\mathit{D}}_{\mathrm{t}}\mathrm{\simeq }\frac{\mathit{\kappa }}{\mathrm{3}}\frac{\mathit{R}{\mathit{i}}_{\mathrm{c}}}{\mathit{Ri}}\mathrm{·}$(1)Since, several attemps have been made to add various physical ingredients to this model, such as μ-gradients (Maeder & Meynet 1996; Maeder 1997) and horizontal diffusion (Talon & Zahn 1997).

As regards massive stars, the use of such models of rotational mixing in stellar evolution codes gives results closer to the observations than the standard model both for CNO products (Heger & Langer 2000; Meynet & Maeder 2000; Maeder & Meynet 2001) and for boron (Fliegner et al. 1996; Mendel et al. 2006; Frischknecht et al. 2010). Thanks to the progress in the computational efficiency of such codes, it is now possible to compute entire synthetic stellar populations from a given distribution of initial conditions that include chemical composition and angular momentum. As shown by Brott et al. (2011) and Potter et al. (2012), these synthetic populations can be compared to observed ones, such as those obtained by large surveys. In particular, the authors recover the main group of stars observed by Hunter et al. (2009).

However, some observational constraints are still not explained by these models. In particular, Hunter et al. (2009) are unable to explain the existence of nitrogen-enriched slow rotators which are actually observed. Besides, current models do not predict enough transport of angular momentum to recover the seismic constraints on internal rotation in red giants and subgiants (Eggenberger et al. 2012; Ceillier et al. 2012; Marques et al. 2013).

Most of these phenomenological models of turbulent transport have never really been tested. Indeed, the physical conditions present in stellar interiors are far too extreme to allow us to perform realistic laboratory experiments. Our approach is to use direct numerical simulations (DNS) as experiments to test the existing prescriptions. Local numerical simulations of homogeneous, sheared, and stably stratified turbulence (such as those performed in a geophysical context by Gerz et al. 1989; Holt et al. 1992; and Jacobitz et al. 1997) have the advantage of being as generic as possible while taking the essential ingredients (shear and stable stratification) into account. There has already been an attempt to test Zahn’s model by Brüggen & Hillebrandt (2001) but because their simulations used numerical viscosity and thermal diffusivity, which depend on the resolution, they were unable to study the effect of thermal diffusion.

In a previous paper (Prat & Lignières 2013, hereafter Paper I), we explored the dependence of the radial turbulent diffusion coefficient on thermal diffusity in the regime of the small Péclet numbers typical of stellar interiors. This was done by considering a parallel-plane flow configuration with forced, uniform, vertical stable stratification and velocity shear. Boundary conditions were periodic in the horizontal directions and impermeability, along with the mean shear, were imposed at the lower and upper boundaries. Since the very high thermal diffusivity introduces a huge gap between the diffusive time scale and the dynamical one, which drastically increases the computational cost of complete Boussinesq simulations, we used an asymptotic development of the Boussinesq equations called the small-Péclet-number approximation (SPNA) described by Lignières (1999). We found that our simulations are qualitatively in good agreement with Zahn’s prescription in this regime and we were able to give a quantitative estimate of the turbulent diffusion coefficient.

The purpose of the present paper is twofold. On the one hand, we extend the previous study to the domain of large Péclet numbers, where Zahn’s model is not valid any more. Although the net effect on global mixing is often assumed to be negligible, some strong-shear situations exist, especially in evolved stars (see for example Hirschi et al. 2004), where the Richardson number is small enough that stellar interiors are unstable with respect to the shear instability even on large scales, thus at large Péclet numbers. On the other hand, we explore the dependence of the turbulent diffusion coefficient with respect to chemical stratification. This effect is likely to be significant at the frontier of the radiative zone and a convective core. Indeed, strong chemical gradients are expected in such a boundary, so that whether mixing occurs or not has a strong influence on stellar evolution. Section 2 presents the formalism and the methods used in our simulations. Section 3 is devoted to the case of a neutral chemical stratification, whereas Sect. 4 deals with the case of stable chemical stratification. Then, the astrophysical consequences of our results are discussed in Sect. 5, along with some prospects.

2. Formalism and methods

The mathematical model and the numerical configuration used to perform the simulations presented in this paper are essentially the same as in Paper I, except for the dynamical effect of stable chemical stratification. The flow is constituted by uniform, forced, vertical mean velocity shear and temperature and concentration gradients. In this configuration, turbulence can be homogeneous and stationary, which implies that turbulent transport of chemical elements can be described as a diffusive process. In Sect. 2.1, we recap the governing equations in the presence of this chemical stratification. Then, we list in Sect. 2.2 physical and numerical parameters on which our simulations depend. Finally, Sect 2.3 is dedicated to the different methods of determining the turbulent diffusion coefficient that we have explored before choosing the current one.

2.1. Governing equations

The flow configuration, including boundary conditions and forcing terms, is detailed in Paper I so will not be reproduced here. In contrast, the governing equations are modified to take the fluctuations of mean molecular weight due to concentration fluctuations into account. The new equations are derived in Sect. 2.1.1 from the Boussinesq equations, and the SPNA, which is valid for high thermal diffusivities, is then applied in Sect. 2.1.2.

2.1.1. Boussinesq approximation

In their classical form, the Boussinesq equations read

$\begin{array}{ccc}& & \mathrm{\nabla }\mathrm{·}{v}\mathrm{=}\mathrm{0}\mathit{,}\\ & & \end{array}$where v is the velocity of the fluid, p the pressure deviation from hydrostatic equilibrium, ρ₀ the mean density, ρ′ the density fluctuations, g = −ge_z the gravitational acceleration, ν the kinematic viscosity, and f_v the forcing. Density fluctuations can be expressed in terms of temperature and concentration (through mean molecular weight) fluctuations θ and c′ with respect to mean temperature and concentration profiles T₀ and c₀, here assumed to be linear: $\frac{{\mathit{\rho }}^{\mathrm{\prime }}}{\mathit{\rho }}\mathrm{=}\mathrm{-}{\mathit{\alpha }}_{\mathrm{T}}\mathit{\theta }\mathrm{-}{\mathit{\alpha }}_{\mathrm{c}}{\mathit{c}}^{\mathrm{\prime }}\mathit{,}$(4)introducing the thermal expansion coefficient α_T and its chemical equivalent α_c. Then, Eq. (3)becomes $\frac{\mathit{\partial }{v}}{\mathit{\partial t}}\mathrm{+}\mathrm{(}{v}\mathrm{·}\mathrm{\nabla }\mathrm{)}{v}\mathrm{=}\mathrm{-}\frac{\mathrm{\nabla }\mathit{p}}{{\mathit{\rho }}_{\mathrm{0}}}\mathrm{-}\mathrm{(}{\mathit{\alpha }}_{\mathrm{T}}\mathit{\theta }\mathrm{+}{\mathit{\alpha }}_{\mathrm{c}}{\mathit{c}}^{\mathrm{\prime }}\mathrm{)}{g}\mathrm{+}\mathit{\nu }\mathrm{\Delta }{v}\mathrm{+}{{f}}_{\mathrm{v}}\mathit{,}$(5)in which temperature T = T₀ + θ and concentration c = c₀ + c′ verify the advection/diffusion equations: $\begin{array}{ccc}\frac{\mathit{\partial \theta }}{\mathit{\partial t}}\mathrm{+}{v}\mathrm{·}\mathrm{\nabla }\mathit{T}& \mathrm{=}& \mathit{\kappa }\mathrm{\Delta }\mathit{\theta }\mathrm{+}{\mathit{f}}_{\mathrm{T}}\mathit{,}\\ \frac{\mathit{\partial }{\mathit{c}}^{\mathrm{\prime }}}{\mathit{\partial t}}\mathrm{+}{v}\mathrm{·}\mathrm{\nabla }\mathit{c}& \mathrm{=}& {\mathit{D}}_{\mathrm{m}}\mathrm{\Delta }{\mathit{c}}^{\mathrm{\prime }}\mathrm{+}{\mathit{f}}_{\mathrm{c}}\mathit{,}\end{array}$with κ and D_m the thermal and molecular diffusivities, and f_T and f_c the thermal and chemical forcing terms.

The dimensionless form of the previous equations is obtained by using L, which is the characteristic size of the numerical domain, S^-1, ΔU = SL, ΔT = LdT₀/ dz, Δc = Ldc₀/ dz, and ρ₀S²L² as length, time, velocity, temperature, concentration, and pressure units, respectively. Denoting dimensionless variables by a tilde, we write

$\begin{array}{ccc}& & \begin{array}{c}\mathrm{˜}\\ \mathrm{\nabla }\end{array}\mathrm{·}\begin{array}{c}{˜}\\ {v}\end{array}\mathrm{=}\mathrm{0}\mathit{,}\\ & & \\ & & \\ & & \end{array}$where there are five dimensionless parameters: the Richardson and Péclet numbers that we have already defined, the chemical Richardson number Ri_μ = α_cg/S²dc₀/ dz, the Reynolds number Re = SL²/ν, and the chemical Péclet number Pe_c = SL²/D_m. In this paper, we study the dependence of the turbulent diffusion coefficient on the turbulent Péclet number and the chemical Richardson number.

2.1.2. Small-Péclet number approximation

Two conditions are required so that the stabilising effect of thermal stratification is affected by thermal diffusion. Thermal diffusion and stratification must both play a dynamical role. In other words, their characteristic time scales (the thermal diffusive one ℓ²/κ and the buoyancy one N^-1) must be smaller than the dynamical one, which can be approximated by S^-1. If we suppose that u ≃ Sℓ, the first condition can be simply written Pe_ℓ< 1. The second one straightforwardly implies that Ri> 1. This is the regime we want to study in stars, but when Pe_ℓ ≪ 1, the problem becomes very costly to solve numerically without approximation, because the diffusive time scale becomes much smaller than the dynamical one.

Mathematically, this can be avoided by using the SPNA (Lignières 1999), which consists in a first-order Taylor expansion of the Boussinesq equations with respect to the Péclet number. In this approximation, Eqs. (9)and (10)become

$\begin{array}{ccc}& & \frac{\mathit{\partial }\begin{array}{c}{˜}\\ {v}\end{array}}{\mathit{\partial }\begin{array}{c}\mathrm{˜}\\ \mathit{t}\end{array}}\mathrm{+}\mathrm{(}\begin{array}{c}{˜}\\ {v}\end{array}\mathrm{·}\begin{array}{c}\mathrm{˜}\\ \mathrm{\nabla }\end{array}\mathrm{)}\begin{array}{c}{˜}\\ {v}\end{array}\mathrm{=}\mathrm{-}\begin{array}{c}\mathrm{˜}\\ \mathrm{\nabla }\end{array}\begin{array}{c}\mathrm{˜}\\ \mathit{p}\end{array}\mathrm{+}\mathrm{(}\mathit{RiPe}\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}\mathrm{+}\mathit{R}{\mathit{i}}_{\mathit{\mu }}\stackrel{}{{{\mathit{c}}^{\mathrm{\prime }}}_{\mathrm{˜}}}\mathrm{)}{{e}}_{\mathit{z}}\mathrm{+}\frac{\mathrm{1}}{\mathit{Re}}\begin{array}{c}\mathrm{˜}\\ \mathrm{\Delta }\end{array}\begin{array}{c}{˜}\\ {v}\end{array}\mathrm{+}\begin{array}{c}{˜}\\ {{f}}_{{v}}\end{array}\mathit{,}\\ & & \begin{array}{c}\mathrm{˜}\\ {\mathit{v}}_{\mathit{z}}\end{array}\mathrm{=}\begin{array}{c}\mathrm{˜}\\ \mathrm{\Delta }\end{array}\begin{array}{c}\mathrm{˜}\\ \mathit{\psi }\end{array}\mathit{,}\end{array}$where $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{\psi }\end{array}$ is defined by $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\theta }\end{array}\mathrm{=}\mathit{Pe}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\psi }\end{array}$. It can be noticed that the Richardson and Péclet numbers do not play independent roles in the SPNA. Instead, the effects of both thermal stratification and diffusion become a single physical effect characterised by the number RiPe, which we call the “Richardson-Péclet” number, and its time scale is $\mathit{\tau }\mathrm{=}\frac{\mathit{\kappa }}{{\mathit{N}}^{\mathrm{2}}{\mathit{\ell }}^{\mathrm{2}}}\mathrm{·}$(14)In constrast, chemical stratification is not affected by thermal diffusion.

2.2. Physical and numerical parameters of the simulations

In this study, we do not want to test Zahn’s assumption that turbulent flows tend to reach a statistical steady state. On the contrary, we tune the Richardson number (or the Richardson-Péclet number in the SPNA) to obtain this steady state. We define the stationary Richarson number Ri_s as the highest Richardson number for which turbulent kinetic energy does not decrease statistically in time. For Ri<Ri_s, the flow eventually saturates into a statistical steady state, but at much larger scales, which we suspect to be close to the dimensions of the simulation domain. We are not interested in such steady states because then the flow and the associated transport are expected to strongly depend on the size of the numerical domain. We also define the stationary Richardson-Péclet number (RiPe)_s in the SPNA similarly. In this regime, changing RiPe is more or less equivalent to changing the turbulent scales, the flow evolving towards a steady state characterised by a unique RiPe_ℓ. For smaller RiPe, the turbulent scales are larger, and then get closer to the dimensions of the simulation domain.

Given the parameters we want to vary (the Péclet number and the chemical Richardson number) and those we want to be fixed (the Reynolds number and the chemical Péclet number), the only remaining physical parameters are thus linked to the initial conditions, which we discuss in Sect. 2.2.1. They are also numerical parameters, namely the dimensions of the simulation domain (Sect. 2.2.2) and the resolution of our simulations (Sect. 2.2.3).

2.2.1. Initial conditions

Initial temperature and concentration fluctuations are set to zero. For initial velocity fluctuations, we use an incompressible, isotropic field with a given power spectral density E(k) generated thanks to the method explained by Orszag (1969). Since the resulting velocity field is not a natural solution of the Navier-Stokes equations, it takes some time to reach a statistical steady state. To reduce this relaxation time, Jacobitz et al. (1997) have deduced, from spectra of just relaxed velocity fields, the optimal power spectral density to use, that is, $\mathit{E}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{\propto }{\mathit{k}}^{\mathrm{2}}{\mathrm{e}}^{\mathrm{-}{\mathit{k}}^{\mathrm{2}}\mathit{/}{{\mathit{k}}_{\mathrm{0}}}^{\mathrm{2}}}\mathit{,}$(15)where k is the norm of the wave vector and k₀ the value corresponding to the peak of the spectrum. With this spectrum, initial conditions are now characterised by two parameters: the amplitude of initial fluctuations, which can be described by the quadratic mean of velocity fluctuations q, and the peak wave vector, which is linked to the horizontal integral length scale of turbulence ℓ_i, defined by ${\mathit{\ell }}_{\mathrm{i}}\mathrm{=}\mathrm{2}\mathit{\pi }\frac{\int \begin{array}{c}\mathrm{+}\mathrm{\infty }\\ \mathrm{0}\end{array}\mathit{E}\mathrm{\left(}{\mathit{k}}_{\mathrm{h}}\mathrm{\right)}\mathit{/}{\mathit{k}}_{\mathrm{h}}\mathrm{d}\mathit{k}h}{{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathit{E}\mathrm{\left(}{\mathit{k}}_{\mathrm{h}}\mathrm{\right)}\mathrm{d}{\mathit{k}}_{\mathrm{h}}}\mathit{,}$(16)where k_h is the norm of the horizontal wave vector.

From these parameters, we build two dimensionless parameters, the turbulent Reynolds number based on q and ℓ_i defined as $\begin{array}{ccc}\mathit{R}{\mathit{e}}_{\mathit{\ell }}\mathrm{=}\frac{\mathit{q}{\mathit{\ell }}_{\mathrm{i}}}{\mathit{\nu }}\mathit{,}& & \end{array}$(17)and what we call the “shear number” S′ defined by ${\mathit{S}}^{\mathrm{\prime }}\mathrm{=}\frac{\mathit{S}{\mathit{\ell }}_{\mathrm{i}}}{\mathit{q}}\mathrm{·}$(18)To determine the influence of initial conditions on the final statistical steady state, we have computed final values of these parameters as a function of their initial values. The main result is that the final value of the shear number only slightly depends on its initial value, as illustrated in Fig. 1.

Fig. 1 Evolution of S′ as a function of time for different initial values.

For initial values ranging from 0.25 to 8.5, we obtain final values between 1 and 2. Another comforting point about the generic nature of our simulations is that in those we have used in this paper, we observe a self-sustained mechanism similar to the one observed by Waleffe (1997). This kind of mechanism, in which streamwise vortices transform into longitudinal streaks that disappear in favour of oscillating modes, which finally regenerate vortices, is known to be quite generic in shear flows. In numerical simulations, it is present as long as the numerical domain and the Reynolds number are large enough. Provided these conditions are met, the existence of the self-sustained process depends neither on the size of the simulation domain (Yakhot 2003) nor on the type of forcing (Waleffe 1997).

2.2.2. Dimensions of the numerical domain

Because of the stochastic nature of turbulence, the efficiency of turbulent transport is an averaged quantity, which can be obtained by statistical analyses. It means that to obtain a physically relevant value of the transport coefficient, we need to have enough large structures (those which actually do the transport) in our simulation domain. The more structures there are, the better the precision in determining the coefficient. Nevertheless, increasing the number of large structures also leads to higher computational cost.

The size of large structures can be approximated by the integral scale ℓ_i. For our simulations, we have chosen a numerical domain such that the vertical and horizontal dimensions L_v and L_h of the domain satisfy L_h/L_v = π/ 4. In all the simulations, we have 0.13 <ℓ_i/L_h< 0.21 and 0.10 <ℓ_i/L_v< 0.17, which ensures that between four and eight large structures in each horizontal direction and at least six in the vertical direction are present in the numerical domain.

2.2.3. Numerical resolution

Here, we consider DNS, where all scales down to the dissipation scale must be numerically resolved, as opposed to LES (large eddy simulations), where small scales are modelled with so-called sub-grid models. Since the dissipation scale depends on the Reynolds number, the maximum value of the latter is limited by the resolution. For the present configuration, 128 grid points in every horizontal direction and 257 ones in the vertical direction are enough to correctly resolve our simulations with a typical Reynolds number of 10⁴ (see Paper I for more details). This is confirmed by the form of the spectrum of turbulent kinetic energy, presented in Fig. 2.

Fig. 2 Spectrum of turbulent kinetic energy. The dotted line corresponds to the largest relevant wave vector available in our simulations.

One can clearly see that far away from the dissipation range, the spectrum follows a k^{− 5 / 3} law typical of Kolmogorov turbulence and rapidly drops once the dissipation scale is reached.

2.3. Test of the methods

After the numerical setup is fixed, the transport coefficient can be determined by several methods. We have tested three of them. The first one is based on Lagrangian fluid particles (Sect. 2.3.1), the second one on the widening of a Gaussian mean concentration profile (Sect. 2.3.2), and the last one on a relation between the turbulent flux and the mean gradient of concentration (Sect. 2.3.3).

2.3.1. Particle tracking

Since the considered flow is stationary and homogeneous from an Eulerian point of view away from the lower and upper boundaries (see Paper I for more details), our turbulence is homogeneous from a Lagrangian point of view in this region. It means that the statistical properties of turbulence do not vary with time, seen from a fluid particle moving in the flow. In these conditions, the theory of Taylor (1921) applies and gives us the value of the transport coefficient D_t through the slope of the mean quadratic vertical dispersion of particles as a function of time: $\mathrm{⟨}\mathit{\delta }{\mathit{z}}^{\mathrm{2}}\mathrm{⟩}\mathrm{=}\mathrm{⟨}\mathrm{\left[}\mathit{z}\mathrm{\right(}\mathit{t}\mathrm{)}\mathrm{-}\mathit{z}\mathrm{(}\mathrm{0}\mathrm{)}{\mathrm{]}}^{\mathrm{2}}\mathrm{⟩}\mathrm{=}\mathrm{2}{\mathit{D}}_{\mathrm{t}}\mathit{t,}$(19)where ⟨⟩ denotes an ensemble average or, equivalently, an average over particles. Initially, we homogeneously put particles in the fluid and follow their displacements due to the velocity field. Figure 3 shows the evolution of the mean quadratic vertical dispersion with time for one of our simulations.

Fig. 3 Mean quadratic vertical dispersion against time. The solid line is the result of the simulation, and the dashed one represents the linear regression up to St = 5. For St< 5, the mean quadratic vertical dispersion grows linearly with time and saturates afterwards.

We observe that the transport is diffusive only for short times.

The non-diffusive regime appears when a significant part of particles have reached the lower or higher boundaries of the domain and thus cannot go further because of the impermeability of these boundaries. Then, the turbulence is no longer homogeneous from a Lagrangian point of view. As a consequence, the main fault of this method is that it only allows us to study turbulent transport for a limited time. Because of natural fluctuations of the turbulent flow, determining the turbulent diffusion coefficient is thus subject to a large dispersion (up to 20% of relative difference).

2.3.2. Gaussian initial profile

Another method of estimating the transport coefficient is to study the evolution of the width of a Gaussian mean concentration field. Whereas the previous method was a Lagrangian one, the present one is Eulerian in the sense that we solve Eq. (11)(without the forcing term). The mean vertical concentration profile is assumed to be a Gaussian given by $\mathit{C}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}\mathrm{=}\frac{{\mathit{C}}_{\mathrm{0}}}{\mathit{\sigma }\mathrm{\left(}\mathit{t}\mathrm{\right)}\sqrt{\mathrm{2}\mathit{\pi }}}{\mathrm{e}}^{\mathrm{-}\mathrm{(}\mathit{z}\mathrm{-}{\mathit{z}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}\mathit{/}\left[\mathrm{2}\mathit{\sigma }\mathrm{(}\mathit{t}{\mathrm{)}}^{\mathrm{2}}\right]}\mathit{,}$(20)where z₀ is the vertical position of the peak of the profile and σ its standard deviation. Indeed, if turbulent transport can be considered as a purely diffusive process, the profile remains Gaussian and the effective diffusion coefficient D_eq verifies ${\mathit{D}}_{\mathrm{eq}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathrm{d}{\mathit{\sigma }}^{\mathrm{2}}}{\mathrm{d}\mathit{t}}\mathrm{·}$(21)Initially, the chemical element is mainly concentrated in the plane z = z₀ = L_v/ 2. An example of evolution of σ² with time is represented in Fig. 4.

Fig. 4 Variance of the mean vertical concentration profile against time. The solid line corresponds to the simulation, and the dashed one to the affine regression.

For short times, the variance of the mean profile grows linearly with time, thus indicating that the transport is diffusive, as already observed with the Lagrangian method. Moreover, the values obtain by the two methods are similar (less than 15% of relative difference). For longer times, the Eulerian has a similar problem to the previous one: once the wings of the profile reach the boundaries, it is no longer Gaussian and Eq. (21)can no longer be used to estimate the transport coefficient. Again, the limited averaging time of this method induces strong temporal fluctuations of the computed diffusion coefficient.

2.3.3. Relation between turbulent flux and mean gradient

In the presence of a mean concentration gradient, one can write the following relation between the mean flux of concentration fluctuations and the mean concentration gradient, provided that transport is diffusive, which is the case here: $\mathrm{⟨}{\mathit{c}}^{\mathrm{\prime }}{\mathit{v}}_{\mathit{z}}\mathrm{⟩}\mathrm{=}\mathrm{-}{\mathit{D}}_{\mathrm{t}}\frac{\mathrm{d}{\mathit{c}}_{\mathrm{0}}}{\mathrm{d}\mathit{z}}\mathrm{·}$(22)It is thus possible to determine the turbulent diffusion coefficient by measuring the mean concentration flux in the flow and computing ${\mathit{D}}_{\mathrm{t}}\mathrm{=}\mathrm{-}\frac{\mathrm{⟨}{\mathit{c}}^{\mathrm{\prime }}{\mathit{v}}_{\mathit{z}}\mathrm{⟩}}{\mathrm{d}{\mathit{c}}_{\mathrm{0}}\mathit{/}\mathrm{d}\mathit{z}}\mathrm{·}$(23)Thanks to the forcing term in Eq. (11)combined with boundary conditions similar to those used for temperature, the mean gradient remains statistically constant, so that the transport coefficient can be averaged on much longer times than for the two other methods. It follows that the precision on the estimate of this coefficient is better in this method, leading typically to relative fluctuations of only 5% when averaging on several hundred shear times S^-1. All the results published in Paper I and those presented thereafter have been obtained with this method.

3. Thermal diffusion in the chemically neutral case

This section deals with the effect of thermal diffusion on the vertical turbulent diffusion coefficient when chemical stratification is neutral (Ri_μ = 0). To study this effect, we performed full Boussinesq simulations with Re = 10⁴ and Pe ranging from 10 to 10⁴ (this corresponds to turbulent Péclet numbers Pe_ℓ = qℓ_i/κ, which are more relevant in this context, from 0.34 to 437), and one simulation in the small-Péclet-number approximation. Table 1 displays the average values of some relevant physical parameters measured in the statistical steady state, including the turbulent Reynolds and Péclet numbers, the stationary Richardson number used to obtain the steady state, the product of the latter and the turbulent Péclet number, the ratio β between the turbulent diffusion coefficient and the product of turbulent scales uℓ_i, and the specific dissipation rate of thermal potential energy ε_P = κ ⟨ (∇b)² ⟩, where b = α_Tgθ/N is the “buoyancy”, which is proportional to temperature fluctuations and has the dimension of a velocity.

Thermal potential energy is the potential energy linked to temperature fluctuations gained by a fluid particle moving away from equilibrium.

Firstly, we focus on the regime of small Péclet numbers, where new results confirm those presented in Paper I (Sect. 3.1). Then, we consider the large-Péclet-number regime and compare our numerical results with existing models (Sect. 3.2).

3.1. Small-Péclet-number regime

In this section, we are interested in the regime of small Péclet numbers. We performed numerical simulations both in the Boussinesq approximation at low Péclet number and in the SPNA. We compare the results in Sect. 3.1.1 and try to interpret them in Sect. 3.1.2.

3.1.1. Results

The first result is that Boussinesq simulations tend towards the SPNA one when the Péclet number decreases, thus validating the use of the SPNA in the non-linear regime for Péclet numbers smaller than 0.34. Indeed, if we look at the results of the SPNA simulation and the Boussinesq one at Pe_ℓ = 0.34, we see that the value of Ri_sPe_ℓ in the Boussinesq simulation and that of (RiPe_ℓ)_s in the SPNA one are very close to each other. Moreover, the values of β = D_t/ (qℓ_i) are also very similar in the two simulations.

Furthermore, this proves that in this regime, stability is no longer characterised by the Richardson number alone, but by the product of the Richardson number and the turbulent Péclet number, as assumed by Zahn (1974). The assumption that the turbulent diffusion coefficient is proportional to uℓ_i is also satisfied. This seems to indicate that in the small-Péclet-number regime, the turbulent diffusion coefficient D_t is proportional to κ/Ri, and that this regime is already reached at Pe_ℓ = 0.34, which is not very low. This is confirmed by Fig. 5, which plots D_t/ (κRi^-1) against the turbulent Péclet number and suggests that for small Péclet numbers, D_t/ (κRi^-1) tends to a constant value that is close to what is obtained in the SPNA simulation.

Fig. 5 Dt/ (κRi-1) as a function of the turbulent Péclet number. Dots correspond to full Boussinesq simulations, whereas the solid line represents the value obtained in the SPNA one.

3.1.2. Interpretation

These results qualitatively agree with Zahn’s model, but we also found a quantitative difference on the proportionality constant of about 30%, as discussed in Paper I. This is not surprising because Zahn’s prescription gives only an order of magnitude for the diffusion coefficient. One of the contributions of numerical simulations is precisely to give a quantitative estimate of such coefficients that cannot be deduced directly from theoretical arguments.

Another interesting physical discussion that was not made in Paper I is that thermal diffusion has two opposite effects on vertical motions: (i) decreasing the amplitude of the buoyancy force, which favours such motions; and (ii) contributing to the dissipation of kinetic energy, by dissipating potential energy, which inhibits these motions. The first effect is taken into account in Zahn’s model, but the second one, which is responsible for the damping of gravity modes, is not. That despite this, Zahn’s model seems to be valid can be understood by considering that when thermal diffusivity is very high, all available potential energy will be instantaneously dissipated, so that increasing the diffusivity does not dissipate more, whereas the amplitude of the buoyancy force is still decreasing. The net effect is thus to reduce the stabilising effect of stratification, as accounted for by Zahn. The balance between both effects may be an explanation for the strange behaviour observed in Fig. 5 for Pe_ℓ = 0.90.

We propose now an alternative derivation of Zahn’s model based on the mixing-length theory (MLT) in the SPNA. The main idea is that turbulence can only occur if the dynamical time scale is smaller than that of the process which tends to inhibit it. Here it is stable stratification modified by thermal diffusion, whose time scale τ is given in Eq. (14). The previous condition reads as ${\mathit{S}}^{-1}\mathit{<}\frac{\mathit{\kappa }}{{\mathit{N}}^{\mathrm{2}}{\mathit{\ell }}^{\mathrm{2}}}\mathit{,}$(24)which yields ${\mathit{\ell }}^{\mathrm{2}}\mathit{<}\frac{\mathit{\kappa S}}{{\mathit{N}}^{\mathrm{2}}}\mathrm{·}$(25)We now assume that the largest scales allowed by this condition dominate the transport and that a mixing-length model D_t ~ Sℓ² holds. It follows that ${\mathit{D}}_{\mathrm{t}}\mathrm{~}\frac{\mathit{\kappa }}{\mathit{Ri}}\mathit{,}$(26)which is the form of Zahn’s model. This derivation does not rely on the assumption of marginal stability or stationarity and may thus be applicable in cases where it is not satisfied.

3.2. Large-Péclet-number regime

It is clear from Fig. 5 that D_t/ (κRi^-1) = const. is not valid in the regime of large Péclet numbers. Maeder (1995) derived a prescription by following Zahn’s approach but using a generalised Richardson criterion valid both for small and large Péclet numbers. This problem has also been studied in geophysics (Pr = ν/κ = Pe/Re ~ 1 in the atmosphere or the oceans), notably by Lindborg & Brethouwer (2008). Section 3.2.1, we compare our simulations, already presented in Table. 1, to Maeder’s model and show they are not compatible. Then, we explain how the model proposed by Lindborg & Brethouwer (2008) is supported by our simulations in Sect. 3.2.2.

3.2.1. Comparison with Maeder’s model

The idea of Maeder (1995) is to use a generalised Richardson criterion for marginal stability of the form $\mathit{Ri}\frac{\mathit{P}{\mathit{e}}_{\mathit{\ell }}}{\mathrm{1}\mathrm{+}\mathit{P}{\mathit{e}}_{\mathit{\ell }}}\mathrm{=}\mathit{R}{\mathit{i}}_{\mathrm{c}}\mathit{,}$(27)which reduces to the classical one Ri = Ri_c for large Péclet numbers and to that of Zahn (1974) for small ones. This criterion then leads to the following expression for the turbulent diffusion coefficient: ${\mathit{D}}_{\mathrm{t}}\mathrm{=}\mathit{\beta \kappa }\frac{\mathit{R}{\mathit{i}}_{\mathrm{c}}}{\mathit{Ri}\mathrm{-}\mathit{R}{\mathit{i}}_{\mathrm{c}}}\mathrm{·}$(28)To compare it with our simulations, the previous relation can be written $\frac{\mathit{\kappa }}{{\mathit{D}}_{\mathrm{t}}}\mathrm{=}\frac{\mathit{Ri}\mathrm{-}\mathit{R}{\mathit{i}}_{\mathrm{c}}}{\mathit{\beta R}{\mathit{i}}_{\mathrm{c}}}\mathit{,}$(29)which shows that in Maeder’s model there is an affine relation between κ/D_t and the Richardson number. As illustrated in Fig. 6, such a relation is not supported by our simulations for Ri< 0.15, which corresponds to Pe_ℓ> 25.

Fig. 6 κ/Dt against the Richardson number. Dots correspond to our simulations, whereas the solid line represents Maeder’s model calibrated in the small-Péclet-number regime, where it reduces to Zahn’s model.

Indeed, Maeder’s model predicts that at large Péclet numbers, the stationary Richardson number should tend towards Ri_c, which is supposed to be constant, whereas in our simulations Ri_s reaches lower values.

Maeder’s model was intended to be a general model valid for all Péclet numbers. If this model agrees with our numerical results in the small-Péclet-number regime, where it reduces to Zahn’s, this is no longer the case in the large-Péclet-number regime, which is however the regime it has been created for. One potential reason is that Eq. (27)assumes that thermal diffusion modifies marginal stability even at a high Péclet number, a property that is not found in the linear stability analysis (Lignières et al. 1999). As we see next, thermal diffusion does play an important role in radial transport in the high-Péclet-number regime, but only on motions with small enough length scales, which are produced by the turbulent cascade from the instability injection scale.

3.2.2. Comparison with Lindborg’s model

Lindborg’s model is based on the buoyancy equation of a Lagrangian fluid particle that reads $\frac{\mathrm{d}\mathit{b}}{\mathrm{d}\mathit{t}}\mathrm{+}\mathit{N}{\mathit{v}}_{\mathit{z}}\mathrm{=}\mathit{\kappa }\mathrm{\Delta }\mathit{b}\mathit{.}$(30)This equation can be used to describe the temporal evolution of the mean vertical dispersion of fluid particles. Thanks to certain assumptions on the properties of turbulence, including Lagrangian statistical homogeneity, Lindborg & Brethouwer (2008) are thus able to write $\mathrm{⟨}\mathit{\delta }{\mathit{z}}^{\mathrm{2}}\mathrm{⟩}\mathrm{=}\frac{\mathrm{⟨}\mathit{\delta }{\mathit{b}}^{\mathrm{2}}\mathrm{⟩}}{{\mathit{N}}^{\mathrm{2}}}\mathrm{+}\mathrm{2}\frac{{\mathit{\epsilon }}_{\mathrm{P}}}{{\mathit{N}}^{\mathrm{2}}}\mathit{t,}$(31)where ⟨ δb² ⟩ is the mean buoyancy dispersion and ε_P the specific dissipation rate of potential energy linked to buoyancy. At late times, the first term, which corresponds to the case of zero thermal diffusion, tends to a constant value, whereas the second one increases linearly with time. The corresponding diffusion coefficient is ${\mathit{D}}_{\mathrm{t}}\mathrm{=}\frac{{\mathit{\epsilon }}_{\mathrm{P}}}{{\mathit{N}}^{\mathrm{2}}}\mathrm{·}$(32)One interesting property of this model is that it does not involve any free parameter, which is not the case for many phenomenological models currently used in stellar evolution codes. We tested this model by computing the ratio D_t/ (ε_PN^-2) for the different simulations of Table 1. As illustrated in Fig. 7, one notices that D_t/ (ε_PN^-2) tends towards one for large Péclet numbers.

Fig. 7 Dt/ (εPN-2) against the turbulent Péclet number. Same legend as in Fig. 5.

Lindborg’s model is thus valid for large Péclet numbers, but clearly not for small ones. This is because Lindborg & Brethouwer (2008) neglected in their study a term that is the dominant one for high thermal diffusivities. Based on this observation, one can then wonder if it is possible to derive a general model based on Zahn’s and Lindborg’s models, which are each valid in its own regime, by taking the sum. Figure 8 thus compares the ratios between each model and the total transport observed in our simulations.

Fig. 8 Ratios between Zahn’s and Lindborg’s models and the total transport observed in our simulations. Same legend as in Fig. 5.

The important result is that the ratio between the sum of the two models and the total observed transport is generally close to one. We can thus write ${\mathit{D}}_{\mathrm{t}}\mathrm{=}\mathit{\beta \kappa }\frac{\mathit{R}{\mathit{i}}_{\mathrm{c}}}{\mathit{Ri}}\mathrm{+}\frac{{\mathit{\epsilon }}_{\mathrm{P}}}{{\mathit{N}}^{\mathrm{2}}}\mathrm{=}\frac{\mathrm{1}}{\mathit{Ri}}\left(\mathit{\beta \kappa R}{\mathit{i}}_{\mathrm{c}}\mathrm{+}\frac{{\mathit{\epsilon }}_{\mathrm{P}}}{{\mathit{S}}^{\mathrm{2}}}\right)\mathit{,}$(33)which is approximately valid for all Péclet numbers. This prescription for turbulent transport cannot be used in stellar evolution codes as long as ε_P is not expressed in terms of the mean quantities characterising the flow in these codes. This will be done in a future work.

4. Chemical stratification

In this section we are interested in the dependence of the turbulent diffusion coefficient with respect to chemical stratification. To study this dependence, we performed several simulations with different values of the chemical Richardson number. These simulations and the empirical form of the turbulent diffusion coefficient deduced from them are presented in Sect. 4.1. In Sect. 4.2, we compare the form of the diffusion coefficient with existing prescriptions currently used in stellar evolution codes. The results are discussed in Sect. 4.3.

4.1. Results

Since the SPNA gives satisfying results in the neutral case, we have chosen to use it here, too. For each value of the chemical Richardson number, we have tuned the Richardson-Péclet number to obtain a statistical steady state and determined the turbulent diffusion coefficient in this state. The results of our simulations are summarised in Table 2.

It is important to note that no stationary state could be reached for chemical Richardson numbers greater than 0.12, since turbulence decays regardless of the value of the thermal Richardson number. Figure 9 displays D_t/ (κRi^-1) as a function of the chemical Richardson number.

Fig. 9 Dt/ (κRi-1) against the chemical Richardson number. Dots correspond to our simulations and the solid line to the affine regression.

It appears that this plot is described well by an affine relation that is written $\frac{{\mathit{D}}_{\mathrm{t}}}{\mathit{\kappa R}{\mathit{i}}^{-1}}\mathrm{=}\mathit{\beta }\mathrm{(}\mathit{R}{\mathit{i}}_{\mathrm{cr}}\mathrm{-}\mathit{R}{\mathit{i}}_{\mathit{\mu }}\mathrm{)}\mathit{,}$(34)with β = 0.45 and Ri_cr = 0.12. The turbulent diffusion coefficient is then given by ${\mathit{D}}_{\mathrm{t}}\mathrm{=}\mathit{\beta \kappa }\frac{\mathit{R}{\mathit{i}}_{\mathrm{cr}}\mathrm{-}\mathit{R}{\mathit{i}}_{\mathit{\mu }}}{\mathit{Ri}}\mathrm{·}$(35)Unlike the neutral case in the small-Péclet-number regime, where we found that D_t only depends on one unknown constant, namely the product βRi_c = D_tRi/κ, here it depends on two constants β and Ri_cr that can be determined independently and that indeed play different physical roles.

4.2. Comparison with existing models

Several models of turbulent transport in the presence of a mean gradient of mean molecular weight have been proposed for use in stellar evolution codes, and we have compared our simulations with two of them: Maeder & Meynet (1996) and Maeder (1997). In particular, we present in Sect. 4.2.1 to what extent our prescription agrees with the model of Maeder & Meynet (1996). Then we explain in Sect. 4.2.2 why our simulations are incompatible with the model of Maeder (1997).

4.2.1. Maeder & Meynet’s model

Maeder & Meynet (1996) extended the model of Maeder (1995) to the chemically stratified case. As seen before, we do not expect such a model to be valid for large Péclet numbers. However, for small ones, the authors found a form of the turbulent diffusion coefficient that is similar to Eq. (35). To obtain this result, they assumed that the relevant marginal stability criterion reads as $\mathit{RiP}{\mathit{e}}_{\mathit{\ell }}\mathrm{+}\mathit{R}{\mathit{i}}_{\mathit{\mu }}\mathrm{=}\mathit{R}{\mathit{i}}_{\mathrm{c}}\mathit{,}$(36)where Ri_c can be seen as a critical chemical Richardson number. Assuming further that D_t = βuℓ, one obtains the same form as Eq. (35). Since the authors used the same numerical values for β and Ri_c as Zahn, 1/3 and 1/4, respectively, their model shows quantitative differences with our numerical results: around 25% of relative difference for β and 50% for Ri_c. An important aspect of this criterion is that chemical stratification is not affected by thermal diffusion. Thus, a sufficiently stable chemical stratification is able to totally inhibit turbulent transport.

4.2.2. Maeder’s model

In contrast, the model of Maeder (1997) is based on the assumption that mixing can occur for any value of the chemical Richardson number, provided that the flow is thermally unstable. As a consequence, stability is then given by the criterion of Zahn (1974), and the turbulent diffusion coefficient is ${\mathit{D}}_{\mathrm{t}}\mathrm{=}\mathit{\beta \kappa }\frac{\mathit{R}{\mathit{i}}_{\mathrm{c}}}{\mathit{Ri}\mathrm{+}\mathit{R}{\mathit{i}}_{\mathit{\mu }}}\mathrm{·}$(37)One can easily notice that when κ and Ri are very high and Ri_μ remains finite, this coefficient reduces to that of Zahn (1992), which does not depend on the chemical Richardson number. This is clearly not what we observe in our simulations. Maeder’s model is thus incompatible with our results.

4.3. Interpretation

Looking at the previous results, one may be tempted to conclude that the model of Maeder & Meynet (1996) is more realistic than that of Maeder (1997) and therefore should be used in stellar evolution codes. However, Meynet et al. (2013) find that Maeder (1997) is the prescription that has the best fit to the observation. This is likely to be because this prescription enables mixing in strongly chemically stratified zones, such as the frontier between the convective core and the radiative envelope of massive stars. The model of Maeder & Meynet (1996) was not included in this study, but it predicts a less efficient mixing in such zones than those considered in the study. If we admit that Maeder (1997) overestimates the mixing efficiency in radiative zones of stars, then an additional source of mixing is needed to explain why it obtains such good agreement with the observations.

5. Discussion

We have tested several models of turbulent transport. For large Péclet numbers, our simulations invalidate the model of Maeder (1995) and validate that of Lindborg & Brethouwer (2008). Nevertheless, the latter cannot be applied as such in stellar evolution codes. As regards the dynamical effect of chemical stratification, our simulations agree with a special case of the model of Maeder & Meynet (1996), but not with that of Maeder (1997), which nevertheless gives the better results when compared with observations.

These results were obtained with a flow configuration that depends on several physical and numerical parameters, such as initial conditions and the Reynolds number. As regards initial conditions, we have chosen a parameter domain in which the flow is likely to be the most generic, as discussed in Sect. 2.2.1. For the Reynolds number, a preliminary study shows that transport does not significantly vary for turbulent Reynolds numbers ranging from 225 to 340. According to Michaud & Zahn (1998), these values are consistent with what we expect in stellar radiative zones.

In addition, this configuration relies on certain hypotheses, including stationarity and homogeneity, that are forced in an uncommon way. This configuration is generic, but not necessarily natural. Other numerical methods such as shearing boxes, used notably in geophysics (e.g. Jacobitz et al. 1997) and in astrophysics for the study of accretion discs (e.g. Lesur & Longaretti 2005), also force uniform mean shear flows, but not exactly in the same way as in the present paper. It would be interesting to see if our results can be recovered with such approaches. Moreover, unstationary and unhomogeneous configurations should also be considered. As an example, direct numerical simulations of free shear layers, such as in Brüggen & Hillebrandt (2001), could be considered.

There are still many other physical ingredients that could play a role in the radial transport due to shear instability. Among them is the effect of a very efficient horizontal turbulent diffusion, which is said to reduce the stabilising effect of chemical stratification (Talon & Zahn 1997). One can also wonder what the effects of rotation (by introducing the Coriolis force in the simulations) and magnetic field are on shear instability. An attempt has been made recently by Maeder et al. (2013) to provide some general diffusion coefficients, when simultaneously taking several hydrodynamical instabilities into account. As regards magnetic field, Lecoanet et al. (2010) show that in the linear regime particular magnetic configurations can destabilise a flow, which would be stable without it. Whether this result remains relevant in the non-linear regime is still an open question.

Until now, we have computed transport coefficients for chemical elements, but turbulent viscosity is also crucial to understanding angular momentum transport and should be estimated in the same way. In a forthcoming paper, we intend to determine relations between turbulent viscosity and turbulent thermal and chemical diffusivities.

Acknowledgments

We thank Georges Meynet for his useful comments. This work was granted access to the HPC resources of CALMIP under the allocation 2013–P0507.

References

All Tables

All Figures

	Fig. 1 Evolution of S′ as a function of time for different initial values.
In the text

	Fig. 2 Spectrum of turbulent kinetic energy. The dotted line corresponds to the largest relevant wave vector available in our simulations.
In the text

	Fig. 3 Mean quadratic vertical dispersion against time. The solid line is the result of the simulation, and the dashed one represents the linear regression up to St = 5. For St< 5, the mean quadratic vertical dispersion grows linearly with time and saturates afterwards.
In the text

	Fig. 4 Variance of the mean vertical concentration profile against time. The solid line corresponds to the simulation, and the dashed one to the affine regression.
In the text

	Fig. 5 Dt/ (κRi-1) as a function of the turbulent Péclet number. Dots correspond to full Boussinesq simulations, whereas the solid line represents the value obtained in the SPNA one.
In the text

	Fig. 6 κ/Dt against the Richardson number. Dots correspond to our simulations, whereas the solid line represents Maeder’s model calibrated in the small-Péclet-number regime, where it reduces to Zahn’s model.
In the text

	Fig. 7 Dt/ (εPN-2) against the turbulent Péclet number. Same legend as in Fig. 5.
In the text

	Fig. 8 Ratios between Zahn’s and Lindborg’s models and the total transport observed in our simulations. Same legend as in Fig. 5.
In the text

	Fig. 9 Dt/ (κRi-1) against the chemical Richardson number. Dots correspond to our simulations and the solid line to the affine regression.
In the text