EDP Sciences
Open Access
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/0004-6361/201629187
Published online 23 November 2018

© ESO 2018

Licence Creative Commons
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 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. (2012, 2014, 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).

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. (2012, 2014, 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).

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; 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; Turck-Chiè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; Acevedo-Arreguin et al. 2013; Barnabé et al. 2017). At the same time, the physical description of shear-induced 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 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 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 eddy-viscosities, and the corresponding eddy-diffusivities, Dv, v and Dv,h. We note that in stellar physics modelling though, the simplifying assumption νv, vDv, v and νv,hDv,h has been made until now. Similarly, four coefficients νh,h, Dh,h, νh,v and Dh,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 Dh,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 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.

Table 1.

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 Dh,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).

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, shear-induced 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

(1)

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 “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; 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).

thumbnail 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, θ).

Open with DEXTER

The control parameter of such dynamics is the ratio between the vertical and horizontal turbulent transport coefficients

(2)

and

(3)

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 non-linear regime is to cancel their source (Richard & Zahn 1999; Mathis et al. 2004). Then, their ratio can be written

(4)

where we introduce ε = l/l. From now on to Sect. 2.2.2, we lighten notations by using {Dv,Dh} that stands for {Dv, v,Dh,v} or {Dv,h,Dh,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

(5)

It compares the relative strength of the inertia associated with the horizontal advection of the horizontal turbulent motions and buoyancy1. 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

(6)

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

(7)

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

(8)

Therefore, as expected, the anisotropy becomes larger when stable stratification is stronger because of the corresponding vertical buoyancy restoring force. Thus, Dv/Dh → 0 (or Dh/Dv → ∞) when N → ∞. 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.

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 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

(9)

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 zones2, they find that

(10)

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 non-linear numerical simulations in Cartesian geometry where they explored the variations of l as a function of Roh. This proposed scaling, in conjunction with Eq. (4), leads to

(11)

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

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 pre-existing source of turbulence characterized by a velocity field u0. 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

(12)

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 u0/τ.

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:

(13)

(14)

Here, s denotes the specific entropy fluctuations and ⟨s⟩ the mean specific entropy. We also introduced the unit-vector basis {ej}j = r,θ,φ in spherical coordinates. Following Zahn (1992), here we assume a priori that Dv​ < ​​ < ​Dh 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 = Ser. 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

(15)

where ũ is the Fourier transform of u and = k/||k|| is the normalized wave vector. This can be written as a matrix equation

(16)

with

(17)

where σ = · Ω/Ω, Ω* = 2τΩ, S* = τS, N* = τN, and δij and εilj are the usual Kronecker and Levi-Civita symbols respectively. The matrix equation can be solved, using

(18)

where

(19)

and μ = r. Some terms of ℳ−1 vanish when applied to ũ0, because the latter also satisfies the continuity equation . This finally leads to

(20)

with

(21)

We introduce the spectral tensor 0, which characterizes the properties of the background turbulence. It is defined by

(22)

where ⟨⟩ denotes a statistical average. Then, it is possible to derive the spectral tensor of the generated turbulence by using the relation

(23)

In the case of purely horizontal random motions, which are an idealized view of turbulent displacements in a strongly stratified medium, we have

(24)

where E(k) is the kinetic energy spectrum function defined by

(25)

The full expression for the tensor 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

(26)

(27)

(28)

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

(29)

and Eqs. (26)–(28) yield

(30)

(31)

(32)

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

(33)

We define the two angles α and β such that

(34)

(35)

(36)

The infinitesimal wave vector element can be written

(37)

where k is the norm of k, and we now have

(38)

Equations (30)–(32) lead to

(39)

(40)

In the previously introduced notations, one has and . Using Eq. (4), we finally obtain

(41)

As expected, the strength of the horizontal turbulent transport relatively to the vertical one is thus increasing with stratification (as N4). 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 DhDh,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 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

(42)

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

(43)

The turbulent transport along the vertical direction is then of the order of , which yields

(44)

Choosing τ = S−1 finally implies

(45)

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) Re = (u l)/ν, where ν is the viscosity of the fluid, is larger than a critical value Re;c (Zahn 1992; Talon & Zahn 1997; Prat et al. 2016). In the STAREVOL code, we chose to set Re;c = 10 that corresponds to a value of Dv, 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).

Another frequency that characterizes differentially rotating fluids and their stability, the epicyclic frequency (e.g. Pringle & King 2007)

(46)

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 main-sequence since at this phase. However, for negative values of , 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.

