2D dynamics of the radiative core of low mass stars

Understanding the internal rotation of low mass stars all along their evolution is of primary interest when studying their rotational dynamics, internal mixing and magnetic field generation. In this context, helio- and asteroseismology probe angular velocity gradients deep within solar type stars at different evolutionary stages. Still the rotation close to the center of such stars on the main sequence is hardly detectable and the dynamical interaction of the radiative core with the surface convective envelope is not well understood. For instance, the influence of the differential rotation profile sustained by convection and applied as a boundary condition to the radiation zone is very important in the formation of tachoclines. In this work, we study a 2D hydrodynamical model of a radiative core when an imposed, solar or anti-solar, differential rotation is applied at the upper boundary. This model uses the Boussinesq approximation and we find that the shear induces a cylindrical differential rotation associated with a unique cell of meridional circulation in each hemisphere (counterclockwise when the shear is solar-like and clockwise when it is anti-solar). The results are discussed in the framework of seismic observables (internal rotation rate, core-to-surface rotation ratio) while perspectives to improve our modeling by including magnetic field or transport by internal gravity waves will be discussed.


Introduction
Asteroseismology has revealed the internal rotation profile of numerous stars over the past few years.On one hand, we now have partial access to the dynamical state of low mass stars on the main sequence and during sub giant and red giant evolution phases ( [3], [12], [13], [14], [4]).On the other hand, the internal rotation profile of the Sun deep within the star, until 0.2R , has been inverted ( [11], [20]).Both helioseismology and asteroseismology show that a strong transport of angular momentum occurs in radiative zones.These constrains call for a new generation of stellar models since current 1D stellar models do not recover inverted quantities such as for instance the flat rotation profile of the radiative zone of the Sun (e.g.[46], [48]) and of more evolved stars ( [15], [9], [32]).
Indeed, one who wants to model the evolution of stars gets immediately confronted with the complexity of stellar rotation ( [29], [31]).It is known to have multiple strong consequences on stellar evolution since it drives rotational mixing and transport of angular momentum through the stable radiative zones ( [50], [28], [34]) on secular timescale.In addition, it is connected with the generation of stellar magnetic fields and the presence of tachoclines ( [43], [23], [6], [45]) requesting a deep understanding of the stellar rotation.One of the most challenging aspects is the multidimensional nature of rotating stars ( [38], [16]).e-mail: delphine.hypolite@cea.frIndeed, the centrifugal acceleration leads to the rise of baroclinic flows which can not be fully described by low spherical harmonics expansion and/or latitudinal averages.For these reasons, we propose a 2D approach to describe the secular evolution of the internal dynamics of rotating radiation zones.In this framework, the convective envelope applies a latitudinal shear through differential rotation on the underlying radiative core in low mass stars.For example, in the Sun, the rotation profile of the envelope is conical with fast equatorial regions and slower poles ( [42], [5]).Past studies have looked for the consequences of such a shear in the solar case ( [17], [19]).However, it calls for a generalisation of the dynamical boundary conditions since 3D numerical simulations reveal that solar-like differential rotation is a function of stellar fundamental parameters and mean rotation rate while anti-solar differential rotation can also be expected in cool stars ( [37], [26], [22], [49]).
In this framework, the question we address is: How does the differential rotation of a convective envelope affect the dynamics (differential rotation and meridional circulation) of the radiative core in low mass stars?
In this work, we use 2D Boussinesq modeling to perform a parametric study over the shear applied as a boundary condition on the radiative core and compute the latitudinally averaged core to the surface rotation ratio of the resulting models which can be compared with seismic observables.
EPJ Web of Conferences 160, 02006 (2017) We describe the 2D modeling of the radiative core of low mass stars and its main assumptions in Sect. 2. In Sect.3, we compute the core-to-surface rotation ratio numerically as a function of the shear applied at the top of the domain and discuss the fast rotating solar case.We give our conclusions and perspectives in Sect. 4.

