Anisotropic turbulent transport in stably stratified rotating stellar radiation zones

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. Turbulent transport in differentially rotating stably stratified stellar radiation zones should be carefully modeled and its strength evaluated. Stratification and rotation imply that this turbulent transport is anisotropic. Only phenomenological prescriptions have been proposed for the transport in the horizontal direction, which however constitutes a cornerstone in current theoretical formalisms for stellar hydrodynamics in evolution codes. We derive 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 implement this framework in the stellar evolution code STAREVOL and quantify its impact on the rotational and structural evolution of low-mass stars from the pre-main-sequence to the red giant branch. The anisotropy of the turbulent transport scales as $N^4\tau^2/\left(2\Omega^2\right)$, $N$ and $\Omega$ being the buoyancy and rotation frequencies respectively and $\tau$ 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. As a consequence, a complementary transport mechanism like internal gravity waves or magnetic fields is still needed to explain the observed strong transport of angular momentum along stellar evolution.


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 astero-seismology 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. (2012Deheuvels et al. ( , 2014Deheuvels et al. ( , 2015, Spada et al. (2016), Gehan et al. (2018) found a weak core-to-surface rotation contrast in low-mass subgiant and red giant stars. Next, Benomar et al. (2015) observed 26 solar-type 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 intermediate-mass 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). A&A 620, A22 (2018) (ν v,v ) Vertical turbulent transport induced by the vertical shear instability D h,v (ν h,v ) Horizontal turbulent transport induced by the 3D motions of the vertical shear instability D h,h (ν h,h ) Horizontal turbulent transport induced by the horizontal shear instability D v,h (ν v,h ) Vertical turbulent transport induced by the 3D motions of the horizontal shear instability 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 rotation-driven vertical and horizontal turbulence and meridional flows (see also Maeder & Zahn 1998;. 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 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.
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 high-resolution 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 2013Prat 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; 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 eddy-viscosities, and the corresponding eddy-diffusivities, D v,v and D v,h . We note that in stellar physics modelling though, the sim- 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, , derive a prescription based on results observed for turbulent transport in a non-stratified 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.
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 pre-main sequence and main-sequence 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).  Zahn 1992;Maeder 2009). Among the physical processes at the origin of such secular transport of chemicals and angular momentum, shear-induced turbulence plays a key role (Knobloch & Spruit 1982;Zahn et al. 1983;Zahn 1992;. 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

Turbulent transport in stably
where we consider only the thermal part with g the gravity, H P = |dr/d ln P| the pressure scale height, and P the pressure. We introduce the logarithmic gradients ∇ = d ln T/d ln P with T the temperature, and ∇ ad = (∂ ln T/∂ ln P) ad , and the thermodynamic coefficient δ = − (∂ ln ρ/∂ ln T ) 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 "pancake-like" 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;) that have been successfully implemented in several stellar evolution codes: the Geneva code (e.g. 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). 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 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, θ).
them (Chaboyer & Zahn 1992). Such an assumption is justified when the effect of the instabilities in their non-linear regime is to cancel their source (Richard & Zahn 1999;). Then, their ratio can be written where we introduce ε = l /l ⊥ . From now on to Sect.
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 non-rotating 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 self-similarity that leads to ε = l l ⊥ ≈ u u ⊥ ≈ Fr ⊥ and thus Fr ≈ 1.
These scaling laws have been verified in high-resolution non-linear 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, However, stars are potentially rapidly rotating, for example, young solar-type stars and intermediate-mass 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.

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 eddy-viscosity ν 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 high-resolution 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 large-scale 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 pre-existing 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 high-resolution non-linear numerical simulations and showed that the results agree very well. This strongly motivates a generalization of their work by taking into account simultaneously a large-scale 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

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 pre-existing source of turbulence characterized by a velocity field u 0 . We write fluid equations using the Boussinesq approximation where we assume that space-scales 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 non-linear terms. It consists of assuming that the fluid tends to reach a steady state with a typical timescale τ. Mathematically, it reads dX dt 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 = r sin θ ∇Ω, 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 unit-vector 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 spacẽ whereũ is the Fourier transform of u andk = k/ k is the normalized wave vector. This can be written as a matrix equation with where σ =k · Ω/Ω, Ω * = 2τΩ, S * = τS , N * = τN, and δ i j and ε il j are the usual Kronecker and Levi-Civita symbols respectively. The matrix equation can be solved, using where and µ =k r . Some terms of M −1 vanish when applied toũ 0 , because the latter also satisfies the continuity equationk jũ 0 j = 0. This finally leads tõ with We introduce the spectral tensorQ 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 turbulenceQ by using the relatioñ In the case of purely horizontal random motions, which are an idealized view of turbulent displacements in a strongly stratified medium, we havẽ where E(k) is the kinetic energy spectrum function defined by The full expression for the tensorQ 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 havẽ 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 and Eqs. (26)-(28) yield 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 solid-body 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 A22, page 5 of 16 A&A 620, A22 (2018) The infinitesimal wave vector element can be written where k is the norm of k, and we now have Equations (30) In the previously introduced notations, one has u 2 ⊥ = u 2 θ + u 2 ϕ and u 2 = u 2 r . 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; 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.

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 = r sin θ∂ 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 N Ω = √ 2Ω(2Ω + S sin θ), which is one of the frequencies characterizing differentially rotating flows and their stability, τ = 1/N Ω . First, we consider the case of unstable non-rotating 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 non-rotating 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 S l 2 , 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;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 non-linear interactions (e.g. Sen et al. 2012).
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 N 2 Ω < 0. In practice, this is generally not a problem during the A22, page 6 of 16 S. Mathis et al.: Anisotropic turbulent transport in stably stratified rotating stellar radiation zones main-sequence since N 2 Ω > 0 at this phase. However, for negative values of N 2 Ω , one would need to model the effect of the Rayleigh-Taylor 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 non-linear 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.

