Issue 
A&A
Volume 610, February 2018



Article Number  A35  
Number of page(s)  11  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/201731381  
Published online  26 February 2018 
The 2D dynamics of radiative zones of lowmass stars
^{1}
Laboratoire AIM ParisSaclay, CEA/DRF  CNRS  Université ParisSaclay, IRFU/DAp Centre de Saclay,
91191
GifsurYvette Cedex,
France
email: hypolite.delphine@gmail.com
^{2}
Institut de Recherche en Astrophysique et Planétologie, Observatoire MidiPyrénées, Université de Toulouse,
14 avenue Edouard Belin,
31400
Toulouse,
France
Received:
15
June
2017
Accepted:
22
November
2017
Context. Helioseismology and asteroseismology allow us to probe the differential rotation deep within lowmass stars. In the solar convective envelope, the rotation varies with latitude with an equator rotating faster than the pole, which results in a shear applied on the radiative zone below. However, a polar acceleration of the convective envelope can be obtained through 3D numerical simulations in other lowmass stars and the dynamical interaction of the surface convective envelope with the radiative core needs to be investigated in the general case.
Aim. In the context of secular evolution, we aim to describe the dynamics of the radiative core of lowmass stars to get a deeper understanding of the internal transport of angular momentum in such stars, which results in a solid rotation in the Sun from 0.7R_{⊙} to 0.2R_{⊙} and a weak radial coreenvelope differential rotation in solartype stars. This study requires at least a 2D description to capture the latitudinal variations of the differential rotation.
Methods. We build 2D numerical models of a radiative core on the top of which we impose a latitudinal shear so as to reproduce a conical or cylindrical differential rotation in a convective envelope. We perform a systematic study over the Rossby number R_{o} = ΔΩ/2Ω^{0} measuring the latitudinal differential rotation at the radiative–;convective interface. We provide a 2D description of the differential rotation and the associated meridional circulation in the incompressible and stably stratified cases using the Boussinesq approximation.
Results. The imposed shear generates a geostrophic flow implying a cylindrical differential rotation in the case of an isotropic viscosity. When compared to the baroclinic flow that arises from the stable stratification, we find that the geostrophic flow is dominant when the Rossby number is high enough (R_{o} ≥ 1) with a cylindrical rotation profile. For low Rossby numbers (R_{o} < 1), the baroclinic solution dominates with a quasishellular rotation profile. Using scaling laws from 3D simulations, we show that slow rotators (Ω_{0} < 30Ω_{⊙}) are expected to have a cylindrical rotation profile. Fast rotators (Ω_{0} > 30Ω_{⊙}) may have a shellular profile at the beginning of the main sequence in stellar radiative zones.
Conclusions. This study enables us to predict different types of differential rotation and emphasizes the need for a new generation of 2D rotating stellar models developed in synergy with 3D numerical simulations. The shear induced by a surface convective zone has a strong impact on the dynamics of the underlying radiative zone in lowmass stars. However, it cannot produce a flat internal rotation profile in a solar configuration calling for additional processes for the transport of angular momentum in both radial and latitudinal directions.
Key words: hydrodynamics / stars: interiors / stars: rotation / stars: evolution
© ESO 2018
1 Introduction
Rotation is a key physical mechanism regarding the dynamical, chemical, and magnetic evolution of stars (Maeder 2009; Charbonneau 2010; Brun et al. 2015). Within the context of secular evolution, rotating radiative zones are particularly interesting since they impose the transport timescales of angular momentum and chemicals.
Indeed, differential rotation induces meridional circulation and potential shear instabilities that transport chemical elements and angular momentum (Zahn 1992; Maeder & Zahn 1998; Mathis & Zahn 2004; Rieutord 2006b) and impact the rotation profiles and chemical abundances.
In this context, helio and asteroseismology have been revolutionary, probing for the first time the internal dynamical state of stars. In the Sun, this reveals a uniform rotation profile in the radiative core (see, e.g., the review by Thompson et al. 2003) at least until 0.2R_{⊙} (Couvidat et al. 2003; García et al. 2007); the estimate of differential rotation is more difficult for deep regions such as near the center of the Sun because of the lack of identification of individual gmodes (Appourchaux et al. 2010). In the convective envelope, the differential rotation is found to be conical with a prograde equatorial acceleration and the azimuthal velocity decreasing monotonically towards higher latitude. A strong radial differential rotation is observed at the bottom of the convective zone forming a thin shear zone called the tachocline (Spiegel & Zahn 1992). This calls for strong transport processes operating at the radiative–convective interface that must be investigated in order to understand the interactions between these two regions (e.g., Garaud 2002a; Brun & Zahn 2006; Brun et al. 2011; Strugarek et al. 2011; Varela et al. 2016).
The internal differential rotation of several other types of stars has been revealed by asteroseismology. Indeed, the coretosurface rotation ratio of numerous lowmass, mainsequence, subgiant, and red giant stars (Beck et al. 2012; Deheuvels et al. 2012, 2014, 2016; Benomar et al. 2015) provides new constraints for stellar modeling. Moreover, a surface differential rotation is always observed in the convective envelope of lowmass stars (Barnes et al. 2005; Reinhold et al. 2013) and it has been confirmed by numerical simulations (Brun et al. 2015).
In this context, 3D numerical simulations have been performed to better understand the magnetic activity of solarlike stars (Brun & Toomre 2002; Brun et al. 2004, 2011; Brown et al. 2008; Matt et al. 2011; Augustson et al. 2012) and determine what influences the sign of the latitudinal differential rotation in such stars (Gastine et al. 2014; Brun et al. 2017). These studies, devoted to short timescales, show that solar parameters can produce antisolar rotation profiles (slow equator and fast pole; e.g., Käpylä et al. 2014). To be precise, depending on the convective “fluid” Rossby number , where v is the velocity, R is the stellar radius, and Ω_{0} is the stellar rotation rate, three rotational states have been identified in the convective envelope of lowmass stars: the antisolar differential rotation for , the solarlike profile when the convective Rossby number is between 0.3 and 0.9, and the “Jupiterlike” profile (cylindrical banded profile with alterning fast and slow jets) for (Brun et al. 2017). In the solar case, the isorotational contours within the convective zone have also been fitted with characteristics of the thermal wind equation showing a strong correlation between the entropy and the angular velocity (Balbus 2009; Balbus et al. 2009, 2012; Balbus & Latter 2010; Brun et al. 2011). They show a very good agreement with helioseismology data. However, these works emphasize their difficulty in reproducing the solar tachocline, where the thermal wind balance breaks, and the underlying flows in the radiative zone.
In this work, we aim to understand the impact of the shear of the convective envelope on the underlying radiative core on secular timescales along stellar evolution. Unfortunately, global 3D numerical simulations focus on the dynamical timescales and are not able to explore intensively the entire HR diagram for now. To this end, considerable efforts have been made to solve the rotational dynamics of stellar radiative zones in 1D stellar evolution models (e.g., Talon et al. 1997; Maeder & Meynet 2000; Palacios et al. 2003; Talon & Charbonnel 2005; Eggenberger et al. 2005; Decressin et al. 2009; Marques et al. 2013; Mathis et al. 2013). When ignoring internal gravity waves and magnetic fields and using the formalism by Zahn (1992), Maeder & Zahn (1998), and Mathis & Zahn (2004) that assumes a shellular differential rotation enforced by a strong horizontal turbulence, 1D stellar evolution models fail to reproduce the rotation profile of the solar radiative core and the coreenvelope rotation contrasts revealed by asteroseismology. Moreover, rotating flows are intrinsically bidimensional. In this case, the differential rotation can be radial and latitudinal and a more general 2D approach is needed. In this context, the recent improved angular momentum evolution models (Gallet & Bouvier 2013, 2015; Amard et al. 2016), which follow the rotational evolution of lowmass stars in clusters, highlight the need for fast rotating models, and therefore 2D models, for example during the early evolution phases (Hypolite & Rieutord 2014).
Indeed, in the case of fast rotation, Rieutord (2006b), Espinosa Lara & Rieutord (2013), Rieutord & Beth (2014), and Hypolite & Rieutord (2014) show that the baroclinic flow that pervades rotating radiative envelopes exhibits a meridional circulation and associated differential rotation that require a 2D description. The models they develop are intermediatemass mainsequence stars models (a radiative envelope lying upon a convective core); what is now needed are lowmass star models with the convective envelope on top of the radiative core. Regarding the solar case, Friedlander (1976) studied the spindown of a radiative zone due to applied surface stresses using the Boussinesq approximation. Garaud (2002b) describes the meridional flows in a radiative zone submitted to an imposed solarlike latitudinal shear, based on observations, to reproduce the presence of the convective zone at its top using the anelastic approximation. However, these works are limited to the solar case and a general study of the dynamics of a radiative zone lying below a convection zone with solar or antisolar stress is an interesting and necessary complement to the existing models.
For these reasons, we propose a new 2D study of the dynamics of the radiative zone of lowmass stars. As a first step, with no magnetic field, convective motions, or internal gravity waves, a 2D description is enough to capture the essentials of the physics of a rotating radiative zone submitted to a shear at its upper boundary. We construct a latitudinal shear boundary condition based on the results from an inverted rotation profile within the solar convective envelope and from numerical simulations of lowmass stars’ convective envelopes. The shear is quantified by the Rossby number , which is the normalized latitudinal differential rotation at the convective–radiative interface relative to the stellar rotation rate. We produce solar and antisolar configurations and solve for the flow of an incompressible fluid and then for a stably stratified fluid using the Boussinesq approximation. In Sect. 2, we give a complete description of the incompressible model and the hydrodynamical equations we solve when imposing a latitudinal shear. We introduce the relevant physical parameters of the problem. An analytical study unravels the dynamics of the bulk of the radiative core. We describe the properties of the flow and compare the analytical solutions to numerical simulation solutions. In Sects. 3 and 4, we study the stratified case using the Boussinesq approximation. Using 1D models of solarlike stars to compute the baroclinic torque amplitude (proportional to the Brunt–Väisälä frequency), we provide the 2D differential rotation and meridional circulation varying systematically the Rossby number. We also compute the coretosurface rotation ratio as a function of the Rossby number. In Sect. 5, we summarize our main results, namely we describe theRossby parameter regime where the geostrophic solution arising from the shear dominates the dynamics () or the baroclinic flow dominates the dynamics (). Using scaling laws derived from 3D numerical simulations, we are able to scale the Rossby number and predict the rotational state of a lowmass star’s radiative core induced by hydrodynamical processes as a function of its age and angular velocity.
2 The flow driven by an imposed differential rotation
2.1 Equations of motion
We consideran incompressible viscous fluid enclosed within a sphere. The system is rotating at a constant rate Ω_{0} aligned with the zaxis (Ω_{0} = Ω_{0}e_{z}) and the sphere is assumed, as a first step, not to suffer any deformation due to rotation. The hydrostatic deformation of the sphere is neglected. Therefore, we focus on the effects of the Coriolis acceleration. Solutions are axisymmetric and symmetric with respect to the equator as no mechanism acting in this setting (gravity and rotation) can separate the velocity field from these symmetries. The dynamics of such a fluid are governed by the momentum equation (1)
which we write in a frame rotating at angular velocity Ω_{0}. In this equation we recognize the Coriolis acceleration 2Ω_{0} ∧v and the centrifugal acceleration , where e_{s} is the radial unit vector associated with the radial cylindrical coordinate s, P is the pressure, g the gravity, and μ the dynamical viscosity. Mass conservation implies (2)
for an incompressible fluid of constant density ρ_{0}.
We then gather the centrifugal acceleration, the pressure, and the gravitational potential Φ into a single scalar so that the momentum equation becomes (3)
where is the kinematic viscosity.
2.2 Boundary conditions
Since we wish to describe a radiative core, the bounding sphere materializes the interface with an outer convective envelope. To make this interface simple, we assume that the convective envelope imposes its azimuthal velocity at the top of the radiative region. We neglect any other motion and any mass exchange. Hence at r = R (R is the radius of the radiative–convective interface) we impose (4)
which is the simplest expression we can take inspired by numerical simulations (e.g., Matt et al. 2011; Käpylä et al. 2014). Differential rotation is then called “solarlike” when Δ Ω > 0, that is, when the equatorial regions rotate faster than the pole, and “antisolar” otherwise when Δ Ω < 0.
In the frame corotating with the pole of the model, the boundary condition reads (6)
For the sake of simplicity, we also assume the meridional components of the velocity field at the upper boundary layer to be zero in the rotating frame, namely (7)
We note that these boundary conditions imply that viscous stresses are applied on top of the radiative core, which is different from the boundary conditions chosen by Garaud (2002b), who assumed that the continuity of the viscous stresses affecting the meridional circulation imposes stressfreelike boundary conditions at the interface. Besides, Brun & Zahn (2006) considered impenetrable walls at the bottom of the convective envelope (u_{r} = 0) and stressfree conditions on the latitudinal and azimuthal components of the velocity field.
Here, we choose to impose the velocity as if the convection zone behaved as a solid that can absorb any stress. This is certainly exaggerated and this excludes any mass exchange between the layers. Nevertheless, if we consider that the turbulent convection zone is endowed with a turbulent viscosity, much larger than that of the radiative zone, the shear stress of the radiative zone is likely to be unimportant to modify the flow in the convective zone. Thus imposing the velocity is likely to be more appropriate than imposing the stress. Our boundary conditions are simple and finally just assume that the flow in the radiative zone does not feed back onto the convection zone mean flows.
2.3 Scaled equations
We make the equations adimensional choosing
R as a length scale;
V = RΔΩ as a velocity scale;
and as the timescale.
This adimensionalization yields (8)
where u, τ, and p are respectively the dimensionless velocity, time, and reduced pressure. We introduced the two numbers (9)
namely the Rossby number and the Ekman number E.
We evaluate the order of magnitude of these numbers with the Sun. We take R ~ 0.7R_{⊙} and a rotation rate at the solar tachocline corresponding to a rotation period of 27 days (Brun et al. 2004). From Espinosa Lara & Rieutord (2013), we estimate the kinematic viscosity as either ν ~ 10^{2} m^{2} s^{−1} or ν ~ 10^{4} m^{2} s^{−1} if some turbulence occurs (e.g., Zahn 1992). The resulting Ekman numbers are E ~ 10^{−10} and E ~ 10^{−8} respectively.
As far as the Rossby number is concerned, a typical value for the Sun is ~0.1 (e.g., Gastine et al. 2014). As a further simplification, we shall set this number to zero so as to deal with linear equations focusing mainly on the case of rapid rotation. This assumption is valid as long as the Coriolis term is dominant compared to the nonlinear advection term, which reads as the condition . This condition is satisfied as shown in the next section as long as the Rossby number is less than one.
Taking the curl of the momentum equation, we get a simple vorticity equation for the steady state, (10)
It is completed by the mass conservation equation (11)
and the boundary condition (12)
where the parameter b = ±1 captures the sign of ΔΩ.
2.4 Analytical solution
The inviscid(E = 0) case of this setting is the geostrophic balance (13)
which has for solution an azimuthal geostrophic flow whose amplitude only depends on the radial cylindrical coordinate s, namely (14)
as a consequence of the Taylor–Proudman theorem. The sum of the inviscid solution and its boundary layer correction has to match the boundary conditions (12) and solves Eqs. (10) and (11) when E≠ 0. In this particular case, no boundary layer correction is necessary on the azimuthal flow and we readily write (15)
which satisfies the geostrophic balance (13) and boundary condition (12). The meridional circulation arises from the conservation of angular momentum expressed by the φcomponent of the momentum equation. In cylindrical coordinates it reads (16)
and leads tothe following expression of the meridional flow (17)
where mass conservation has been used.
The meridional stream function defined as u_{merid} = ∇× [ψ(r, θ)e_{φ}] and leading to (18)
We note that the foregoing meridional circulation does not meet the imposed boundary conditions at r = 1 for the radial and latitudinal velocities. Indeed, u_{r}(1)≠0. This implies the existence of a thin Ekman layer that absorbs this mass flux and generates an correction to the interior ū_{θ} of (17). The numerical solutions of the next section validate the foregoing predictions on the flow.
2.5 Numerical solutions
To prepare for the study of a more complex situation (the stably stratified case), we now solve Eqs. (10) and (11) with a spectral numerical method (Rieutord 1987). Briefly, we expand the velocity fields on the vectorial spherical harmonics, and discretize the radial functions on the Gauss–Lobatto grid associated with Chebyshev polynomials. A more detailed description of the method is given in Appendix A.
Solving numerically Eqs. (10) and (11) lead to velocity fields that can be described by a differential rotation Δ Ω = u_{φ}∕s and the stream function ψ of the meridional flow. As we can see in Fig. 1, the differential rotation is cylindrical as the Taylor–Proudman theorem predicts. The meridional circulation has a unique cell in each hemisphere whose circulation sign depends also on the sign of Δ Ω. Namely, the circulation is clockwise when ΔΩ < 0 and counter clockwise when ΔΩ > 0.
We study the dependence of the velocity field amplitude with the Ekman number and summarize the obtained results in Fig. 2. The amplitude of the azimuthal velocity does not depend on the Ekman number value as expected. The amplitude of the meridional circulation does depend on the Ekman number value and is proportional to E as also expected.
Inside the Ekman boundary layer that develops at the outer boundary near , one can notice the boundary layer corrections. The radial component is , while the latitudinal one is closer to .
In Fig. 3, we show the velocity components in the radial (orange), the latitudinal (green), and the azimuthal (purple) directions obtained numerically (solid lines). We overplot the analytical expressions that we derived in the last section with squares. The numerical and analytical solutions are in perfect agreement except near the outer boundary where boundary layer corrections apply.
Fig. 1 Differential rotation δΩ and meridional circulation stream function ψ (red: direct sense, blue: clockwise sense) shown in the meridional plane for E = 10^{−6} and b = {1, −1} (top and bottom respectively). The stellar rotation axis is vertical. 
Fig. 2 Logarithm of the absolute value of the radial (orange lines) and latitudinal (green lines) components of the velocity field as a function of E at , θ = π∕2 and for b= +1. Scaling laws are {u_{r}, u_{θ}}∝ E^{d}, where the index d is indicated for each case on the legend panel. 
Fig. 3 Comparison of the analytical expressions of the velocity components derived in the last section (squares) and the numerical results (solid lines) as a function of the radius at the colatitude θ = 0.5 for E = 10^{−6} and b = +1 in the radial(orange), latitudinal (green), and azimuthal (purple) directions. 
3 The stably stratified case
We now move to the stably stratified case. For that we insert into the fluid cold sources (following Rieutord 2006b, hereafter R06) that set a stable density stratification and trigger a baroclinic flow. We wish to determine the parameter regime where the baroclinic differential rotation and the foregoing viscositydriven differential rotation respectively dominate.
3.1 The flow equations
To keep the model simple, we use the Boussinesq approximation. Thus density fluctuations only appear in the buoyancy term here driven by the effective gravity (i.e., the combined effect of centrifugal and gravitational accelerations).
Thus, we introduce heat sinks δQ that will perturb the previous flow V_{0} = bΔΩRs^{3}e_{φ}. The heat equation for a steadystate solution thus reads (20)
where δT is the temperature perturbation generated by the heat sinks and κ is the heat diffusivity of the fluid. In this equation v = V_{0} + δv where δv is the velocity perturbation arising from the introduction of the heat sinks. We note that δv is not necessarily small compared to V_{0}. However, we shall neglect the advection heat term v ⋅∇δT altogether on the ground that we expect that δT be axisymmetric and v be dominated by its azimuthal component. We already know that the meridional circulation is smaller than the azimuthal flow in the unstratified case. The baroclinic solutions derived by Rieutord (2006b) share the same property, so we can confidently expect that this nonlinear term is small (see below for the discussion). We are thus led to a simple equation for the steady temperature field introduced by the heat sinks, namely (21)
which is solved by δT(r), where we assume that the heat sinks δQ have a spherically symmetric distribution. Let us now move to the equations for the velocity field. Mass conservation still demands (22)
because of the use of the Boussinesq approximation. The momentum equation now reads (23)
which is written in a frame rotating at angular velocity Ω_{0}, that is, the angular velocity of the pole. We have now included the associated centrifugal acceleration into the effective gravity g_{eff}.
The density is perturbed by the temperature variations so that (24)
Here ρ_{0} is constant and associated with the reference temperature T_{0}. So we also use the simple equation of state (usually associated with the Boussinesq approximation) (25)
where α is the dilation coefficient at constant pressure.
As the Boussinesq approximation commands it, we neglect δρ everywhere except in the buoyancy term. Thus, we rewrite Eq. (23) as (26)
where Π is a reduced pressure that now includes the barotropic term ρ_{0}g_{eff}. Finally, Eq. (26) is rewritten (27)
This equation can be further simplified by remarking that δT ≡ δT(r) and that . The spherically symmetric part of the buoyancy force can be incorporated into the reduced pressure, so that only the centrifugal force term needs to be kept. Finally, taking the curl of this equation we eliminate the reduced pressure gradient and obtain (28)
where the prime indicates a radial derivative. In this equation the new driving by the baroclinic torque appears explicitly.
We may now introduce the Brunt–Väisälä frequency, which quantifies the stratification of the fluid, namely (29)
where g(r) = rg_{s} for a constant density fluidand with g_{s} the gravity at the surface of the sphere. Equation (28) now reads (30)
3.2 Scaled equations
For the unstratified case, we now rescale the equations using Δ Ω R as the velocity scale and R as the length scale. We thus find (31)
where we introduced the scaled Brunt–Väisälä frequency (32)
and the relative amplitude of the centrifugal force, namely (33)
In the Sun, ε_{⊙}~ 10^{−5}. If, as in Sect. 2, we consider only the limit of small Rossby numbers, we can linearize system (31) and solve (34)
completed with boundary conditions (7) and (12). Hence, we assume , but also εn^{2} ≲ 1, which is actually possible (see below).
At this stage we may remark that unlike in the case treated in Rieutord (2006b), the heat equation has disappeared and the Prandtl number does not appear in the problem. The reason is that we neglected at the outset the heat advection by meridional currents. As shown in Rieutord (2006b), this is strictly valid in the limit of the vanishing λparameter
where is the Prandtl number defined as the ratio of the kinematic viscosity to the thermal diffusivity. In the Sun, λ_{⊙} ~ 10^{−2}. Moreover, λ ≪ 1 is obtained for stars rotating sufficiently rapidly. However, if the star is a slow rotator, steady state flows are not relevant since we know that initial conditions then also control the actual flow, because baroclinic modes are damped on the Eddington–Sweet time, which tends to infinity as rotation vanishes (e.g., Busse 1981; Rieutord 2006a). This setting generates a thermal wind solution arising from the baroclinic torque and the geostrophic solution generated by the shear imposed by the boundary conditions described in the previous section. It allows us to evaluate their competition.
3.3 Stars matching the model
The foregoing model is rather simple but uses a number of hypothesis. We now need to identify stars that match these conditions.
The righthand side amplitude of Eq. (34) can be written (35)
where M is the stellar mass and G is the gravitational constant. We compute mainsequence Modules for Experiments in Stellar Astrophysics (MESA) models (Paxton et al. 2011) with a metallicity of Z = 0.02 and the mixinglength parameter α_{MLT} = 2 (the defaultvalue used by MESA) when the central hydrogen fraction reaches 0.5 for masses in the range [0.4–1.4]M_{⊙}, constituting F, G, K, and M stars. As shown in Fig. 4, the amplitude of the term R^{3}N^{2}∕4GM on the righthand side of Eq. (34) is always less than unity. With small Rossby numbers, the term εn^{2} (r) can be of order unity, which means that the baroclinic flow has an amplitude of the same order of magnitude as the geostrophic one driven by the shear. For young stars, such as ZAMS stars, the amplitude of the term R^{3} N^{2}∕4GM is of the same order of magnitude leading to the same competition between the baroclinic and the geostrophic flows. The dynamics set by the shear may be highly modified by the baroclinicity. For these reason, we need to resort to numerical solutions to determine which flow dominates, and the corresponding parameter regime. We also determine the features of the dynamics in each identified regime.
Fig. 4 Top: Brunt–Väisälä frequency as a function of the radius, in the radiative core of lowmass stars during the main sequence (X_{c} = 0.5) for stellar masses in the range [0.4 − 1.4]M_{⊙} from MESA models. Bottom: amplitude of the righthand side of Eq. (34) (to be multiplied by ) as a function of the radius. 
4 Numerical solutions
4.1 Differential rotation and meridional circulation
Using the spectral numerical method described in Appendix A, we numerically solve Eq. (34). We use the 1M_{⊙} MESA model as an input for the righthand side of the vorticity equation and we vary the Rossby number systematically.
We show in Fig. 5, the differential rotation and the associated meridional circulation for Rossby numbers between 10^{−3} and 10. When (i.e., for weak imposed differential rotation), the dynamics is typical of the flow that arises when a baroclinic torque is applied as describedby R06. The thermal wind balance breaks the Taylor–Proudman one. The differential rotation is roughly shellular and the number of cells of the meridional circulation is equal to the number of inflection points of the Brunt–Väisälä frequency profile plus one, here two. These patterns are aligned with the cylindrical zdirection because of rotation.
Differences with Fig. 4 of R06 come from the upper boundary conditions we set. As shown by Eq. (7), we impose a surface shear on the azimuthal velocity coupled with no penetrative boundary conditions on the meridional components of the velocity field, while R06 imposes regular stressfree conditions.
For Rossby numbers higher than the threshold (i.e., for strong imposed differential rotation), the differential rotation profile tends to be cylindrical. The thermal wind from the baroclinic torque is weaker in comparison with the geostrophic flow and the Taylor–Proudman balance is restored. The meridional circulation is dominated by a single, global circulation pattern in each hemisphere. For positive b (solarlike differential rotation case), the meridional circulation is counterclockwise and the differential rotation shows an equatorial acceleration as the imposed differential rotation at the boundary. For negative b, that is, when the equator is slower than the pole, we observe the same behavior but the direction of the meridional circulation is clockwise and the field of differential rotation is reversed with a polar acceleration. At the surface, the fluid moves toward the equator, which rotates slower than the pole.
With MESA models we explored the influence of the stellar mass on the foregoing results, but remaining with the case of solartype stars. In all cases, the differential rotation and associated meridional circulation turn out to have the same properties as previously described since the shape of the Brunt–Väisälä frequency profile is the same for masses in [0.4–1.4]M_{⊙}. The amplitudes of the flows are slightly different but remain of the same order of magnitude as in the solar case.
Fig. 5 Differential rotation δΩ and meridional circulation stream function ψ (red: direct sense, blue: clockwise sense) shown in the meridional plane for E = 2 × 10^{−6} and (top to bottom). The two left columns are for b = +1 and the two on the right for b = −1. The stellar rotation axis is vertical. 
4.2 Coretosurface rotation ratio
Asteroseismic analysis provides the coretosurface rotation ratio of numerous lowmass stars (Benomar et al. 2015). Being the only internal rotation diagnosis as of today because of the low angular resolution provided by asteroseismology, we therefore compute this ratio as the latitudinal average of the angular velocity over a radius close to the center, namely r_{core} ~ 0.15, and near the surface, r_{surf} ~ 1 (36)
We do not compute at the center r = 0 since individual gmodes are hardly identifiable in this region (García et al. 2007; Appourchaux et al. 2010; García 2010). When plotting the differential rotation in Fig. 5, we subtracted the rotation rate of the pole at the surface to the rotation rate so that negative value of the rotation rate means that the examined zone rotates retrogradely in the reference frame rotating with the pole. The coretosurface rotation rate, computed in this framework, is shown in Fig. 6 as a function of the Rossby number (Δ Ω∕2Ω_{0}).
The surface shear drives completely the averaged surface rotation rate, which has the sign of b and whose amplitude is independent of the Rossby number. The averaged core rotation rate is strongly impacted by the baroclinic torque amplitude, which induces a slow shellularlike rotation when compared to the outer latitudinally averaged rotation. Thus, the averaged core rotation rate is negative in most cases and tends towards zero when increasing (i.e., decreasing the baroclinic torque amplitude). In conclusion, the coretosurface rotation ratio (in absolute value) decreases as increases. The only departure from this behavior (the small bump for positive b and high ) is due to a change of sign of the averaged core rotation rate, which becomes positive for Rossby numbers higher than 5, that is, when the geostrophic solution completely dominates. In the retrograde case, the baroclinic flow has to have a higher amplitude to overcome the (retrograde) shearinduced flow than in the prograde case. The scaling law index is m = −1 when writing . At low Rossby numbers, when the baroclinic solution is dominant, the coretosurface rotation ratio is high, which means that the differential rotation between the core and the surface is large in absolute value but, when considering the pole rotation, the core is slower than the surface according to our simulation. Conversely, for large Rossby numbers, the geostrophic solution is dominant and the coretosurface rotation ratio is small, with a rapidly rotating core (respectively slowly) when b = −1 (respectively b = +1), as illustrated in Fig. 5.
Fig. 6 Logarithm of the coretosurface rotation ratio vs. the logarithm of the Rossby number for E = 2 × 10^{−6}. 
5 Discussion and conclusion
5.1 Rossby number and timescales
Regarding the 3D numerical simulation results for the differentially rotating convective envelopes of solarlike stars, the relative latitudinal shear has been scaled as (37)
with m = −0.11 in the hydrodynamical case or m = −0.56 according toMagneto Hydro Dynamical numerical simulations (Varela et al. 2016). This scaling law is closer to the scaling law derived from the observations ΔΩ ∝Ω^{0.15} (Barnes et al. 2005; Reinhold et al. 2013) than the pure hydrodynamical one. The relative shear, hence the Rossby number in our study, decreases as the stellar rotation rate increases. For angular velocities between 1 and 100 Ω_{⊙}, which is the range of observed angular velocities in solarlike stars reported in Gallet & Bouvier (2013, 2015), the corresponding Rossby number is between 7 × 10^{−2} and 5 × 10^{−3} using the MHD scaling law and n = 0.44 from Varela et al. (2016). For such values of the Rossby number, our numerical results predict that the dynamics are driven by the baroclinicity with a quasishellular differential rotation as illustrated in Fig. 5. Indeed, we recall that with our setting, the baroclinic torque amplitude increases when the Rossby number decreases.
However, our study compares the steady solution of the baroclinic and the geostrophic flows. Meanwhile it is not certain whether these solutions reach their steady states during the mainsequence lifetime (10^{10} yr for a solarlike star).
The characteristic timescale for the settling of the geostrophic solution is the one of a spinup (38)
where P is the rotation period and E the Ekman number. In the stellar case, the Ekman number is around 10^{−10} leading to τ_{SU} = 10^{5}P. If P is of the order of tens of days, as is the Sun, the transient solution does not last more than ~10^{4} yr, which is short in comparison to the time a star spends on the main sequence.
The baroclinic modes dampen on the Eddington–Sweet timescale, namely on (39)
where τ_{KH} = R^{2}∕κ ≈ 10^{8} yr in the solar case, is the thermal Kelvin–Helmholtz timescale. When the stellar rotation rate increases, the Eddington–Sweet timescale gets shorter and tends to the Kelvin–Helmholtz timescale. This is often the case on the premain sequence.
Hence, in rapid or young rotators, the baroclinic steady state is likely to be reached. For slower (older) rotators like the present Sun, the ratio is around 10^{4} leading to τ_{ES} ~ 10^{12} yr. Since it is very large, it means that there are likely residuals of the baroclinic modes (initial conditions). Therefore, for solarlike stars, the geostrophic solution is most certainly steady, while it is not clear for the baroclinic one depending on the stellar rotation rate and thus the age.
We compute the Eddington–Sweet timescale as a function of the angular velocity and find that the limit for the steady state to be reached on a solar mainsequence lifetime is Ω_{0} ~ 30Ω_{⊙}. Stars rotating faster than Ω_{0} ~ 30Ω_{⊙} can thus reach a steady baroclinic state and have baroclinic dynamics within their radiative zone. For stars rotating slower than Ω_{0} ~ 30Ω_{⊙}, the dynamics are not determined directly because the baroclinic flow is unsteady and depends on initial conditions.
If we consider that the geostrophic flow is steady and that the baroclinic flow amplitude superposes on it according to the following transient time evolution (40)
its amplitude becomes comparable to that of the geostrophic flow at (41)
which is within the range of Ω_{0} ∈ [Ω_{⊙}, 10^{2}Ω_{⊙}], always longer than or comparable to the mainsequence lifetime because of the amplitude of the Rossby number. Therefore, the dynamics of the radiative core of solartype stars rotating slower than Ω_{0} ~ 30Ω_{⊙} may be dominated by the shear with a cylindrical differential rotation imposed by the geostrophic balance according to the initial conditions. These would be the differential rotation profiles at the end of the premain sequence when the gravitational contraction ends. Hypolite & Rieutord (2014) have shown that it would be also cylindrical.
We summarize these results in Fig. 7, where we use the results from Gallet & Bouvier (2013). We plot the schematical rotational evolution of the envelope of solarlike stars Gallet et Bouvier obtained in a rapid, median, and slow rotating case. We delimit the region where we can expect the steady baroclinic solution to drive a shellular differential rotation in red, and the region where the differential rotation may be cylindrical due to the imposed shear (yellow region). The steady state analysis is clearly limited for this study even though we make such an approximation in order to keep the equations linear since the nonlinear term is also proportional to . A time evolution analysis of the baroclinic torque over initial conditions of the geostrophic flow would insure this prediction.
If we use the scaling law derived from the observations by Barnes et al. (2005) and Reinhold et al. (2013), we obtain high Rossby numbers in such a way that, comparing only the steady state amplitudes, the baroclinic torque would be very small and the cylindrical rotation profile induced by the shear would be dominant for all rotation rates. Also, on the mainsequence, solartype stars, age is correlated to the rotation rate through a gyrochronological law (Skumanich 1972; Kawaler 1988; Réville et al. 2015; van Saders et al. 2016) since stars undergo a spindown through stellar wind. Indeed, the wind is responsible for a mass and angular momentum loss expected to generate a spindown geostrophic solution in the bulk of the star (Rieutord & Beth 2014) on a timescale similar to the spinup case given by Eq. (38), that is, a short timescale regarding the duration of the main sequence. Therefore, a complete study would require us to also consider the spindown flow from Rieutord & Beth (2014) in competition with the baroclinic flow and the flow driven by the imposed convective differential rotation.
Our simple model predicts that the Sun’s radiative core should have a cylindrical differential rotation with an angular velocity gradient that is rather mild and similar to the one deduced from asteroseismic observations of solarlike (M ∈ [1, 1.6]M_{⊙}) mainsequence stars (Benomar et al. 2015). However, the cylindrical differential rotation is in contradiction with the observed solidbody rotation of the present solar radiative zone (at least for r ≥ 0.2R_{⊙}). Other processes possibly therefore come into play, such as internal gravity waves (Zahn et al. 1997), hydrodynamic shearinduced anisotropic turbulence (Zahn 1992), or magnetic fields (Mestel & Weiss 1987; Gough & McIntyre 1998).
Fig. 7 Schematical trends of angular velocity Ω_{⋆}∕Ω_{⊙} of the envelope of fast (blue line), median (green line), and slow (red line) solartype rotators as a function of age from Gallet & Bouvier (2013). The open circle is the angular velocity of the present Sun. Using the scaling law from Varela et al. (2016), the red area delimits where the dynamics should be driven by the baroclinicity leading to a shellular rotation. The yellow area shows where our modeling allows us to expect a cylindrical differential rotation in the radiative core of solarlike stars if taking only into account rotation and induced largescale flows. 
5.2 Conclusion
In order to study the dynamics induced by the shear that a surface convective envelope imposes on an internal radiative core in mainsequence lowmass stars, we built a simplified model of a central radiative zone on the top of which we impose a latitudinal shear boundary condition. We considered a simple latitudinal shear at the surface of the model inspired by heliosismology inversion profiles and global numerical simulations of the convective envelope of differentially rotating lowmass stars. This allowed us to quantify the impact of the shear between the pole and the equator of the convective–radiative interface on the steady flow within the radiative core.
Analytically, we find that the imposed shear drives a geostrophic solution in addition to the thermal wind rising from the stable stratification. The baroclinic flow is characterized by a quasishellular differential rotation with multiple cells of meridional circulation, while the geostrophic flow tends to sustain a cylindrical differential rotation. When the geostrophic solution induced by the shear is dominant (high Rossby number regime ), there is a unique cell in each hemisphere, which is quite similar to the previous results obtained by Friedlander (1976) and Garaud (2002b) in the solar case. The high Rossby number case must be interpreted very carefully since it is formally out of the linear regime. Indeed, for , the nonlinear advection term may be important and we did not include it in this work.
Since the baroclinic torque amplitude is proportional to the inverse of the Rossby number, we evaluate its value according to the scaling laws found in 3D numerical simulations by Varela et al. (2016) to determine if the spinup flow from the shear overtakes the baroclinic flow. The Rossby number is found to decrease as the global stellar rotation rate increases. For such values of Rossby numbers (), the baroclinic flow is expected to be dominant for stars rotating faster than 30Ω_{⊙}. For stars rotating more slowly, the baroclinic flow is probably still transient and less important in amplitude than the geostrophic flow. Therefore, we expect these stars to have a cylindrical differential rotation. Scaling laws for Rossby numbers directly derived from observations (Reinhold et al. 2013) suggest that the cylindrical differential rotation profile is dominant for all rotators. Regarding our 2D models with parameters closest to the Sun, we do not reproduce the flat rotation profile observed both in radius and latitude in the solar radiative zone. As in previous results obtained in 1D, this strengthens the need to take into account an efficient process responsible for extra transport of angular momentum. Such a process could be internal gravity waves (Talon & Charbonnel 2005), magnetic fields (Mathis & Zahn 2005; Strugarek et al. 2011; AcevedoArreguin et al. 2013), and the anisotropic turbulent transport (Mathis et al. 2017). The effects of these processes will be studied in forthcoming works.
In addition, the difference in the sign of the differential rotation between Rieutord (2006b) models and the compressible Evolution STellaire en Rotation (ESTER) mainsequence models (Espinosa Lara & Rieutord 2013) also suggests that the compressibility may play an important role and that the Boussinesq approximation is just a first step before using a more detailed modeling such as one using the anelastic approximation. The general method presented here will be also applied to other types of stars such as intermediatemass and massive stars with a differentially rotating convective core and Ftype stars with both a differentially rotating convective core and envelope.
Acknowledgements
The authors thank the referee for their comments that allowed us to improve the article. D. Hypolite and S. Mathis acknowledge funding by the European Research Council through ERC SPIRE grant 647383 and the CNES PLATO grant at CEA Saclay. The authors also wish to thank the MESA website. M. Rieutord acknowledges funding by the French Agence Nationale de la Recherche (ANR), under grant ESRR (ANR16CE31000701). This work was also supported by the Centre National de la Recherche Scientifique (CNRS, UMR 5277), through the Programme National de Physique Stellaire (PNPS), and by the CNES PLATO grant at IRAP.
Appendix A Spectral expansion of hydrodynamical equations
The radial grid points, corresponding to the expansion onto a Gauss–Lobatto grid derived with Chebyshev polynomials, are defined as (A.1)
The colatitude dependence of the solution is described with the vectorial spherical harmonics (e.g. Rieutord 1987). We represent the velocity field as follows (A.2)
In this formula, are the normalized spherical harmonics (e.g. CohenTannoudji et al. 1997), e_{r} is the unit radial vector, and the horizontal gradient operator is defined on the unit sphere. The axisymmetry of the solutions imposes m = 0 (∂_{φ} = 0). For this reason we will not write the index m in the following.
A.1 Hydrodynamical equations
The continuity Eq. (22) reads on this expansion (A.4)
The vorticity Eq. (34) projected onto the R_{l} function is written (A.5)
and onto the T_{l} direction (A.6)
where δ_{ij} is the Kronecker symbol and . The projection onto the S_{l} function is redundant with the first one and the component of the velocity field v^{l} is computed with the continuity equation projection, Eq. (A.4). The coupling coefficients read (A.7)
A.2 Boundary conditions at the upper boundary
The azimuthal velocity written on the (R_{l}, S_{l}, T_{l}) basis reads (A.8)
Because the model is symmetric with respect to the equator, the azimuthal velocity has the property to be fully described by odd l. Using only the two first harmonics l = 1 and l = 3, the azimuthal velocity reads (A.9)
In orderto set the shear described by expression (6), we then set at the surface (A.10)
All other lcomponents of the azimuthal velocity are zero at the upper boundary. Components with even l numbers (u^{l}, v^{l}) are also zero because we set to zero the meridional velocity components at r = 1.
References
 AcevedoArreguin, L. A., Garaud, P., & Wood, T. S. 2013, MNRAS, 434, 720 [NASA ADS] [CrossRef] [Google Scholar]
 Amard, L., Palacios, A., Charbonnel, C., Gallet, F., & Bouvier, J. 2016, A&A, 587, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Appourchaux, T., Belkacem, K., Broomhall, A.M., et al. 2010, A&ARv, 18, 197 [NASA ADS] [CrossRef] [Google Scholar]
 Augustson, K. C., Brown, B. P., Brun, A. S., Miesch, M. S., & Toomre, J. 2012, ApJ, 756, 169 [NASA ADS] [CrossRef] [Google Scholar]
 Balbus, S. A. 2009, MNRAS, 395, 2056 [NASA ADS] [CrossRef] [Google Scholar]
 Balbus, S. A.,& Latter, H. N. 2010, MNRAS, 407, 2565 [NASA ADS] [CrossRef] [Google Scholar]
 Balbus, S. A., Bonart, J., Latter, H. N., & Weiss, N. O. 2009, MNRAS, 400, 176 [NASA ADS] [CrossRef] [Google Scholar]
 Balbus, S. A., Latter, H., & Weiss, N. 2012, MNRAS, 420, 2457 [NASA ADS] [CrossRef] [Google Scholar]
 Barnes, J. R., Collier Cameron, A., Donati, J.F., et al. 2005, MNRAS, 357, L1 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Beck, P. G., Montalban, J., Kallinger, T., et al. 2012, Nature, 481, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Benomar, O., Takata, M., Shibahashi, H., Ceillier, T., & García, R. A. 2015, MNRAS, 452, 2654 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, B. P., Browning, M. K., Brun, A. S., Miesch, M. S., & Toomre, J. 2008, ApJ, 689, 1354 [CrossRef] [Google Scholar]
 Brun, A. S., & Toomre, J. 2002, ApJ, 570, 865 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brun, A. S., & Zahn, J.P. 2006, A&A, 457, 665 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brun, A. S., Miesch, M. S., & Toomre, J. 2004, ApJ, 614, 1073 [NASA ADS] [CrossRef] [Google Scholar]
 Brun, A. S., Miesch, M. S., & Toomre, J. 2011, ApJ, 742, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Brun, A. S., García, R. A., Houdek, G., Nandy, D., & Pinsonneault, M. 2015, Space Sci. Rev., 196, 303 [NASA ADS] [CrossRef] [Google Scholar]
 Brun, A. S., Strugarek, A., Varela, J., et al. 2017, ApJ, 836, 192 [NASA ADS] [CrossRef] [Google Scholar]
 Busse, F. H. 1981, Geophysical and Astrophysical Fluid Dynamics, 17, 215 [Google Scholar]
 Charbonneau, P. 2010, Living Reviews in Solar Physics, 7, 3 [Google Scholar]
 CohenTannoudji, C., & Diu, B., & Laloë, F. 1997, Mécanique quantique (Ed. Hermann) [Google Scholar]
 Couvidat, S., García, R. A., TurckChièze, S., et al. 2003, ApJ, 597, L77 [NASA ADS] [CrossRef] [Google Scholar]
 Decressin, T., Mathis, S., Palacios, A., et al. 2009, A&A, 495, 271 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Deheuvels, S., García, R. A., Chaplin, W. J., et al. 2012, ApJ, 756, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Deheuvels, S., Doğan, G., Goupil, M. J., et al. 2014, A&A, 564, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Deheuvels, S., Brandão, I., Silva Aguirre, V., et al. 2016, A&A, 589, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Eggenberger, P., Maeder, A., & Meynet, G. 2005, A&A, 440, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Espinosa Lara, F., & Rieutord, M. 2013, A&A, 552, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Friedlander, S. 1976, J. Fluid Mech., 76, 209 [NASA ADS] [CrossRef] [Google Scholar]
 Gallet, F., & Bouvier, J. 2013, A&A, 556, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gallet, F., & Bouvier, J. 2015, A&A, 577, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Garaud, P. 2002a, MNRAS, 329, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Garaud, P. 2002b, MNRAS, 335, 707 [NASA ADS] [CrossRef] [Google Scholar]
 García, R. A. 2010, Highlights of Astronomy, 15, 345 [NASA ADS] [Google Scholar]
 García, R. A., TurckChièze, S., JiménezReyes, S. J., et al. 2007, Science, 316, 1591 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Gastine, T., Yadav, R. K., Morin, J., Reiners, A., & Wicht, J. 2014, MNRAS, 438, L76 [NASA ADS] [CrossRef] [Google Scholar]
 Gough, D. O., & McIntyre, M. E. 1998, Nature, 394, 755 [NASA ADS] [CrossRef] [Google Scholar]
 Hypolite, D., & Rieutord, M. 2014, A&A, 572, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Käpylä, P. J., Käpylä, M. J., & Brandenburg, A. 2014, A&A, 570, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kawaler, S. D. 1988, ApJ, 333, 236 [NASA ADS] [CrossRef] [Google Scholar]
 Maeder, A. 2009, Physics, Formation and Evolution of Rotating Stars (Berlin, Heidelberg: SpringerVerlag) [Google Scholar]
 Maeder, A., & Meynet, G. 2000, ARA&A, 38, 143 [NASA ADS] [CrossRef] [Google Scholar]
 Maeder, A., & Zahn, J.P. 1998, A&A, 334, 1000 [NASA ADS] [Google Scholar]
 Marques, J. P., Goupil, M. J., Lebreton, Y., et al. 2013, A&A, 549, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., & Zahn, J.P. 2004, A&A, 425, 229 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., & Zahn, J.P. 2005, A&A, 440, 653 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., Decressin, T., Eggenberger, P., & Charbonnel, C. 2013, A&A, 558, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., Prat, V., Amard, L., Charbonnel, C., & Palacios, A. 2017, A&A, submitted [Google Scholar]
 Matt, S. P., Do Cao, O., Brown, B. P., & Brun, A. S. 2011, Astron. Nachr., 332, 897 [NASA ADS] [CrossRef] [Google Scholar]
 Mestel, L., & Weiss, N. O. 1987, MNRAS, 226, 123 [NASA ADS] [CrossRef] [Google Scholar]
 Palacios, A., Talon, S., Charbonnel, C., & Forestini, M. 2003, A&A, 399, 603 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
 Reinhold, T., Reiners, A., & Basri, G. 2013, A&A, 560, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Réville, V., Brun, A. S., Matt, S. P., Strugarek, A., & Pinto, R. F. 2015, ApJ, 798, 116 [NASA ADS] [CrossRef] [Google Scholar]
 Rieutord, M. 1987, Geophysical and Astrophysical Fluid Dynamics, 39, 163 [NASA ADS] [CrossRef] [Google Scholar]
 Rieutord, M. 2006a, in EAS Pub. Ser., eds. M. Rieutord, & B. Dubrulle, 21, 275 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rieutord, M. 2006b, A&A, 451, 1025 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rieutord, M., & Beth, A. 2014, A&A, 570, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Skumanich, A. 1972, ApJ, 171, 565 [NASA ADS] [CrossRef] [Google Scholar]
 Spiegel, E. A., & Zahn, J.P. 1992, A&A, 265, 106 [NASA ADS] [Google Scholar]
 Strugarek, A., Brun, A. S., & Zahn, J.P. 2011, A&A, 532, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Talon, S., & Charbonnel, C. 2005, A&A, 440, 981 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Talon, S., Zahn, J.P., Maeder, A., & Meynet, G. 1997, A&A, 322, 209 [NASA ADS] [Google Scholar]
 Thompson, M. J., ChristensenDalsgaard, J., Miesch, M. S., & Toomre, J. 2003, ARA&A, 41, 599 [NASA ADS] [CrossRef] [Google Scholar]
 van Saders, J. L., Ceillier, T., Metcalfe, T. S., et al. 2016, Nature, 529, 181 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Varela, J., Strugarek, A., & Brun, A. S. 2016, Advances in Space Research, 58, 1507 [NASA ADS] [CrossRef] [Google Scholar]
 Zahn, J.P. 1992, A&A, 265, 115 [NASA ADS] [Google Scholar]
 Zahn, J.P., Talon, S., & Matias, J. 1997, A&A, 322, 320 [Google Scholar]