Modeling low mass stars radiative core in 2D
We consider a viscous flow enclosed within a rotating spherical shell.We suppose the centrifugal acceleration does not deform it.We study the interaction of the combined effect of the inertia with rotation in a stably stratified environment, leaving aside magnetic fields ( [18], [35], [45]), internal gravity waves ( [10], [47]) and anisotropic turbulence ( [50], [33], [27]) as a first step.Because we are interested in the evolution of the fluid on secular timescale, we solve as a first step the steady equation of the vorticity combined with the equations of energy and continuity using the Boussinesq approximation ( [38], [40], [24]).This approximation retains density variations only in the buoyancy term taking into account both the gravity and centrifugal accelerations and implies that the stratification is only taken into account through the Brunt-Väisälä frequency profile (positive in a radiative zone).It is given as an input of the simulations using the 1D MESA stellar evolution code ( [41]) to generate ZAMS low mass stars radiative core models (we use a metallicity of Z = 0.02, and a mixing length theory parameter of α MLT = 2) and we seek the resulting velocity field, i.e. the differential rotation and the meridional circulation for a 1M model.As developed by [38], hereafter R06, and [25], the scaled equations we solve are where u is the scaled velocity field in the corotating frame, θ T is the scaled temperature profile and n is the scaled Brunt-Väisälä frequency profile.We remove dimensions using R c , the radius of the radiative core, as a length scale, V, the baroclinic velocity as a velocity scale, T as a temperature scale and N the scale of the Brunt-Väisälä frequency profile both from R06 (We refer the reader to R06 for the detailed description of the adimensionalization).
The physical parameters that characterize the system are the Ekman number E which quantifies the importance of the viscosity (we note the kinematic viscosity ν) over the Coriolis effects where Ω 0 is the rotation rate of the corotating frame, and the Prandtl number Pr which is the ratio of the kinematic viscosity over the thermal diffusivity κ The Ekman number is evaluated and found to be very small in radiative core of low mass stars (around 10 −11 ).
We use here the smallest value that numerical resolution allows, i.e.E = 10 −6 .Finally, we use a solar value for the Prandtl number, namely Pr = 2.10 −6 ([6]).In this framework, 3D simulations ( [37], [26], [22], [49]) found that the differential rotation is conical in convective envelopes of low mass stars.We therefore impose at the top of the radiative core In the corotating frame, the azimuthal dimensionless velocity reads where the dimensionless number b = R c ∆Ω V quantifies the shear at the top of the radiative core.The differential rotation is solar-like when the equatorial regions rotate faster than the pole (b > 0) and anti-solar when they rotate slower (b < 0).
At the radiative-convective envelope, the meridional components of the velocity are set to zero illustrating the non penetrative case u r = 0 and assuming that the presence of a convective envelope induces u θ = 0 in the corotating frame.
The shear we impose by this dynamical boundary condition generates a geostrophic flow, i.e. a flow respecting the geostrophic balance where the pressure gradient compensates the Coriolis acceleration, characterized by a columnar differential rotation which depends only on the distance to the rotation axis as predicted by the Taylor-Proudman theorem.Without any shear or for a very weak shear |b| < 10 −2 , the flow is driven by the baroclinic torque −n 2 (r) sin θ cos θ e ϕ in the vorticity equation.The solution generated and described by R06 in Eq. ( 7) is the thermal wind characterized here by a quasi shellular differential rotation as shown in Fig. 1 (left panel).The difference between the flow obtained by R06 and the one we show in Fig. 1 comes from the fact that the Brunt-Väisälä frequency profile used is different in R06 (polynomial interpolation) than the one we use (input of a 1M ZAMS MESA model) and the boundary conditions are stress-free in R06.
When |b| > 10 −2 , the amplitude of the geostrophic flow arising from the shear applied at the upper boundary is larger than the baroclinic solution, which would dominate otherwise, and the Taylor-Proudman balance tends to be restored leading to a quasi-columnar structure.In Fig. 2, we show the resulting differential rotation δΩ for b = 10 and b = −10.The field is normalized by its extremum value and we subtract the rotation of the pole so it is zero there and the differential rotation in the bulk is relative to it.Therefore, the differential rotation tends towards a cylindrical profile which is very different from the shellular profile assumed by 1D stellar models.The associated meridional circulation is dominated by a single, global circulation cell in each hemisphere.For b positive (solar-like differential rotation), the meridional circulation is counter clockwise.At the surface of the model of the radiative core, the fluid moves toward the pole which rotates slower than the equator.Conversely, for negative b, the meridional circulation is clockwise and towards the equatorial regions at the radiative-convective interface.

