© ESO, 2014

1. Introduction

The influence of rotation has long been known to be crucial for understanding mixing in radiative regions of stars and in interpreting the observed surface abundances (Strittmatter 1969). Many (if not all) stellar evolution codes now include some modelling of the transport of angular momentum and chemical elements by the flows induced by rotation in the stably stratified radiative zones. The difficulty is that stellar evolution codes are one-dimensional while fluid flows are generally multi-dimensional. Basically, two types of modelling are currently used: one is based on a simple (turbulent) diffusion process (e.g. Pinsonneault 1997) and the other includes a first-order modelling of meridional advection, distinguishing the transport of angular momentum and that of chemicals (Zahn 1992). However, both of them need adjustment of diffusion coefficients with respect to observations. While this modelling has succeeded in explaining various features of abundance patterns or evolutionary effects like the surface abundance of lithium as a function of mass (Charbonnel & Talon 1999), the (relative) high number of red super-giants in low-metallicity galaxies (Maeder & Meynet 2001), or the ratio of type Ibc to type II supernovae (Meynet & Maeder 2005), many recent results challenge our understanding of this so-called rotational mixing. One of the most famous is the distribution of Large Magellanic Cloud B-stars in a diagram plotting the nitrogen abundance versus the rotational velocity (the so-called Hunter diagram). As shown by Brott et al. (2011), many slowly rotating stars show an over abundance of nitrogen compared to the predictions of models, while some fast rotating stars show much less nitrogen than expected.

In order to progress in this difficult problem, it is clear that a better view of the dynamics of rotating stars is needed. Many processes contribute to rotational mixing. Let us recall that the timescale of element transport inside stars is essentially controlled by the radiative zones. Indeed, convective regions mix almost instantaneously. When the star is non-rotating, elements may migrate through radiative regions thanks to microscopic diffusive processes or propagating waves. In rotating stars, radiative regions are no longer at rest: beyond the global rotation baroclinic flows arise as a differential rotation and an associated meridional circulation (Rieutord 2006b). While meridional currents can obviously carry elements from deep to surface layers, differential rotation can also contribute to the transport through shear driven instabilities and associated turbulence. As pointed out by Zahn (1992), baroclinicity may be helped by angular momentum losses resulting from mass loss. To be complete, a sequence of gravitational contraction may also drive a redistribution of angular momentum and elements during the pre-main sequence or just at the end of the main sequence.

Going beyond the above mentioned modelling of rotational mixing requires relaxing the spherical symmetry of the models. A first step in this direction is to simplify the structure of the stars so as to focus on their interior hydrodynamics. This line of research was followed in Rieutord (2006a). The star was reduced to a spherical ball filled with an incompressible fluid to study the properties of the baroclinic flows. Besides giving a simplified view of the dynamic processes controlling a radiative zone, this work has led to a simplified set-up for boundary conditions in complete (compressible) two-dimensional models of fast rotating stars (Espinosa Lara & Rieutord 2013). We address the effect of mass contraction/expansion of a star or a giant planet on the global flows that affect the envelope of the rotating object. Although this driving results from the compressibility of the fluid, we decided to investigate its consequences using an incompressible fluid. The contraction/expansion of the envelope is mimicked by a core that absorbs or injects matter at a constant rate in a surrounding stably or neutrally stratified envelope. We wish to determine the resulting flows and the circumstances when it overwhelms the eventual baroclinic flows. This will complete the work of Rieutord & Beth (2014) who studied this same competition in the case of a spin-down driven by mass losses.

The paper is organized as follows: in Sect. 2 we describe the model and the questions of fluid dynamics that are addressed. We next consider the case where outer boundary conditions are no-slip and use the ensuing analytical solutions to enlight the dynamics (Sect. 3). In Sect. 4, we discuss the presumably more realistic case with stress-free boundary conditions. Conclusions and a discussion of extrapolations of the results to models with a compressible fluid follow.

2. The model

2.1. Description

To model simply the contracting/expanding star or giant planet, we consider a self-gravitating incompressible viscous fluid enclosed between two spherical shells. These shells are assumed not to be distorted by rotation. The inner shell may represents a core-envelope boundary through which a flux of matter is imposed. This flux is described by a uniform radial velocity V_s at this interface. For a contracting envelope, V_s is negative. Note that as the fluid of the envelope enters the core and as its radius is assumed fixed, the core’s density linearly increases with time, namely $\begin{array}{ccc}{\mathit{\rho }}_{\mathrm{core}}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}\left(\mathrm{1}\mathrm{+}\frac{\mathit{t}}{{\mathit{t}}_{\mathrm{s}}}\right)& & \end{array}$where we introduced the initial density of the core ρ₀ and the “suction time”, which is the time the core needs to increase its density by a factor 2, namely $\begin{array}{ccc}{\mathit{t}}_{\mathrm{s}}\mathrm{=}\frac{{\mathit{R}}_{\mathrm{core}}}{\mathrm{3}{\mathit{V}}_{\mathrm{s}}}\begin{array}{c}\mathrm{˜}\\ \mathit{\rho }\end{array}& & \end{array}$(1)where $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\rho }\end{array}\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{0}}}{\mathit{\rho }}$. This time is of same order of magnitude as the Kelvin-Helmholtz time. Here, ρ is the density of the envelope. We note that changing the sign of V_s can readily describe an expanding envelope due to a wind. The outer bounding sphere is of constant radius and lets matter through so as to ensure mass conservation in the envelope.

2.2. A digest of the following fluid dynamics

The remainder of the paper is fluid dynamics. A summary of the various steps are described below and the astrophysical implications of the results are addressed at the end of the paper.

As well known, the contraction of the star with the conservation of angular momentum induces a global acceleration of the rotation rate, namely a spin-up, of the star. The problem of spin-up flows has been considered many times in fluid dynamics literature (see the review of Duck & Foster 2001). However, these studies have been mostly motivated by engineering applications and therefore have considered a driving by boundaries and not, as in our case, by a radial flow. The influence of a stable stratification, that we need to know to deal with radiative regions of stars, has been more seldom considered in astrophysically relevant geometries. The most relevant study is certainly the work of Friedlander (1976), who considered the spin-down of the radiative core of the Sun driven by Reynolds stresses at the interface of the convective and radiative regions. So, here too, the driving of the flows is by boundaries. Friedlander’s model is much simplified (as is ours) as it uses the Boussinesq approximation (i.e. neglecting the compressibility of the fluid). It also neglects baroclinic flows associated with this set-up. This work however establishes that the spin-down (or spin-up) timescale of a stably stratified fluid is the classical Eddington-Sweet timescale, namely the product of the heat diffusion timescale (also known as Kelvin-Helmholtz time scale) and the factor , where is the Brunt-Väisälä frequency in the radiative region and Ω is the rotation rate. In slowly rotating stars, the Eddington-Sweet timescale is much longer than the Kelvin-Helmholtz, but in fast rotators that are considered here, these timescales are similar. Hence, the spin-up driven by a gravitational contraction has never been considered in the literature as far as we know.

The first step of our analysis considers the case of a neutrally stratified envelope like that found in a star with an outer convection zone. The translation of this constant entropy medium into the incompressible fluid model is the simple constant density fluid. This simple model allows us to derive an asymptotic solution for small Ekman numbers (i.e. small viscosity as appropriate for stellar applications). This solution shows that the spinning up core controls the flow inside its tangent cylinder and gives the amplitude of the quasi-steady flow there. Outside the tangent cylinder we find that the solution depends very much on the outer boundary conditions. We analyse the rather artificial case where no-slip conditions are imposed on the outer boundary, since this case also allows for the derivation of an analytic asymptotic solution. Moreover, with these boundary conditions, we can also solve the case where the envelope is stably stratified (but not contracting) and derive the associated differential rotation forced by the baroclinic torque.

The competition between the two forcings (spin-up and baroclinicity) necessarily arises in the pre-main-sequence phase of an intermediate mass star. We indeed expect that for such a mass range, an outer radiative envelope sets in after the birth line, likely after the disappearance of convective flows (e.g. Maeder 2009).

With our simplified model we can appreciate the result of this competition. As a first step, still using the artificial no-slip outer boundary condition, we compare the amplitudes of the baroclinic flow and the contraction-driven spin-up flow when they are taken separately. This comparison gives us a criterion on the parameters of the star. Using numerical solutions of the full problem, including the two forcing mechanisms simultaneously, we confirm the validity of the criterion.

Our next step is to leave aside the no-slip outer boundary condition and concentrate on the more realistic stress-free condition. In such a case, the derivation of an asymptotic analytic solution is much more involved and we had to resort to numerics. Numerical solutions show that outside the tangent cylinder of the core, a steady spin-up flow of an unstratified fluid driven by contraction has an amplitude that scales with the ratio of two small parameters, namely Ro/E, where Ro is the Rossby number measuring the driving and E is the Ekman number measuring the viscosity. With stellar parameters this ratio is expected to be larger than unity showing that a steady solution is necessarily of large amplitude. In addition, the time needed to reach such a steady state is roughly the contraction time, thus suggesting that this steady state is actually never reached. This result prompted us to study the transient state of the spin-up flow first without stratification and then account for a stable stratification in the envelope. Before the transient state reaches a significant amplitude, it can be studied with linear equations that are easier to solve. This transient flow shows the emergence of a quasi-self-similar solution that simply grows in amplitude in the volume outside the core’s tangent cylinder. Since our original question is to discover whether such flow is able to supersede the baroclinic flow driven by the combination of rotation and stratification, we compared the two flows. The easy way is to compare the amplitudes of the flow taken separately and quite clearly we find that the contraction-induced spin-up rapidly overtakes the amplitude of a baroclinic flow. The numerical solution of the transient starting from an established baroclinic flow confirms this result.

