Issue 
A&A
Volume 620, December 2018



Article Number  A22  
Number of page(s)  16  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/201629187  
Published online  23 November 2018 
Anisotropic turbulent transport in stably stratified rotating stellar radiation zones
^{1} IRFU, CEA, Université ParisSaclay, 91191 GifsurYvette, France
email: stephane.mathis@cea.fr
^{2} Université Paris Diderot, AIM, Sorbonne Paris Cité, CEA, CNRS, 91191 GifsurYvette, France
^{3} LUPM, Université de Montpellier, CNRS, Place E. Bataillon  cc 072, 34095 Montpellier Cedex 05, France
^{4} Department of Astronomy, University of Geneva, Chemin des Maillettes 51, 1290 Versoix, Switzerland
^{5} University of Exeter, Department of Physics & Astronomy, Stoker Road, Devon, Exeter, EX4 4QL, UK
^{6} IRAP, UMR 5277, CNRS and Université de Toulouse, 14 Av. E. Belin, 31400 Toulouse, France
^{7} Institut UTINAM, CNRS UMR 6213, Univ. Bourgogne FrancheComté, OSU THETA FrancheComtéBourgogne, Observatoire de Besançon, BP 1615, 25010 Besançon Cedex, France
Received:
26
June
2016
Accepted:
3
August
2016
Context. Rotation is one of the key physical mechanisms that deeply impact the evolution of stars. Helio and asteroseismology reveal a strong extraction of angular momentum from stellar radiation zones over the whole Hertzsprung–Russell diagram.
Aims. Turbulent transport in differentially rotating, stably stratified stellar radiation zones should be carefully modelled and its strength evaluated. Stratification and rotation imply that this turbulent transport is anisotropic. So far only phenomenological prescriptions have been proposed for the transport in the horizontal direction. This, however, constitutes a cornerstone in current theoretical formalisms for stellar hydrodynamics in evolution codes. We aim to improve its modelling.
Methods. We derived a new theoretical prescription for the anisotropy of the turbulent transport in radiation zones using a spectral formalism for turbulence that takes simultaneously stable stratification, rotation, and a radial shear into account. Then, the horizontal turbulent transport resulting from 3D turbulent motions sustained by the instability of the radial differential rotation is derived. We implemented this framework in the stellar evolution code STAREVOL and quantified its impact on the rotational and structural evolution of solar metallicity lowmass stars from the premainsequence to the red giant branch.
Results. The anisotropy of the turbulent transport scales as N^{4}τ^{2}/(2Ω^{2}), N and Ω being the buoyancy and rotation frequencies respectively and τ a time characterizing the source of turbulence. This leads to a horizontal turbulent transport of similar strength in average that those obtained with previously proposed prescriptions even if it can be locally larger below the convective envelope. Hence the models computed with the new formalism still build up too steep internal rotation gradients compared to helioseismic and asteroseismic constraints. As a consequence, a complementary transport mechanism such as internal gravity waves or magnetic fields is still needed to explain the observed strong transport of angular momentum along stellar evolution.
Conclusions. The new prescription links for the first time the anisotropy of the turbulent transport in radiation zones to their stratification and rotation. This constitutes important theoretical progress and demonstrates how turbulent closure models should be improved to get firm conclusions on the potential importance of other processes that transport angular momentum and chemicals inside stars along their evolution.
Key words: hydrodynamics / turbulence / stars: rotation / stars: evolution
© ESO 2018
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
Rotation is one of the key physical mechanisms that deeply modify the dynamics and evolution of stars (e.g. Maeder 2009). The transport of angular momentum and chemicals it induces in their stably stratified radiation zones drives their secular rotational and chemical evolution. Rotation modifies the evolutionary path of stars in the Hertzsprung–Russell diagram (hereafter HRD), their lifetime, their nucleosynthesis, chemical stratification and yields, and their magnetism.
In this context, helio and asteroseismology provide key information through the insight they give on the internal rotation profiles of the Sun and stars. On one hand, helioseismic data show that the radiative core of the Sun is rotating as an almost solid body down to r = 0.2 R_{⊙} with a potential central acceleration (Brown et al. 1989; Thompson et al. 2003; García et al. 2007; Fossat et al. 2017). On the other hand, asteroseismic data reveal a strong extraction of angular momentum over the whole HRD. First, Beck et al. (2012), Mosser et al. (2012), Deheuvels et al. (2012, 2014, 2015), Spada et al. (2016), Gehan et al. (2018) found a weak coretosurface rotation contrast in lowmass subgiant and red giant stars. Next, Benomar et al. (2015) observed 26 solartype stars with a small differential rotation between the base of their convective envelope and the upper part of their radiative core. In addition, weak differential rotation rates are found in the radiative envelope of intermediatemass and massive stars (Kurtz et al. 2014; Saio et al. 2015; Triana et al. 2015; Murphy et al. 2016; Aerts et al. 2017). Finally, a strong extraction of angular momentum is required to explain the rotation rates of white dwarfs (e.g. Suijs et al. 2008; Hermes et al. 2017) and neutron stars (e.g. Heger et al. 2005; Hirschi & Maeder 2010).
In this context, helio and asteroseismology provide key information through the insight they give on the internal rotation profiles of the Sun and stars. On one hand, helioseismic data show that the radiative core of the Sun is rotating as an almost solid body down to r = 0.2 R_{⊙} with a potential central acceleration (Brown et al. 1989; Thompson et al. 2003; García et al. 2007; Fossat et al. 2017). On the other hand, asteroseismic data reveal a strong extraction of angular momentum over the whole HRD. First, Beck et al. (2012), Mosser et al. (2012); Deheuvels et al. (2012, 2014, 2015), Spada et al. (2016); Gehan et al. (2018) found a weak coretosurface rotation contrast in lowmass subgiant and red giant stars. Next, Benomar et al. (2015) observed 26 solartype stars with a small differential rotation between the base of their convective envelope and the upper part of their radiative core. In addition, weak differential rotation rates are found in the radiative envelope of intermediatemass and massive stars (Kurtz et al. 2014; Saio et al. 2015; Triana et al. 2015; Murphy et al. 2016; Aerts et al. 2017). Finally, a strong extraction of angular momentum is required to explain the rotation rates of white dwarfs (e.g. Suijs et al. 2008; Hermes et al. 2017) and neutron stars (e.g. Heger et al. 2005; Hirschi & Maeder 2010).
In his seminal paper, Zahn (1992) proposed for the first time a consistent and complete formalism to describe the secular transport of angular momentum and chemicals under the combined action of rotationdriven vertical and horizontal turbulence and meridional flows (see also Maeder & Zahn 1998; Mathis & Zahn 2004). This theoretical treatment of rotation relies on a key physical assumption: an anisotropic turbulent transport in stellar radiation zones stronger in the horizontal (latitudinal) direction than in the vertical one because of the restoring buoyancy force along the radial entropy (and chemical) stratification. This strong horizontal transport erases horizontal gradients of physical quantities, among which angular velocity, thus enforcing a “shellular” rotation depending only on the radial coordinate. This framework allowed for a successful implementation in 1D stellar evolution codes (Talon et al. 1997; Meynet & Maeder 2000; Palacios et al. 2003; Decressin et al. 2009; Ekström et al. 2012; Marques et al. 2013; Chieffi & Limongi 2013) with numerous applications over a broad range of stellar types and evolutionary stages.
However, disagreements between the predictions of this formalism and seismic data on internal stellar rotation at various stages of the evolution (Talon & Zahn 1998; TurckChièze et al. 2010; Eggenberger et al. 2012; Marques et al. 2013; Ceillier et al. 2013) led the community to examine the role of other transport mechanisms for angular momentum, such as internal gravity waves (e.g. Talon et al. 2002; Talon & Charbonnel 2005; Charbonnel & Talon 2005; Charbonnel et al. 2013; Rogers 2015; Pinçon et al. 2017) and magnetic fields (e.g. Gough & McIntyre 1998; Spruit 1999; Mathis & Zahn 2005; Eggenberger et al. 2005; Denissenkov & Pinsonneault 2007; Strugarek et al. 2011; AcevedoArreguin et al. 2013; Barnabé et al. 2017). At the same time, the physical description of shearinduced vertical and horizontal turbulence in stably stratified stellar layers was not questioned.
The aim of the present paper is to shine new light on the hydrodynamical processes induced by rotation and on their role in the whole picture, based on the most recent numerical and theoretical developments. In this context, recent highresolution numerical simulations in local Cartesian configurations that were performed to estimate the turbulent transport induced by the instability of a vertical shear (Prat & Lignières 2013, 2014; Prat et al. 2016; Garaud & Kulenthirarajah 2016; Garaud et al. 2017) can be used for guidance. Their comparison with former phenomenological prescriptions for related turbulent transport coefficients in the radial direction (Zahn 1992; Talon & Zahn 1997) indeed invites stellar astrophysicists to reconsider this status.
This is particularly true for turbulence in the horizontal direction. Indeed, while strong stratification is invoked as the source of the strong anisotropic transport (Zahn 1992), none of the prescriptions currently in use for horizontal turbulent transport coefficients in stellar evolution codes depend explicitly on the two restoring forces in stellar radiation zones: stratification and rotation. Two mechanisms can be identified to sustain the horizontal turbulent transport: (i) the instability of the shear of the latitudinal differential rotation, but also (ii) the 3D turbulent motions induced by the instability of the radial differential rotation that transport momentum and chemicals both in the vertical and horizontal directions. Symmetrically, Zahn (1992) already identified the two sources for the vertical turbulent transport, that is, (a) the instability of the shear of the radial differential rotation and (b) the transport induced by the 3D turbulent motions induced by the instability of the latitudinal differential rotation. He introduced the corresponding vertical turbulent transport coefficients, ν_{v, v} and ν_{v,h}, modelled as eddyviscosities, and the corresponding eddydiffusivities, D_{v, v} and D_{v,h}. We note that in stellar physics modelling though, the simplifying assumption ν_{v, v} ≡ D_{v, v} and ν_{v,h} ≡ D_{v,h} has been made until now. Similarly, four coefficients ν_{h,h}, D_{h,h}, ν_{h,v} and D_{h,v} can also be defined for the horizontal transport (Table 1 is recapitulating the different turbulent diffusivities (viscosities) and their source). On one hand, the first developments focused on the instability of the latitudinal differential rotation (i.e. on ν_{h,h} and D_{h,h}). Zahn (1992) initially proposed a prescription based only on phenomenological arguments. Next, Maeder (2003) derived a prescription based on the evaluation of the dissipation of the energy contained in a horizontal shear. Finally, Mathis et al. (2004), derive a prescription based on results observed for turbulent transport in a nonstratified Taylor–Couette experiment (Richard & Zahn 1999). On the other hand, theoretical works (e.g. Billant & Chomaz 2001; Kitchatinov & Brandenburg 2012) and numerical simulations (Waite & Bartello 2006; Kitchatinov & Brandenburg 2012) in fundamental and astrophysical fluid dynamics have been devoted to characterize key properties of anisotropic turbulent flows in rotating stably stratified media, such as velocities and length scales in the vertical and horizontal directions. This is a good motivation to also consider the horizontal transport induced by 3D turbulent motions sustained by the instability of the radial differential rotation, which has been ignored up to now.
The turbulent diffusivities (viscosities) and their source.
In this work, we derive a new prescription for the corresponding turbulent transport coefficients (i.e. ν_{h,v} and D_{h,v}). This allows us to propose for the first time an expression that depends explicitly on stratification, rotation, and their ratio. In Sect. 2, we generalize the spectral model introduced by Kitchatinov & Brandenburg (2012) to study the anisotropy of turbulent flows in differentially rotating stably stratified layers and we derive scaling laws allowing us to establish our new prescriptions. In Sect. 3, we implement these scaling laws in the stellar evolution code STAREVOL. We study their impact on the rotational and structural evolution of a 1.0 M_{⊙} star during its premain sequence and mainsequence and on the subgiant and giant phases of a 1.25 M_{⊙} star. Finally, we discuss results and perspectives of this work in the conclusion (Sect. 4).
2. Turbulent transport in stably stratified differentially rotating layers
2.1. Theoretical background
2.1.1. Stably stratified fluids
Stably stratified stellar radiation zones are supposed to be the place of a mild mixing induced by rotation that deeply impacts the evolution of stars (e.g. Zahn 1992; Maeder 2009). Among the physical processes at the origin of such secular transport of chemicals and angular momentum, shearinduced turbulence plays a key role (Knobloch & Spruit 1982; Zahn et al. 1983; Zahn 1992; Talon & Zahn 1997; Mathis et al. 2004).
In this context, the impact of stable stratification on turbulent motions that develop in radiative regions must be understood in detail and carefully evaluated. Indeed, the motions become highly anisotropic because of the potentially strong buoyancy force along the direction of the entropy and mean molecular weight stratifications (e.g. Riley & Lelong 2000; Billant & Chomaz 2001; Brethouwer et al. 2007; Davidson 2013; Marino et al. 2014, and references therein). It is quantified by the Brunt–Väisälä frequency N, defined by
where we consider only the thermal part with g the gravity, H_{ P } = dr/dlnP the pressure scale height, and P the pressure. We introduce the logarithmic gradients ∇ = dlnT/dlnP with T the temperature, and ∇_{ad} = (∂lnT/∂lnP)_{ad}, and the thermodynamic coefficient δ = −(∂lnρ/∂lnT)_{ P, μ }, where ρ is the density and μ the mean molecular weight. The component of the buoyancy frequency related to μgradients is here ignored as a first step. Buoyancy force reduces the amplitude of the vertical displacement of turbulent eddies along the direction parallel to the stratification (hereafter denoted ∥) while no restoring force is applied along the perpendicular (generally horizontal) direction (hereafter denoted ⊥), if ignoring the action of the Coriolis acceleration (see Sect. 2.2.2). This anisotropy induces turbulent velocities u_{∥} and u_{⊥} in the directions parallel and perpendicular to the stratification respectively. Moreover, flows become organized as horizontal “pancakelike” structures with characteristic vertical and horizontal length scales, l_{∥} and l_{⊥} (see Fig. 1). When occurring in stellar radiation zones, such anisotropic turbulence can have important consequences (Zahn 1992). Indeed, because of the action of the buoyancy force in the vertical direction and of the lack of restoring force in the horizontal one, it would lead to a stronger transport of angular momentum and chemicals horizontally than vertically. This may lead to only weak horizontal gradients of angular velocity (Ω), entropy (and temperature) and mean molecular weight (Chaboyer & Zahn 1992). As proposed by Zahn (1992), rotation then can become “shellular” with related simplifications of 2D transport equations (Zahn 1992; Maeder & Zahn 1998; Mathis & Zahn 2004) that have been successfully implemented in several stellar evolution codes: the Geneva code (e.g. Talon et al. 1997; Meynet & Maeder 2000; Eggenberger et al. 2008; Ekström et al. 2012), the code STAREVOL (e.g. Palacios et al. 2003; Decressin et al. 2009; Charbonnel & Lagarde 2010; Amard et al. 2016), the code CESTAM (e.g. Marques et al. 2013), and the code FRANEC (e.g. Chieffi & Limongi 2013).
Fig. 1. Anisotropic turbulent transport in a stably stratified stellar radiation zone with a buoyancy frequency N rotating with an angular velocity Ω. We introduce the characteristic velocities (u_{∥}, u_{⊥}) and length scales (l_{∥}, l_{⊥}) of a turbulent “pancake” in the direction of the entropy stratification and in the horizontal one respectively and usual spherical coordinates (r, θ). 
The control parameter of such dynamics is the ratio between the vertical and horizontal turbulent transport coefficients
and
where τ is a dynamical time scale characterizing the turbulence and its source, when assuming a diffusive description for each of them (Chaboyer & Zahn 1992). Such an assumption is justified when the effect of the instabilities in their nonlinear regime is to cancel their source (Richard & Zahn 1999; Mathis et al. 2004). Then, their ratio can be written
where we introduce ε = l_{∥}/l_{⊥}. From now on to Sect. 2.2.2, we lighten notations by using {D_{v},D_{h}} that stands for {D_{v, v},D_{h,v}} or {D_{v,h},D_{h,h}}, depending on the source of turbulence studied, namely the vertical and horizontal shear, respectively.
To go further, scaling laws can be obtained using different approaches. Following Billant & Chomaz (2001) (see also Davidson 2013), we can introduce the horizontal Froude number of the flow
It compares the relative strength of the inertia associated with the horizontal advection of the horizontal turbulent motions and buoyancy^{1}. It can also be interpreted as the ratio between the time associated to internal gravity waves (τ_{IGWs} = 1/N) and the turbulent characteristic time τ. In their approach, Billant & Chomaz (2001) consider u_{⊥} and l_{⊥} as fixed quantities. The vertical characteristic velocity (u_{∥}) and length scale (l_{∥}) result from the dynamics of the stably stratified flow. Using dimensionless quantities in the Navier–Stockes, continuity, and heat transport equations (in the inviscid and adiabatic limits) describing the dynamics of nonrotating stably stratified flows (we refer the reader to the Eqs. (7)–(10) of their article), they showed that a second vertical Froude number
is also a characteristic number of the system. It compares the inertia of the vertical advection of the horizontal turbulent motions and buoyancy. Making an asymptotic expansion of these equations in the regime of strong stratification for which Fr_{⊥} < < 1, they identified a selfsimilarity that leads to
These scaling laws have been verified in highresolution nonlinear Cartesian numerical simulations computed by Brethouwer et al. (2007) for high Reynolds numbers.
This implies that
Therefore, as expected, the anisotropy becomes larger when stable stratification is stronger because of the corresponding vertical buoyancy restoring force. Thus, D_{v}/D_{h} → 0 (or D_{h}/D_{v} → ∞) when N → ∞. However, stars are potentially rapidly rotating, for example, young solartype stars and intermediatemass and massive stars (e.g. Meynet & Maeder 2000; Gallet & Bouvier 2013, 2015). Moreover, their angular velocities vary over several orders of magnitude along their evolution. It is thus mandatory to understand and quantify the effect of the Coriolis acceleration, which will also constrain turbulent motions, particularly in the horizontal direction.
2.1.2. Effects of rotation
The effects of the stable stratification in the vertical direction on turbulent flows being now introduced, the balance in the horizontal direction and the effects of rotation through the Coriolis acceleration should be examined. Indeed, if turbulent motions induced by the vertical shear instability are considered, they are intrinsically 3D and they induce turbulent transport and mixing both in the vertical direction, where they are submitted to the buoyancy force, and in the horizontal direction, where the Coriolis acceleration is the restoring force.
In his seminal paper, Zahn (1992) already pointed out the key role of rotation to control the anisotropy of turbulent transport in stably stratified stellar radiation zones. He discussed the transition between 3D isotropic to 2D rotationally controled turbulent flows (see e.g. Hopfinger et al. 1982, for experimental studies) and its consequence for the vertical turbulent transport. In this context, he also showed how the instability of an horizontal shear contributes to this transport in the radial direction, and he derived the prescriptions for the eddyviscosity ν_{v,h} (we refer the reader to Eqs. (2.19), (2.20), (2.22), and (2.23) in Zahn 1992). This already points out how the 3D turbulence induced by the instability of a shear varying along a given direction (radial or latitudinal) sustains turbulent transport both in the vertical and horizontal directions. Therefore, one must take into account the contributions of both vertical and horizontal shear instabilities when studying the horizontal (and vertical) turbulent transport in stellar radiation zones.
In addition, several works in fundamental fluid mechanics examined the strength of the anisotropy of turbulence in stably stratified rotating flows as a function of their buoyancy frequency (N) and angular velocity (Ω). Depending on the rotation rate, Billant & Chomaz (2001) propose that the anisotropy is driven either mainly by stratification (in the regime of “slow” rotation) or by the combined action of stratification and rotation (in the regime of “rapid” rotation). The physical parameter that controls the transition from one regime to the other is the dimensionless Rossby number
which quantifies the relative strength of advection in the horizontal direction to the Coriolis acceleration. In the slowly rotating regime (Ro ^{h} → ∞), their formalism recovers results discussed in the previous section when ignoring rotation. In the rapidly rotating regime (Ro ^{h} < < 1), which is of interest for stellar radiation zones^{2}, they find that
We here recognize the internal Rossby deformation radius in the case where the full Coriolis acceleration is taken into account (Dellar & Salmon 2005, see also Pedlosky 1982). It characterizes the restoring action of the Coriolis acceleration in the horizontal direction, along which rotation limits the size of turbulent stratified pancakes with l_{⊥} decreasing with increasing Ω. Therefore, rotation competes with stratification, which sustains the anisotropy of turbulent flows with l_{⊥} that increases with N. These theoretical scaling laws have been confirmed by Waite & Bartello (2006) who computed highresolution nonlinear numerical simulations in Cartesian geometry where they explored the variations of l_{∥} as a function of Ro^{h}. This proposed scaling, in conjunction with Eq. (4), leads to
However, these studies do not take into account a largescale shear, the key physical ingredient of stellar radiation zones we are studying here.
In this framework, the results obtained by Kitchatinov & Brandenburg (2012) are of great interest. Indeed, although they focused in their work only on the case of uniform rotation as in Billant & Chomaz (2001), they derived a general analytical spectral formalism that allows one to derive the response of a stratified rotating fluid to a preexisting source of turbulence and the anisotropy of the induced turbulent transport. In the cases they studied, they compared the predictions of this turbulent model to direct highresolution nonlinear numerical simulations and showed that the results agree very well. This strongly motivates a generalization of their work by taking into account simultaneously a largescale shear, stable stratification and rotation to study the problem of anisotropic turbulent transport induced by shear instabilities in stellar radiation zones. This is what we propose below.
2.2. Anisotropic turbulent transport in stably stratified rotating radially sheared flows
2.2.1. Spectral formalism
In this section, we generalize the formalism of Kitchatinov & Brandenburg (2012) in the case of a sheared, rotating, stably stratified flow. We intend to determine the response of the fluid to a preexisting source of turbulence characterized by a velocity field u^{0}. We write fluid equations using the Boussinesq approximation where we assume that spacescales characterizing turbulent motions are smaller than those of the variations of background quantities.
Moreover, following Kitchatinov & Brandenburg (2012), a relaxation approximation is used to remove nonlinear terms. It consists of assuming that the fluid tends to reach a steady state with a typical timescale τ. Mathematically, it reads
for any quantity X. Another simplification is that the background turbulence forces the flow with the same time scale, so that the forcing term can be written u_{0}/τ.
In the presence of a shear S = rsinθ ∇Ω, the flow is thus described by the continuity equation ∇ ⋅ u = 0 and the momentum and heat transport equations for adiabatic motions, where diffusion and viscosity are neglected as a first step:
Here, s^{′} denotes the specific entropy fluctuations and ⟨s⟩ the mean specific entropy. We also introduced the unitvector basis {e_{j}}_{j = r,θ,φ} in spherical coordinates. Following Zahn (1992), here we assume a priori that D_{v} < < D_{h} leading to a shellular differential rotation and corresponding shear which mainly depends on the radial coordinate. We thus consider the case of a purely radial shear S = Se_{r}. This hypothesis will be verified in Appendix B.
After combining the momentum and heat transport equations and eliminating the pressure fluctuations thanks to the continuity equation, one obtains in the Fourier space
where ũ is the Fourier transform of u and k̂ = k/k is the normalized wave vector. This can be written as a matrix equation
with
where σ = k̂ · Ω/Ω, Ω^{*} = 2τΩ, S^{*} = τS, N^{*} = τN, and δ_{ij} and ε_{ilj} are the usual Kronecker and LeviCivita symbols respectively. The matrix equation can be solved, using
where
and μ = k̂_{r}. Some terms of ℳ^{−1} vanish when applied to ũ^{0}, because the latter also satisfies the continuity equation . This finally leads to
with
We introduce the spectral tensor Q̃^{0}, which characterizes the properties of the background turbulence. It is defined by
where ⟨⟩ denotes a statistical average. Then, it is possible to derive the spectral tensor of the generated turbulence Q̃ by using the relation
In the case of purely horizontal random motions, which are an idealized view of turbulent displacements in a strongly stratified medium, we have
where E(k) is the kinetic energy spectrum function defined by
The full expression for the tensor Q̃ is too complex to be shown here. Instead, we have focussed on the components needed to determine the anisotropy of the transport using Eq. (4). We have
Assuming that N^{*} ≫ 1 (i.e. N ≫ τ^{−1}) and N^{*} ≫ Ω^{*} (i.e. N ≫ 2Ω as expected in the bulk of stellar radiation zones), these expressions can be simplified by retaining only the dominant terms. In this regime, Eq. (19) becomes
which show no explicit dependence on the shear number S^{*}. Therefore, we recover in the case of a radial differential rotation the expressions found by Kitchatinov & Brandenburg (2012) in the case of solidbody rotation.
Turbulent velocity correlations, which allow us to evaluate transport properties, are finally obtained thanks to the relation
We define the two angles α and β such that
The infinitesimal wave vector element can be written
where k is the norm of k, and we now have
In the previously introduced notations, one has and . Using Eq. (4), we finally obtain
As expected, the strength of the horizontal turbulent transport relatively to the vertical one is thus increasing with stratification (as N^{4}). On the other hand, it decreases with rotation (as Ω^{−2}) since rotation is limiting horizontal turbulent motions through the Coriolis acceleration.
This result constitutes a new key to understanding the anisotropic turbulent transport in stellar radiation zones. If we are able to compute the turbulent transport in the vertical (horizontal) direction induced by a given instability, it allows us to deduce the induced transport in the horizontal (vertical) one as a function of stratification and rotation. This strongly improves the current state of the art with respect to prescriptions previously proposed for D_{h} ≡ D_{h,h} (Zahn 1992; Maeder 2003; Mathis et al. 2004) for the turbulent transport induced by the instability of the horizontal shear, that were based on phenomenological laws and which did not take the simultaneous action of the stratification and the Coriolis acceleration into account.
However, the time τ that characterizes the turbulence is still a free parameter in Eq. (41), and should be specified. The only constraint proposed by Kitchatinov & Brandenburg (2012) is that it should be much larger than the time characterizing the buoyancy τ_{N} = 1/N.
2.2.2. Characterization of the turbulent timescale
We propose here three physically motivated possibilities to estimate the turbulent timescale: (i) the time characterizing the radial shear, τ = 1/S, where S = rsinθ∂_{r}Ω; (ii) a time characterizing the Coriolis acceleration of a rotating radially sheared flow, τ = 1/(2Ω + S); and (iii) the time associated to the epicyclic frequency , which is one of the frequencies characterizing differentially rotating flows and their stability, τ = 1/N_{Ω}.
First, we consider the case of unstable nonrotating vertically sheared flows. The results of numerical simulations of shear instability without rotation (Prat & Lignières 2014) suggest that the typical turbulent timescale, which corresponds to the turnover time of large eddies, is of the order of 1/S in the nonrotating strongly stratified regime. Therefore, τ = 1/S seems to be a reasonable choice when the dynamics is dominated by the shear.
Moreover, this choice allows one to recover the model of vertical turbulent transport by Zahn (1992) using arguments on timescales. Indeed, the fact that turbulence can be sustained only if the dynamical timescale τ is smaller than the timescale of the process that inhibits turbulence in the vertical direction, here the stratification modified by thermal diffusion, can be written
where K is the thermal diffusivity and l_{∥} the typical turbulent length scale along the stratification direction. The scales that induce the most important transport are the largest unstable scales given by
The turbulent transport along the vertical direction is then of the order of , which yields
Choosing τ = S^{−1} finally implies
which is exactly the same model as the one by Zahn (1992). We recall that this turbulent transport occurs if and only if the turbulent Reynolds number (i.e. the ratio between the viscous time and the dynamical time of turbulence) R_{e} = (u_{∥} l_{∥})/ν, where ν is the viscosity of the fluid, is larger than a critical value R_{e;c} (Zahn 1992; Talon & Zahn 1997; Prat et al. 2016). In the STAREVOL code, we chose to set R_{e;c} = 10 that corresponds to a value of D_{v, v} of the order of ten times the molecular viscosity.
When the dynamics is no longer dominated by the shear because of rotation, τ = 1/(2Ω + S) becomes a plausible choice. Indeed, the frequency 2Ω + S characterizes the Coriolis acceleration along the azimuthal direction in a rotating sheared flow (e.g. Mathis 2009). Moreover, it allows one to recover τ = 1/S in the asymptotic limit where Ω ≪ S. In contrast, when the shear is negligible, τ ≃ 1/(2Ω), which is the characteristic time scale of inertial waves. This can be expected since inertial waves are known to structure turbulence in rapidly rotating flows because of their nonlinear interactions (e.g. Sen et al. 2012).
Another frequency that characterizes differentially rotating fluids and their stability, the epicyclic frequency (e.g. Pringle & King 2007)
with ϖ = rsinθ, can also be considered with the corresponding characteristic time τ = 1/N_{Ω}. In the limit where the shear is negligible, one obtains τ ≃ 1/(2Ω), as with τ = 1/(2Ω+S). In the opposite limit where rotation is negligible, , which is not compatible with the two other choices (τ = 1/S, 1/(2Ω+S)), since τ still depends on Ω.
We should here point out that one will be able to make a choice between τ = 1/(2Ω+S) and 1/N_{Ω} only when direct numerical simulations of the vertical shear instability in rotating stably stratified fluids will be available for stellar regimes. This is the reason why here we study and discuss both cases. Another important point to keep in mind when choosing one of these characteristic times is to verify that the rotation profile is not subject to the Rayleigh–Taylor instability for which . In practice, this is generally not a problem during the mainsequence since at this phase. However, for negative values of , one would need to model the effect of the RayleighTaylor instability itself and the turbulent transport it induces (Maeder et al. 2013).
Finally, if one considers the choice τ = 1/N, even though Kitchatinov & Brandenburg (2012) assumed it was not the case, one obtains an anisotropy which is consistent with the work by Billant & Chomaz (2001). However, 1/N characterizes the stratification that inhibits the turbulence while the nonlinear timescale should rather be given by the process that drives the turbulence, that is, the shear. Therefore, we rule out this possibility in the following study.
2.2.3. Horizontal turbulent transport
Using Eq. (41), it is finally possible to derive D_{h} when knowing D_{v} (or viceversa):
for a given choice of turbulent source and characteristic time τ.
In stellar radiation zones, the main important source of turbulence in the vertical direction is the secular vertical shear instability (Zahn 1992; Talon & Zahn 1997). This instability induces 3D turbulent motions both in the vertical and in the horizontal directions (e.g. Prat & Lignières 2013; Garaud & Kulenthirarajah 2016). The horizontal turbulent diffusion coefficient computed using Eq. (47), when assuming that D_{v} is associated with this secular radial shear instability (i.e. D_{v} ≡ D_{v, v} and thus D_{h} ≡ D_{h,v}) and nonzero (for R_{e} > R_{e;c}), thus corresponds to the horizontal turbulent transport induced by 3D anisotropic turbulent motions this vertical instability sustains because of the stable stratification and rotation.
It is important to point out that if the horizontal shear becomes large enough to be simultaneously unstable, it will induce a horizontal turbulent transport with the corresponding turbulent diffusion coefficient (D_{h,h}, see Mathis et al. 2004) as well as a vertical transport (and corresponding vertical turbulent diffusion coefficient D_{v,h}). This is because of the 3D anisotropic turbulent motions this horizontal instability would sustain. This was already pointed out by Zahn (1992) (in his Sect. 2.4.2), even if it has never been considered in stellar evolution studies.
In this work, we have chosen to focus on the first of these two cases. To describe the turbulent transport along the radial direction due to the secular radial shear instability, we adopted the model derived by Zahn (1992):
which has been recently examined and validated by independent highresolution Cartesian nonlinear numerical simulations (Prat & Lignières 2013; Prat et al. 2016; Garaud & Kulenthirarajah 2016). The horizontal turbulent diffusion coefficient associated with the horizontal turbulent transport induced by 3D motions that the vertical instability sustains when R_{e} > R_{e;c}, D_{h,v}, is derived using Eq. (47).
The first possible choice for τ, namely τ = S^{−1}, leads to
In this case, we thus obtain a discontinuous variation for D_{h,v} if the Reynolds criterion is not met and the vertical shearinduced turbulence vanishes. This may be regularized once a prescription for D_{v, v} that will take into account the action of the Coriolis acceleration will be derived. Since we expect the Coriolis acceleration to stabilize the flow, the transport predicted by Eq. (49) can be considered as an upper limit of the actual transport efficiency.
We then examined the two choices for characteristic turbulent timescales that already include rotation. First, assuming a shellular rotation (see Appendix B), the expression
leads to
where X = r∂_{r}Ω/(2Ω). For the implementation in 1D stellar evolution codes, we take a latitudinal average of this expression
where .
Finally, considering the epicyclic characteristic time
we obtain
The horizontal averages in Eqs. (52) and (54) are analytically computed in the Appendix A to allow a robust numerical implementation.
Interestingly, all the choices for τ lead to the same dependence of D_{h,v} on the ratio N/(2Ω). Including rotation in the prescription for τ also allows us to regularize the behaviour of D_{h,v} for a vanishing vertical shearinduced turbulence. This is a significant improvement towards a consistent description of anisotropic turbulence in stellar radiation zones. As already pointed out above, the next step requires a model for D_{v, v} that takes the action of the Coriolis acceleration into account. This should be achieved with direct numerical simulations of the vertical shear instability in rotating stably stratified fluids in stellar regimes, which are beyond the scope of this work, and which will allow us to test our analytical prescriptions. In this framework, it is important to take also the action of the chemical stratification (through the μgradients) and of the viscosity into account as in Prat et al. (2016). Finally, it is necessary to take into account the possible feedback of the horizontal turbulence on the vertical turbulent transport (Talon & Zahn 1997). In the following section, we investigate the impact of the different expressions we derived for D_{h,v} on the transport of angular momentum along the evolution of lowmass main sequence (MS), subgiant and red giant stars.
3. Application to stellar interiors and evolution
3.1. Method
We implemented Eqs. (48), (49), (52), and (54) in the stellar evolution code STAREVOL as an additional source of turbulence in the horizontal direction. The total horizontal turbulent diffusion coefficient D_{h} is then defined as the sum of D_{h,v} as given by these equations (we test the different expressions for τ) and of D_{h,h} from Mathis et al. (2004; Eq. (19)). The treatment of rotationinduced processes in STAREVOL is described in detail, for example in Decressin et al. (2009) and Amard et al. (2016), and we refer to this latest paper for the information on the basic input physics (i.e., equation of state, nuclear reactions, opacities, etc.) used here (see also Siess et al. 2000; Palacios et al. 2003; Lagarde et al. 2012). The angular velocity profile in the radiative interior derives from the redistribution of angular momentum by shearinduced turbulence and meridional flows, which are treated using the formalism derived by Zahn (1992), Maeder & Zahn (1998), and Mathis & Zahn (2004); uniform rotation is assumed in convective regions. Angular momentum losses at the stellar surface due to pressuredriven magnetized stellar winds are taken into account using the prescription by Matt et al. (2015). Convective regions are modelled according to the mixing length theory. The adopted mixinglength parameter is defined by the calibration of a standard solar model and is fixed to α_{MLT} = 1.973. We adopted Z = 0.013446 as the initial metallicity that corresponds to the solar metallicity according to Asplund et al. (2009).
We explored the impact of our new modelling of the shearinduced horizontal turbulence on the transport of angular momentum in lowmass stars along the MS and the red giant branch. We start with the case of a 1 M_{⊙}, Z_{⊙} model for which we have constraints both on the internal rotation rate at the age of the Sun and on the evolution of the surface rotation with time (Sect. 3.2). We then moved to the case of a 1.25 M_{⊙}, Z_{⊙} model which internal and surface rotation states are constrained by asteroseismic data from mixed modes of subgiant and giant stars (Sect. 3.3).
3.2. Case of a 1.0 M_{⊙}, Z_{⊙} main sequence star
We focus on a 1.0 M_{⊙} star at solar metallicity and present a comparative analysis of the evolution of the different diffusion coefficients and meaningful quantities as a function of time along the MS evolution of this lowmass star. For this, we compared several models computed from the premain sequence (PMS) up to 9 Gyr with the different expressions for τ (i.e., with the different expressions for the horizontal eddy diffusivity D_{h,v} given in Eqs. (49), (52), and (54), added to the horizontal turbulent diffusivity D_{h,h} derived in Mathis et al. 2004), keeping all the other ingredients unchanged. These models are also compared with a model that ignores D_{h,v} and includes only the prescription for D_{h,h} given by Mathis et al. (2004).
3.2.1. Model setup
We used similar assumptions for the initial rotation rate and prescription for magnetic braking as in Amard et al. (2016). A phase of disc coupling, during which the surface rotation rate is maintained to its original value of Ω_{ini} = 1.6 × 10^{−5} s^{−1} (corresponding to an initial rotation period of 4.55 days), is accounted for during the first 5 Myr of the evolution on the PMS. This corresponds to a median rotator, as observed in statistical samples of PMS stars in young clusters (see Amard et al. 2016 and Gallet & Bouvier 2015 for details). Four median rotator models are thus computed, three with the different expressions for D_{h,v} added to D_{h,h} from Mathis et al. (2004), and the last one with only D_{h,h}. For the case where D_{h,v} is computed with τ = 1/S (Eq. (49)), we also compute two models with different initial velocities to investigate the effect of the global velocity on the turbulent transport; these fast and slow rotators start with an initial period of 1.4 and 9.0 days, respectively. For all the models angular momentum losses at the stellar surface due to pressuredriven magnetized stellar winds are taken into account using the prescription by Matt et al. (2015) with K_{wind} = 7 × 10^{30} to fit the solar rotation rate at 4.57 Gyr, m = 0.22 as expected from a purely dipolar magnetic field and the parameter p = 2.4 to fit the surface rotation rates in open clusters at best. We then cover the statistical distributions of rotation periods in open clusters and associations from 1 Myr to 2.5 Gyr (Gallet & Bouvier 2013; Gallet & Bouvier 2015, Amard et al. 2016). We summarize in Table 2 the specifications of the different 1 M_{⊙}, Z_{⊙} models that we compare below.
Specifications of the different 1 M_{⊙}, Z_{⊙} models.
3.2.2. Comparing different prescriptions for the horizontal turbulent diffusivity
In Fig. 2, we compare the profiles of the total horizontal diffusivity D_{h} and of the vertical diffusivity D_{v, v} within all the 1.0 M_{⊙}, Z_{⊙} models described in Table 2 at three different epochs on the MS (zero age main sequence, age of the Sun, and when the central hydrogen mass fraction is reduced to 10% of its original value). We note first that below ∼0.5 R_{⋆}, where R_{⋆} is the radius of the star, the novel expressions for the horizontal turbulent diffusion (D_{h,v}) lead to values for the total D_{h} = D_{h,v} + D_{h,h} that are very close to the ones given by the expression of D_{h,h} by Mathis et al. (2004) alone. The main reason for this behaviour is that the differential rotation needs to be high enough to trigger on the vertical shear instability (i.e. R_{e} > R_{e;c}), and this is not the case in the central part of the star for most of the MS. Therefore, the horizontal component of D_{h} (D_{h,h}) given by Mathis et al. (2004) dominates the transport of angular momentum along the isobars, all along the MS in the deep interior. This can be seen clearly by looking at the behaviour of D_{h,v}, which becomes important only in the outer part of the radiative region above ∼0.5 R_{⋆} where differential rotation is stronger (and able to trigger the vertical shear instability), as imposed by the extraction of angular momentum by the magnetized wind. At the age of the Sun and up to the turnoff, D_{h,v} can be several orders of magnitude higher than D_{h,h} in these outer layers below the convective envelope.
Fig. 2. Top row: total horizontal turbulent diffusivity D_{h} (full line) and its D_{h,v} component (when active if R_{e} > R_{e;c}; dotted lines) for the 1.0 M_{⊙}, Z_{⊙} models described in Table 2 at three evolutionary stages (from left to right: ZAMS, age of the Sun, and when the central H mass fraction is reduced to 10% of its original value). The hatched areas correspond to the convective regions. Bottom row: vertical turbulent diffusion coefficient in the same models. 
Figure 2 also shows the profiles of the vertical turbulent diffusivity D_{v} = D_{v, v} at the same evolutionary stages for all the prescriptions. The strong coupling at the ZAMS obtained in the fast rotator (light green curve, model A2) is due to a very efficient transport by the meridional circulation, and results in a weak vertical turbulent transport as can be seen in the corresponding bottom left panel. As the models evolve on the MS a flattening of the D_{v} profile appears below the convective envelope (above r ≈ 0.5 R_{⋆}) for models A1, A2, A3, B and C. The vertical shearinduced turbulent transport seems to be weaker in average around log(D_{v})≈2.5 in this region when the radial differential rotation is accounted in addition to the horizontal shear as a source for the horizontal turbulent transport (i.e. when adding D_{h,v} to the D_{h,h} of Mathis et al. 2004). The radial variations of D_{v} can also be understood by the radial differential rotation obtained for r > 0.5 R_{⋆} reported in Fig. 5.
On the upper panel of Fig. 3 we show the anisotropy ratio of the turbulence given by Eq. (47). We note that this representation actually corresponds to the potential anisotropy of turbulence since we plot it even in the vertical shearstable regions. We can distinguish two behaviours with the new prescriptions: the case of models B and C where the turbulent timescale depends on the rotation rate on one hand, and the case of models A where such a dependence is not present (τ = 1/S) on the other hand. While both groups show an anisotropy ratio that increases moderately during the evolution, in the first group of models the anisotropy is relatively low in the central region due to the higher angular velocity there, while in the second case the anisotropy is strong in the central region. Quantitatively during the MS, the region below the envelope shows an anisotropy ratio that is about 10^{9}. This means that whenever the vertical shear is high enough to overcome the Reynolds criterion (i.e. R_{e} > R_{e;c}) below the convective envelope for r > 0.5 R_{⋆}, the total horizontal diffusion coefficient log(D_{h}) rises from 7.5 up to 11.5 ( ≃ 10^{9} × D_{v}) that leads to a stronger meridional circulation and transport of angular momentum as explained in Appendix C.
Fig. 3. Top row: anisotropic ratio given by Eq. (47) for the different choices of τ. It represents the potential anisotropy of turbulence sustained by the vertical shear instability in the stably stratified differentially rotating medium. It is thus equivalent to the ratio between the horizontal and vertical turbulent diffusivities (D_{h,v}/D_{v, v}) when the latter is active. Bottom row: product of the thermal part of the Brunt–Väisälä frequency (N) by the computed turbulent characteristic time (τ). Both are shown at the same evolutionary stages and with the same colour code as in Fig 2. 
On the lower panel of Fig. 3 we show the product of the Brunt–Väisälä frequency N (its thermal part given in Eq. (1)) by the turbulent characteristic timescale τ for models A, B and C at the same evolution points on the MS. As initially assumed by Kitchatinov & Brandenburg (2012), this quantity turns out to be greater than unity all along the MS evolution. It stays quite constant during this latter and shows very similar values in models B and C. Indeed, rotation tends to compensate for the disappearance of the shear, and even leads to a decrease of the turbulent timescale, thus limiting the increase of the product Nτ. For models A though, the product Nτ is high in the central region where the shear is weaker at early stages (see below, Fig. 5) and the Brunt–Väisälä frequency is important.
Finally, we show the effective diffusivity (Fig. 4) that is also affected by the change of horizontal turbulent diffusion coefficient and the potentially enhanced meridional circulation (see Appendix C), consistently with what was already pointed out in Mathis et al. (2007), the erratum associated with Mathis et al. (2004). By definition, , where U_{2}(r) is the radial function of the expansion of the radial component of the meridional circulation on Legendre polynomials (Chaboyer & Zahn 1992). In the case of a strong horizontal turbulence, one can see in Decressin et al. (2009; we refer the reader to Eqs. (32) and (40) in this article) that entropy advection by the meridional circulation balances the horizontal turbulent heat diffusion. Then, we have U_{2} ∝ D_{h} and D_{eff} ∝ D_{h}. This explains both the amplitude and the radial variations of D_{eff} observed for r > 0.5R_{⋆} in models A, B, and C.
Fig. 4. Effective diffusivity D_{eff} in the vertical direction induced by the advective transport at the same evolutionary stages and with the same colour code as in Fig 2. 
3.2.3. Impact on the evolution of the internal and surface rotation
We now discuss the impact of introducing the new formalism for D_{h} on the evolution of internal rotation and on the surface rotation rates of the 1.0 M_{⊙}, Z_{⊙} models. First, we verify and confirm in Appendix B that the hypothesis of a shellular rotation are respected, that is, D_{v} < < D_{h} and . Next, Fig. 5 shows the internal rotation profiles and the differential rotation rate (defined as −rΩ^{′}/2Ω, where ^{′} denotes the radial derivative d/dr) for all the models at the same ages as previously (upper and lower panels respectively). Given the results on the total horizontal diffusivity D_{h} presented in Sect. 3.2.2, adding the vertical source for horizontal turbulence D_{h,v}, has globally little impact on the rotation profiles, which remain strongly differential during the overall MS evolution. The ZAMS almost flat profile exhibited by model A2 (fast rotator, light green) reflects a strong coupling at the ZAMS for the fast rotators, that is nonetheless rapidly wiped out by wind braking later on the MS evolution. At the ZAMS, the only distinction that can be made between the different models is related to their initial spin rate. This is related to the shallow differential rotation at that time. The models keep the rotational properties they have been given initially, meaning that a fast (slow) rotator is still fast (slow) rotator. Moreover, the PMS contraction and the relatively high angular velocity in the early phases where magnetic braking is not yet at its maximum efficiency lead to a very high D_{h,h}. Therefore, the impact of the new component D_{h,v} is negligible at this phase and the four models initialised with the same rotation rate are superimposed. At the age of the Sun and later at the end of the MS, we can no longer distinguish between the different rotators in terms of surface rotation. However, the effects of the different new turbulence prescriptions are now noticeable on the rotation profiles below the convective envelope. When including D_{h,v} (if active for R_{e} > R_{e;c} for R > 0.5 R_{⋆}), a slightly larger extraction of angular momentum is observed below the envelope because of the more efficient transport from the inside of the star to the surface driven by the enhanced meridional circulation sustained by the additional horizontal turbulent diffusion (Appendix C). As seen on the bottom panel of Fig. 5, the external layers of the radiative core are the place where the radial differential rotation rate is the more important. This is due both to the radius increasing and to the torque applied at the surface by the magnetised stellar wind. The vertical shear instability is then more able to develop and contributes to the horizontal turbulent diffusion coefficient (see Fig. 2). Also, one can note that the rotation profile is independent of the chosen characteristic timescale for turbulence (τ). Indeed, models A, B, and C lead to very similar rotation profiles at the age of the Sun and at the end of the MS.
Fig. 5. Top row: angular velocity profile represented at the same evolutionary stages and with the same colour code as in Fig 2. Bottom row: ratio between the differential rotation and the rotation rate, equivalent to the dimensionless parameter X in Eqs. (49), (51), and (54). 
The rotation profile flatness obtained with Mathis et al. (2004; model D) at the ZAMS is associated to a very strong meridional circulation triggered by the contraction during the PMS. This is especially true in the case of the fast rotating model since in their expression, D_{h,h} depends at the first order on the rotation rate.
The differential rotation obtained with the new prescriptions can also be seen in Fig. 6 where we draw on the top panel the evolution of the normalised core to the surface differential rotation (i.e. the difference between the mean rotation rate of the radiative core, Ω_{rad}, and the surface rotation, Ω_{s}, as defined in Amard et al. 2016 normalised by their sum) as a function of time. A value close to one indicates a very strong gradient of rotation between the radiative core and the convective envelope of the star, while ΔΩ = 0 means that the model rotates as a solid body. This diagram has to be read in parallel to the bottom panel which shows the evolution with time of both the surface angular velocity (gyrotracks) and the mean radiative core angular velocity. Predictions of the models are compared with rotation periods observed in stars with masses between 0.9 and 1.1 M_{⊙} in open clusters over a large age range collected by Bouvier et al. (2014, and references therein). We note that the three models associated to different initial rotation rates reproduce the statistical data points for surface rotation with time fairly well, in a similar way to Amard et al. (2016). Even though the vertical shear turbulence prescription is different, it has little effects on the global transport of angular momentum, hence on the surface rotation.
Fig. 6. Bottom panel: evolution of the angular velocity as a function of time for models computed with the prescriptions for D_{h} from Mathis et al. (2004), Zahn (1992), and that of the present paper with the three different characteristic time scales for turbulence (black, red, green (τ = 1/S), blue (τ = 1/(2Ω + S)) and magenta (τ = 1/N_{Ω}), respectively). Scaled on the left is the rotational period in days and on the right the angular velocity in solar units with Ω_{⊙} = 2.86 × 10^{−6} s^{−1}. Each cross represents a star belonging to an open cluster which age and mass (between 1.4 and 1.6 M_{⊙}) have been taken from literature. All observational data are taken from Gallet & Bouvier (2015) and Bouvier et al. (2014) and references therein. The solid lines show the evolution of the surface rotation, and dashed lines the evolution of the mean radiative core rotation rate. Top panel: evolution of the normalised core to surface differential rotation rate for the six 1.0 M_{⊙}, Z_{⊙} models as colourcoded in Fig. 2. 
The same prescription for the extraction of angular momentum by the wind is used for all the models (see Sect. 3.2). Therefore, the surface and the core rotation rates evolve differently only depending on the way we describe the horizontal turbulence. We see that the global evolutions of both the surface and the core rotation rates are almost unchanged from one model to the other. Indeed, even though the transport generated by the added term to the horizontal turbulent diffusion can be very high in amplitude locally, it is active only when the vertical shear instability is triggered. For a MS solarmass star, this is the case for approximately 10% of the lifetime, and this fraction decreases with the strength of the stellarwind torque applied at the surface that generates radial differential rotation.
While the horizontal turbulent transport induced by 3D motions triggered by the vertical shear instability can have a very high amplitude locally (for instance, below the convective envelope), it is a selfregulated mechanism. Indeed, if the vertical shear is unstable and leads to important D_{h}, an efficient meridional circulation is triggered that weakens the radial differential rotation (see Appendix C) because of the advection of angular momentum towards the surface to values below the threshold needed to sustain the vertical shear instability and the source of the associated strong horizontal turbulent transport thus vanishes. This selfregulation is the cause of the oscillations of the surface angular velocity that appear in the lower panel of Fig. 6 for ages larger that 600 Myrs approximately.
In addition, while the potentially strong horizontal turbulent transport and the induced efficient advection of angular momentum locally flattens the rotation profile, it also steepens the gradient of angular velocity on each side of the turbulent region because of the spatial discretisation that is inherent to any stellar evolution numerical code. This enhances the shear with the neighbouring layers that in turn triggers the turbulence during the next time step. That transient phenomenon is responsible for the radial oscillations that are observed on Figs. 2–5 below the convective envelope on the main sequence when D_{h,v} is locally several orders of magnitude higher than D_{h,h}. A smoother transition would likely be possible with a nonlocal instability criterion (see e.g. Gagnier & Garaud 2018).
As a partial conclusion, additional transport mechanisms such as internal gravity waves, magnetic fields or other instabilities (we refer the reader to the introduction for relevant references) are therefore still necessary to explain the efficient extraction of angular momentum from the core of lowmass stars during their PMS and MS. The case of subgiant and red giant stars, which have been successfully sounded thanks to spacebased asteroseismology, should also be examined.
3.3. Case of red giant stars illustrated with a 1.25 M_{⊙} model
The asteroseismic data obtained for the stellar core rotation from the mixed modes stochastically excited in subgiants and red giants provide very good constraints for the transport of angular momentum in evolved stars (e.g. Beck et al. 2012; Mosser et al. 2012; Deheuvels et al. 2012, 2014, 2015; Spada et al. 2016; Gehan et al. 2018). The observed stars have masses ranging between 1.1 and 1.4 M_{⊙}. We thus focus now on the intermediate mass of 1.25 M_{⊙} for which we evolve models from the PMS to the RGB phase with the same assumptions and prescription for the loss of angular momentum through stellar winds as for the previous 1.0 M_{⊙} models.
We see in Sect. 3.2.3 that the timescale associated with the turbulence (τ) does not change drastically the results obtained for the rotation evolution. Therefore, to simplify the comparison, we only show here two 1.25 M_{⊙} models with D_{h} computed following Mathis et al. (2004) only and the new one assuming that τ = 1/S (same transport as model A)^{3}.
Figure 7 shows the evolution of the surface rotation of these two models compared to the data given by Deheuvels et al. (2014) for a few subgiant stars. The braking is not calibrated to reproduce the surface velocity of such stars, therefore we do not expect to reproduce quantitatively the data. As such, the surface angular velocity of the models is too high. Two effects may be responsible for this: a lack of angular momentum extraction by the wind as well as a lack of transport in the radiative region of the star. While the first seems obvious, the latter needs some explanation: a more efficient transport of angular momentum on the MS would lead to a smaller reservoir of angular momentum when the star begins to expand on the subgiant branch, thus to a smaller angular velocity.
Fig. 7. Evolution of the surface and core rotation rates (solid and dashed lines respectively) as a function of the stellar radius that increases with time (MS, subgiant, and RGB correspond respectively to R_{eff}/R_{⊙} between ∼1.0 and 1.8, 1.8 and 3, and higher than 3). We show theoretical predictions for two 1.25 M_{⊙} models at solar metallicity computed with the prescription for D_{h} = D_{h,h} by Mathis et al. (2004; black) and the new prescription including both components of D_{h} = D_{h,v} + D_{h,h} for the case τ = 1/S (green). We compare the model predictions with relevant data for stars with asteroseismic masses between 1.1 and 1.4 M_{⊙}. For SGB stars we show surface and core rotation rates (black and orange dots) derived by Deheuvels et al. (2014). For RGB stars, we plot the core rotation rates derived by Mosser et al. (2012; red dots). 
In Fig. 7, we also show the angular velocity of the rotating core of the same models. To be able to compare them with asteroseismic data, we have computed the core rotation rate as the mass average of the angular velocity in the resonant cavity of gravitydominated modes, that is in the region below the peak of the chemical Brunt–Väisälä frequency (, with H_{P} the pressure scale height, ∇_{μ} the mean molecular weight gradient, and ϕ = (∂lnρ/∂lnμ)_{P, T}). As for the MS evolution, the rotational tracks are very close. The model with the new prescription shows a core rotating slightly slower as it is expected from a model with a more efficient coupling. However, this difference is far from sufficient to explain the asteroseismic data derived for the subgiants and red giants cores by Mosser et al. (2012). Hence again, other mechanisms such as propagative internal gravity waves (Talon & Charbonnel 2008; Pinçon et al. 2017), mixed modes (Belkacem et al. 2015a,b), or MHD instabilities (Rüdiger et al. 2015) should be invoked to potentially reproduce observations.
4. Conclusion
In this work, we derive new theoretical prescriptions for the anisotropy of the turbulent transport in differentially rotating stellar radiation zones. Extending the theoretical formalism derived by Kitchatinov & Brandenburg (2012) to the case of rotating stably stratified flows with a vertical shear, we find that the ratio between the horizontal and the vertical turbulent transport scales as N^{4}τ^{2}/(2Ω^{2}), where we recall that N and Ω are the buoyancy and rotation frequencies, respectively, and τ is the time characterizing the source of the turbulence. This shows that if anisotropy increases with stratification, the Coriolis acceleration constitutes the restoring force in the horizontal direction, which was ignored in previously derived prescriptions.
We propose here three physicallymotivated expressions for τ: τ = 1/S, 1/(2Ω+S), and 1/N_{Ω}, where S is the shear and N_{Ω} the epicyclic frequency. The first choice corresponds to the model for the turbulent transport induced along the radial direction by the vertical shear instability proposed by Zahn (1992) and validated by recent direct numerical simulations (Prat & Lignières 2013; Garaud & Kulenthirarajah 2016). The second and third choices correspond to the introduction of the influence of rotation on the shearinduced turbulence. Their robustness should be validated in a near future by direct numerical simulations of rotating stably stratified flows with an unstable vertical shear in the range of parameters that corresponds to stellar regimes. These required simulations, which are out of the scope of this work, will also allow us to improve the Zahn (1992)’s turbulent model by taking into account the action of the Coriolis acceleration, of the chemical stratification and viscosity as in Prat et al. (2016), and of the interactions with the horizontal turbulence (Talon & Zahn 1997).
Then, the turbulent transport coefficient in the horizontal direction, induced by the 3D turbulent motions sustained by the vertical sheardriven instability, scales as Ri_{c}(N/2Ω)^{2}Kf(S, Ω), where Ri_{c} and K are the critical Richardson number and thermal diffusivity, respectively, and f is a function that depends on the prescription chosen for τ, when assuming the model proposed by Zahn (1992) for the vertical turbulent transport coefficient. We identify that f = 1 for τ = 1/S. We recall that all the previous prescriptions (Zahn 1992; Maeder 2003; Mathis et al. 2004) have no explicit dependence on the entropy (and chemicals) stratification and do not take the action of the Coriolis acceleration into account, although these are the two restoring forces for turbulent flows in differentially rotating stellar radiation zones.
When applied to complete stellar evolution models of rotating lowmass, solarmetallicity stars, accounting for the feedback of the vertical shear in the horizontal direction (D_{h} = D_{h,h} + D_{h,v}), the new prescriptions do not modify greatly the results previously obtained using the Mathis et al. (2004) prescription (D_{h} = D_{h,h}) for the internal and surface rotation rates. Although the horizontal turbulent transport induced by 3D motions triggered by the vertical shear instability can have a very high amplitude locally (for instance below the convective envelope of lowmass mainsequence stars), it is a selfregulated mechanism. Indeed, if the vertical shear is unstable and leads to important D_{h}, an efficient meridional circulation is triggered that weakens the radial differential rotation to values below the threshold needed to sustain the vertical shear instability and the source of the associated strong horizontal turbulent transport thus vanishes. This also leads to a quenching of the vertical turbulent transport. As a consequence, this mechanism is not able to provide the extraction of angular momentum from stellar cores needed to explain their rotation rates in the Sun and in mainsequence, subgiant and red giant lowmass and intermediatemass stars. A supplementary mechanism such as internal gravity waves or/and magnetic fields and their instabilities should thus be invoked. In this framework, this work demonstrates the great importance of pursuing the efforts to consistently model turbulent transport induced by rotationdriven hydrodynamical instabilities, that have a major impact on the structural and rotational evolution of stars. This strongly motivates future works to improve: (i) the physical description of the radial turbulent transport induced by the vertical shear instability with taking the action of the Coriolis acceleration, the effects of the chemical stratification and of viscosity (Prat et al. 2016), and the interactions with the horizontal turbulence (Talon & Zahn 1997) into account; (ii) the prescriptions for the horizontal and vertical turbulent transports induced by instabilities of the horizontal differential rotation. Their strengths have to be carefully evaluated in order to be able to discuss in a consistent framework the effects of other transport mechanisms of angular momentum such as internal gravity waves and magnetic fields.
See Appendix B.
Acknowledgments
We dedicate this work to the memory of our friend and mentor JeanPaul Zahn, in the steps of whom we walk. Pursuing the path he opened is a great honor. We warmly thank the referees for their careful reading and constructive reports that allowed us to improve our work. S. M. and V. P. acknowledge support by ERC SPIRE 647383 and FP7 SpaceInn grants. S. M., V. P. and A. P. acknowledge support by CNES PLATO funding. C. C. and L. A. thank the Equal Opportunity Office of the University of Geneva. The authors acknowledge financial support from the French Programme National de Physique Stellaire PNPS of CNRS/INSU, the grant ANR 2011 Blanc SIMI56 020 01 Toupies (Towards understanding the spin evolution of stars) and the Swiss National Science Foundation (SNF).
References
 AcevedoArreguin, L. A., Garaud, P., & Wood, T. S. 2013, MNRAS, 434, 720 [NASA ADS] [CrossRef] [Google Scholar]
 Aerts, C., Van Reeth, T., & Tkachenko, A. 2017, ApJ, 847, L7 [NASA ADS] [CrossRef] [Google Scholar]
 Amard, L., Palacios, A., Charbonnel, C., Gallet, F., & Bouvier, J. 2016, A&A, 587, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Barnabé, R., Strugarek, A., Charbonneau, P., Brun, A. S., & Zahn, J.P. 2017, A&A, 601, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Beck, P. G., Montalban, J., Kallinger, T., et al. 2012, Nature, 481, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Belkacem, K., Marques, J. P., Goupil, M. J., et al. 2015a, A&A, 579, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belkacem, K., Marques, J. P., Goupil, M. J., et al. 2015b, A&A, 579, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Benomar, O., Takata, M., Shibahashi, H., Ceillier, T., & García, R. A. 2015, MNRAS, 452, 2654 [NASA ADS] [CrossRef] [Google Scholar]
 Billant, P., & Chomaz, J.M. 2001, Phys. Fluids, 13, 1645 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Bouvier, J., Matt, S. P., Mohanty, S., et al. 2014, Protostars and Planets VI, 433 [Google Scholar]
 Brethouwer, G., Billant, P., Lindborg, E., & Chomaz, J.M. 2007, J. Fluid Mech., 585, 343 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, T. M., ChristensenDalsgaard, J., Dziembowski, W. A., et al. 1989, ApJ, 343, 526 [NASA ADS] [CrossRef] [Google Scholar]
 Ceillier, T., Eggenberger, P., García, R. A., & Mathis, S. 2013, A&A, 555, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chaboyer, B., & Zahn, J.P. 1992, A&A, 253, 173 [Google Scholar]
 Charbonnel, C., & Lagarde, N. 2010, A&A, 522, A10 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Charbonnel, C., & Talon, S. 2005, Science, 309, 2189 [Google Scholar]
 Charbonnel, C., Decressin, T., Amard, L., Palacios, A., & Talon, S. 2013, A&A, 554, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chieffi, A., & Limongi, M. 2013, ApJ, 764, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Davidson, P. A. 2013, Turbulence in rotating, stratified, and electrically conducting fluids (Cambridge University Press) [CrossRef] [Google Scholar]
 Decressin, T., Mathis, S., Palacios, A., et al. 2009, A&A, 495, 271 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Deheuvels, S., García, R. A., Chaplin, W. J., et al. 2012, ApJ, 756, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Deheuvels, S., Doğan, G., Goupil, M. J., et al. 2014, A&A, 564, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Deheuvels, S., Ballot, J., Beck, P. G., et al. 2015, A&A, 580, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dellar, P. J., & Salmon, R. 2005, Phys. Fluids, 17, 106601 [NASA ADS] [CrossRef] [Google Scholar]
 Denissenkov, P. A., & Pinsonneault, M. 2007, ApJ, 655, 1157 [NASA ADS] [CrossRef] [Google Scholar]
 Eggenberger, P., Maeder, A., & Meynet, G. 2005, A&A, 440, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Eggenberger, P., Meynet, G., & Maeder, A. 2008, Ap&SS, 316, 43 [NASA ADS] [CrossRef] [Google Scholar]
 Eggenberger, P., Montalbán, J., & Miglio, A. 2012, A&A, 544, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fossat, E., Boumier, P., Corbard, T., et al. 2017, A&A, 604, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gagnier, D., & Garaud, P. 2018, ApJ, 862, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Gallet, F., & Bouvier, J. 2013, A&A, 556, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gallet, F., & Bouvier, J. 2015, A&A, 577, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Garaud, P., & Kulenthirarajah, L. 2016, ApJ, 821, 49 [Google Scholar]
 Garaud, P., Gagnier, D., & Verhoeven, J. 2017, ApJ, 837, 133 [NASA ADS] [CrossRef] [Google Scholar]
 García, R. A., TurckChièze, S., JiménezReyes, S. J., et al. 2007, Science, 316, 1591 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Gehan, C., Mosser, B., Michel, E., et al. 2018, A&A, 616, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gough, D. O., & McIntyre, M. E. 1998, Nature, 394, 755 [NASA ADS] [CrossRef] [Google Scholar]
 Heger, A., Woosley, S. E., & Spruit, H. C. 2005, ApJ, 626, 350 [NASA ADS] [CrossRef] [Google Scholar]
 Hermes, J. J., Gänsicke, B. T., Kawaler, S. D., et al. 2017, ApJS, 232, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Hirschi, R., & Maeder, A. 2010, A&A, 519, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hopfinger, E. J., Gagne, Y., & Browand, F. K. 1982, J. Fluid Mech., 125, 505 [NASA ADS] [CrossRef] [Google Scholar]
 Kitchatinov, L. L., & Brandenburg, A. 2012, Astron. Nachr., 333, 230 [NASA ADS] [CrossRef] [Google Scholar]
 Knobloch, E., & Spruit, H. C. 1982, A&A, 113, 261 [NASA ADS] [Google Scholar]
 Kurtz, D. W., Saio, H., Takata, M., et al. 2014, MNRAS, 444, 102 [NASA ADS] [CrossRef] [Google Scholar]
 Lagarde, N., Decressin, T., Charbonnel, C., et al. 2012, A&A, 543, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Maeder, A. 2003, A&A, 399, 263 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Maeder, A. 2009, Physics, Formation and Evolution of Rotating Stars (Berlin Heidelberg: Springer) [Google Scholar]
 Maeder, A., & Zahn, J.P. 1998, A&A, 334, 1000 [NASA ADS] [Google Scholar]
 Maeder, A., Meynet, G., Lagarde, N., & Charbonnel, C. 2013, A&A, 553, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marino, R., Mininni, P. D., Rosenberg, D. L., & Pouquet, A. 2014, Phys. Rev. E, 90, 023018 [NASA ADS] [CrossRef] [Google Scholar]
 Marques, J. P., Goupil, M. J., Lebreton, Y., et al. 2013, A&A, 549, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S. 2009, A&A, 506, 811 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., & Zahn, J.P. 2004, A&A, 425, 229 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., & Zahn, J.P. 2005, A&A, 440, 653 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., Palacios, A., & Zahn, J.P. 2004, A&A, 425, 243 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., Palacios, A., & Zahn, J.P. 2007, A&A, 462, 1063 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Matt, S. P., Brun, A. S., Baraffe, I., Bouvier, J., & Chabrier, G. 2015, ApJ, 799, L23 [NASA ADS] [CrossRef] [Google Scholar]
 Meynet, G., & Maeder, A. 2000, A&A, 361, 101 [NASA ADS] [Google Scholar]
 Mosser, B., Goupil, M. J., Belkacem, K., et al. 2012, A&A, 548, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Murphy, S. J., Fossati, L., Bedding, T. R., et al. 2016, MNRAS, 459, 1201 [NASA ADS] [CrossRef] [Google Scholar]
 Palacios, A., Talon, S., Charbonnel, C., & Forestini, M. 2003, A&A, 399, 603 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pedlosky, J. 1982, Geophysical Fluid Dynamics (New York and Berlin: SpringerVerlag) [Google Scholar]
 Pinçon, C., Belkacem, K., Goupil, M. J., & Marques, J. P. 2017, A&A, 605, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Prat, V., & Lignières, F. 2013, A&A, 551, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Prat, V., & Lignières, F. 2014, A&A, 566, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Prat, V., Guilet, J., Viallet, M., & Müller, E. 2016, A&A, 592, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pringle, J. E., & King, A. 2007, Astrophysical Flows (Cambridge University Press) [CrossRef] [Google Scholar]
 Richard, D., & Zahn, J.P. 1999, A&A, 347, 734 [Google Scholar]
 Rieutord, M. 2006, A&A, 451, 1025 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Riley, J. J., & Lelong, M.P. 2000, Ann. Rev. Fluid Mech., 32, 613 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, T. M. 2015, ApJ, 815, L30 [NASA ADS] [CrossRef] [Google Scholar]
 Rüdiger, G., Gellert, M., Spada, F., & Tereshin, I. 2015, A&A, 573, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Saio, H., Kurtz, D. W., Takata, M., et al. 2015, MNRAS, 447, 3264 [NASA ADS] [CrossRef] [Google Scholar]
 Sen, A., Mininni, P. D., Rosenberg, D., & Pouquet, A. 2012, Phys. Rev. E, 86, 036319 [NASA ADS] [CrossRef] [Google Scholar]
 Siess, L., Dufour, E., & Forestini, M. 2000, A&A, 358, 593 [Google Scholar]
 Spada, F., Gellert, M., Arlt, R., & Deheuvels, S. 2016, A&A, 589, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Spiegel, E. A., & Zahn, J.P. 1992, A&A, 265, 106 [Google Scholar]
 Spruit, H. C. 1999, A&A, 349, 189 [NASA ADS] [Google Scholar]
 Strugarek, A., Brun, A. S., & Zahn, J.P. 2011, A&A, 532, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Suijs, M. P. L., Langer, N., Poelarends, A.J., et al. 2008, A&A, 481, L87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Talon, S., & Charbonnel, C. 2005, A&A, 440, 981 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Talon, S., & Charbonnel, C. 2008, A&A, 482, 597 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Talon, S., & Zahn, J.P. 1997, A&A, 317, 749 [Google Scholar]
 Talon, S., & Zahn, J.P. 1998, A&A, 329, 315 [NASA ADS] [Google Scholar]
 Talon, S., Zahn, J.P., Maeder, A., & Meynet, G. 1997, A&A, 322, 209 [Google Scholar]
 Talon, S., Kumar, P., & Zahn, J.P. 2002, ApJ, 574, L175 [Google Scholar]
 Thompson, M. J., ChristensenDalsgaard, J., Miesch, M. S., & Toomre, J. 2003, ARA&A, 41, 599 [NASA ADS] [CrossRef] [Google Scholar]
 Triana, S. A., Moravveji, E., Pápics, P. I., et al. 2015, ApJ, 810, 16 [NASA ADS] [CrossRef] [Google Scholar]
 TurckChièze, S., Palacios, A., Marques, J. P., & Nghiem, P. A. P. 2010, ApJ, 715, 1539 [NASA ADS] [CrossRef] [Google Scholar]
 Waite, M. L., & Bartello, P. 2006, J. Fluid Mech., 568, 89 [NASA ADS] [CrossRef] [Google Scholar]
 Zahn, J.P. 1983, in SaasFee Advanced Course 13: Astrophysical Processes in Upper Main Sequence Stars, eds. A. N. Cox, S. Vauclair, & J. P. Zahn, 253 [Google Scholar]
 Zahn, J.P. 1992, A&A, 265, 115 [NASA ADS] [Google Scholar]