Seismic observables
The core-to-surface rotation ratio is the latitudinally averaged rotation rate at the top of the radiative core divided by the latitudinally averaged central rotation rate δΩ C θ / δΩ S θ where . . .θ = π/2 0 . . .sin θdθ.This ratio has been inverted in asteroseismology on several low mass stars.This only seismic diagnosis as of today is computed to characterize its behavior according to the shear strength we impose.
In our set-up, at the radiative-convective interface, an Ekman layer arises due to rotation and viscosity.This layer permits the velocity field (the sum of both the baroclinic flow and the geostrophic one) to match the boundary conditions by adding corrections whose amplitudes are important within the layer and weak outside of it.The Ekman layer has a characteristic thickness of √ E, where we recall that E is the dimensionless Ekman number.We therefore compute the averaged surface rotation rate δΩ S θ just outside of it to avoid the boundary layer corrections.Moreover, we do not compute it at r = 1 since it would tend toward zero for small b and consequently the ratio δΩ C θ / δΩ S θ would diverge.The averaged core rotation rate δΩ C θ is computed at the radius r = 0.15R c because g-modes, that would allow to probe central regions, are hardly identifiable closer to the center ( [2], [21]).
We display the ratio δΩ C θ / δΩ S θ as a function of b in Fig. 3 in logarithmic values.When |b| tends towards zero, the ratio tends to a finite value.This is because, when the shear is very small, the dynamics is dominated by the baroclinic flow described by R06 visible in Fig. 1.With the Brunt-Väisälä frequency profile from MESA, the differential rotation is shellular and since we subtract the rotation of the pole, the differential rotation is negative in the entire zone.Therefore the ratio is positive that corre-sponds to negative averaged core rotation rate divided by a negative averaged surface rotation rate.
When b is negative, the ratio stays positive even when the shear becomes important and the geostrophic flow dominates.Moreover, we observe that when the absolute value of the shear increases, the ratio decreases.Indeed, the shear implies a transport of angular momentum, proportional to b both in the radial and latitudinal directions, which tends to accelerate the superficial regions more than the central regions.
If b is positive, the ratio displays a singularity around b ∼ 0.053.This singularity is due to the change of sign of δΩ S θ .Indeed for b ≥ 0.053, the acceleration of the equator we impose through the shear boundary condition induces a positive δΩ S θ .The ratio is then negative and has the same behavior than the branch of negative b. Namely, when b increases, the absolute value of the ratio decreases due to the enhanced transport of angular momentum.Roughly on the right side of the singularity, the solution is dominated by the geostrophic flow.
The position of the singularity depends on the zero we choose as a reference for the differential rotation; here, we recall that we choose the rotation of the pole to be zero.It also depends on the radius at which we compute δΩ S θ .Indeed, if we compute it at r = 1, the ratio δΩ C θ / δΩ S θ would diverge for small |b| since the boundary conditions would reduce to no-slip ones.The singularity would disappear as it would coincide with the divergence at small b and the sign of the ratio δΩ C θ / δΩ S θ would be the opposite of the sign of b for all b.
Since we work in the corotating frame at an unknown rotation rate and we set the rotation of the pole to be zero, this study provides only a qualitative behavior of the coreto-surface rotation ratio and the sign has a relative physical meaning.The important behavior to remember is the decrease of the ratio as the shear is strong.
In the derivation of the Eqs.(1), we assumed that the Brunt-Väisälä frequency N 2 is close to 4Ω 2 which is the fast rotating case where the transient baroclinic timescale is short enough for the steady state study to be relevant (R06, [39]).This means that a fast rotating solar case (a young Sun) is reachable for b = 10.For this value, the averaged core-to-surface rotation rate ratio is expected to be rather small, less than one in agreement with observations ( [4]).But the internal rotation profile is far from flat.As we can see in Fig. 4, the rotation profiles for different colatitudes within the radiative core do not coincide as observed by heliosismology.Moreover, we do not produce a tachocline, pointing out that our current modeling is lacking a physical process responsible for the current rotation field of the solar radiative core (for reviews, see [36], [7]).

Conclusions
We build a 2D model of fast rotating radiative cores in low mass stars undergoing the shear of a differentially rotating convective envelope.
This shear is responsible for a redistribution of angular momentum deep within the star and gives rise to  a geostrophic solution with a cylindrical differential rotation.Such a profile is very different from the shellular rotation profile assumed in 1D models.These models should therefore be used carefully in the radiative core of low mass stars when the shear applied by a convective envelope is important, i.e. b > 10 −2 .Regarding the seismic observables such as the averaged core-to-surface rotation rate or the internal rotation profile, we conclude that we have to take into account other processes responsible for efficient transport of angular momentum both in the radial and in the latitudinal directions.Such processes may be internal gravity waves ( [51]), magnetic fields ( [23], [44], [45], [1]) or an anisotropic turbulence ( [50], [30], [33]).These candidates are considered for implementation in our current 2D modeling of radiative zones.

Figure 1 :
Figure 1: Differential rotation δΩ and meridional circulation stream function ψ (red : direct sens, blue: clockwise sens) shown in the meridional plane for an Ekman number E = 10 −6 , a Prandtl number Pr = 2.10 −6 and b = 10 −2 .This is the pure baroclinic dynamics.The stellar rotation axis is vertical.

Figure 2 :
Figure 2: Differential rotation δΩ and meridional circulation stream function ψ (red : direct sens, blue: clockwise sens) shown in the meridional plane for an Ekman number E = 10 −6 , a Prandtl number Pr = 2.10 −6 with b = 10 (solar like rotation) on the left and b = −10 (anti-solar rotation) on the right.The stellar rotation axis is vertical.

Figure 3 :
Figure 3: Logarithm of the averaged core-to-surface rotation ratio as a function of the logarithm of the shear parameter b for an Ekman number E = 10 −6 and a Prandtl number Pr = 2.10 −6 .Orange dots are for positive b, green dots for negative b.The circled dots are negative core-to-surface rotation ratio (non circled are positive).When b 1 the dynamics is dominated by the baroclinic flow and the core-to-surface rotation ratio is constant (horizontal dashed line).For b positive, δΩ S θ changes sign inducing a singularity located at the vertical dashed line.

Figure 4 :
Figure 4: Relative radial differential rotation profiles within the radiative zone for b = 10 and E = 10 −6 at different colatitudes.