We now present the detailed derivation of these fluid dynamics results. The astrophysical side of the problem is addressed in Sect. 5.

2.3. Equations of motion

In an inertial frame, the dynamics of an incompressible fluid enclosed within the two shells is governed by the equations of momentum and mass conservation: $\begin{array}{ccc}\{\begin{array}{c}\\ \\ \end{array}& & \end{array}$(2)where v is the velocity field, P the pressure, and ν the kinematic viscosity of the envelope. Let us now remove the bulk rotation and set $\begin{array}{ccc}{v}\mathrm{=}{\Omega }\mathrm{\wedge }{r}\mathrm{+}{w}\mathit{.}& & \end{array}$(3)We define Ω as the angular velocity of the core. Since the field v is axisymmetric, $\begin{array}{ccc}\mathrm{(}{v}\mathrm{·}{\nabla }\mathrm{)}{v}\mathrm{=}{\Omega }\mathrm{\wedge }\mathrm{(}{\Omega }\mathrm{\wedge }{r}\mathrm{)}\mathrm{+}\mathrm{2}{\Omega }\mathrm{\wedge }{w}\mathrm{+}\mathrm{(}{w}\mathrm{·}{\nabla }\mathrm{)}{w}& & \end{array}$where we recognise the centrifugal and Coriolis accelerations. The centrifugal potential is gathered with the pressure into Π. Substituting (3) in the set of Eqs. (2), we find that the relative velocity w verifies $\begin{array}{ccc}\{\begin{array}{c}\\ \\ \end{array}& & \end{array}$(4)Note that the RHS-term now depends on the acceleration of the rotation rate of the core where needs to be determined.

2.4. Scaled equations and linearization

We scale the equations using the Kelvin-Helmholtz timescale $\frac{\mathit{R}}{{\mathit{V}}_{\mathrm{s}}}$, which is the timescale of the complete absorption of the envelope by the core. The length scale R is the outer radius of the envelope, and the velocity scale is the suction velocity w = V_su. The system (5) now reads $\begin{array}{ccc}\{\begin{array}{c}\\ \\ \end{array}& & \end{array}$(5)where p is the reduced pressure. The following dimensionless numbers appear: $\begin{array}{ccc}\mathit{\omega ̇}\mathrm{=}\frac{\mathrm{\Omega ̇}\mathit{R}}{\mathrm{2}\mathrm{\Omega }{\mathit{V}}_{\mathrm{s}}}\mathit{,} \mathit{Ro}\mathrm{=}\frac{{\mathit{V}}_{\mathrm{s}}}{\mathrm{2}\mathrm{\Omega }\mathit{R}}\mathit{,} \mathit{E}\mathrm{=}\frac{\mathit{\nu }}{\mathrm{2}\mathrm{\Omega }{\mathit{R}}^{\mathrm{2}}}\mathrm{·}& & \end{array}$(6)Here we introduced the non-dimensional acceleration of the core rotation rate , a Rossby number Ro and the Ekman number E, which measures the viscosity of the envelope.

Obviously, the suction velocity is very small compared to the rotation velocity. Hence, we expect Ro ≪ 1. As a first step setting Ro = 0 seems reasonable as long as Rou is less than unity (so as to be able to neglect quadratic terms). Thus, we are left with a steady state problem that describes the quasi-steady evolution of the system as long as the non-dimensional time τ verifies: $\begin{array}{ccc}\mathit{Ro}\mathrm{\ll }\mathit{\tau }\mathrm{\ll }{\mathit{\tau }}_{\mathrm{s}}& & \end{array}$(7)where ${\mathit{\tau }}_{\mathrm{s}}\mathrm{=}\frac{\mathit{\eta }\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\rho }\end{array}}{\mathrm{3}}$ is the scaled suction time ($\mathit{\eta }\mathrm{=}\frac{{\mathit{R}}_{\mathrm{core}}}{\mathit{R}}$ is the scaled inner radius). The left part of the inequality τ ≫ Ro means that we neglect the transients corresponding to a few rotation periods where boundary layers form. Likewise, τ ≪ τ_s means that the rotation rate has not been changed, namely that .

2.5. The acceleration of the core

2.5.1. General equation

Equation (5) need the expression of . By absorbing the envelope, the mass of the core grows, as its angular momentum. Evolution of the angular momentum L_z of the core is governed by $\begin{array}{ccc}\frac{\mathrm{d}{\mathit{L}}_{\mathit{z}}}{\mathrm{d}\mathit{t}}\mathrm{=}{{e}}_{\mathit{z}}\mathrm{·}\left\{\mathrm{-}{\mathrm{\int }}_{\mathrm{\left(}\mathit{S}\mathrm{\right)}}\mathrm{(}{r}\mathrm{\wedge }\mathit{\rho }{v}\mathrm{)}{v}\mathrm{·}{d}{S}\mathrm{+}{\mathrm{\int }}_{\mathrm{\left(}\mathit{S}\mathrm{\right)}}{r}\mathrm{\wedge }\mathrm{\left[}\mathit{\sigma }\mathrm{\right]}{d}{S}\right\}\mathit{.}& & \end{array}$(8)The first integral is the flux of incoming angular momentum and the second integral is the viscous torque applied on the core surface. The tensor [ σ ] is the stress tensor. We find $\begin{array}{ccc}{{e}}_{\mathit{z}}\mathrm{·}\mathrm{(}{r}\mathrm{\wedge }\mathrm{[}\mathit{\sigma }\mathrm{\left]}{{e}}_{{r}}\mathrm{\right)}\mathrm{=}\mathit{r}\sin \mathit{\theta }{\mathit{\sigma }}_{\mathit{r\phi }}\mathrm{=}\mathit{r}\sin \mathit{\theta \rho \nu }{\frac{\mathit{\partial }{\mathit{w}}_{\mathit{\phi }}}{\mathit{\partial r}}|}_{\mathit{r}\mathrm{=}{\mathit{R}}_{\mathrm{core}}}\mathit{.}& & \end{array}$(9)Besides, for a sphere of mass M_core and radius ηR: $\begin{array}{ccc}{L}\mathrm{=}\mathit{I}{\Omega }\mathrm{=}\frac{\mathrm{2}}{\mathrm{5}}{\mathit{M}}_{\mathrm{core}}\mathrm{(}\mathit{\eta R}{\mathrm{)}}^{\mathrm{2}}{\Omega }\mathrm{·}& & \end{array}$(10)With the previous scaling, we obtain $\begin{array}{ccc}\begin{array}{c}\\ & & \\ \end{array}& & \end{array}$(11)The evaluation of the remaining integral needs the expression of the azimuthal flow u_φ at the core-envelope boundary r = η, namely in the Ekman boundary layers.

2.5.2. Boundary layer analysis

First of all, we change the boundary conditions with mass flux to ordinary boundary condition by making the substitution $\begin{array}{ccc}{u}\mathrm{=}{{u}}^{{\prime }}\mathrm{-}\frac{{\mathit{\eta }}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{2}}}{{e}}_{\mathit{r}}\mathit{.}& & \end{array}$The new velocity field u^′ verifies $\begin{array}{ccc}\{\begin{array}{c}\\ & & \\ & & \\ & & \end{array}& & \end{array}$(12)with the boundary conditions $\begin{array}{ccc}{{u}}^{\mathrm{\prime }}\mathrm{=}{0} \mathrm{on} \mathit{r}\mathrm{=}\mathit{\eta }& & \end{array}$and $\begin{array}{ccc}{{e}}_{\mathit{r}}\mathrm{\times }\mathrm{\left[}\mathit{\sigma }\mathrm{\right]}{{e}}_{\mathit{r}}\mathrm{=}{0}\mathit{,} \mathrm{and} {\mathit{u}}_{\mathit{r}}^{\mathrm{\prime }}\mathrm{=}\mathrm{0} \mathrm{on} \mathit{r}\mathrm{=}\mathrm{1.}& & \end{array}$These latter conditions describe stress-free boundary conditions. Indeed, we assume that the upper layers outside the outer shell have unimportant dynamic effects and are just providing/absorbing some mass. In the following, we drop the prime on the velocity field.

As shown in Espinosa Lara & Rieutord (2013), the Ekman number in stars is always very small thus leading to the formation of boundary layers near the boundaries. To derive the expression of the flow, we therefore examine the asymptotic case E ≪ 1. We first note that if we consider the inviscid case E = 0, the set of Eqs. (12) is solved by $\begin{array}{ccc}\{\begin{array}{c}\\ & & \\ & & \end{array}& & \end{array}$(13)In the azimuthal direction, we look for a flow such as as dictated by the Taylor-Proudman theorem (e.g. Greenspan 1969). Such a flow does not verify the no-slip boundary conditions at the interface r = η. It needs boundary layer corrections so as to satisfy $\mathit{u̅}{}_{\mathrm{0}}\mathrm{+}\begin{array}{c}{}^{\mathrm{􏽥}}\\ {\mathit{u}}_{\mathrm{0}}\end{array}\mathrm{=}\mathrm{0}$. The bar refers to the solution within the envelope i.e. outside the boundary layer and tilded quantities are for boundary layer corrections.

Following Rieutord (2006a), we write the boundary layer corrections as: $\begin{array}{ccc}\mathrm{(}{n}\mathrm{\wedge }\begin{array}{c}\mathrm{􏽥}\\ {\mathit{u}}_{\mathrm{0}}\end{array}\mathrm{+}\mathrm{i}\begin{array}{c}\mathrm{􏽥}\\ {\mathit{u}}_{\mathrm{0}}\end{array}\mathrm{)}\mathrm{=}\mathrm{-}\mathrm{(}{n}\mathrm{\wedge }\mathit{u̅}\mathrm{0}\mathrm{+}\mathrm{i}\mathit{u̅}\mathrm{0}{\mathrm{)}}_{\mathit{\alpha }\mathrm{=}\mathrm{0}}\exp \mathrm{(}\mathrm{-}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{i}\mathrm{\left)}\mathit{\alpha }\mathrm{\right)}& & \end{array}$(14)where $\mathit{\alpha }\mathrm{=}\mathrm{(}\mathit{r}\mathrm{-}\mathit{\eta }\mathrm{)}\sqrt{\frac{\mathrm{\left|}\cos \mathit{\theta }\mathrm{\right|}}{\mathrm{2}\mathit{E}}}$ and n = −e_r. We keep only the decreasing solution. The corrections thus read $\begin{array}{ccc}\{\begin{array}{c}\\ & & \\ & & \end{array}& & \end{array}$(15)As in Rieutord (2006a), mass conservation gives the relation between the Ekman pumping $\begin{array}{c}{\mathrm{􏽥}}^{}\\ {\mathit{u}}_{\mathit{r}}\end{array}$ and the geostrophic flow U. It yields$\begin{array}{ccc}\mathrm{1}\mathrm{+}\mathrm{(}\mathrm{2}\mathrm{-}\mathrm{3}{\sin }^{\mathrm{2}}\mathit{\theta }\mathrm{)}\mathit{\eta }\mathit{\omega ̇}\mathrm{=}\frac{\mathrm{1}}{\mathit{\eta }\sin \mathit{\theta }}\sqrt{\frac{\mathit{E}}{\mathrm{2}}}\frac{\mathit{\partial }}{\mathit{\partial \theta }}\left(\frac{\sin \mathit{\theta }\mathit{U}\mathrm{\left(}\mathit{\eta }\sin \mathit{\theta }\mathrm{\right)}}{\sqrt{\mathrm{\left|}\cos \mathit{\theta }\mathrm{\right|}}}\right)& & \end{array}$(16)where we keep only $\mathrm{𝒪}\mathrm{\left(}\sqrt{\mathit{E}}\mathrm{\right)}$ terms. Integrating with respect to θ, we find $\begin{array}{ccc}\mathit{U}\mathrm{\left(}\mathit{\eta }\sin \mathit{\theta }\mathrm{\right)}\mathrm{=}\mathit{\eta }\sqrt{\frac{\mathrm{2}}{\mathit{E}}}\frac{\sqrt{\mathrm{\left|}\cos \mathit{\theta }\mathrm{\right|}}}{\sin \mathit{\theta }}\mathrm{(}\mathrm{1}\mathrm{-}\cos \mathit{\theta }\mathrm{+}\mathit{\eta }\mathit{\omega ̇}\cos \mathit{\theta }{\sin }^{\mathrm{2}}{\mathit{\theta }}^{\mathrm{)}}\mathit{.}& & \end{array}$(17)Note that this expression defines U only within the tangent cylinder of the core defined as s = rsinθ = η (s is the radial cylindrical coordinate).

Near the surface r = η i.e. within the boundary layer, the shear is dominated by the boundary layer correction. It simplifies the computation of the radial derivative of the azimuthal flow u_φ which reads $\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{u}}_{\mathit{\phi }}}{\mathit{\partial r}}\mathrm{=}\frac{\mathit{\eta }}{\mathit{E}}\mathrm{\left|}\cot \mathit{\theta }\mathrm{\right|}\mathrm{(}\mathrm{1}\mathrm{-}\cos \mathit{\theta }\mathrm{+}\mathit{\eta }\mathit{\omega ̇}\cos \mathit{\theta }{\sin }^{\mathrm{2}}{\mathit{\theta }}^{\mathrm{)}}\mathit{.}& & \end{array}$(18)Integral in (11) can now be evaluated: $\begin{array}{ccc}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }}{\sin }^{\mathrm{2}}\mathit{\theta }{\frac{\mathit{\partial }{\mathit{u}}_{\mathit{\phi }}}{\mathit{\partial r}}|}_{\mathit{r}\mathrm{=}\mathit{\eta }}\mathrm{d}\mathit{\theta }\mathrm{=}\mathit{\eta }\frac{\mathrm{2}}{\mathrm{3}\mathit{E}}\left(\frac{\mathrm{1}}{\mathrm{2}}\mathrm{+}\frac{\mathrm{2}}{\mathrm{5}}\mathit{\eta }\mathit{\omega ̇}\right)\mathit{.}& & \end{array}$(19)We can now derive the acceleration of the angular velocity of the core. Considering quasi-steady solutions that arise when Ro ≪ τ ≪ τ_s, (11) leads to $\begin{array}{ccc}\mathit{\omega ̇}\mathrm{=}\frac{\mathrm{9}}{\mathrm{4}\begin{array}{c}\mathrm{˜}\\ \mathit{\rho }\end{array}\mathit{\eta }}\mathit{,}& & \end{array}$(20)which completes the Eqs. (12).