All Figures
Fig. 1 Differential rotation δΩ and meridional circulation stream function ψ (red: direct sense, blue: clockwise sense) shown in the meridional plane for E = 10^{−6} and b = {1, −1} (top and bottom respectively). The stellar rotation axis is vertical. 

In the text 
Fig. 2 Logarithm of the absolute value of the radial (orange lines) and latitudinal (green lines) components of the velocity field as a function of E at , θ = π∕2 and for b= +1. Scaling laws are {u_{r}, u_{θ}}∝ E^{d}, where the index d is indicated for each case on the legend panel. 

In the text 
Fig. 3 Comparison of the analytical expressions of the velocity components derived in the last section (squares) and the numerical results (solid lines) as a function of the radius at the colatitude θ = 0.5 for E = 10^{−6} and b = +1 in the radial(orange), latitudinal (green), and azimuthal (purple) directions. 

In the text 
Fig. 4 Top: Brunt–Väisälä frequency as a function of the radius, in the radiative core of lowmass stars during the main sequence (X_{c} = 0.5) for stellar masses in the range [0.4 − 1.4]M_{⊙} from MESA models. Bottom: amplitude of the righthand side of Eq. (34) (to be multiplied by ) as a function of the radius. 

In the text 
Fig. 5 Differential rotation δΩ and meridional circulation stream function ψ (red: direct sense, blue: clockwise sense) shown in the meridional plane for E = 2 × 10^{−6} and (top to bottom). The two left columns are for b = +1 and the two on the right for b = −1. The stellar rotation axis is vertical. 

In the text 
Fig. 6 Logarithm of the coretosurface rotation ratio vs. the logarithm of the Rossby number for E = 2 × 10^{−6}. 

In the text 
Fig. 7 Schematical trends of angular velocity Ω_{⋆}∕Ω_{⊙} of the envelope of fast (blue line), median (green line), and slow (red line) solartype rotators as a function of age from Gallet & Bouvier (2013). The open circle is the angular velocity of the present Sun. Using the scaling law from Varela et al. (2016), the red area delimits where the dynamics should be driven by the baroclinicity leading to a shellular rotation. The yellow area shows where our modeling allows us to expect a cylindrical differential rotation in the radiative core of solarlike stars if taking only into account rotation and induced largescale flows. 

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