2.2.3. Horizontal turbulent transport

Using Eq. (41), it is finally possible to derive Dh when knowing Dv (or vice-versa):

(47)

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 Dv is associated with this secular radial shear instability (i.e. DvDv, v and thus DhDh,v) and non-zero (for Re > Re;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 (Dh,h, see Mathis et al. 2004) as well as a vertical transport (and corresponding vertical turbulent diffusion coefficient Dv,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):

(48)

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 Re > Re;c, Dh,v, is derived using Eq. (47).

The first possible choice for τ, namely τ = S−1, leads to

(49)

In this case, we thus obtain a discontinuous variation for Dh,v if the Reynolds criterion is not met and the vertical shear-induced turbulence vanishes. This may be regularized once a prescription for Dv, 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

(50)

leads to

(51)

where X = rrΩ/(2Ω). For the implementation in 1D stellar evolution codes, we take a latitudinal average of this expression

(52)

where .

Finally, considering the epicyclic characteristic time

(53)

we obtain

(54)

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 Dh,v on the ratio N/(2Ω). Including rotation in the prescription for τ also allows us to regularize the behaviour of Dh,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 Dv, 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 (Talon & Zahn 1997). In the following section, we investigate the impact of the different expressions we derived for Dh,v on the transport of angular momentum along the evolution of low-mass 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 Dh is then defined as the sum of Dh,v as given by these equations (we test the different expressions for τ) and of Dh,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 Mathis & Zahn (2004); 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 mixing-length 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).

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 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 Dh,v given in Eqs. (49), (52), and (54), added to the horizontal turbulent diffusivity Dh,h derived in Mathis et al. 2004), keeping all the other ingredients unchanged. These models are also compared with a model that ignores Dh,v and includes only the prescription for Dh,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 Dh,v added to Dh,h from Mathis et al. (2004), and the last one with only Dh,h. For the case where Dh,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 pressure-driven magnetized stellar winds are taken into account using the prescription by Matt et al. (2015) with Kwind = 7 × 1030 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.