The foregoing solution (17) shows that the differential rotation driven by the mass contraction scales as . It means that the linear regime that we solved is valid only when $\mathit{Ro}\mathrm{\ll }\sqrt{\mathit{E}}$, which is actually the case (see below). The foregoing solution however describes the fluid flow only within the tangent cylinder of the core.

Outside the cylinder, the solution depends on the outer boundary conditions and on the Stewartson layer that lies upon the cylinder. This makes the global solution quite involved, all the more that we should also account for a possible stable stratification of the envelope. Indeed, during the PMS phase of intermediate mass stars, the envelope is completely radiative. Therefore the contraction-induced differential rotation competes with the differential rotation induced by the baroclinicity of the envelope.

Before getting any further, we need to evaluate stellar numbers that have appeared.

2.6. Orders of magnitude

As a test case of the foregoing problem, we first consider the contraction of a fully radiative 3 M_⊙ PMS star. On the birthline, the star’s surface temperature is around T_∗ ~ 5600 K, its luminosity L_∗ ~ 10²L_⊙ and its radius is R_∗ ~ 10¹⁰ m. The young star contracts on the Kelvin-Helmholtz time upon the PMS (Henyey) track, namely: $\begin{array}{ccc}{\mathit{t}}_{\mathrm{KH}}\mathrm{=}\frac{\mathit{G}{\mathit{M}}^{\mathrm{2}}}{\mathit{RL}}\mathit{,}& & \end{array}$according to Maeder (2009). This leads to t_KH ~ 2.6 × 10⁵ yr.

Setting arbitrarily R_core = 0.15 R_∗, we find V_s ~ 6 × 10^-5 m s^-1. Considering a rotation velocity of 10 km s^-1 (such that near the end of the PMS after a gravitational contraction at constant angular momentum we obtain a star like HD 37806 studied by Boehm & Catala 1995), we find a Rossby number $\begin{array}{ccc}\mathit{Ro}\mathrm{~}\mathrm{3}\mathrm{\times }{\mathrm{10}}^{-9}& & \end{array}$that is very small as expected.

The estimate of the Ekman number requires a value of the kinematic viscosity. If we use Zahn’s prescription (Zahn 1992) for a turbulent viscosity, we obtain ν ~ 10⁴ m² s^-1 and thus $\begin{array}{ccc}\mathit{E}\mathrm{~}{\mathrm{10}}^{-10}\mathit{.}& & \end{array}$With the radiative viscosity (e.g. Espinosa Lara & Rieutord 2013), ν ~ 10² m² s^-1, we get $\begin{array}{ccc}\mathit{E}\mathrm{~}{\mathrm{10}}^{-12}\mathit{.}& & \end{array}$For both values the condition $\mathit{Ro}\mathrm{\ll }\sqrt{\mathit{E}}$ is satisfied.

During the contraction, the quasi-steady state within the cylinder is reached after a spin-up time, namely after $\frac{\mathrm{(}\mathrm{2}\mathrm{\Omega }{\mathrm{)}}^{-1}}{\sqrt{\mathit{E}}}$. Using the previous numbers, we find that this state occurs after ~10³ or 10⁴ years so rather soon after the start of contraction. Therefore, within the tangent cylinder, we can neglect the transient phase.

If we consider the contracting envelope of a Jupiter-like giant planet, Ekman numbers are also very small although with larger uncertainties: Ogilvie & Lin (2004) give E ~ 10^-7 while Wu (2005) suggest E ~ 10^-13. The typical contraction time of giant planets is over 1 Gyrs (Fortney & Nettelmann 2010) so that Ro ~ P_rot/t_contraction ~ 10^-12. The condition $\mathit{Ro}\mathrm{\ll }\sqrt{\mathit{E}}$ is also easily satisfied.