Appendix A: Horizontal averages
Case of τ = 1/(2Ω+S)
The use of this prescription for the turbulent characteristic time requires to compute integrals of the form
which are also defined when x > −1. The substitution gives
The partial fraction decomposition of the integrand is
When 1 − x^{2} > 0 (i.e., when −1 < x < 1), this leads to
When 1 − x^{2} < 0 (for x > 1), the last two terms of Eq. (A.3) can be further decomposed, and one finally obtains
These analytical results allow us to avoid to compute latitudinal averages in Eqs. (52) and (54) numerically and to speed up their evaluation.
Case of τ = 1/N_{Ω}
The use of this prescription requires to compute integrals of the form
which are defined when x > −1. The classical substitution t = cos θ gives
When 1 + 1/x > 0, that corresponds to x > 0, we use the partial fraction decomposition
with α^{2} = 1 + 1/x to obtain
When 1 + 1/x < 0, that corresponds to −1 < x < 0, the partial fraction decomposition
leads to
Appendix B: Validation of the shellular rotation assumption
Fig. B.1. Ratio between the horizontal variation of the angular velocity and the mean angular velocity on an isobar, expressed as a horizontal Rossby number, for the same models and colour code as in Fig. 2. 
We recall the firstorder expansion of the angular velocity introduced by Mathis & Zahn (2004):
with , where P_{2}(θ) is the secondorder Legendre polynomial, and
the socalled shellular rotation introduced by Zahn (1992), which is the average of the angular velocity on an isobar.
According to Zahn (1992) and Mathis & Zahn (2004), the asymptotic solution for Ω_{2} can be expressed as
when solving the equation for the transport of the angular momentum in the latitudinal direction assuming a short turbulent diffusion timescale r^{2}/D_{h}. Here, U_{2} and V_{2} are the l = 2 radial functions of the vertical and latitudinal components of the rotationdriven meridional circulation respectively
which is expanded on Legendre polynomials. This subsonic largescale flow verifies the anelastic approximation, i.e. ∇ ⋅ (ρ 𝒰_{M}) = 0, where acoustic waves are filtered out, that leads to
With D_{h}, U_{2} and known, V_{2} hence Ω_{2} can be computed at each time step and at each radius in the process of solving the angular momentum transport equation (Mathis & Zahn 2004).
Then, we introduce the horizontal Rossby number (see also Billant & Chomaz 2001)
To follow this quantity using an 1D stellar evolution code, we introduce its horizontal average over an hemisphere only
because of the antisymmetry of sin^{2}θ ∂_{θ}Q_{2}(θ) around the equator. The horizontal Rossby number is thus given by the ratio between the horizontal differential rotation (Ω_{2}) and the mean shellular rotation () similarly to the one introduced by Spiegel & Zahn (1992).
To verify the shellular rotation assumption, we expect in the strongly stratified regime D_{v} ≪ D_{h} and thus (e.g. Zahn 1992), i.e. . This is verified in our simulations where our new prescriptions are implemented as illustrated in Fig. B.1.
appendix C: Transport loop and strong horizontal turbulence
A larger D_{h} leads to an increase of the transport of angular momentum by the meridional circulation and of the horizontal diffusion of heat over an isobar (we refer the reader to the heat transport equation derived in Maeder & Zahn (1998; Eq. (4.36)) and in Mathis & Zahn (2004; Eqs. (101) and (102)). As can be seen in Fig. 5, which shows the rotation profile for the different D_{h} prescriptions, the amplified transport leads to a locally weaker differential rotation (its evolution is discussed in details in Sect. 3.2.3) than with smaller horizontal turbulent diffusion coefficients. To understand the details of what is going on with a stronger horizontal transport, it is interesting to consider the hydrodynamical transport loop in stellar radiation zones as introduced by Rieutord (2006) and Decressin et al. (2009). First, the meridional circulation is mechanically driven by the wind and the shearinduced vertical transport (see Eq. (16) in Decressin et al. 2009). Then, the circulation advects heat and because of the strong horizontal diffusion of entropy over the isobar, the temperature relaxes to a state where its fluctuation on an isobar is smaller than in the case with a weaker D_{h}. This can be understood by looking at Eq. (19) in Decressin et al. (2009) and considering a simplified balance between heat advection and horizontal diffusion. Next, the radial gradient of rotation is deduced from the thermal wind balance given by Eq. (17) in Decressin et al. (2009). A smaller temperature fluctuation leads to a smaller radial gradient of angular velocity. The vertical viscous shearinduced turbulent transport thus diminishes as well as the corresponding term driving the meridional circulation (Eq. (16) in Decressin et al. 2009). This term (noted U_{v} in Decressin et al. 2009) has an opposite sign compared to the term driven by the angular momentum extraction due to the stellar wind (noted U_{Γ} in Decressin et al. 2009). The amplitude of the meridional circulation thus increases in absolute value with a stronger advected flux of angular momentum towards the stellar surface and the transport loop is closed.
All Tables
All Figures
Fig. 1. Anisotropic turbulent transport in a stably stratified stellar radiation zone with a buoyancy frequency N rotating with an angular velocity Ω. We introduce the characteristic velocities (u_{∥}, u_{⊥}) and length scales (l_{∥}, l_{⊥}) of a turbulent “pancake” in the direction of the entropy stratification and in the horizontal one respectively and usual spherical coordinates (r, θ). 

In the text 
Fig. 2. Top row: total horizontal turbulent diffusivity D_{h} (full line) and its D_{h,v} component (when active if R_{e} > R_{e;c}; dotted lines) for the 1.0 M_{⊙}, Z_{⊙} models described in Table 2 at three evolutionary stages (from left to right: ZAMS, age of the Sun, and when the central H mass fraction is reduced to 10% of its original value). The hatched areas correspond to the convective regions. Bottom row: vertical turbulent diffusion coefficient in the same models. 

In the text 
Fig. 3. Top row: anisotropic ratio given by Eq. (47) for the different choices of τ. It represents the potential anisotropy of turbulence sustained by the vertical shear instability in the stably stratified differentially rotating medium. It is thus equivalent to the ratio between the horizontal and vertical turbulent diffusivities (D_{h,v}/D_{v, v}) when the latter is active. Bottom row: product of the thermal part of the Brunt–Väisälä frequency (N) by the computed turbulent characteristic time (τ). Both are shown at the same evolutionary stages and with the same colour code as in Fig 2. 

In the text 
Fig. 4. Effective diffusivity D_{eff} in the vertical direction induced by the advective transport at the same evolutionary stages and with the same colour code as in Fig 2. 

In the text 
Fig. 5. Top row: angular velocity profile represented at the same evolutionary stages and with the same colour code as in Fig 2. Bottom row: ratio between the differential rotation and the rotation rate, equivalent to the dimensionless parameter X in Eqs. (49), (51), and (54). 

In the text 
Fig. 6. Bottom panel: evolution of the angular velocity as a function of time for models computed with the prescriptions for D_{h} from Mathis et al. (2004), Zahn (1992), and that of the present paper with the three different characteristic time scales for turbulence (black, red, green (τ = 1/S), blue (τ = 1/(2Ω + S)) and magenta (τ = 1/N_{Ω}), respectively). Scaled on the left is the rotational period in days and on the right the angular velocity in solar units with Ω_{⊙} = 2.86 × 10^{−6} s^{−1}. Each cross represents a star belonging to an open cluster which age and mass (between 1.4 and 1.6 M_{⊙}) have been taken from literature. All observational data are taken from Gallet & Bouvier (2015) and Bouvier et al. (2014) and references therein. The solid lines show the evolution of the surface rotation, and dashed lines the evolution of the mean radiative core rotation rate. Top panel: evolution of the normalised core to surface differential rotation rate for the six 1.0 M_{⊙}, Z_{⊙} models as colourcoded in Fig. 2. 

In the text 
Fig. 7. Evolution of the surface and core rotation rates (solid and dashed lines respectively) as a function of the stellar radius that increases with time (MS, subgiant, and RGB correspond respectively to R_{eff}/R_{⊙} between ∼1.0 and 1.8, 1.8 and 3, and higher than 3). We show theoretical predictions for two 1.25 M_{⊙} models at solar metallicity computed with the prescription for D_{h} = D_{h,h} by Mathis et al. (2004; black) and the new prescription including both components of D_{h} = D_{h,v} + D_{h,h} for the case τ = 1/S (green). We compare the model predictions with relevant data for stars with asteroseismic masses between 1.1 and 1.4 M_{⊙}. For SGB stars we show surface and core rotation rates (black and orange dots) derived by Deheuvels et al. (2014). For RGB stars, we plot the core rotation rates derived by Mosser et al. (2012; red dots). 

In the text 
Fig. B.1. Ratio between the horizontal variation of the angular velocity and the mean angular velocity on an isobar, expressed as a horizontal Rossby number, for the same models and colour code as in Fig. 2. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.