Table 2.

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 Dh and of the vertical diffusivity Dv, 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 (Dh,v) lead to values for the total Dh = Dh,v + Dh,h that are very close to the ones given by the expression of Dh,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. Re > Re;c), and this is not the case in the central part of the star for most of the MS. Therefore, the horizontal component of Dh (Dh,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 Dh,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 turn-off, Dh,v can be several orders of magnitude higher than Dh,h in these outer layers below the convective envelope.

thumbnail Fig. 2.

Top row: total horizontal turbulent diffusivity Dh (full line) and its Dh,v component (when active if Re > Re;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.

Open with DEXTER

Figure 2 also shows the profiles of the vertical turbulent diffusivity Dv = Dv, 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 Dv 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(Dv)≈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 Dh,v to the Dh,h of Mathis et al. 2004). The radial variations of Dv 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 109. This means that whenever the vertical shear is high enough to overcome the Reynolds criterion (i.e. Re > Re;c) below the convective envelope for r > 0.5 R, the total horizontal diffusion coefficient log(Dh) rises from 7.5 up to 11.5 ( ≃ 109 × Dv) that leads to a stronger meridional circulation and transport of angular momentum as explained in Appendix C.

thumbnail 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 (Dh,v/Dv, 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.

Open with DEXTER

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 U2(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 U2Dh and DeffDh. This explains both the amplitude and the radial variations of Deff observed for r > 0.5R in models A, B, and C.

thumbnail Fig. 4.

Effective diffusivity Deff in the vertical direction induced by the advective transport at the same evolutionary stages and with the same colour code as in Fig 2.

Open with DEXTER

3.2.3. Impact on the evolution of the internal and surface rotation

We now discuss the impact of introducing the new formalism for Dh 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, Dv​ < ​​ < ​Dh 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 Dh presented in Sect. 3.2.2, adding the vertical source for horizontal turbulence Dh,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 Dh,h. Therefore, the impact of the new component Dh,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 Dh,v (if active for Re > Re;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.

thumbnail 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).

Open with DEXTER

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, Dh,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.

thumbnail Fig. 6.

Bottom panel: evolution of the angular velocity as a function of time for models computed with the prescriptions for Dh 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 colour-coded in Fig. 2.

Open with DEXTER

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 Dh, 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. 25 below the convective envelope on the main sequence when Dh,v is locally several orders of magnitude higher than Dh,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 space-based 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 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. 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 Dh 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 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 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 sub-giant branch, thus to a smaller angular velocity.

thumbnail 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 Reff/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 Dh = Dh,h by Mathis et al. (2004; black) and the new prescription including both components of Dh = Dh,v + Dh,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).

Open with DEXTER

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 (, with HP 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 N4τ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 (Talon & Zahn 1997).

Then, the turbulent transport coefficient in the horizontal direction, induced by the 3D turbulent motions sustained by the vertical shear-driven instability, scales as Ric(N/2Ω)2Kf(S, Ω), where Ric 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 low-mass, solar-metallicity stars, accounting for the feedback of the vertical shear in the horizontal direction (Dh = Dh,h + Dh,v), the new prescriptions do not modify greatly the results previously obtained using the Mathis et al. (2004) prescription (Dh = Dh,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 Dh, 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 main-sequence, subgiant and red giant low-mass and intermediate-mass 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 (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.


1

The Froude number is the equivalent for stably stratified flows of the Rossby number for rotating flows.

2

See Appendix B.

3

Models with τ = 1/(2Ω + S) and τ = 1/NΩ have been computed and are indistinguishable from the model computed with τ = 1/S.

Acknowledgments

We dedicate this work to the memory of our friend and mentor Jean-Paul 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 SIMI5-6 020 01 Toupies (Towards understanding the spin evolution of stars) and the Swiss National Science Foundation (SNF).

References

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

(A.1)

which are also defined when x > −1. The substitution gives

(A.2)

The partial fraction decomposition of the integrand is

(A.3)

When 1 − x2 > 0 (i.e., when −1 < x < 1), this leads to

(A.4)

When 1 − x2 < 0 (for x > 1), the last two terms of Eq. (A.3) can be further decomposed, and one finally obtains

(A.5)

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

(A.6)

which are defined when x > −1. The classical substitution t = cos θ gives

(A.7)

When 1 + 1/x > 0, that corresponds to x > 0, we use the partial fraction decomposition

(A.8)

with α2 = 1 + 1/x to obtain

(A.9)

When 1 + 1/x < 0, that corresponds to −1 < x < 0, the partial fraction decomposition

(A.10)

leads to

(A.11)

Appendix B: Validation of the shellular rotation assumption

thumbnail 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.

Open with DEXTER

We recall the first-order expansion of the angular velocity introduced by Mathis & Zahn (2004):

(B.1)

with , where P2(θ) is the second-order Legendre polynomial, and

(B.2)

the so-called 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

(B.3)

when solving the equation for the transport of the angular momentum in the latitudinal direction assuming a short turbulent diffusion timescale r2/Dh. Here, U2 and V2 are the l = 2 radial functions of the vertical and latitudinal components of the rotation-driven meridional circulation respectively

(B.4)

which is expanded on Legendre polynomials. This sub-sonic large-scale flow verifies the anelastic approximation, i.e. ⋅ (ρ 𝒰M) = 0, where acoustic waves are filtered out, that leads to

(B.5)

With Dh, U2 and known, V2 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)

(B.6)

To follow this quantity using an 1D stellar evolution code, we introduce its horizontal average over an hemisphere only

(B.7)

because of the anti-symmetry of sin2θθQ2(θ) 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 DvDh 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 Dh 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 Dh 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 shear-induced 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 Dh. 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 shear-induced turbulent transport thus diminishes as well as the corresponding term driving the meridional circulation (Eq. (16) in Decressin et al. 2009). This term (noted Uv 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

Table 1.

The turbulent diffusivities (viscosities) and their source.

Table 2.

Specifications of the different 1 M, Z models.

All Figures

thumbnail 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, θ).

Open with DEXTER
In the text
thumbnail Fig. 2.

Top row: total horizontal turbulent diffusivity Dh (full line) and its Dh,v component (when active if Re > Re;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.

Open with DEXTER
In the text
thumbnail 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 (Dh,v/Dv, 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.

Open with DEXTER
In the text
thumbnail Fig. 4.

Effective diffusivity Deff in the vertical direction induced by the advective transport at the same evolutionary stages and with the same colour code as in Fig 2.

Open with DEXTER
In the text
thumbnail 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).

Open with DEXTER
In the text
thumbnail Fig. 6.

Bottom panel: evolution of the angular velocity as a function of time for models computed with the prescriptions for Dh 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 colour-coded in Fig. 2.

Open with DEXTER
In the text
thumbnail 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 Reff/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 Dh = Dh,h by Mathis et al. (2004; black) and the new prescription including both components of Dh = Dh,v + Dh,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).

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.