2.7. Adding stratification

2.7.1. Scaled equations

To account for a stable stratification in the envelope, we now generalize the set of Eqs. (12) by taking the buoyancy force and the equation for temperature fluctuations into account. In PMS stars, the stable stratification of the envelope may come after a convective episode and thus may be evolving with time. To simplify and get an upper bound on the effects of stratification we impose the Brunt-Väisälä frequency as constant in time.

To be consistent with the foregoing model that uses an incompressible fluid we use the Boussinesq approximation. Following Rieutord (2006a) and combining with (12) we find $\begin{array}{ccc}\{\begin{array}{c}\\ \\ \\ \\ \\ \\ \end{array}& & \end{array}$(21)We use the same scales and notations as Rieutord (2006a). The temperature perturbation is scaled as δT = ϵT_∗θ_T where $\mathit{ϵ}\mathrm{=}\frac{{\mathrm{\Omega }}^{\mathrm{2}}\mathit{R}}{{\mathit{g}}_{\mathrm{s}}}$ is the ratio of centrifugal acceleration to surface gravity. Recall that the centrifugal acceleration is driving the baroclinic flow. The Brunt-Väisälä profile ${\mathit{n}}_{\mathit{T}}^{\mathrm{2}}$ is scaled with ${\mathrm{𝒩}}^{\mathrm{2}}\mathrm{=}\frac{\mathit{\alpha }{\mathit{T}}_{\mathrm{\ast }}{\mathit{g}}_{\mathit{S}}}{\mathit{R}}$ where α is the dilation coefficient and g_s the surface gravity.

The dimensionless number B monitors the ratio of the forcings. From the expression of the scaling of the baroclinic flow (see Rieutord 2006a), we have $\mathit{B}\mathrm{=}\frac{\mathit{ϵ}{\mathrm{𝒩}}^{\mathrm{2}}\mathit{R}}{\mathrm{2}\mathrm{\Omega }{\mathit{V}}_{\mathrm{s}}}\mathrm{·}$(22)Finally, the dimensionless number $\begin{array}{c}\mathrm{˜}\\ {\mathit{E}}_{\mathit{T}}\end{array}\mathrm{=}\frac{\mathit{E}}{\mathit{\lambda }} \mathrm{with} \mathit{\lambda }\mathrm{=}\mathrm{𝒫}\frac{{\mathrm{𝒩}}^{\mathrm{2}}}{\mathrm{4}{\mathrm{\Omega }}^{\mathrm{2}}}$(23)measures heat diffusion. The dimensionless number is the Prandtl number. For fast rotators, λ is a small parameter that we set to 10^-4, following the estimate of Rieutord & Beth (2014).

Based on stellar models, a typical profile of the Brunt-Väisälä frequency is shown in Fig. 1. We use the polynomial expression $\begin{array}{ccc}\{\begin{array}{c}\\ & & \\ \\ & & \end{array}& & \end{array}$(24)to represent this function. The coefficients α, β and γ result from the polynomial fit.

Fig. 1 Typical and scaled Brunt-Väisälä frequency profile as a function of the normalized radius with η = 0.15.

3. An interesting solution with rigid outer boundary conditions

Before getting into the full numerics, it is interesting to consider the case where the outer bounding sphere of the envelope is rigidly rotating at the same rate as the inner core. Outer no-slip conditions can be expected if a turbulent layer threaded by magnetic fields covers the stellar surface (see Rieutord & Beth 2014), however the synchronism between this layer and the surface is here an ad hoc assumption (which can be easily removed). The interesting point that is addressed below comes from the simple analytical solution that can be derived for the flow outside the tangent cylinder and offers an interesting view of the properties of the system.

3.1. The steady mass contraction induced flow

With no-slip conditions on the outer boundary r = 1, we may easily derive the expression of the geostrophic flow¹ out of the tangent cylinder of the core. When no stratification is present, the azimuthal velocity reads $\begin{array}{ccc}\mathit{U}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}\sqrt{\frac{\mathrm{2}}{\mathit{E}}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{s}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{4}}\left(\frac{{\mathit{\eta }}^{\mathrm{2}}}{\mathit{s}}\mathrm{-}\mathit{\omega ̇}\mathit{s}\right) \mathit{s}\mathrm{\ge }\mathit{\eta }\mathit{.}& & \end{array}$(25)As shown in Fig. 2, this analytical solution nicely matches the numerical solution.

Fig. 2 Comparison between numerical (solid line) and analytical (star line) solutions of the geostrophic flow at the equator ${\mathit{u}}_{\mathit{\phi }}\mathrm{(}\mathit{r,\theta }\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}\mathrm{)}$ for E = 10-7, η = 0.15, $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\rho }\end{array}\mathrm{=}\mathrm{10}$ without stratification when the envelope is rigidly rotating at the same rate as the inner one.

3.2. The steady baroclinic flow

Let us now consider the opposite case where a pure baroclinic flow (no contraction) meets no-slip boundary conditions. It verifies $\begin{array}{ccc}\{\begin{array}{c}\\ \\ \\ \end{array}& & \end{array}$(26)As shown in Rieutord (2006a), neglecting temperature perturbations and viscosity, the φ-component of the vorticity Eq. (26) leads to $\begin{array}{ccc}\{\begin{array}{c}\\ & & \\ & & \end{array}& & \end{array}$(27)where F(s) is a geostrophic solution determined by the boundary conditions. In order to get the expression of F(s), we write $\begin{array}{ccc}\mathit{u̅}\mathrm{0}\mathrm{=}\left(\mathrm{-}\mathit{sB}\begin{array}{c}\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{{\mathrm{\int }}_{\mathit{r}}^{\mathrm{1}}\frac{{\mathit{n}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathit{r}}\mathrm{d}\mathit{r}}\\ \mathit{g}\mathrm{\left(}\mathit{r}\mathrm{\right)}\end{array}\mathrm{+}\mathit{F}\mathrm{\left(}\mathit{s}\mathrm{\right)}\right){{e}}_{\mathit{\phi }}& & \end{array}$(28)and look for the boundary layer corrections at r = 1. These are $\begin{array}{ccc}\begin{array}{c}\mathrm{􏽥}\\ {\mathit{u}}_{\mathrm{0}}\end{array}\mathrm{=}\mathrm{ℐ}\mathit{m}\mathrm{\{}\mathrm{-}\mathrm{(}{{e}}_{\mathit{r}}\mathrm{\wedge }\mathit{u̅}\mathrm{0}\mathrm{+}\mathrm{i}\mathit{u̅}\mathrm{0}{\mathrm{)}}_{\mathit{r}\mathrm{=}\mathrm{1}}{\mathit{e}}^{\mathrm{-}\mathit{\alpha }\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{i}\mathrm{)}}\mathrm{\}}& & \end{array}$where $\mathit{\alpha }\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{r}\mathrm{)}\sqrt{\frac{\mathrm{\left|}\cos \mathit{\theta }\mathrm{\right|}}{\mathrm{2}\mathit{E}}}\mathrm{=}\mathit{\zeta }\sqrt{\frac{\mathrm{\left|}\cos \mathit{\theta }\mathrm{\right|}}{\mathrm{2}}}$. Using (27) with the foregoing expression yields $\begin{array}{ccc}\{\begin{array}{c}\\ \\ \end{array}& & \end{array}$(29)We note that g(1) = 0. Mass conservation implies $\begin{array}{ccc}\frac{\mathrm{1}}{\sqrt{\mathit{E}}}\frac{\mathit{\partial }\begin{array}{c}\mathrm{􏽥}\\ {\mathit{u}}_{\mathit{r}}\end{array}}{\mathit{\partial \zeta }}\mathrm{=}\frac{\mathrm{1}}{\sin \mathit{\theta }}\frac{\mathit{\partial }\mathrm{\left(}\sin \mathit{\theta }\begin{array}{c}\mathrm{􏽥}\\ {\mathit{u}}_{\mathit{\theta }}\end{array}\mathrm{\right)}}{\mathit{\partial \theta }}& & \end{array}$(30)so that $\begin{array}{ccc}\begin{array}{c}\mathrm{􏽥}\\ {\mathit{u}}_{\mathit{r}}\end{array}\mathrm{\left(}\mathrm{1}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\sqrt{\mathit{E}}}{\sin \mathit{\theta }}\frac{\mathit{\partial }}{\mathit{\partial \theta }}\left(\sin \mathit{\theta }{\mathrm{\int }}_{\mathit{\zeta }}^{\mathrm{\infty }}\begin{array}{c}\mathrm{􏽥}\\ {\mathit{u}}_{\mathit{\theta }}\end{array}\mathrm{d}\mathit{\zeta }\right)\mathit{.}& & \end{array}$(31)Since $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{u}}_{\mathit{r}}\end{array}\mathrm{+}{\overline{\mathit{u}}}_{\mathit{r}}\mathrm{=}\mathrm{0}$ at r = 1, we finally find $\begin{array}{ccc}\mathrm{-}\mathit{u̅}\mathit{r}\mathrm{(}\mathit{r}\mathrm{=}\mathrm{1}\mathrm{)}\mathrm{=}\sqrt{\frac{\mathit{E}}{\mathrm{2}}}\frac{\mathrm{1}}{\sin \mathit{\theta }}\frac{\mathit{\partial }}{\mathit{\partial \theta }}\left(\frac{\sin \mathit{\theta }}{\sqrt{\cos \mathit{\theta }}}\mathit{F}\mathrm{\left(}\sin \mathit{\theta }\mathrm{\right)}\right)\mathit{.}& & \end{array}$(32)However, the radial component of the vorticity equation (or the angular momentum equation) in the interior leads to $\begin{array}{ccc}\sin \mathit{\theta }{\overline{\mathit{u}}}_{\mathit{r}}\mathrm{+}\mathrm{cth}{\overline{\mathit{u}}}_{\mathit{\theta }}\mathrm{=}\mathit{E}{\mathrm{\Delta }}^{\mathrm{\prime }}{\overline{\mathit{u}}}_{\mathit{\varphi }}\mathit{.}& & \end{array}$Here ${\overline{\mathit{u}}}_{\mathit{\varphi }}\mathrm{=}\mathrm{-}\mathit{sBg}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{+}\mathit{F}\mathrm{\left(}\mathit{s}\mathrm{\right)}$ so that consistency of the solution with (32) requires that $\begin{array}{ccc}\mathit{F}\mathrm{\equiv }\mathrm{𝒪}\mathrm{\left(}\mathit{B}\sqrt{\mathit{E}}\mathrm{\right)}\mathit{.}& & \end{array}$It means that in the limit of vanishing Ekman numbers the function F can be neglected. Therefore, at leading order, the envelope differential rotation is dominated by the shellular flow: $\begin{array}{ccc}\mathit{u̅}\mathit{\phi }\mathrm{=}\mathrm{-}\mathit{sB}{\mathrm{\int }}_{\mathit{r}}^{\mathrm{1}}\frac{{\mathit{n}}^{\mathrm{2}}\mathrm{\left(}\mathit{x}\mathrm{\right)}}{\mathit{x}}\mathrm{d}\mathit{x}\mathrm{=}\mathrm{-}\mathit{sBg}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{.}& & \end{array}$Figure 3 shows the comparison of the analytical solution − sBg(r) with the numerical solution at the equator $\mathit{\theta }\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}$. The difference on the left edge comes from the fact that g(η) is not zero and would require an additional boundary layer correction.

