© 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 g-modes (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 core-to-surface rotation ratio of numerous low-mass, main-sequence, 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 low-mass 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 solar-like 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 anti-solar rotation profiles (slow equator and fast pole; e.g., Käpylä et al. 2014). To be precise, depending on the convective “fluid” Rossby number $\mathcal{R}{o}_f=v/2{\Omega }_{0}R$, where v is the velocity, R is the stellar radius, and Ω₀ is the stellar rotation rate, three rotational states have been identified in the convective envelope of low-mass stars: the anti-solar differential rotation for $\mathcal{R}{o}_f>1$, the solar-like profile when the convective Rossby number is between 0.3 and 0.9, and the “Jupiter-like” profile (cylindrical banded profile with alterning fast and slow jets) for $\mathcal{R}{o}_f<0.3$ (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 H-R 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 core-envelope 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 low-mass 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 intermediate-mass main-sequence stars models (a radiative envelope lying upon a convective core); what is now needed are low-mass star models with the convective envelope on top of the radiative core. Regarding the solar case, Friedlander (1976) studied the spin-down 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 solar-like 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 anti-solar 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 low-mass 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 low-mass stars’ convective envelopes. The shear is quantified by the Rossby number $\mathcal{R}o=\Delta \Omega /2{\Omega }_0$, which is the normalized latitudinal differential rotation at the convective–radiative interface relative to the stellar rotation rate. We produce solar and anti-solar 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 solar-like 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 core-to-surface 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 ($\mathcal{R}o\ge 1$) or the baroclinic flow dominates the dynamics ($\mathcal{R}o<1$). Using scaling laws derived from 3D numerical simulations, we are able to scale the Rossby number and predict the rotational state of a low-mass 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 Ω₀ aligned with the z-axis (Ω₀ = Ω₀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 $$\begin{array}{lll}\hfill & \hfill & {\rho }_0[{\partial }_{t}v+(v\cdot \nabla )v+2{\Omega }_0\wedge v-{\Omega }_0^{2}s{e}_s]=-\nabla P+{\rho }_{0}g+\mu \Delta v,\hfill \end{array}$$(1)

which we write in a frame rotating at angular velocity Ω₀. In this equation we recognize the Coriolis acceleration 2Ω₀ ∧v and the centrifugal acceleration $-{\Omega }_0^{2}s{e}_s$, 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 $$\nabla \cdot v=0$$(2)

for an incompressible fluid of constant density ρ₀.

We then gather the centrifugal acceleration, the pressure, and the gravitational potential Φ into a single scalar $\Pi =\frac{P}{{\rho }_0}+\Phi -\frac{1}{2}{s}^2{\Omega }_0^2$ so that the momentum equation becomes $${\partial }_{t}v+(v\cdot \nabla )v+2{\Omega }_0\wedge v=-\nabla \Pi +\nu \Delta v,$$(3)

where $\nu =\frac{\mu }{{\rho }_0}$ 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 $$v=R\sin \theta {\Omega }_{cz}(r=R,\theta )e_{\phi },$$(4)

with $${\Omega }_{cz}(r=R,\theta )={\Omega }_0+\Delta \Omega {\sin }^2\theta ,$$(5)

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 “solar-like” when Δ Ω > 0, that is, when the equatorial regions rotate faster than the pole, and “anti-solar” otherwise when Δ Ω < 0.

In the frame corotating with the pole of the model, the boundary condition reads $$v_{\phi }(r=R,\theta )=R{\sin }^3\theta \Delta \Omega .$$(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 $$v_r(r=R,\theta )=v_{\theta }(r=R,\theta )=0.$$(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 stress-free-like boundary conditions at the interface. Besides, Brun & Zahn (2006) considered impenetrable walls at the bottom of the convective envelope (u_r = 0) and stress-free 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 ${(2{\Omega }_0)}^{-1}$ as the timescale.

This adimensionalization yields $$\{\begin{array}{c}{\partial }_{\tau }u+\mathcal{R}o (u\cdot \nabla )u+e_z\wedge u=-\nabla p+E\Delta u,\\ \\ \nabla \cdot u=0,\end{array}$$(8)

where u, τ, and p are respectively the dimensionless velocity, time, and reduced pressure. We introduced the two numbers $$\mathcal{R}o=\frac{\Delta \Omega }{2{\Omega }_0},\text{\hspace{1em}and\hspace{1em}\hspace{1em}}E=\frac{\nu }{2{\Omega }_0{R}^2},$$(9)

namely the Rossby number $\mathcal{R}o$ 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² m² s⁻¹ or ν ~ 10⁴ m² s⁻¹ if some turbulence occurs (e.g., Zahn 1992). The resulting Ekman numbers are E ~ 10⁻¹⁰ and E ~ 10⁻⁸ 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 $u_{\phi }\ll \mathcal{R}{o}^{-1}$. 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, $$\nabla \wedge (e_z\wedge u-E\Delta u)=0.$$(10)

It is completed by the mass conservation equation $$\nabla \cdot u=0,$$(11)

and the boundary condition $$u=b{\sin }^3\theta e_{\phi }\text{\hspace{1em}at\hspace{1em}}r=1,$$(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 $$e_z\wedge u=-\nabla p,$$(13)

which has for solution an azimuthal geostrophic flow whose amplitude only depends on the radial cylindrical coordinate s, namely $$\overline{u}=F(s)e_{\phi },$$(14)

as a consequence of the Taylor–Proudman theorem. The sum of the inviscid solution and its boundary layer correction $\overline{u}+\tilde{u}$ 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 $$F(s)=b{s}^3,$$(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 $$u_s=E\left(\Delta -\frac{1}{s^2}\right)u_{\phi },$$(16)

and leads tothe following expression of the meridional flow $$u_{\text{merid}}=8bE(s{e}_s-2z{e}_z),$$(17)

where mass conservation has been used.

The meridional stream function defined as u_merid = ∇× [ψ(r, θ)e_φ] and leading to $$u_r=\frac{1}{r\sin \theta }{\partial }_{\theta }(\sin \theta \psi )\text{\hspace{1em}and\hspace{1em}}{u}_{\theta }=-\frac{1}{r}{\partial }_r(r\psi ),$$(18)

reads $$\psi =-8bEsz.$$(19)

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 $O(\sqrt{E})-{\tilde{u}}_{\theta }$ 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 $r=1-\frac{\sqrt{E}}{2}$, one can notice the boundary layer corrections. The radial component is $\mathcal{O}(E)$, while the latitudinal one is closer to $\mathcal{O}(\sqrt{E})$.

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 $r=\{0.5,0.75,1-\frac{\sqrt{E}}{2}\}$, θ = π∕2 and for b= +1. Scaling laws are {ur,  uθ}∝ Ed, 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 viscosity-driven 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₀ = b|ΔΩ|Rs³e_φ. The heat equation for a steady-state solution thus reads $$v\cdot \nabla \delta T=\kappa \Delta \delta T+\delta Q,$$(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₀ + δ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₀. 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 $O(E)$ 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 $$\kappa \Delta \delta T+\delta Q=0,$$(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 $$\nabla \cdot v=0,$$(22)

because of the use of the Boussinesq approximation. The momentum equation now reads $$\rho (v\cdot \nabla v+2{\Omega }_0\wedge v)=-\nabla P+\rho g_{\text{eff}}+\mu \Delta v,$$(23)

which is written in a frame rotating at angular velocity Ω₀, 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 $$\rho ={\rho }_0+\delta \rho .$$(24)

Here ρ₀ is constant and associated with the reference temperature T₀. So we also use the simple equation of state (usually associated with the Boussinesq approximation) $$\frac{\delta \rho }{{\rho }_0}=-\alpha (T-T_0)=-\alpha \delta T,$$(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 $${\rho }_0[(v\cdot \nabla )v+2{\Omega }_0\wedge v]=-\nabla \Pi +\delta \rho g_{\text{eff}}+\mu \Delta g,$$(26)

where Π is a reduced pressure that now includes the barotropic term ρ₀g_eff. Finally, Eq. (26) is rewritten $$(v\cdot \nabla )v+2{\Omega }_0\wedge v=-\nabla \Pi /{\rho }_0-\alpha \delta T{g}_{\text{eff}}+\nu \Delta v.$$(27)

This equation can be further simplified by remarking that δT ≡ δT(r) and that $g_{\text{eff}}=g_0(r)+{\Omega }_0^{2}s{e}_s$. 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 $$\nabla \wedge [(v\cdot \nabla )v+2{\Omega }_0\wedge v-\nu \Delta v]=-\alpha {\Omega }_0^2\delta T^{\prime }(r)r\sin \theta \cos \theta e_{\phi },$$(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 $$N^2(r)=\alpha \delta T^{\prime }g(r),$$(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 $$\nabla \wedge [(v\cdot \nabla )v+2{\Omega }_0\wedge v-\nu \Delta v]=-N^2(r)\frac{{\Omega }_0^{2}r}{g(r)}\sin \theta \cos \theta e_{\phi }.$$(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 $$\{\begin{array}{l}\nabla \wedge [\mathcal{R}o(u\cdot \nabla )u+e_z\wedge u-E\Delta u]=-\epsilon n^2(r)\sin \theta \cos \theta e_{\phi },\hfill \\ \hfill \\ \nabla \cdot u=0,\hfill \end{array}$$(31)

where we introduced the scaled Brunt–Väisälä frequency $$n^2(r)=\frac{N^2(r)}{2{\Omega }_0|\Delta \Omega |}=\frac{N^2(r)}{4{\Omega }_0^2}\mathcal{R}{o}^{-1},$$(32)

and the relative amplitude of the centrifugal force, namely $$\epsilon =\frac{{\Omega }_0^{2}R}{g_s}.$$(33)

In the Sun, ε_⊙~ 10⁻⁵. If, as in Sect. 2, we consider only the limit of small Rossby numbers, we can linearize system (31) and solve $$\{\begin{array}{l}\nabla \wedge (e_z\wedge u-E\Delta u)=-\epsilon n^2(r)\sin \theta \cos \theta e_{\phi },\hfill \\ \hfill \\ \nabla \cdot u=0,\hfill \end{array}$$(34)

completed with boundary conditions (7) and (12). Hence, we assume $\mathcal{R}o\ll 1$, but also εn² ≲ 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 $$\lambda =P\frac{N^2}{4{\Omega }_0^2},$$

where $P=\nu /\kappa$ is the Prandtl number defined as the ratio of the kinematic viscosity to the thermal diffusivity. In the Sun, λ_⊙ ~ 10⁻². 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 right-hand side amplitude of Eq. (34) can be written $$\epsilon n^2(r)=\frac{1}{\mathcal{R}o}\frac{R^3{N}^2(r)}{4GM},$$(35)

where M is the stellar mass and G is the gravitational constant. We compute main-sequence Modules for Experiments in Stellar Astrophysics (MESA) models (Paxton et al. 2011) with a metallicity of Z = 0.02 and the mixing-length 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³N²∕4GM on the right-hand side of Eq. (34) is always less than unity. With small Rossby numbers, the term εn² (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³ N²∕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 low-mass stars during the main sequence (Xc = 0.5) for stellar masses in the range [0.4 − 1.4]M⊙ from MESA models. Bottom: amplitude of the right-hand side of Eq. (34) (to be multiplied by $\mathcal{R}{o}^{-1}$) 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 right-hand 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⁻³ and 10. When $\mathcal{R}o\le {10}^{-2}$ (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 z-direction 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 stress-free conditions.

For Rossby numbers higher than the threshold $\mathcal{R}o\approx 1$ (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 (solar-like differential rotation case), the meridional circulation is counter-clockwise 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 solar-type 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 $\mathcal{R}o=\{10,1,{10}^{-1},{10}^{-2},{10}^{-3}\}$ (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 Core-to-surface rotation ratio

Asteroseismic analysis provides the core-to-surface rotation ratio of numerous low-mass 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 $$\frac{{\langle \delta {\Omega }_C\rangle }_{\theta }}{{\langle \delta {\Omega }_S\rangle }_{\theta }}=\frac{{\displaystyle {\int }_0^{\pi /2}\delta }\Omega (r_{\text{core}},\theta )\sin \theta \text{d}\theta }{{\displaystyle {\int }_0^{\pi /2}\delta }\Omega (r_{\text{surf}},\theta )\sin \theta \text{d}\theta }.$$(36)

We do not compute ${\langle \delta {\Omega }_C\rangle }_{\theta }$ at the center r = 0 since individual g-modes 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 core-to-surface rotation rate, computed in this framework, is shown in Fig. 6 as a function of the Rossby number (Δ Ω∕2Ω₀).

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 shellular-like 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 $\mathcal{R}o$ (i.e., decreasing the baroclinic torque amplitude). In conclusion, the core-to-surface rotation ratio (in absolute value) decreases as $\mathcal{R}o$ increases. The only departure from this behavior (the small bump for positive b and high $\mathcal{R}o$) 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) shear-induced flow than in the prograde case. The scaling law index is m = −1 when writing $|{\langle \delta {\Omega }_C\rangle }_{\theta }/{\langle \delta {\Omega }_S\rangle }_{\theta }|\propto \mathcal{R}{o}^m$. At low Rossby numbers, when the baroclinic solution is dominant, the core-to-surface 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 core-to-surface 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 core-to-surface 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 solar-like stars, the relative latitudinal shear has been scaled as $$\frac{\Delta \Omega }{{\Omega }_0} \propto {\Omega }_0^m,$$(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 solar-like stars reported in Gallet & Bouvier (2013, 2015), the corresponding Rossby number is between 7 × 10⁻² and 5 × 10⁻³ using the MHD scaling law $\Delta \Omega =6.5\times {10}^{-5}{\Omega }_0^n$ 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 quasi-shellular 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 main-sequence lifetime (10¹⁰ yr for a solar-like star).

The characteristic timescale for the settling of the geostrophic solution is the one of a spin-up $${\tau }_{\text{SU}}=\frac{P}{\sqrt{E}},$$(38)

where P is the rotation period and E the Ekman number. In the stellar case, the Ekman number is around 10⁻¹⁰ leading to τ_SU = 10⁵P. If P is of the order of tens of days, as is the Sun, the transient solution does not last more than ~10⁴ 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 $${\tau }_{\text{ES}}={\tau }_{\text{KH}}\frac{N^2}{4{\Omega }_0^2},$$(39)

where τ_KH = R²∕κ ≈ 10⁸ 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 pre-main 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 $\frac{N^2}{4{\Omega }_0^2}$ is around 10⁴ leading to τ_ES ~ 10¹² yr. Since it is very large, it means that there are likely residuals of the baroclinic modes (initial conditions). Therefore, for solar-like 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 main-sequence lifetime is Ω₀ ~ 30Ω_⊙. Stars rotating faster than Ω₀ ~ 30Ω_⊙ can thus reach a steady baroclinic state and have baroclinic dynamics within their radiative zone. For stars rotating slower than Ω₀ ~ 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 $$\mathcal{R}{o}^{-1}\left(\frac{1-e^{-\tau /{\tau }_{\text{ES}}}}{1-e^{-1}}\right),$$(40)

its amplitude becomes comparable to that of the geostrophic flow at $$\tau =-{\tau }_{\text{ES}}\ln [1-\mathcal{R}o(1-e^{-1})],$$(41)

which is within the range of Ω₀ ∈ [Ω_⊙, 10²Ω_⊙], always longer than or comparable to the main-sequence lifetime because of the amplitude of the Rossby number. Therefore, the dynamics of the radiative core of solar-type stars rotating slower than Ω₀ ~ 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 pre-main 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 solar-like 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 $\mathcal{R}o$. 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 main-sequence, solar-type 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 spin-down through stellar wind. Indeed, the wind is responsible for a mass and angular momentum loss expected to generate a spin-down geostrophic solution in the bulk of the star (Rieutord & Beth 2014) on a timescale similar to the spin-up 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 spin-down 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 solar-like (M ∈ [1, 1.6]M_⊙) main-sequence stars (Benomar et al. 2015). However, the cylindrical differential rotation is in contradiction with the observed solid-body 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 shear-induced 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) solar-type 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 solar-like stars if taking only into account rotation and induced large-scale 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 main-sequence low-mass 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 low-mass 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 $\mathcal{O}(1)$ geostrophic solution in addition to the thermal wind rising from the stable stratification. The baroclinic flow is characterized by a quasi-shellular 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 $\mathcal{R}o=\Delta \Omega /2{\Omega }_0>1$), 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 $\mathcal{R}o>1$, 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 spin-up 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 ($\mathcal{R}o~{10}^{-2}$), 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; Acevedo-Arreguin 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) main-sequence 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 intermediate-mass and massive stars with a differentially rotating convective core and F-type 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 (ANR-16-CE31-0007-01). 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 $$\{\begin{array}{l}{r}_j=\frac{1}{2}\left(1-\cos \left(\frac{j\pi }{N_r-1}\right)\right),\hfill \\ 0\le j\le N_r-1,\hfill \\ r_j\in [0;1].\hfill \end{array}$$(A.1)

The colatitude dependence of the solution is described with the vectorial spherical harmonics $(R_l^m,S_l^m,T_l^m)$ (e.g. Rieutord 1987). We represent the velocity field as follows $$\begin{array}{lll}u\hfill & =\hfill & {\displaystyle \sum _{l=0}^{+\infty }{\displaystyle \sum _{m=-l}^{+l}\{u_m^l(}}r)R_l^m(\theta ,\phi )\hfill \\ \hfill & \hfill & +v_m^l(r)S_l^m(\theta ,\phi )+w_m^l(r)T_l^m(\theta ,\phi )\},\hfill \end{array}$$(A.2)

where $$R_l^m=Y_l^m(\theta ,\phi )e_r,S_l^m={\nabla }_H{Y}_l^m,T_l^m={\nabla }_H\wedge R_l^m.$$(A.3)

In this formula, $Y_l^m$ are the normalized spherical harmonics (e.g. Cohen-Tannoudji et al. 1997), e_r is the unit radial vector, and the horizontal gradient operator ${\nabla }_H={\partial }_{\theta }{e}_{\theta }+\frac{1}{\sin \theta }{\partial }_{\phi }{e}_{\phi }$ 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 $$v^l=\frac{1}{l(l+1)}\frac{1}{r}\frac{\partial }{\partial r}(r^2{u}^l).$$(A.4)

The vorticity Eq. (34) projected onto the R_l function is written $$A_{l-1}^l{r}^{l-1}\frac{\partial }{\partial r}\left(\frac{u^{l-1}}{r^{l-2}}\right)+A_{l+1}^l{r}^{-l-2}\frac{\partial }{\partial r}\left(r^{l+3}{u}^{l+1}\right)+E{\Delta }_l{w}^l=0,$$(A.5)

and onto the T_l direction $$\begin{array}{lll}\hfill & \hfill & B_{l-1}^l{r}^{l-1}\frac{\partial }{\partial r}\left(\frac{w^l}{r^{l-1}}\right)+B_{l+1}^l{r}^{-l-2}\frac{\partial }{\partial r}\left(r^{l+2}{w}^{l+1}\right)\hfill \\ \hfill & \hfill & -E{\Delta }_l{\Delta }_l(r{u}^l)=\sqrt{\frac{16\pi }{5}}\epsilon n^2(r){\delta }_{l2},\hfill \end{array}$$(A.6)

where δ_ij is the Kronecker symbol and ${\Delta }_l=\frac{1}{r}\frac{{\partial }^2}{\partial r^2}r-\frac{l(l+1)}{r^2}$. 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 $$\{\begin{array}{l}{A}_{l+1}^l=\frac{1}{(l+1)}\frac{1}{\sqrt{(2l+1)(2l+3)}},\hfill \\ A_{l-1}^l=\frac{1}{l}\frac{1}{\sqrt{(2l-1)(2l+1)}},\hfill \\ B_{l+1}^l=\frac{l(l+1)(l+2)}{\sqrt{(2l+1)(2l+3)}},\hfill \\ B_{l-1}^l=\frac{l(l^2-1)}{\sqrt{(2l-1)(2l+1)}}.\hfill \end{array}$$(A.7)

A.2 Boundary conditions at the upper boundary

The azimuthal velocity written on the (R_l, S_l, T_l) basis reads $$u_{\phi }(r,\theta )=-{\displaystyle \sum _{l=1}^{+\infty }{w}^l}(r){\partial }_{\theta }{Y}_l(\theta ).$$(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 $$\begin{array}{lll}{u}_{\phi }(r,\theta )\hfill & =\hfill & w^{l=1}(r)\sin \theta \sqrt{\frac{3}{4\pi }}\hfill \\ \hfill & \hfill & -w^{l=3}(r)\left(4\sin \theta +5{\sin }^3\theta \right)\sqrt{\frac{7}{4\pi }}\frac{3}{2}.\hfill \end{array}$$(A.9)

In orderto set the shear described by expression (6), we then set at the surface $$w^{l=1}(r=1)=-b\frac{4}{5}\sqrt{\frac{4\pi }{3}},$$(A.10)

and $$w^{l=3}(r=1)=-\frac{b}{5}\frac{2}{3}\sqrt{\frac{4\pi }{7}}.$$(A.11)

All other l-components 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

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 $r=\{0.5,0.75,1-\frac{\sqrt{E}}{2}\}$, θ = π∕2 and for b= +1. Scaling laws are {ur,  uθ}∝ Ed, 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 low-mass stars during the main sequence (Xc = 0.5) for stellar masses in the range [0.4 − 1.4]M⊙ from MESA models. Bottom: amplitude of the right-hand side of Eq. (34) (to be multiplied by $\mathcal{R}{o}^{-1}$) 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 $\mathcal{R}o=\{10,1,{10}^{-1},{10}^{-2},{10}^{-3}\}$ (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 core-to-surface 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) solar-type 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 solar-like stars if taking only into account rotation and induced large-scale flows.

In the text