Horizontal turbulent transport
Using Eq. (41), it is finally possible to derive D h when knowing D v (or vice-versa): 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;. 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 non-zero (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 ) 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 high-resolution Cartesian non-linear 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 ... θ = π 0 ...sin θ dθ/ π 0 sin θ dθ. 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 shear-induced 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 feed-back of the horizontal turbulence on the vertical turbulent transport ). 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 low-mass main sequence (MS), subgiant and red giant stars.

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 rotation-induced 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 shear-induced turbulence and meridional flows, which are treated using the formalism derived by Zahn (1992), Maeder & Zahn (1998), and ; uniform rotation is assumed in convective regions. Angular momentum losses at the stellar surface due to pressure-driven 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 shear-induced horizontal turbulence on the transport of angular momentum in low-mass 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).

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 low-mass star. For this, we compared several models computed from the pre-main 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 , 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 .

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  (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, 2015Amard et al. 2016). We summarize in Table 2 the specifications of the different 1 M , Z models that we compare below.

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  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  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. 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 A22, page 8 of 16  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. 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 shear-induced 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 . 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 shear-stable 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.
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 . By definition, D eff = 1 30 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.

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 A22, page 9 of 16 A&A 620, A22 (2018)  respected, that is, D v < < D h and Ω (r, θ) ≈ Ω (r). 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 A22, page 10 of 16 S. Mathis et al.: Anisotropic turbulent transport in stably stratified rotating stellar radiation zones 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. The rotation profile flatness obtained with 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.
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 solar-mass star, this is the case for approximately 10% of the lifetime, and this fraction decreases with the strength of the stellar-wind 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 self-regulated 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 self-regulation 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 non-local 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 low-mass 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.  , 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 colour-coded in Fig. 2. 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 sub-giants 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. 2012Deheuvels et al. , 2014Deheuvels et al. , 2015Spada 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  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 sub-giant 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 sur- 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 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). face 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 sub-giant branch, thus to a smaller angular velocity.
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 gravity-dominated modes, that is in the region below the peak of the chemical Brunt-Väisälä frequency (N 2 µ = gϕ H P ∇ µ , 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.

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 physically-motivated 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 shear-induced 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 .
Then, the turbulent transport coefficient in the horizontal direction, induced by the 3D turbulent motions sustained by the vertical shear-driven instability, scales as Ri c (N/2Ω) 2 K f (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; 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 low-mass, solar-metallicity 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  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 low-mass main-sequence stars), it is a self-regulated 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 low-mass 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 rotation-driven 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 ) 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.