3.3. Transition between the two steady flows: the baroclinic flow versus that induced by mass contraction

Fig. 3 Comparison between numerical (solid line) and analytical (stars) solutions of the geostrophic flow at the equator ${\mathit{u}}_{\mathit{\phi }}\mathrm{(}\mathit{r,\theta }\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}\mathrm{)}$ for E = 10-7, η = 0.15, $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\rho }\end{array}\mathrm{=}\mathrm{10}$, B = 1, and λ = 10-4 without mass contraction.

Fig. 4 Evolution of the azimuthal velocity ${\mathit{u}}_{\mathit{\phi }}\mathrm{(}\mathit{r,\theta }\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}\mathrm{)}$ for η = 0.15, $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\rho }\end{array}\mathrm{=}\mathrm{10}$, E = 10-7, λ = 10-4, and B = 0,104,105,107 (downward). The solid line is the numerical solution, the dashed line is the full analytical solution, and the two dot lines are the solutions of each flow. The transition is at B = 104, the differential rotation is then governed by baroclinicity as shown in Fig. 5. Boundary conditions are no-slip on both sides r = η and r = 1.

When gravitational contraction occurs in a baroclinic envelope, two drivings compete: the baroclinic torque and that induced by mass contraction. For a global view of this competition we resort to numerical solutions. The numerical method is detailed in appendix. The global problem is the superposition of the two flows: the flow induced by mass contraction (25) and the baroclinic flow (27). Since the system is linear, the full solution is a linear combination of both. At the equator it reads $\begin{array}{ccc}{\mathit{u}}_{\mathit{\phi }}\left(\mathit{r,\theta }\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}\right)\mathrm{=}\sqrt{\frac{\mathrm{2}}{\mathit{E}}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{r}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{4}}\left(\frac{{\mathit{\eta }}^{\mathrm{2}}}{\mathit{r}}\mathrm{-}\frac{\mathrm{9}}{\mathrm{4}\begin{array}{c}\mathrm{˜}\\ \mathit{\rho }\end{array}\mathit{\eta }}\mathit{r}\right)\mathrm{-}\mathit{rB}{\mathrm{\int }}_{\mathit{r}}^{\mathrm{1}}\frac{{\mathit{n}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathit{r}}\mathrm{d}\mathit{r}\mathit{.}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}& & \end{array}$(33)From the foregoing equation, we see that the differential rotation is governed by baroclinicity when $\begin{array}{ccc}\mathit{B}\mathrm{\gg }{\mathit{E}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}& & \end{array}$(34)It is shown in Fig. 4. However, we also know (from the angular momentum flux balance of a steady flow, see Rieutord 2006a) that the meridional circulation associated with baroclinic flows is , while the meridional circulation of the contraction-induced spin-up is . Hence, a baroclinic flow completely controls the dynamics as long as BE ≫ 1. Note that this inequality implies (34) since E ≪ 1. When the spin-up strengthens, the foregoing inequalities predict an intermediate regime $\begin{array}{ccc}{\mathit{E}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{\ll }\mathit{B}\mathrm{\ll }{\mathit{E}}^{-1}& & \end{array}$(35)where the differential rotation is of baroclinic origin (BE^1/2) but the meridional circulation is driven by contraction (BE ≪ 1). Finally, when BE^1/2 ≪ 1, the flow is fully controlled by the contraction induced spin-up.

This two step transition is confirmed by numerical solutions as illustrated in Figs. 5 and 6. There, for a given E, B is progressively increased and we clearly see the intermediate regime where differential rotation and meridional circulation are of different origin (third row).

4. The case of stress-free boundary conditions

The outer layers of the envelope are not rigidly attached to the core. Therefore the use of outer stress-free boundary conditions is more realistic. In this case, however, we no longer have access to an analytical expression of the flow in the envelope and have to resort to numerical solutions.

4.1. Scaling the steady mass contraction induced flow

Let us first study the steady solution of the spin-up flow. As shown in Fig. 7, it exhibits the typical cylindrical differential rotation of a dominating mass contraction flow. The equatorial surface region rotates faster than the core and the pole is slower. The meridional circulation displays two cells with a strong Stewartson layer at s = η (compared with the previous with no-slip conditions).

Fig. 5 Evolution of the differential rotation δΩ, the meridional circulation ψ for η = 0.15, $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\rho }\end{array}\mathrm{=}\mathrm{10}$, E = 10-7, λ = 10-4, and B = 0,104,105,107 (downwards). In the first column, difference in differential rotation is shown with contours: solid contours are faster than the core and dashed contours are slower than the core. When the spin-up flow dominates the differential rotation we obtain a fast equator while baroclinicity induces a fast pole. In the second column, meridional circulation is described with dotted lines for clockwise circulation (solid lines for counter-clockwise circulation). Boundary conditions are no-slip on both sides r = η and r = 1.

In Fig. 8, we show the amplitudes of the mass contraction flow at two positions: inside and outside of the tangent cylinder. Since the boundary conditions on the core are no-slip, the flow within the tangent cylinder is expected to be according to expression (17). For small cores η = 0.15 or 0.25, numerical solutions show that the Stewartson layer impacts the interior of the tangent cylinder and that the asymptotic state is reached only at extremely small values of the Ekman number E ≤ 10^-8. When the core is bigger, for instance η = 0.5, the E^{− 1 / 2} scaling inside the tangent cylinder clearly appears for all Ekman numbers less than 10^-5.

Outside the tangent cylinder, numerics show that the differential rotation is always when the outer boundary conditions are stress-free. This important amplitude indicates that the steady state may not be reached during the contraction phase and may not be studied with linear equations since quadratic terms are expected to be important, namely $\mathrm{𝒪}\mathrm{\left(}\frac{\mathit{Ro}}{{\mathit{E}}^{\mathrm{2}}}\mathrm{\right)}$.

4.2. The transient phase

The large amplitude of the steady state outside the tangent cylinder forces us to consider the time evolution of the solution of the mass contraction induced flow. To do so, we solve the set of Eqs. (5) with an order one implicit scheme (Euler’s method) so as to eliminate inertial waves and concentrate on secular evolution.

In Fig. 9, we plot a proxy of the amplitude of the differential rotation for various Ekman numbers and for no-slip and stress-free outer boundary conditions. With no-slip conditions, we see that the steady spin-up flow is quickly established and justifies the use of a steady solution. On the contrary, the use of outer stress-free boundary conditions leads to a much longer transient flow that lasts more than the typical timescale of the driving by gravitational contraction.

To go further it is interesting to characterize this transient flow with respect to the parameters of the problem. From the numerical solution we find that the transient duration τ_sf scales like ${\mathit{\tau }}_{\mathrm{sf}}\mathrm{\propto }\mathit{Ro}{\mathit{E}}^{-0.86}\mathit{.}$(36)This scaling of the Ekman number is very close to E^{− 6 / 7}, which is reminiscent of one of the scalings of the Stewarston layer in the spherical Couette flow (see Stewartson 1966). In these layers that surround the core along the tangent cylinder, a typical thickness is E^{2 / 7}. This might control the amplitude of the flow outside the tangent cylinder when stress-free outer conditions are met. The analysis of the Stewartson layer associated with this transient flow is difficult and beyond the scope of the present work.

