Issue 
A&A
Volume 572, December 2014



Article Number  A15  
Number of page(s)  12  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/201423378  
Published online  20 November 2014 
Dynamics of the envelope of a rapidly rotating star or giant planet in gravitational contraction
^{1}
Université de Toulouse, UPSOMP, IRAP, 31400
Toulouse, France
email: hypolite.delphine@gmail.com
^{2}
CNRS, IRAP, 14
avenue Édouard Belin, 31400
Toulouse,
France
Received:
7
January
2014
Accepted:
20
September
2014
Aims. We wish to understand the processes that control the fluid flows of a gravitationally contracting and rotating star or giant planet.
Methods. We consider a spherical shell containing an incompressible fluid that is slowly absorbed by the core so as to mimic gravitational contraction. We also consider the effects of a stable stratification that may modify the dynamics of a premainsequence star of intermediate mass.
Results. This simple model reveals the importance of both the Stewartson layer attached to the core and the boundary conditions met by the fluid at the surface of the object. In the case of a premainsequence star of intermediate mass where the envelope is stably stratified, shortly after the birth line, the spinup flow driven by contraction overwhelms the baroclinic flow that would take place otherwise. This model also shows that for a contracting envelope, a selfsimilar flow of growing amplitude controls the dynamics. It suggests that initial conditions on the birth line are most probably forgotten. Finally, the model shows that the shear (Stewartson) layer that lies on the tangent cylinder of the core is likely a key feature of the dynamics that is missing in 1D models. This layer can explain the core and envelope rotational coupling that is required to explain the slow rotation of cores in giant and subgiant stars.
Key words: stars: rotation / hydrodynamics
© 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 onedimensional while fluid flows are generally multidimensional. 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 firstorder 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 supergiants in lowmetallicity 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 socalled rotational mixing. One of the most famous is the distribution of Large Magellanic Cloud Bstars in a diagram plotting the nitrogen abundance versus the rotational velocity (the socalled 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 nonrotating, 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 premain 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 setup for boundary conditions in complete (compressible) twodimensional 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 spindown 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 noslip and use the ensuing analytical solutions to enlight the dynamics (Sect. 3). In Sect. 4, we discuss the presumably more realistic case with stressfree 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 selfgravitating incompressible viscous fluid enclosed between two spherical shells. These shells are assumed not to be distorted by rotation. The inner shell may represents a coreenvelope 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 where we introduced the initial density of the core ρ_{0} and the “suction time”, which is the time the core needs to increase its density by a factor 2, namely (1)where . This time is of same order of magnitude as the KelvinHelmholtz 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 spinup, of the star. The problem of spinup 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 spindown 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 setup. This work however establishes that the spindown (or spinup) timescale of a stably stratified fluid is the classical EddingtonSweet timescale, namely the product of the heat diffusion timescale (also known as KelvinHelmholtz time scale) and the factor , where is the BruntVäisälä frequency in the radiative region and Ω is the rotation rate. In slowly rotating stars, the EddingtonSweet timescale is much longer than the KelvinHelmholtz, but in fast rotators that are considered here, these timescales are similar. Hence, the spinup 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 quasisteady 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 noslip 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 (spinup and baroclinicity) necessarily arises in the premainsequence 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 noslip outer boundary condition, we compare the amplitudes of the baroclinic flow and the contractiondriven spinup 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 noslip outer boundary condition and concentrate on the more realistic stressfree 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 spinup 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 spinup 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 quasiselfsimilar 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 contractioninduced spinup 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: (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 (3)We define Ω as the angular velocity of the core. Since the field v is axisymmetric, 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 (4)Note that the RHSterm 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 KelvinHelmholtz timescale , 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_{s}u. The system (5) now reads (5)where p is the reduced pressure. The following dimensionless numbers appear: (6)Here we introduced the nondimensional 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 quasisteady evolution of the system as long as the nondimensional time τ verifies: (7)where is the scaled suction time ( 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 (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 (9)Besides, for a sphere of mass M_{core} and radius ηR: (10)With the previous scaling, we obtain (11)The evaluation of the remaining integral needs the expression of the azimuthal flow u_{φ} at the coreenvelope 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 The new velocity field u^{′} verifies (12)with the boundary conditions and These latter conditions describe stressfree 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 (13)In the azimuthal direction, we look for a flow such as as dictated by the TaylorProudman theorem (e.g. Greenspan 1969). Such a flow does not verify the noslip boundary conditions at the interface r = η. It needs boundary layer corrections so as to satisfy . 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: (14)where and n = −e_{r}. We keep only the decreasing solution. The corrections thus read (15)As in Rieutord (2006a), mass conservation gives the relation between the Ekman pumping and the geostrophic flow U. It yields(16)where we keep only terms. Integrating with respect to θ, we find (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 (18)Integral in (11) can now be evaluated: (19)We can now derive the acceleration of the angular velocity of the core. Considering quasisteady solutions that arise when Ro ≪ τ ≪ τ_{s}, (11) leads to (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 , 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 contractioninduced 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^{2}L_{⊙} and its radius is R_{∗} ~ 10^{10} m. The young star contracts on the KelvinHelmholtz time upon the PMS (Henyey) track, namely: according to Maeder (2009). This leads to t_{KH} ~ 2.6 × 10^{5} 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 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^{4} m^{2} s^{1} and thus With the radiative viscosity (e.g. Espinosa Lara & Rieutord 2013), ν ~ 10^{2} m^{2} s^{1}, we get For both values the condition is satisfied.
During the contraction, the quasisteady state within the cylinder is reached after a spinup time, namely after . Using the previous numbers, we find that this state occurs after ~10^{3} or 10^{4} 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 Jupiterlike 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 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 BruntVä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 (21)We use the same scales and notations as Rieutord (2006a). The temperature perturbation is scaled as δT = ϵT_{∗}θ_{T} where is the ratio of centrifugal acceleration to surface gravity. Recall that the centrifugal acceleration is driving the baroclinic flow. The BruntVäisälä profile is scaled with 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 (22)Finally, the dimensionless number (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 BruntVäisälä frequency is shown in Fig. 1. We use the polynomial expression (24)to represent this function. The coefficients α, β and γ result from the polynomial fit.
Fig. 1 Typical and scaled BruntVä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 noslip 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 noslip conditions on the outer boundary r = 1, we may easily derive the expression of the geostrophic flow^{1} out of the tangent cylinder of the core. When no stratification is present, the azimuthal velocity reads (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 for E = 10^{7}, η = 0.15, 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 noslip boundary conditions. It verifies (26)As shown in Rieutord (2006a), neglecting temperature perturbations and viscosity, the φcomponent of the vorticity Eq. (26) leads to (27)where F(s) is a geostrophic solution determined by the boundary conditions. In order to get the expression of F(s), we write (28)and look for the boundary layer corrections at r = 1. These are where . Using (27) with the foregoing expression yields (29)We note that g(1) = 0. Mass conservation implies (30)so that (31)Since at r = 1, we finally find (32)However, the radial component of the vorticity equation (or the angular momentum equation) in the interior leads to Here so that consistency of the solution with (32) requires that 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: Figure 3 shows the comparison of the analytical solution − sBg(r) with the numerical solution at the equator . 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 for E = 10^{7}, η = 0.15, , B = 1, and λ = 10^{4} without mass contraction. 
Fig. 4 Evolution of the azimuthal velocity for η = 0.15, , E = 10^{7}, λ = 10^{4}, and B = 0,10^{4},10^{5},10^{7} (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 = 10^{4}, the differential rotation is then governed by baroclinicity as shown in Fig. 5. Boundary conditions are noslip 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 (33)From the foregoing equation, we see that the differential rotation is governed by baroclinicity when (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 contractioninduced spinup 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 spinup strengthens, the foregoing inequalities predict an intermediate regime (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 spinup.
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 stressfree boundary conditions
The outer layers of the envelope are not rigidly attached to the core. Therefore the use of outer stressfree 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 spinup 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 noslip conditions).
Fig. 5 Evolution of the differential rotation δΩ, the meridional circulation ψ for η = 0.15, , E = 10^{7}, λ = 10^{4}, and B = 0,10^{4},10^{5},10^{7} (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 spinup 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 counterclockwise circulation). Boundary conditions are noslip 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 noslip, 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 stressfree. 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 .
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 noslip and stressfree outer boundary conditions. With noslip conditions, we see that the steady spinup flow is quickly established and justifies the use of a steady solution. On the contrary, the use of outer stressfree 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 (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 stressfree 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 . The solid line is the numerical solution, dotted lines are the analytical solutions of the spinup 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 = 10^{3},3 × 10^{3},10^{4}, respectively. Beyond this threshold, the baroclinicity governs the differential rotation. Boundary conditions are noslip 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 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 (counterclockwise) circulation. Boundary conditions are noslip at r = η and stressfree 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 w^{l = 1} normalized by its steady solution as a function of nondimensional 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 stressfree 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 selfsimilarity. Its spatial shape remains almost unchanged, while its amplitude grows as time passes. The associated differential rotation is parallel to the zaxis as shown in Fig. 10. Its amplitude grows according to the time profile displayed in Fig. 9. This time dependence can be approximated as (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 selfsimilar spinup 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 extraforcing 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 spinup 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 spinup 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 where k ≲ 0.2. During the phase of gravitational contraction the spinup flow grows according to (37). At some time τ_{c}, the spinup flow overwhelms the baroclinic flow. At such time we have with V_{s} ~ R/t_{KH}, t_{KH} being the KelvinHelmholtz time of the star. Hence, (38)Using numbers of a typical 3 M_{⊙} star on the birthline, we find that This small ratio indicates that τ_{c} ≪ τ_{sf} so that the contractioninduced flow overwhelms the baroclinic flow during the linear growth of the transient flow (37). Hence, from (36) we obtain With the definition of Ro and V_{b} it turns out that Because of the very small value of the Ekman number, τ_{c} is clearly less than unity showing that the spinup 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 spinup flow from a preexisting 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 spinup 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 spinup to overwhelm baroclinic flows. Let us finally note that, once the spinup flow is settled, the flow remains approximately selfsimilar for at least 80% of the KelvinHelmholtz 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 = 10^{2}, ϵ = 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 spinup 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 bellshaped 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 intermediatemass 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 preexisting spinup 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 spinup controls the largescale 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 stressfree conditions we should expect an unsteady flow. In addition, it shows that this transient flow keeps a selfsimilar 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 EddingtonSweet timescale. As far as intermediatemass stars are concerned, because of their fast rotation, the EddingtonSweet timescale is close to the KelvinHelmholtz timescale and the transition to the quasisteady state of the main sequence is quite short (for instance for 7 M_{⊙} star, the EddingtonSweet timescale is 2.8 Myr for a rotation near breakup, 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 coreenvelope interface the buildup 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 TaylorProudman 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 (ANR09BLAN0140). 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
 Boehm, T., & Catala, C. 1995, A&A, 301, 155 [NASA ADS] [Google Scholar]
 Brott, I., Evans, C. J., Hunter, I., et al. 2011, A&A, 530, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Charbonnel, C., & Talon, S. 1999, A&A, 351, 635 [NASA ADS] [Google Scholar]
 Deheuvels, S., García, R. A., Chaplin, W. J., et al. 2012, ApJ, 756, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Duck, P., & Foster, M. 2001, Annu. Rev. Fluid Mech., 33, 231 [NASA ADS] [CrossRef] [Google Scholar]
 Eggenberger, P., Montalbán, J., & Miglio, A. 2012, A&A, 544, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Espinosa Lara, F., & Rieutord, M. 2007, A&A, 470, 1013 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Espinosa Lara, F., & Rieutord, M. 2013, A&A, 552, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fortney, J. J., & Nettelmann, N. 2010, Space Sci. Rev., 152, 423 [NASA ADS] [CrossRef] [Google Scholar]
 Friedlander, S. 1976, J. Fluid Mech., 76, 209 [NASA ADS] [CrossRef] [Google Scholar]
 Greenspan, H. P. 1969, The Theory of Rotating Fluids (Cambridge University Press) [Google Scholar]
 Maeder, A. 2009, Physics, Formation and Evolution of Rotating stars (Springer) [Google Scholar]
 Maeder, A., & Meynet, G. 2001, A&A, 373, 555 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Meynet, G., & Maeder, A. 2005, A&A, 429, 581 [Google Scholar]
 Mosser, B., Goupil, M. J., Belkacem, K., et al. 2012, A&A, 548, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ogilvie, G. I., & Lin, D. N. C. 2004, ApJ, 610, 477 [NASA ADS] [CrossRef] [Google Scholar]
 Pinsonneault, M. 1997, Annu. Rev. Astron. Astrophys., 35, 557 [Google Scholar]
 Rieutord, M. 2006a, A&A, 451, 1025 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rieutord, M. 2006b, in Stellar fluid dynamics and numerical simulations: From the sun to neutron stars, eds. M. Rieutord, & B. Dubrulle (EAS Pub.) 21, 275 [Google Scholar]
 Rieutord, M., & Beth, A. 2014, A&A, 570, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stewartson, K. 1966, J. Fluid Mech., 26, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Strittmatter, P. A. 1969, Annu. Rev. Astron. Astrophys., 7, 665 [NASA ADS] [CrossRef] [Google Scholar]
 Wu, Y. 2005, ApJ, 635, 688 [NASA ADS] [CrossRef] [Google Scholar]
 Zahn, J.P. 1992, A&A, 265, 115 [NASA ADS] [Google Scholar]
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(A.1)where(A.2) 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:(A.3)We finally add the boundary conditions on this field: (A.4)We discretize the radial direction for r ∈ [ η;1 ] onto the GaussLobatto 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(A.5)where Λ = l(l + 1).
The energy equation reads (A.6)The equation of motion is projected onto two directions because equations on and on are redundant.
On , it reads (A.7)δ_{ij} is the Kronecker symbol.
On , it reads
(A.8)We note ,,, and the coupling coefficients: (A.9)\newpage\noindentNoting that the forcing implies equatorially symmetric solutions (the resulting differential rotation is equatorially symmetric), the radial functions w^{l} are nonzero only for odd l while u^{l} and t^{l} only for even l. The series is therefore w^{1},u^{2},t^{2},w^{3},...
All Figures
Fig. 1 Typical and scaled BruntVä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 for E = 10^{7}, η = 0.15, 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 for E = 10^{7}, η = 0.15, , B = 1, and λ = 10^{4} without mass contraction. 

In the text 
Fig. 4 Evolution of the azimuthal velocity for η = 0.15, , E = 10^{7}, λ = 10^{4}, and B = 0,10^{4},10^{5},10^{7} (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 = 10^{4}, the differential rotation is then governed by baroclinicity as shown in Fig. 5. Boundary conditions are noslip on both sides r = η and r = 1. 

In the text 
Fig. 5 Evolution of the differential rotation δΩ, the meridional circulation ψ for η = 0.15, , E = 10^{7}, λ = 10^{4}, and B = 0,10^{4},10^{5},10^{7} (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 spinup 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 counterclockwise circulation). Boundary conditions are noslip 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 . The solid line is the numerical solution, dotted lines are the analytical solutions of the spinup 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 = 10^{3},3 × 10^{3},10^{4}, respectively. Beyond this threshold, the baroclinicity governs the differential rotation. Boundary conditions are noslip 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 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 (counterclockwise) circulation. Boundary conditions are noslip at r = η and stressfree 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 w^{l = 1} normalized by its steady solution as a function of nondimensional 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 stressfree 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 spinup 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 = 10^{2}, ϵ = 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 spinup dominated dynamics. The dashed bold curve shows the initial profile and the solid bold curve shows the last profile. 

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