Fig. 6 Azimuthal velocity uφ as a function of the normalized radius for λ = 10-4, η = 0.15, and $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\rho }\end{array}\mathrm{=}\mathrm{10}$. The solid line is the numerical solution, dotted lines are the analytical solutions of the spin-up flow and the baroclinic flow, the dashed line is the sum of both. Downwards: E = 10-5,10-6,10-7 and the transition value on the differential rotation is B = 103,3 × 103,104, respectively. Beyond this threshold, the baroclinicity governs the differential rotation. Boundary conditions are no-slip on both sides r = η and r = 1. Note that as the Ekman number decreases, the discrepancy between the numerical and the analytical solutions decreases as well.

Fig. 7 Differential rotation and meridional circulation for E = 10-7, η = 0.15, and $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\rho }\end{array}\mathrm{=}\mathrm{10}$ for an unstratified configuration. For δΩ, dashed (resp. solid) lines represent rotation slower (resp. faster) than the core. For ψ, dotted (resp. solid) lines are for clockwise (counter-clockwise) circulation. Boundary conditions are no-slip at r = η and stress-free at r = 1.

Fig. 8 Logarithm of the absolute value of the amplitude of the numerical azimuthal velocity as a function of the logarithm of the Ekman number at z = 0.5. The left panel is for η = 0.15, the middle one is for η = 0.25 and the right one for η = 0.5. Inside the tangent cylinder, measures were taken respectively at the points s = 0.1,0.2,0.25, the azimuthal velocity scales around E− 1 / 2. It is proportional to E-1 within the envelope s = 0.65. Note that the azimuthal velocity is mainly negative for small cores and positive for bigger cores like η = 0.5.

Fig. 9 Time evolution of the major component of the azimuthal velocity wl = 1 normalized by its steady solution ${\mathit{w}}_{\mathrm{st}}^{\mathit{l}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}$ as a function of non-dimensional time with E = 10-5,10-6,10-7 (the dot, the dashed, and the solid line, respectively) and Ro = 10-6 at the radius r = 0.88. On the left side are the solutions with rigid boundary conditions on both sides and on the right side the solutions with stress-free boundary conditions on the outer. In the first case, the steady state is reached in a time smaller than the dimensionless time of contraction, while in the second case this time is longer and scales as E-0.86.

Another remarkable property of the transient flow is its approximate self-similarity. Its spatial shape remains almost unchanged, while its amplitude grows as time passes. The associated differential rotation is parallel to the z-axis as shown in Fig. 10. Its amplitude grows according to the time profile displayed in Fig. 9. This time dependence can be approximated as $\begin{array}{ccc}\frac{\mathit{A}}{\mathit{E}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathit{\tau }\ln \mathrm{10}\mathit{/}{\mathit{\tau }}_{\mathrm{sf}}}\mathrm{)}& & \end{array}$(37)where A is a constant of order unity.

The foregoing result may be translated in the stellar case. It shows that a contracting, fully convective star, may reach a self-similar spin-up flow with cylindrical rotation. Neglecting viscous force (in fact Reynolds stresses), we may expect that substituting ρv to the foregoing incompressible velocity field, we can get an good representation of the actual flow in a compressible envelope. This is supported by the fact that the geostrophic balance is unchanged in this case. Of course, this conjecture has to be verified by the study of the compressible case.

4.3. Transient phase and stratification

During the contraction of a star on the PMS, a radiative envelope progressively takes the place of an initial convective envelope when the star is massive enough. This radiative envelope is stably stratified and without any extra-forcing from gravitational contraction it would relax to the steady baroclinic state that we mentioned before.

Fig. 10 Left: differential rotation associated with the transient phase of a spin-up induced by mass contraction for E = 10-6 and Ro = 10-6 at time τ = 0.02. Dashed (solid) lines represent rotation slower (faster) than the core. Right: profile of the corresponding azimuthal velocity in the equatorial plane as a function of the normalized radius for various times. The curves are scaled by the absolute value of the amplitude at the minimum.

One question is therefore whether the contraction induced spin-up is strong enough to overwhelm the foregoing baroclinic flows that are themselves transient flows. To get an idea of the result, we may use the amplitude of a steady baroclinic flows as an upper limit of the actual flows. From this perspective, we can use the results of Espinosa Lara & Rieutord (2007) and Espinosa Lara & Rieutord (2013) who showed that the baroclinic flow in a radiative envelope is characterized by a differential rotation that is typically 15% of the bulk rotation. Let us assume that the maximum amplitude of the baroclinic flow reads $\begin{array}{ccc}{\mathit{V}}_{\mathrm{b}}\mathrm{=}\mathit{k}{\mathit{V}}_{\mathrm{eq}}& & \end{array}$where k ≲ 0.2. During the phase of gravitational contraction the spin-up flow grows according to (37). At some time τ_c, the spin-up flow overwhelms the baroclinic flow. At such time we have $\begin{array}{ccc}\frac{\mathit{A}}{\mathit{E}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{c}}\ln \mathrm{10}\mathit{/}{\mathit{\tau }}_{\mathrm{sf}}}\mathrm{)}\mathrm{=}\frac{{\mathit{V}}_{\mathrm{b}}}{{\mathit{V}}_{\mathrm{s}}}& & \end{array}$with V_s ~ R/t_KH, t_KH being the Kelvin-Helmholtz time of the star. Hence, $\frac{{\mathit{\tau }}_{\mathrm{c}}}{{\mathit{\tau }}_{\mathrm{sf}}}\mathrm{~}\mathrm{-}\log \left(\mathrm{1}\mathrm{-}\frac{\mathit{E}{\mathit{V}}_{\mathit{b}}{\mathit{t}}_{\mathrm{KH}}}{\mathit{AR}}\right)\mathrm{·}$(38)Using numbers of a typical 3 M_⊙ star on the birthline, we find that $\begin{array}{ccc}\frac{\mathit{E}{\mathit{V}}_{\mathrm{b}}{\mathit{t}}_{\mathrm{KH}}}{\mathit{AR}}\mathrm{=}\frac{\mathit{k}}{\mathrm{2}\mathit{A}}\frac{\mathit{\nu G}{\mathit{M}}^{\mathrm{2}}}{{\mathit{R}}^{\mathrm{3}}\mathit{L}}\lesssim {\mathrm{10}}^{-4}\mathit{.}& & \end{array}$This small ratio indicates that τ_c ≪ τ_sf so that the contraction-induced flow overwhelms the baroclinic flow during the linear growth of the transient flow (37). Hence, from (36) we obtain $\begin{array}{ccc}{\mathit{\tau }}_{\mathrm{c}}\mathrm{~}\frac{{\mathit{E}}^{\mathrm{0.14}}\mathit{Ro}}{\mathit{A}}\frac{{\mathit{V}}_{\mathrm{b}}{\mathit{t}}_{\mathrm{KH}}}{\mathit{R}}\mathrm{·}& & \end{array}$With the definition of Ro and V_b it turns out that $\begin{array}{ccc}{\mathit{\tau }}_{\mathrm{c}}\mathrm{~}\frac{\mathit{k}}{\mathrm{2}\mathit{A}}{\mathit{E}}^{\mathrm{0.14}}\mathit{.}& & \end{array}$Because of the very small value of the Ekman number, τ_c is clearly less than unity showing that the spin-up flow will in the end take over the baroclinic flow, likely much before this latter flow can be established.

We verified this conclusion with a numerical simulation integrating the spin-up flow from a pre-existing baroclinic flow driven by a fixed stable stratification. In Fig. 11, we show the time evolution of the azimuthal velocity in the equatorial plane of the star. This describes the transient phase from a steady baroclinic flow to a growing spin-up flow. The transition between the two flows happens between τ = 10^-3 and τ = 10^-2. This is less than τ_c ~ 0.2 (at E = 10^-5), namely less than our first evaluation obtained by a comparison of the amplitudes of the flows (see Eq. (38)). The parameters have been chosen such that BE ≪ 1, as expected in real situations. The time τ_c therefore appears as a good indicator of the time needed by spin-up to overwhelm baroclinic flows. Let us finally note that, once the spin-up flow is settled, the flow remains approximately self-similar for at least 80% of the Kelvin-Helmholtz time.

Fig. 11 Radial profiles taken at various times of the azimuthal velocity in the equatorial plane of the star for E = 10-5, B = 102, ϵ = 10-7, λ = 10-4, and Ro = 10-5. Profiles are scaled as in Fig. 10. The dimensionless time τ is specified for each curve. Dashed lines are profiles for baroclinicly dominated dynamics while solid lines are for spin-up dominated dynamics. The dashed bold curve shows the initial profile and the solid bold curve shows the last profile.

5. Discussion and conclusions

The gravitational contraction that occurs before or after the main sequence strongly influences the rotation rate of the stars and their internal dynamics. In the foregoing study, we have investigated the consequences of the combination of rotation and gravitational contraction with a very simplified model to decipher this complicated dynamics and be prepared for the construction of more elaborated models of rotating stars like the ESTER (for the French: Evolution STellaire En Rotation, in English: stellar evolution in rotation) models of Espinosa Lara & Rieutord (2013).

Thus, to the compressible gas of the star, we substituted an incompressible fluid that may also be stably stratified. Schematically, the bell-shaped profile of the star density is replaced by a step function delineating a central core that absorbs matter from the envelope as the star contracts. The size of the core is small and arbitrary. The envelope is either neutrally or stably stratified.

During PMS contraction, the stellar envelope of an intermediate-mass star usually passes from a convective to a radiative state. But the radiative state is stably stratified. The combined effect of rotation and stable stratification drives a baroclinic flow that may contest a pre-existing spin-up flow built up during a previous convective phase of the envelope. The question is which of these flows will govern the dynamics of the contracting star and finally determine the initial conditions of the dynamics on the main sequence.

Using our simplified model we compared the strength of these two flows and found that although contraction is slow, the induced spin-up controls the large-scale flows of the outer envelope, namely its differential rotation and meridional circulation. Moreover, our model underlines the role of the outer boundary conditions and shows that with realistic stress-free conditions we should expect an unsteady flow. In addition, it shows that this transient flow keeps a self-similar shape during its growth (if we omit boundary layers).

When the star reaches the main sequence, the contraction turns off and the flows in the envelope relax towards the steady baroclinic flow on the Eddington-Sweet timescale. As far as intermediate-mass stars are concerned, because of their fast rotation, the Eddington-Sweet timescale is close to the Kelvin-Helmholtz timescale and the transition to the quasi-steady state of the main sequence is quite short (for instance for 7 M_⊙ star, the Eddington-Sweet timescale is 2.8 Myr for a rotation near break-up, to be compared to the 46 Myr of the lifetime of such stars). On the other hand, if for some reason (such as the combination of magnetic fields and mass loss) the star loses much angular momentum so that 2Ω ≪ N at the beginning of the main sequence, then the dynamic state of an outer envelope will be controlled by slowly decaying baroclinic modes excited by the contraction phase. The slow decay may occupy a significant fraction of the main sequence and affect the mixing processes.

Back to fluid dynamics, the simple model that we used shows other details about the dynamics of this system, like for instance the shear layer (the Stewartson layer) that circumvent the core on its tangent cylinder. This feature is clearly an artefact of the model for stars with no convective cores like PMS stars, but it is certainly an important feature for stars leaving the main sequence where the core contracts and the envelope expands. At the core-envelope interface the build-up of a density jump due to nuclear evolution, combined with rotation, triggers a Stewartson layer on the tangent cylinder of the discontinuity. This layer may indeed explain the efficient transport of angular momentum between the core and envelope of giant or subgiant stars that is needed to explain the rather mild radial differential rotation observed in these stars (Deheuvels et al. 2012; Mosser et al. 2012; Eggenberger et al. 2012). Indeed, our model shows that there is a tight coupling between the inner part of the tangent cylinder and the core itself. This coupling is essentially a consequence of the Taylor-Proudman theorem that imposes no velocity gradient along the rotation axis (columnar flow). Therefore, the transport of angular momentum between the core and the envelope is much enhanced by the Stewartson layer. Since such a layer has a surface that is R_star/R_core larger than the surface of the core, we expect that the flux of angular momentum between the core and the envelope will be enhanced by a similar factor (viscosity and velocity gradient being assumed similar) with respect to a 1D shellular profile. Noting that the moment of inertia of the core and of the matter inside its tangent cylinder are not much different, we expect that the Stewartson layer plays a crucial role in the angular momentum exchange between the core and envelope and might well be the key feature that reconcile models and observations. It is clear that present 1D models do not take this fluid dynamics feature into account and that the final answer will be given by 2D models incorporating this flow. A dedicated study is clearly needed to give a quantitative estimate of this effect and to offer a new comparison with observations.

Hence, more than the numbers and the applicability to a given object, the foregoing Boussinesq model underlines the main features of the dynamics of a contracting and rotating envelope. It stresses the key role of outer boundary conditions and the various flows that might govern a contracting phase depending on the strength of the stratification. The side effect of the core in this model underlines the role of a Stewartson layer that may appear either after a rapid change in density or in viscosity. The model also stresses the fact that no steady state can be expected as for the interior flows, but that these flows may converge towards a universal flow when gravitational contraction ceases. The next step of these investigations focusing specifically on stars will be developed with full physics using the ESTER code of Espinosa Lara & Rieutord (2013).

Acknowledgments

We would like to thank the referee for his constructive remarks on the first version of the manuscript. We also acknowledge the support of the French Agence Nationale de la Recherche (ANR), under grant ESTER (ANR-09-BLAN-0140). The numerical calculations have been carried out on the CalMip machine of the “Centre Interuniversitaire de Calcul de Toulouse” (CICT), which is gratefully acknowledged.

References

Appendix A: Numerical method

To solve the set of linear Eqs. (12), we discretize the equations using a spectral method. We project the velocity field onto the harmonics spherical base$\begin{array}{ccc}{u}\mathrm{=}\sum _{\mathit{l}\mathrm{=}\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\sum _{\mathit{m}\mathrm{=}\mathrm{-}\mathit{l}}^{\mathrm{+}\mathit{l}}{\mathit{u}}_{\mathit{m}}^{\mathit{l}}{{R}}_{\mathit{l}}^{\mathit{m}}\mathrm{+}{\mathit{v}}_{\mathit{m}}^{\mathit{l}}{{S}}_{\mathit{l}}^{\mathit{m}}\mathrm{+}{\mathit{w}}_{\mathit{m}}^{\mathit{l}}{{T}}_{\mathit{l}}^{\mathit{m}}& & \end{array}$(A.1)where$\begin{array}{ccc}{{R}}_{\mathit{l}}^{\mathit{m}}\mathrm{=}{\mathit{Y}}_{\mathit{l}}^{\mathit{m}}{{e}}_{{r}}\mathit{,} {{S}}_{\mathit{l}}^{\mathit{m}}\mathrm{=}{\nabla }{\mathit{Y}}_{\mathit{l}}^{\mathit{m}}\mathit{,} {{T}}_{\mathit{l}}^{\mathit{m}}\mathrm{=}{\nabla }\mathrm{\times }{{R}}_{\mathit{l}}^{\mathit{m}}& & \end{array}$(A.2)${\mathit{Y}}_{\mathit{l}}^{\mathit{m}}$ are the normalized spherical harmonics, e_r is the radial unity vector, and ∇ is defined on the unity sphere. We write the temperature perturbation onto the spherical harmonics base too:$\begin{array}{ccc}{\mathit{\theta }}_{\mathit{T}}\mathrm{=}\sum _{\mathit{l}\mathrm{=}\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\sum _{\mathit{m}\mathrm{=}\mathrm{-}\mathit{l}}^{\mathrm{+}\mathit{l}}{\mathit{t}}_{\mathit{m}}^{\mathit{l}}{\mathit{Y}}_{\mathit{l}}^{\mathit{m}}\mathit{.}& & \end{array}$(A.3)We finally add the boundary conditions on this field: $\begin{array}{ccc}\{\begin{array}{c}\\ \\ \end{array}& & \end{array}$(A.4)We discretize the radial direction for r ∈ [ η;1 ] onto the Gauss-Lobatto grid associated with the Chebyshev polynomials. Thereby, equations are solved in two dimensions (r,θ). The system is axisymmetric, which implies m = 0.

The equation of continuity reads$\begin{array}{ccc}{\mathit{v}}_{\mathit{m}}^{\mathit{l}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{\Lambda }}\frac{\mathrm{1}}{\mathit{r}}\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{(}{\mathit{r}}^{\mathrm{2}}{{\mathit{u}}_{\mathit{m}}^{\mathit{l}}}^{\mathrm{)}}& & \end{array}$(A.5)where Λ = l(l + 1).

The energy equation reads $\begin{array}{ccc}\mathit{B}\begin{array}{c}\mathrm{˜}\\ {\mathit{E}}_{\mathit{T}}\end{array}\left({\mathit{r}}^{\mathrm{2}}\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{r}}^{\mathrm{2}}}{\mathit{t}}_{\mathit{m}}^{\mathit{l}}\mathrm{+}\mathrm{2}\mathit{r}\frac{\mathit{\partial }}{\mathit{\partial r}}{\mathit{t}}_{\mathit{m}}^{\mathit{l}}\mathrm{-}\mathrm{\Lambda }{\mathit{t}}_{\mathit{m}}^{\mathit{l}}\right)\mathrm{=}{\mathit{n}}_{\mathit{T}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathit{r}}^{\mathrm{2}}{\mathit{u}}_{\mathit{m}}^{\mathit{l}}\mathit{.}& & \end{array}$(A.6)The equation of motion is projected onto two directions because equations on ${{R}}_{\mathit{l}}^{\mathit{m}}$ and on ${{S}}_{\mathit{l}}^{\mathit{m}}$ are redundant.

On ${{R}}_{\mathit{l}}^{\mathit{m}}$, it reads $\begin{array}{ccc}{\mathit{A}}_{\mathit{l}\mathrm{-}\mathrm{1}}^{\mathit{l}}{\mathit{r}}^{\mathit{l}\mathrm{-}\mathrm{1}}\frac{\mathit{\partial }}{\mathit{\partial r}}\left(\frac{{\mathit{u}}_{\mathit{m}}^{\mathit{l}\mathrm{-}\mathrm{1}}}{{\mathit{r}}^{\mathit{l}\mathrm{-}\mathrm{2}}}\right)\mathrm{+}{\mathit{A}}_{\mathit{l}\mathrm{+}\mathrm{1}}^{\mathit{l}}{\mathit{r}}^{\mathrm{-}\mathit{l}\mathrm{-}\mathrm{2}}\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{(}{\mathit{r}}^{\mathit{l}\mathrm{+}\mathrm{3}}{{\mathit{u}}_{\mathit{m}}^{\mathit{l}\mathrm{+}\mathrm{1}}}^{\mathrm{)}}& & \\ \mathrm{+}\mathit{E}{\mathrm{\Delta }}_{\mathit{l}}{\mathit{w}}_{\mathit{m}}^{\mathit{l}}\mathrm{=}\mathrm{-}\sqrt{\frac{\mathrm{4}\mathit{\pi }}{\mathrm{3}}}{\mathit{\delta }}_{\mathit{l}\mathrm{1}}\mathrm{(}\frac{{\mathit{\eta }}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{2}}}\mathrm{-}\mathit{r}\mathit{\omega ̇}\mathrm{)}& & \end{array}$(A.7)δ_ij is the Kronecker symbol.

On ${{T}}_{\mathit{l}}^{\mathit{m}}$, it reads

$\begin{array}{ccc}\mathrm{-}{\mathit{B}}_{\mathit{l}\mathrm{-}\mathrm{1}}^{\mathit{l}}{\mathit{r}}^{\mathit{l}\mathrm{-}\mathrm{1}}\frac{\mathit{\partial }}{\mathit{\partial r}}\left(\frac{{\mathit{w}}_{\mathit{m}}^{\mathit{l}}}{{\mathit{r}}^{\mathit{l}\mathrm{-}\mathrm{1}}}\right)\mathrm{-}{\mathit{B}}_{\mathit{l}\mathrm{+}\mathrm{1}}^{\mathit{l}}{\mathit{r}}^{\mathrm{-}\mathit{l}\mathrm{-}\mathrm{2}}\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{(}{\mathit{r}}^{\mathit{l}\mathrm{+}\mathrm{2}}{{\mathit{w}}_{\mathit{m}}^{\mathit{l}\mathrm{+}\mathrm{1}}}^{\mathrm{)}}& & \\ \mathrm{-}\mathit{l}\mathrm{(}\mathit{l}\mathrm{+}\mathrm{1}\mathrm{)}{\mathit{t}}_{\mathit{m}}^{\mathit{l}}\mathrm{+}\mathit{E}{\mathrm{\Delta }}_{\mathit{l}}{\mathrm{\Delta }}_{\mathit{l}}\mathrm{\left(}\mathit{r}{\mathit{u}}_{\mathit{m}}^{\mathit{l}}\mathrm{\right)}\mathrm{=}\mathrm{-}\sqrt{\frac{\mathrm{16}\mathit{\pi }}{\mathrm{5}}}{\mathit{n}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathit{\delta }}_{\mathit{l}\mathrm{2}}\mathit{.}& & \end{array}$(A.8)We note ${\mathit{A}}_{\mathit{l}\mathrm{-}\mathrm{1}}^{\mathit{l}}$,${\mathit{A}}_{\mathit{l}\mathrm{+}\mathrm{1}}^{\mathit{l}}$,${\mathit{B}}_{\mathit{l}\mathrm{-}\mathrm{1}}^{\mathit{l}}$, and ${\mathit{B}}_{\mathit{l}\mathrm{+}\mathrm{1}}^{\mathit{l}}$ the coupling coefficients: $\begin{array}{ccc}\{\begin{array}{c}\\ \\ \end{array}& & \end{array}$(A.9)\newpage\noindentNoting that the forcing implies equatorially symmetric solutions (the resulting differential rotation is equatorially symmetric), the radial functions w^l are non-zero only for odd l while u^l and t^l only for even l. The series is therefore w¹,u²,t²,w³,...

All Figures

	Fig. 1 Typical and scaled Brunt-Väisälä frequency profile as a function of the normalized radius with η = 0.15.
In the text

Fig. 2 Comparison between numerical (solid line) and analytical (star line) solutions of the geostrophic flow at the equator ${\mathit{u}}_{\mathit{\phi }}\mathrm{(}\mathit{r,\theta }\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}\mathrm{)}$ for E = 10-7, η = 0.15, $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\rho }\end{array}\mathrm{=}\mathrm{10}$ without stratification when the envelope is rigidly rotating at the same rate as the inner one.

In the text

	Fig. 3 Comparison between numerical (solid line) and analytical (stars) solutions of the geostrophic flow at the equator ${\mathit{u}}_{\mathit{\phi }}\mathrm{(}\mathit{r,\theta }\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}\mathrm{)}$ for E = 10-7, η = 0.15, $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\rho }\end{array}\mathrm{=}\mathrm{10}$, B = 1, and λ = 10-4 without mass contraction.
In the text

Fig. 4 Evolution of the azimuthal velocity ${\mathit{u}}_{\mathit{\phi }}\mathrm{(}\mathit{r,\theta }\mathrm{=}\frac{\mathit{\pi }}{\mathrm{2}}\mathrm{)}$ for η = 0.15, $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\rho }\end{array}\mathrm{=}\mathrm{10}$, E = 10-7, λ = 10-4, and B = 0,104,105,107 (downward). The solid line is the numerical solution, the dashed line is the full analytical solution, and the two dot lines are the solutions of each flow. The transition is at B = 104, the differential rotation is then governed by baroclinicity as shown in Fig. 5. Boundary conditions are no-slip on both sides r = η and r = 1.

In the text

Fig. 5 Evolution of the differential rotation δΩ, the meridional circulation ψ for η = 0.15, $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\rho }\end{array}\mathrm{=}\mathrm{10}$, E = 10-7, λ = 10-4, and B = 0,104,105,107 (downwards). In the first column, difference in differential rotation is shown with contours: solid contours are faster than the core and dashed contours are slower than the core. When the spin-up flow dominates the differential rotation we obtain a fast equator while baroclinicity induces a fast pole. In the second column, meridional circulation is described with dotted lines for clockwise circulation (solid lines for counter-clockwise circulation). Boundary conditions are no-slip on both sides r = η and r = 1.

In the text

Fig. 6 Azimuthal velocity uφ as a function of the normalized radius for λ = 10-4, η = 0.15, and $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\rho }\end{array}\mathrm{=}\mathrm{10}$. The solid line is the numerical solution, dotted lines are the analytical solutions of the spin-up flow and the baroclinic flow, the dashed line is the sum of both. Downwards: E = 10-5,10-6,10-7 and the transition value on the differential rotation is B = 103,3 × 103,104, respectively. Beyond this threshold, the baroclinicity governs the differential rotation. Boundary conditions are no-slip on both sides r = η and r = 1. Note that as the Ekman number decreases, the discrepancy between the numerical and the analytical solutions decreases as well.

In the text

Fig. 7 Differential rotation and meridional circulation for E = 10-7, η = 0.15, and $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\rho }\end{array}\mathrm{=}\mathrm{10}$ for an unstratified configuration. For δΩ, dashed (resp. solid) lines represent rotation slower (resp. faster) than the core. For ψ, dotted (resp. solid) lines are for clockwise (counter-clockwise) circulation. Boundary conditions are no-slip at r = η and stress-free at r = 1.

In the text

Fig. 8 Logarithm of the absolute value of the amplitude of the numerical azimuthal velocity as a function of the logarithm of the Ekman number at z = 0.5. The left panel is for η = 0.15, the middle one is for η = 0.25 and the right one for η = 0.5. Inside the tangent cylinder, measures were taken respectively at the points s = 0.1,0.2,0.25, the azimuthal velocity scales around E− 1 / 2. It is proportional to E-1 within the envelope s = 0.65. Note that the azimuthal velocity is mainly negative for small cores and positive for bigger cores like η = 0.5.

In the text

Fig. 9 Time evolution of the major component of the azimuthal velocity wl = 1 normalized by its steady solution ${\mathit{w}}_{\mathrm{st}}^{\mathit{l}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}$ as a function of non-dimensional time with E = 10-5,10-6,10-7 (the dot, the dashed, and the solid line, respectively) and Ro = 10-6 at the radius r = 0.88. On the left side are the solutions with rigid boundary conditions on both sides and on the right side the solutions with stress-free boundary conditions on the outer. In the first case, the steady state is reached in a time smaller than the dimensionless time of contraction, while in the second case this time is longer and scales as E-0.86.

In the text

Fig. 10 Left: differential rotation associated with the transient phase of a spin-up induced by mass contraction for E = 10-6 and Ro = 10-6 at time τ = 0.02. Dashed (solid) lines represent rotation slower (faster) than the core. Right: profile of the corresponding azimuthal velocity in the equatorial plane as a function of the normalized radius for various times. The curves are scaled by the absolute value of the amplitude at the minimum.

In the text

Fig. 11 Radial profiles taken at various times of the azimuthal velocity in the equatorial plane of the star for E = 10-5, B = 102, ϵ = 10-7, λ = 10-4, and Ro = 10-5. Profiles are scaled as in Fig. 10. The dimensionless time τ is specified for each curve. Dashed lines are profiles for baroclinicly dominated dynamics while solid lines are for spin-up dominated dynamics. The dashed bold curve shows the initial profile and the solid bold curve shows the last profile.

In the text