Issue 
A&A
Volume 648, April 2021



Article Number  A109  
Number of page(s)  22  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202039248  
Published online  05 May 2021 
Axisymmetric investigation of differential rotation in contracting stellar radiative zones
Institut de Recherche en Astrophysique et Planétologie (IRAP), Université de Toulouse, 14 Avenue Edouard Belin, 31400 Toulouse, France
email: bgouhier@irap.omp.eu; flignieres@irap.omp.eu; ljouve@irap.omp.eu
Received:
24
August
2020
Accepted:
11
November
2020
Context. Stars experience rapid contraction or expansion at different phases of their evolution. Modelling the transport of angular momentum and the transport of chemical elements occurring during these phases remains an unsolved problem.
Aims. We study a stellar radiative zone undergoing radial contraction and investigate the induced differential rotation and meridional circulation.
Methods. We consider a rotating spherical layer crossed by an imposed radial velocity field that mimics the contraction, and numerically solve the axisymmetric hydrodynamical equations in both the Boussinesq and anelastic approximations. An extensive parametric study is conducted to cover regimes of contraction, rotation, stable stratification, and density stratification that are relevant for stars.
Results. The differential rotation and the meridional circulation result from a competition between the contractiondriven inward transport of angular momentum and an outward transport dominated by either viscosity or an Eddington–Sweettype circulation, depending on the value of the P_{r}(N_{0}/Ω_{0})^{2} parameter, where P_{r} is the Prandtl number, N_{0} the Brunt–Väisäilä frequency, and Ω_{0} the rotation rate. Taking the density stratification into account is important to study more realistic radial contraction fields, and also because the resulting flow is less affected by unwanted effects of the boundary conditions. In these different regimes and for a weak differential rotation we derive scaling laws that relate the amplitude of the differential rotation to the contraction timescale.
Key words: hydrodynamics / methods: numerical / stars: rotation
© B. Gouhier et al. 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
The transport of angular momentum (AM) and chemical elements in stars strongly affects their evolution, from premain sequence (PMS) to evolved stages. These processes are particularly crucial in the stellar radiative zones, but their modelling remains an open question. In standard evolution models, these stably stratified zones are assumed to be motionless despite early works pointing out the lack of a static solution in uniformly rotating stars (Von Zeipel 1924) and a related largescale meridional circulation driven by the centrifugal acceleration (Eddington 1925; Sweet 1950). A more complete formalism was then introduced by Zahn (1992) including a selfconsistent meridional flow and models of the turbulent transport driven by hydrodynamical instabilities. The importance of additional processes such as internal waves (e.g., Talon & Charbonnel 2005) and magnetic fields (Spruit 2002) has also been investigated. Zahn’s formalism successfully explained a number of observed stellar properties, such as the nitrogen abundances at the surface of red supergiants or the observed bluetored supergiant ratio in the Small Magellanic Cloud (Maeder & Meynet 2001). The dynamics of internal gravity waves could possibly explain the flat rotation profile of the solar radiative zone (Talon et al. 2002 and Charbonnel & Talon 2005) inferred through helioseismology (Schou et al. 1998), as well as the lithium dip in solarlike stars (Charbonnel & Talon 2005). Despite these early encouraging results many stellar observations remain unexplained, especially the internal rotation rates revealed by asteroseismology in evolved stars and in intermediatemass main sequence stars (see Aerts et al. 2019 for a review). The current theoretical understanding of the structure of the differential rotation in stellar radiative zones, crucial for the development of instabilities and thus for the related turbulent transport, is still largely incomplete. In particular, the assumption that the differential rotation is mostly radial, rather than radial and latitudinal, is at the base of Zahn’s formalism, but the validity of this assumption has never been thoroughly tested. In this paper we are particularly interested in the differential rotation and the largescale meridional flows generated in periods of the stellar life when contraction, expansion, or both processes occur. This is the case for example for PMS stars which are gravitationally contracting before starting their core nuclear reactions, or for subgiant and giant stars which undergo contraction of their core and expansion of their envelope. Within these stars a strong redistribution of AM is expected to happen, producing differential rotation and thus potentially unstable shear flows.
Solarlike PMS stars such as T Tauri stars contract during their evolution up to the main sequence (MS). The increase in the rotation rate due to this contraction can be estimated using AM conservation and assuming uniform rotation. For example, taking typical rotation rates from Gallet & Bouvier (2013), a star rotating at 8% of the breakup velocity at 1 Myr that experiences a decrease in radius from 2.58 R_{⊙} to 0.9 R_{⊙} (while the moment of inertia goes from 19.6 I_{⊙} to 1.03 I_{⊙}, where I_{⊙} = 6.41 ⋅ 10^{53} g⋅cm^{2} is the moment of inertia of the present Sun) would rotate at 32% of the breakup velocity at the zero age main sequence (ZAMS). The observation of stars in young open clusters instead indicates that their surface typically rotates at only 10% of their breakup velocity, showing the need to add other processes like AM exchanges with a circumstellar disk and AM extraction through a magnetised stellar wind (Bouvier et al. 1986). Interestingly, simple models of the internal AM transport, where a coupling time between a uniformly rotating outer convective zone and a uniformly rotating radiative core determines the surface rotation, strongly suggest that a certain level of radial differential rotation is present in the interior of these contracting young solartype stars (Gallet & Bouvier 2013).
Another stellar class of interest for our purpose are the more massive Herbig Ae/Be stars. Just like the T Tauri stars, these objects experience a contraction during the PMS phase which should lead them to spin up to the MS. Using typical values from Böhm & Catala (1995) and assuming uniform rotation, we can estimate that a 2 M_{⊙} star rotates beyond its breakup velocity on the ZAMS while a 3 M_{⊙} star should rotate at around 65% of its breakup velocity and a 5 M_{⊙} at around 52% breakup. Again, observations indicate lower surface rotation rates on the ZAMS, corresponding to 30% breakup for a 2 M_{⊙} star and around 40% for the 3 and 5 M_{⊙} stars (Alecian et al. 2013a). These data would be compatible either with a ∝1/r^{2} inner differential rotation and no AM losses or with winddriven AM losses and smoother radial gradients. An additional interest of the Herbig Ae/Be stellar class is the small fraction of Herbig Ae/Be stars that are known to be strongly magnetic (Wade et al. 2005; Alecian et al. 2013b) because the contractiondriven AM transport in their radiative envelope would surely be affected by these magnetic fields.
When studying AM transport, red giant stars and particularly early red giants and subgiants are objects of particular interest. Asteroseismology puts constraints on the rotation rates of both the core and the envelope of these stars, and was able to show for the first time that a significant differential rotation exists in a stellar radiative zone (Deheuvels et al. 2012; Beck et al. 2012). By studying a sample of six Kepler stars at different evolutionary stages, Deheuvels et al. (2014) also showed that the level of differential rotation increases with the stellar age. During the subgiant phase the combination of core contraction and envelope expansion indeed tends to produce a strong radial differential rotation. However, uptodate hydrodynamical models based on Zahn’s formalism predict core rotation rates that are two or three orders of magnitude higher than observed (Eggenberger et al. 2012; Marques et al. 2013; Ceillier et al. 2013), thus showing that some model assumptions are not valid or that additional processes extract the AM from the core. Internal gravity waves (Pinçon et al. 2017) or magnetic fields (Fuller et al. 2019; Eggenberger et al. 2019; den Hartogh et al. 2020) have been invoked as possible candidates.
In this work we investigate the largescale flows induced by a steady mass flux through a stably stratified radiative zone, thus modelling the contraction of the star either during the PMS or the subgiant and giant phases. More specifically, we describe the structure of the steady axisymmetric solutions (i.e. the differential rotation and the meridional circulation) in the different physical regimes relevant for stars undergoing contraction. To our knowledge, such studies are scarce. In Hypolite & Rieutord (2014), the effect of the gravitational contraction was studied in a linear regime without stratification or assuming an approximate treatment of the stable stratification proposed by Rieutord (2006). The expansion was considered with a similar approach by Rieutord & Beth (2014). The differential rotation induced by the contraction can in principle give rise to various kinds of instabilities that would in turn affect the AM distribution. In this context our present work on steady axisymmetric solutions can be viewed as a first step towards a future study on the instabilities of contractiondriven flows and their role on the AM transport in stars.
The paper is organised as follows. In Sect. 2 the mathematical model is presented. In Sect. 3 we discuss the different physical timescales. After shedding light on the values of the relevant parameters in stars in Sect. 4, the numerical method is presented in Sect. 5. The results under the Boussinesq and anelastic approximations are then respectively described in Sects. 6 and 7. Finally, a summary and a conclusion are given in Sect. 8.
2. Mathematical formulation
In this work we model the contraction of a portion of stellar radiative zone. To this end, the fluid dynamics equations are solved numerically using one of two different approaches, either the Boussinesq approximation or the anelastic approximation. These two approximations are based on a soundproof approach where the acoustic waves are filtered out. The Boussinesq case is studied to analyse the response to contraction of an incompressible fluid, without losing the effects of buoyancy. Then the anelastic approximation is introduced to study the effects of the density stratification. In this section we first present the dimensionless Boussinesq and anelastic equations then discuss the reference state and the choice of the forcing velocity field that mimics the contraction of the stellar radiative zone.
2.1. Governing equations
To model a stellar radiative zone we consider a spherical shell of inner radius r_{i} and outer radius r_{0} filled with a selfgravitating ideal gas subject to a forced radial velocity field intended to mimic the stellar contraction. In the following the nondimensional form of the governing equations in a rotating frame is obtained using the radius of the outer sphere r_{0} as the reference length, the value of the contraction velocity field at the outer sphere V_{0} as a characteristic velocity, and the contraction timescale τ_{c} = r_{0}/V_{0} as the reference timescale. The frame rotates at Ω_{0}, the rotation rate of the outer sphere. All the thermodynamics variables are expanded as background plus fluctuations, the background variables being identified with an overbar and the fluctuations with a prime. The background density and temperature field are nondimensionalised respectively by the outer sphere density ρ_{0} and temperature T_{0}, while the background gradients of temperature and entropy fields are adimensionalised using the temperature and entropy difference and between the two spheres. The pressure fluctuations are nondimensionalised by ρ_{0}r_{0}Ω_{0}V_{0} and the entropy fluctuations by C_{p}Ω_{0}V_{0}/g_{0}, where C_{p} is the heat capacity and g_{0} the gravity at the outer sphere. Hereafter, all dimensionless variables are identified with a tilde. Under these conditions, the dimensionless form of the governing equations for conservation of mass, momentum, and entropy is given by (Jones et al. 2011; Gastine & Wicht 2012; Wicht et al. 2018; Dietrich & Wicht 2018)
where V_{f}(r) is the massconserving contraction velocity field whose expression is given in Sect. 2.3. The dimensionless momentum equation Eq. (2) was derived using the LantzBraginskyRoberts (LBR) approximation (Lantz 1992; Braginsky & Roberts 1995) and assuming that the kinematic viscosity ν is uniform. The dimensionless equation of entropy evolution Eq. (3) was obtained using the formalism described by Braginsky & Roberts (1995) and Clune et al. (1999) to introduce the diffusion of entropy instead of the diffusion of temperature. The thermal diffusivity κ is assumed to be uniform. In the above equations is the dimensionless stress tensor and is the dimensionless viscous heating. Furthermore, the gravity model assumes that the mass of the star is concentrated in the centre:
Here G is the gravitational constant and M_{c} is the mass at the centre.
In Eqs. (1)–(3) R_{o} is a Rossby number based on the amplitude of the contraction velocity field R_{o} = V_{0}/Ω_{0}r_{0}, E is the Ekman number , P_{r} is the Prandtl number P_{r} = ν/κ, and N_{0} is a reference value of the Brunt–Väisälä frequency defined in Sect. 2.2. It gives a mean measure of the magnitude of the stable stratification. The ratio of the Rossby number to the Ekman number allows us to define the contraction Reynolds number Re_{c}, from which we define the Péclet number associated with the contraction Pe_{c} = P_{r} Re_{c}. In Eq. (3) a parameter related to the density stratification appears, namely the dissipation number:
The Boussinesq approximation is then found in the limit where the effects of compressibility are neglected, except in the buoyancy term of the momentum equation. In this case the temperature fluctuations are adimensionalised using Ω_{0}V_{0}T_{0}/g_{0} as the temperature scale. The dimensionless governing equations for continuity, momentum evolution, and temperature evolution become
In Eq. (7) we use a gravity profile corresponding to a uniform fluid density ρ_{c}:
To conclude this section we note that when the anelastic approximation is used the parameter space is defined by six dimensionless numbers: ϵ_{s} (defined in Sect. 2.2), D_{i}, Re_{c}, E, P_{r}, and . Instead, only the last four are necessary in the Boussinesq approximation.
2.2. Background state
In the anelastic approximation we use a nonadiabatic reference state which enables us to prescribe a profile for the background entropy gradient. We simply use a uniform positive entropy gradient to produce a stable stratification. The uniform background entropy gradient is related to the Brunt–Väisälä frequency through
from which we define a reference Brunt–Väisälä frequency:
For the anelastic approximation to remain valid, we nonetheless need to ensure that the deviation to the isentropic state is small. The magnitude of this deviation is controlled by the parameter ϵ_{s}:
This parameter, together with the dissipation number, sets the temperature profile and the magnitude of the density stratification.
Finally, the prescribed background gradient of entropy , the dissipation number D_{i}, and the ϵ_{s} parameter give the background temperature profile through the equation
as well as the background density profile
where γ is the adiabatic index that is equal to 5/3 for a monatomic gas. When the Boussinesq approximation is used, the background temperature is the solution of the conduction equation
for prescribed temperatures at the spherical boundaries:
Here , again ensuring stable stratification. The background temperature gradient is then related to the Brunt–Väisälä frequency through the relationship
from which the reference Brunt–Väisälä frequency is defined as
2.3. Contraction velocity field
The specificity of our model is to introduce a contraction process that will force largescale flows inside the stably stratified stellar layer. We assume here that the contraction is produced by a steady mass flux oriented towards the centre of the star. In other words, a fixed radial velocity field V_{f} is added to the flow velocity, and this forcing velocity field is chosen to be radial and directed towards the centre of the star. In order for this additional velocity field to conserve mass, we need to ensure that . In the anelastic case the contraction velocity field thus simply reads:
Under the Boussinesq approximation, since the background density is uniform, this velocity field simplifies into
We would like to point out that to model stellar expansion, this forcing velocity field V_{f} (r) can be chosen to be positive.
3. Timescales of physical processes
In this section we use the equation of AM conservation to derive the timescales of the physical processes acting on the differential rotation. This allows us to distinguish a priori the parameter regimes where each process dominates.
In its dimensional form, the equation of AM conservation reads
where D^{2} = (∇^{2}−1/r^{2}sin^{2}θ) is the azimuthal component of the vector Laplacian operator, U_{s} = cos θ U_{θ} + sin θ U_{r} is the cylindrical velocity field, and NL denotes the nonlinear advection term. By considering the balance between the partial time derivative of azimuthal velocity field and the contraction term, we recover the contraction timescale defined in Sect. 2.1, namely τ_{c} = r_{0}/V_{0}, when ΔΩ/Ω_{0} ≳ 1, or its linear version , when ΔΩ/Ω_{0} ≪ 1 so that rU_{ϕ} ≪ r^{2} sin θ Ω_{0}.
From the balance with the viscous and Coriolis terms, we now define three additional timescales for the AM transport. The first is the viscous timescale, which characterises a redistribution of AM by viscous processes:
The two other timescales are related to the AM transport by meridional flows. In the linear regime the timescale of AM redistribution by the meridional motions is
and thus it depends on the nature of the circulation through its velocity scale U_{s}. We now consider two wellknown types of circulation, the Eddington–Sweet and the Ekman circulations, that will be relevant to analyse our simulation results. In this paper the steadystate Eddington–Sweet circulation debated by Busse (1981) or Busse (1982), is consistent with the conservation of AM (see Sect. 8).
In a stably stratified medium, largescale meridional motions are efficiently prevented as the increasing buoyancy content of rising (or sinking) fluid elements produces strong restoring buoyancy forces. However, by diminishing the amplitude of temperature fluctuations, and thus of the fluid element buoyancy, thermal diffusion enables the Eddington–Sweettype circulation. To derive its characteristic timescale we use the thermal wind balance, obtained by taking the curl of Eq. (A.17) projected in the azimuthal direction:
We then use the steady linear dimensional form of the thermal balance Eq. (A.16)
where and the temperature deviations δΘ(r, θ) = Θ′(r, θ)−T′(r) are respectively the mean spherically symmetric temperature field and the associated fluctuations obtained by taking the additional radial thermal background T′(r) induced by the contraction field into account (for details, see Appendix A).
Accordingly, the characteristic amplitudes of the differential rotation ΔΩ, the temperature deviations δΘ, and the radial velocity field U_{r} are related as
from which we get the following relationship between the radial velocity field and the level of differential rotation:
Then, from the continuity equation, we have U_{θ} ≈ U_{r} ≈ U_{s}, and by injecting the expected order of magnitude of U_{s} in Eq. (23), we obtain the Eddington–Sweet timescale:
In an unstratified medium, the Ekman pumping mechanism produces a largescale circulation controlled by the boundary layers. The mismatch between boundary conditions applied at spherical surfaces and the inviscid interior flow solutions leads to thin viscous boundary layers that in practice modify the boundary conditions seen by the interior flow. For example, the noslip boundary conditions at the outer sphere induce an Ekman layer that imposes the following relationship between the azimuthal and the radial velocity fields,
as the boundary condition satisfied by the interior flow at the outer sphere. A wellknown example of Ekman circulation is the classical unstratified spherical Couette flow where the necessary jumps in rotation rate between the interior and the spherical boundaries, denoted ΔΩ_{J}, force a meridional circulation of the order of that extends over the whole domain (L is the characteristic length in the Couette problem). From Eq. (23), the associated timescale for the AM transport, also called the Ekman spinup timescale, is
For a stably stratified medium, Ekman layers can still exist at the boundaries because the buoyancy force has a negligible effect on motions of very small radial extent, although they are not expected to drive a largescale global circulation (Barcilon & Pedlosky 1967). As we discuss later, an Ekman layer is present at the outer sphere of our numerical simulations to connect the interior circulation to the prescribed boundary condition, and because the Ekman number is small but finite this has a nonnegligible effect on the resulting differential rotation.
The ratios that compare the relative importance of these three AM transport processes,
are controlled by two dimensionless numbers, the Ekman number E and the P_{r}(N_{0}/Ω_{0})^{2} parameter. As we are interested in cases with low Ekman numbers, , we can thus define three regimes of interest, namely
The actual differential rotation will result from a competition between the forcing associated with the contraction and these three AM transport processes. The corresponding timescale ratios are as follows:
This introduces an additional parameter, either R_{o} or Re_{c}, which will have a direct impact on the level of differential rotation. In our numerical study this parameter is varied to explore the linear regime, ΔΩ/Ω_{0} ≪ 1, and the beginning of the nonlinear regime, ΔΩ/Ω_{0} ≳ 1.
A nonlinear regime is indeed expected if the contraction is the dominant transport process. In this case, the steady state of the AM conservation Eq. (21) reads
which can be solved to obtain
Thus, for our setup with r_{i} = 0.3 r_{0}, the maximum value of the differential rotation is
which is clearly in the nonlinear regime. We thus expect to reach this regime when the timescale ratios Eq. (34) are higher than one, while a linear regime ΔΩ/Ω_{0} ≪ 1 should exist when these timescale ratios are low enough.
4. Estimate of the relevant physical parameters in stars
In this section we estimate the Ekman number E, the P_{r}(N_{0}/Ω_{0})^{2} parameter, and the contraction Reynolds number Re_{c} in stars undergoing contraction.
4.1. Premain sequence stars
Among PMS stars, lowmass stars are characterised by a bimodal distribution of their rotational period (P_{rot} ≈ 8 days for those experiencing a disklocking phase and P_{rot} ≈ 1 day for the others; Maeder 2008), while intermediate and highmass stars are rather rapid rotators (P_{rot} ≈ 1 day) except when they are magnetised (P_{rot} ≈ 1.6 − 4.3 days) (Alecian et al. 2013a).
Stellar models of a 3 M_{⊙} Population I PMS star (Talon & Charbonnel 2008) show that the thermal part of the Brunt–Väisälä frequency varies from at the beginning of the PMS to at the ZAMS. The Prandtl number is very low in nondegenerate stellar interiors, typically P_{r} ≈ 10^{−6} (Garaud et al. 2015). We thus obtain the orders of magnitude of P_{r}(N_{0}/Ω_{0})^{2} in PMS stars listed in Table 1.
Estimates of P_{r}(N_{0}/Ω_{0})^{2} in PMS stars.
Another relevant parameter is the contraction Reynolds number Re_{c}. Since this number can be expressed as the ratio of the contraction timescale to the viscous timescale (Re_{c} = τ_{ν}/τ_{c}), a rough estimate can be derived by replacing the contraction timescale by the Kelvin–Helmholtz timescale . In this case , thus τ_{ED}/τ_{c} = Re_{c} P_{r}(N_{0}/Ω_{0})^{2} ≈ (N_{0}/Ω_{0})^{2} is always greater than one for PMS stars, regardless of their evolutionary status or rotation period.
The Ekman number is always very small in stellar interiors. According to Lara & Rieutord (2013), at the ZAMS the radiative viscosity in stellar radiative zones ranges from ν ≈ 8 ⋅ 10^{1} cm^{2} ⋅ s^{−1} in the core to ν ≈ 5 ⋅ 10^{2} cm^{2} ⋅ s^{−1} in the middle of the envelope for a 3 M_{⊙} star, while its value remains approximately constant between the core and the middle of the envelope for a 7 M_{⊙} or a 10 M_{⊙} star (around ν ≈ 10^{3} cm^{2} ⋅ s^{−1}). During the PMS, the radius of a star varies a great deal. As an example, for a solarmass star, the radius ranges from 2.6 R_{⊙} during the early PMS phase to 0.9 R_{⊙} for more advanced stages of the PMS (Maeder 2008; Gallet & Bouvier 2013). During the PMS evolution of more massive stars, the radius can decrease from 8 R_{⊙} to 3.4 R_{⊙} for a 7 M_{⊙} star, and from 4.5 R_{⊙} to 2.3 R_{⊙} for a 3 M_{⊙} star (Alecian et al. 2013b). For such values, the Ekman number is indeed always very small, varying from E = 10^{−15} to E = 10^{−17}. We thus conclude that the dominant parameter regime for PMS stars is the nonlinear Eddington–Sweet regime defined by
4.2. Subgiants
In the degenerate cores of subgiants, the Prandtl number increases to ∼10^{−3} due to the degeneracy of electrons (Garaud et al. 2015). Meanwhile, the Brunt–Väisälä frequency in these dense cores is quite large, its value varying between at the beginning of the subgiant phase to at a more advanced stage (Talon & Charbonnel 2008). Asteroseismology puts strong constraints on the core rotation rate of these stars (Deheuvels et al. 2014). Considering a typical value of Ω_{core} = 3.34 ⋅ 10^{−6} rad ⋅ s^{−1}, we find that the degenerate core is characterised by very large values of the parameter P_{r}(N_{0}/Ω_{0})^{2} up to ∼10^{5}.
Outside the degenerate core, typical P_{r}(N_{0}/Ω_{0})^{2} are 10^{−3}, much smaller than one, as in PMS stars. Given the small Ekman numbers, we find that the nonlinear viscous regime is the relevant one inside the degenerate core of the subgiants,
whereas outside the degenerate core the nonlinear Eddington–Sweet regime dominates:
5. Numerical method
In order to study numerically the hydrodynamical steady states obtained in the various regimes described above, we carry out numerical simulations using the pseudospectral code MagIC. MagIC is based on a Chebyshev discretisation in the radial direction and spherical harmonic discretisation for the latitudinal and azimuthal directions. MagIC is a fully documented, publicly available code^{1} which solves the (magneto)hydrodynamical equations in a spherical shell. An axisymmetric version of Wicht (2002) for the Boussinesq approximation and of Gastine & Wicht (2012) for the anelastic approximation is used to solve Eqs. (2) and (3) as well as Eqs. (7) and (8) numerically. The solenoidal condition of Eqs. (1) and (6) is ensured by a poloidal–toroidal decomposition of the mass flux. Our domain extent is r ∈ [r_{i} = 0.3;r_{0} = 1.0] and θ ∈ [0;π].
For all simulations we impose stressfree conditions at the inner sphere for the latitudinal and azimuthal velocity fields, and impermeability condition for the radial velocity field:
At the outer sphere we impose an impermeability condition on the radial velocity field, a noslip condition on the latitudinal velocity field, and uniform rotation:
The boundary conditions for Θ′ (respectively S′) are set to zero at the top and bottom of the domain in the Boussinesq (respectively anelastic) case. Initially we impose a uniform rotation in the whole domain Ω(r, θ, t = 0) = Ω_{0} in the Boussinesq case and a fixed level of differential rotation
in the anelastic case, where σ controls the magnitude of the rotation contrast.
For most of the simulations the grid resolution is N_{r} × N_{θ} = 193 × 256. The large number of points in the radial direction is necessary because we also aim to characterise the different boundary layers and their potential roles. When the Ekman number becomes very small (typically 10^{−7}), the resolution is sometimes increased to N_{r} × N_{θ} = 217 × 288. We always ensure that the number of radial grid cells in the Ekman layer is nonnegligible. Thus, for the best resolved cases the number of radial grid cells in the Ekman layer ranges from 5 grid cells in the E = 10^{−7} case to 12 grid cells in the E = 10^{−5} case.
6. Results in the Boussinesq case
The analysis of the different transport mechanisms performed in Sect. 3 led to define three different regimes depending on the Ekman number and the P_{r}(N_{0}/Ω_{0})^{2} parameter. A third parameter, Re_{c} or R_{o}, appears to control the level of differential rotation, and thus whether the nonlinear terms play a role. The numerical investigation presented below has been carried out by varying these parameters and the results are analysed with the guidance of the timescale analysis of Sect. 3.
6.1. Taylor–Proudman regime
We find that the differential rotation is near cylindrical when . According to the thermal wind balance Eq. (24), such a regime is possible if the buoyancy force is weak enough. In the following the regimes of cylindrical rotation are called Taylor–Proudman regimes, although strictly speaking the Taylor–Proudman theorem states that all the velocity components must be cylindrical.
The meridional cut in the left panel of Fig. 1 shows the structure of the flow for a simulation in this regime (E = 10^{−4}, P_{r}(N_{0}/Ω_{0})^{2} = 10^{−4}, Re_{c} = 10^{−2}). The differential rotation exhibits a cylindrical profile, the rotation rate being highest near the tangent cylinder. For this particular simulation, the maximum amplitude of the normalised differential rotation (Ω(r,θ)−Ω_{0})/Ω_{0} is barely 0.001, which thus corresponds to a linear case. The right panel of the same figure displays the norm and the streamlines of the total meridional flow . We see that a fluid particle initially located at the outer sphere is deflected into the Ekman layer towards the pole before it enters into the interior flow. From then on, if the particle is outside the tangent cylinder, it goes back towards the outer sphere thus forming a cell. On the contrary, if the particle is inside the tangent cylinder, it goes towards the inner sphere. A remarkable feature is the strong vertical jet directed towards the equator of the inner sphere.
Fig. 1. Steady axisymmetric flow in the Taylor–Proudman regime. In the left panel, the coloured contours represent the differential rotation δΩ(r, θ) = Ω(r, θ)−Ω_{0} normalised to the rotation rate of the outer sphere Ω_{0} while in the right panel, we show the norm of the total meridional velocity field . The parameters of the simulation are E = 10^{−4}, P_{r}(N_{0}/Ω_{0})^{2} = 10^{−4} and Re_{c} = 10^{−2} (simulation 1.3 in Table 2). The black lines show the streamlines of the total meridional circulation i.e. taking the effect of the contraction into account. Outside the tangent cylinder, the dashed lines correspond to an anticlockwise circulation and the solid lines to a clockwise one while inside the tangent cylinder, the dashed lines represent a downward circulation and the full ones, an upward circulation. 
In order to characterise this regime we performed several simulations for a range of Ekman and contraction Reynolds numbers, E = 10^{−2} − 10^{−6}, Re_{c} = 10^{−2} − 10, as listed in Table 2. These numerical results are compared with an analytical solution that can be obtained in the linear and low Ekman number regime. Considering the balance between the contraction term and the Coriolis term in the linear inviscid steady equation of AM evolution Eq. (21), we first get the cylindrical radial component of the velocity field:
Parameters of the simulations carried out in the Taylor–Proudman regime.
The full inviscid meridional flow is then determined using the U_{r} = 0 condition at the inner sphere. This solution does not fulfil the boundary conditions at the outer sphere; this is done through an Ekman boundary layer, which comes with a jump in the rotation rate. Then, using the Ekman pumping condition Eq. (30) as a boundary condition for the interior flow, it is possible to express a mass budget inside a cylinder of radius s = r sin θ, which involves the inward flow from the Ekman layer and the outward flow through the cylinder. The inward and the outward mass fluxes are obtained by integrating Eqs. (30) and (44), respectively. After some algebra (see details in Appendix B) we obtain the following^{2}:
This solution can be interpreted as follows. At leading order the contraction forces a global meridional circulation that is unable to match the boundary conditions at the outer sphere. This is done in an Ekman boundary layer, which in turn produces a certain level of differential rotation determined by the Ekman pumping relation for the prescribed interior radial field Eq. (30). This differential rotation is communicated to the bulk of the flow through the Taylor–Proudman constraint ∂U_{ϕ}/∂z = 0.
The analytical solution scaled by is plotted with a black dashed line in the left panel of Fig. 2 and compared with numerical solutions obtained for different Ekman and Rossby numbers. It shows that the flow behaves differently inside and outside the tangent cylinder. The rotation rate is highest on the tangent cylinder. At this location the analytical inviscid solution is not differentiable and the vertical velocity is also discontinuous because it tends to minus infinity when the tangent cylinder is approached from the left. In low Ekman number flows such discontinuities are expected to produce boundary layers whose thickness decreases with the Ekman number. This is visible in Fig. 2, which also shows, as expected, that the agreement with the analytical solution improves for smaller Ekman numbers. The right panel further shows that for low enough Ekman numbers and in the linear regime the global differential rotation between the inner and the outer spheres taken along a given radius scales as
Fig. 2. Axisymmetric steady solutions of the differential rotation as a function of radius at the latitude θ = π/4 in the Taylor–Proudman regime. Left panel: analytical solution Eq. (45) is plotted in black dashed lines and compared with the numerical simulations obtained for fixed P_{r}(N_{0}/Ω_{0})^{2} parameter (10^{−4}) and contraction Reynolds number Re_{c} (10^{−2}) with various Ekman numbers: 10^{−2} (in grey), 10^{−3} (in cyan), 10^{−4} (in light blue), 10^{−5} (in dark blue), and 10^{−6} (in purple) respectively for runs 1.1–1.5 in Table 2. For E = 10^{−5} and E = 10^{−6}, additional contraction Reynolds numbers are studied, namely 10^{−1} (runs 2.1 and 3.1 in Table 2), 1 (runs 2.2 and 3.2 in Table 2), and 10 (runs 2.3 and 3.3 in Table 2). All curves are rescaled by . Right panel: differential rotation between the inner and outer spheres ΔΩ = Ω(r_{i})−Ω(r_{0}) at the same latitude and normalised with the rotation rate taken at the outer sphere Ω_{0}, now plotted as a function of . The previous numerical solutions obtained for E = 10^{−5} (runs 1.4 and 2.1–2.3 in Table 2) are plotted in dark blue, while the simulations performed for E = 10^{−6} (runs 1.5 and 3.1–3.3 in Table 2) are plotted in purple. 
Finally, the existence of the Taylor–Proudman regime in a stably stratified atmosphere where can be explained as follows. From Eq. (44) we see that U_{r} ≈ V_{0}. Thus, according to the thermal balance Eq. (25), δΘ ≈ P_{r}(N_{0}/Ω_{0})^{2} E^{−1} V_{0}, while according to Eq. (30), U_{ϕ} ≈ E^{−1/2} V_{0}. The ratio of these two quantities is . As a consequence, in the regime , the effect of the buoyancy forces in the thermal wind equation Eq. (24) is too weak and the Taylor–Proudman constraint ∂U_{ϕ}/∂z = 0 holds.
6.2. Viscous regime
According to the timescale analysis of Sect. 3, the viscous transport of AM dominates when 1 ≪ P_{r}(N_{0}/Ω_{0})^{2}. To investigate this regime, numerical simulations have been performed for values of P_{r}(N_{0}/Ω_{0})^{2} varying between 10^{2} and 10^{ 4}. The different runs characterised by the three parameters, P_{r}(N_{0}/Ω_{0})^{2}, E, Re_{c}, are listed in Table 3.
Parameters of performed simulations in the viscous regime.
In this parameter range we find that the stationary differential rotation is mostly radial and that the level of differential rotation increases with Re_{c}. This is illustrated in the first and third panels in Fig. 3 where the differential rotation is shown for two runs with equally high P_{r}(N_{0}/Ω_{0})^{2} = 10^{ 4} and different Re_{c}. The differential rotation is closely radial in both cases, while the rotation contrast between the outer and the inner spheres goes from ∼10^{−2} for Re_{c} = 10^{−2} to ∼8 as the contraction rate is increased to Re_{c} = 10. The other two panels display the total meridional velocity field, through its norm and streamlines, showing that it is clearly dominated by the imposed radial contraction flow. Subtracting this contribution, we would see a meridional circulation confined into boundary layers at the outer and inner spheres, the confinement being stronger for smaller Ekman numbers.
Fig. 3. Differential rotation normalised to the rotation rate of the outer sphere δΩ(r, θ)/Ω_{0} in the viscous regime (first and third panels) and norm of the total meridional velocity field (second and fourth panels) with the streamlines of the associated total meridional circulation (in black). In the two left panels the parameters of the simulation are P_{r}(N_{0}/Ω_{0})^{2} = 10^{ 4}, Re_{c} = 10^{−2}, and E = 10^{−5} (run 3.4 in Table 3), and in the two right panels the parameters are P_{r}(N_{0}/Ω_{0})^{2} = 10^{4}, Re_{c} = 10, and E = 10^{−4} (run 4.3 in Table 3). In the second and fourth panels the black solid lines represent the streamlines of an ascending and inwardly directed circulation, while the black dashed lines correspond to a descending and inwardly directed one. 
The results of the numerical simulations are now compared with the analytical solution that can be derived by neglecting the Coriolis and the nonlinear advection terms in the AM conservation Eq. (21). The steady balance between the contraction and viscous terms indeed reads
whose solution is given by the following expression:
If the differential rotation δΩ_{ ν}(r)/Ω_{0} is less than one, the forcing term in Eq. (47) is reduced to , and the analytical solution for the differential rotation simplifies into
which implies that this approximate solution should be valid for low Reynolds numbers Re_{c}.
The analytical solution is compared with the numerical solutions in Fig. 4. The agreement is very good showing that the viscous term indeed dominates the transport in the regime P_{r}(N_{0}/Ω_{0})^{2} ≥ 10^{2} considered. We also find that the linear relation ΔΩ/Ω_{0} ∝ Re_{c} = τ_{ν}/τ_{c} valid for low Re_{c} still provides a good order of magnitude estimate of the differential rotation between the inner and the outer spheres up to the highest Re_{c} considered.
Fig. 4. Differential rotation in the viscous regime scaled by Re_{c} as a function of radius at the latitude θ = π/8. The analytical solution Eq. (48) (black dashed lines) is compared with the numerical results. Low Reynolds number runs (1.1–4.1 in Table 3, with Re_{c} = 10^{−2} − 10^{−1}) are plotted with violet dotted lines and are all very close to the low Reynolds number version of the analytical solution, Eq. (49). Simulations at Re_{c} = 1 (runs 4.2 and 5.1 in Table 3) and Re_{c} = 10 (run 4.3 Table 3) are also in agreement with the analytical solution. 
In this subsection we show that for high enough P_{r}(N_{0}/Ω_{0})^{2} ≳ 10 stable stratification is so strong that the meridional circulation is inhibited and the effect of the contraction is balanced by the viscous transport of AM. An analytical solution of the resulting radial differential rotation is available.
6.3. Eddington–Sweet regime
When the stratification is not too strong, , we expect that the transport of AM is dominated by an Eddington–Sweet type meridional circulation.
The meridional cut displayed in the left panel of Fig. 5 shows the typical flow structure obtained in this regime. The differential rotation exhibits a dependence on both latitude and radius. Thus, contrary to the Taylor–Proudman and viscous regimes, it is neither cylindrical nor spherical. Meanwhile the total meridional circulation (right panel) is similar to the one already described in the Taylor–Proudman case with an inward flow concentrated inside the tangent cylinder and more particularly in a vertical jet towards the equator of the inner sphere. To understand what determines the level and distribution of the differential rotation in the Eddington–Sweet regime, we performed a parametric study where E, P_{r}(N_{0}/Ω_{0})^{2}, and Re_{c} were varied; the values of these parameters are listed in Table 4.
Fig. 5. Steady differential rotation (left panel) normalised to the top value, and norm and streamlines of the total meridional velocity field (right panel). The parameters of the simulation (run 1.4 in Table 4) are E = 10^{−6}, P_{r}(N_{0}/Ω_{0})^{2} = 10^{−1}, and Re_{c} = 10^{−2}. Outside the tangent cylinder the black dashed lines represent the streamlines of an anticlockwise circulation, and the solid lines those of a clockwise circulation. Inside the tangent cylinder the dashed lines simply correspond to the streamlines of a downward circulation, while the solid lines are for an upward circulation. 
Parameters of the simulations in the Eddington–Sweet regime.
Owing to the flow complexity, we were not able to find an analytical solution for the azimuthal velocity field. We instead compare our numerical results with a scaling relation for the amplitude of differential rotation, which we derive as follows. Assuming a linear regime, the balance between the Coriolis term and the contraction term in the AM evolution equation Eq. (21) implies that the circulation velocities scale as the imposed contraction radial velocity field:
Then the thermal wind balance Eq. (24) and the thermal balance Eq. (25) provide an estimate for the radial velocity amplitude given by Eq. (28). Injecting Eq. (50) in Eq. (28), we obtain
which is equivalent to
Physically, it says that the steady differential rotation is determined by the balance between the inward transport of AM forced by the contraction with timescale , and the outward transport of AM by the Eddington–Sweet circulation with timescale τ_{ED}.
Coming back to the numerical results, the left panel in Fig. 5 shows that the differential rotation taken between the inner and the outer spheres depends on the latitude considered at the inner sphere, and that it is greatest at the equator. Figure 6 displays the differential rotation between the equator of the inner sphere and the outer sphere for simulations performed at various R_{o}, considering two different P_{r}(N_{0}/Ω_{0})^{2} = 10^{−1} and 10^{−2}, and a fixed low Ekman number E = 10^{−6}.
Fig. 6. Differential rotation between the inner equator and the outer sphere normalised to the outer value of the rotation rate Ω_{0} in the Eddington–Sweet regime, plotted as a function of P_{r}(N_{0}/Ω_{0})^{2}R_{o}/E. The Rossby number R_{o} is varied to encompass the linear and nonlinear regimes. The Ekman number is fixed to E = 10^{−6}. The green curve corresponds to simulations with P_{r}(N_{0}/Ω_{0})^{2} = 10^{−1} (runs 1.4, 3.1, 4.1, 5.1, and 6.1 in Table 4) and the blue curve to simulations with P_{r}(N_{0}/Ω_{0})^{2} = 10^{−2} (runs 2.2, 3.2, 4.2, and 5.2 in the same table). 
This figure first shows that a linear regime exists as long as ΔΩ/Ω_{0} is low enough. The differential rotation is clearly linearly dependent on R_{o} until ΔΩ/Ω_{0} ≳ 10^{−1}, where it starts to deviate. The differential rotation then falls below the linear relation because the nonlinear advection terms, both in the momentum and temperature equations, increase the efficiency of the AM redistribution.
This figure also shows that the scaling Eq. (52) does not reproduce the differential rotation amplitude because the curves corresponding to the two different P_{r}(N_{0}/Ω_{0})^{2} do not collapse. The numerical results would be better reproduced by a relation of the form .
However, we now argue that Eq. (52) should still be approximately valid for the very low stellar Ekman numbers by showing that the observed discrepancy is due to the relatively high Ekman numbers of our simulations. We find that our results are best described by a model where the differential rotation is the superposition of two terms,
the first term being a byproduct of the Ekman layer induced by the mismatch between the inviscid contractiondriven interior circulation (Eqs. (B.8) and (B.9)) and the boundary condition U_{r} = 0 at the outer sphere. This term is actually close to the unstratified Taylor–Proudman solution described in Sect. 6.1 that shows a jump in differential rotation across the outer Ekman layer. The second term scales with P_{r}(N_{0}/Ω_{0})^{2}/E and is of the Eddington–Sweet type. It connects smoothly to the outer sphere boundary condition, that is δΩ_{ED}(r_{0}, θ) = 0. Therefore, the Eddington–Sweet term is dominant in the interior for small enough , whereas the first term is never fully negligible close to the outer sphere. In this model the interior meridional circulation velocities and temperature fluctuations scale as Eqs. (50) and (26), respectively (for details, see Appendix D).
As illustrated in Fig. 7, our numerical results are consistent with this model. Considering flows in the linear regime for two different P_{r}(N_{0}/Ω_{0})^{2} = 10^{−1} and 10^{−2}, and the same Ekman number E = 10^{−6}, we test the relevance of the Eddington–Sweet scaling given by Eq. (52) when applied to the full differential rotation on the one hand, and to the field obtained after subtracting the unstratified solution Eq. (45) on the other. In this second case we observe that the distribution and amplitude of the rescaled differential rotations are very similar, whereas significant amplitude differences persist when the unstratified solution is not subtracted before rescaling. The fact that the scaling Eq. (52) is not accurate for the full differential rotational thus appears to be due to the nonnegligible contribution of the term. This is confirmed by the comparison of the two terms of Eq. (53) displayed in the right column of Fig. 7. This comparison clearly shows that has to be as low as 10^{−2} and that one has to consider layers away from the outer sphere for the Eddington–Sweet term to clearly dominate over the Taylor–Proudman term.
Fig. 7. Differential rotation in the Eddington–Sweet linear regime for P_{r}(N_{0}/Ω_{0})^{2} = 10^{−1} (top panel) and 10^{−2} (bottom panel), the other two parameters E = 10^{−6} and Re_{c} = 10^{−2} being identical (runs 1.4 and 2.2 in Table 4). The coloured contours are normalised to the rotation rate of the outer sphere Ω_{0}. Left column: differential rotation δΩ(r, θ)/Ω_{0} rescaled by E (P_{r}(N_{0}/Ω_{0})^{2}R_{o})^{−1}. Middle column: differential rotation obtained after subtracting the unstratified Taylor–Proudman solution given by Eq. (45) and rescaling by E (P_{r}(N_{0}/Ω_{0})^{2}R_{o})^{−1} (see Eq. (53)). Right column: ratio . 
Runs performed at other Ekman numbers confirm the relevance of this description. In particular, as expected from the expression Eq. (53), simulations with the same ratio (such as runs 1.3 and 2.3 in Table 4) have similar scaled differential rotation .
While the Eddington–Sweet scaling and its Taylor–Proudman correction provide a global description of the differential rotation, additional physics is required to describe regions close to boundary layers. This involves the 𝒪(E^{1/2}) Ekman layer at the outer sphere, the free shear layer that smooths the vertical velocity jump across the tangent cylinder, and also the Ekman layer at the inner sphere induced by the mismatch between the interior flow and the boundary conditions there. As for the classical spherical Couette flow, the Ekman layer at the inner sphere is stronger (i.e. more extended and with higher amplitudes) near the equator than at higher latitudes. As can be observed in Fig. 7, we find that the thicknesses of the layer along the tangent cylinder and of the equatorial layer decrease with E and P_{r}(N_{0}/Ω_{0})^{2}, although determining their exact scaling with these nondimensional numbers might be tricky and is beyond the scope of the present study.
7. Effect of the density stratification: The anelastic case
In this section we analyse the effects of the density stratification using the anelastic approximation. To this end, the magnitude of the background density contrast between the inner and the outer spheres, ρ_{i}/ρ_{0}, is varied from 1.75 to 60.2. Table 5 lists the different density contrasts considered together with the associated parameters of the background state, the dissipation number D_{i}, and the ϵ_{s} parameter (see Eqs. (13) and (14)).
Density contrasts between the inner and outer spheres investigated in the anelastic approximation, and their associated parameters.
The first effect of the density stratification is to modify the profile of the contraction velocity field (see Eqs. (19) and (20)). Figure 8 shows that for a density contrast of 1.75 the absolute value of the contraction velocity at the inner sphere is twice as small as in the Boussinesq case. For stronger contrasts (i.e. when ρ_{i}/ρ_{0} > (r_{0}/r_{i})^{2} ≈ 11.11) the absolute value of the contraction velocity field even decreases with depth.
Fig. 8. Contraction radial velocity field V_{f}, Eq. (19), normalised to the top value V_{0} and plotted as a function of radius r. The lines (from light to dark) show the weakest to the strongest density contrasts between the inner and outer spheres, which are respectively ρ_{i}/ρ_{0} = 1.75, 10.8, 20.9, and 60.2. The contraction radial velocity field in the Boussinesq case given by Eq. (20) (i.e. corresponding to a density ratio of 1) is plotted as a dashed line as a reference. 
7.1. Anelastic Taylor–Proudman regime
We first study the effects of the density stratification in the regime where the Boussinesq simulations have shown that the rotation is cylindrical because the buoyancy force has a negligible effect on the dynamics. The Ekman number, the P_{r}(N_{0}/Ω_{0})^{2} parameter, and the contraction Reynolds number are fixed to 10^{−5}, 10^{−4}, and 10^{−2}, respectively, while the density contrast between the inner and outer spheres is varied from 1.75 to 60.2. Two additional simulations have been also performed at E = 10^{−4} and E = 10^{−6}, the other parameters being fixed. The different runs and their associated parameters are summarised in Table 6.
Parameters of the simulations in the anelastic Taylor–Proudman regime.
Figure 9 shows the steady flow obtained for a density contrast ρ_{i}/ρ_{0} = 20.9. It can be compared to the Boussinesq simulation displayed in Fig. 1, a simulation performed with the same nondimensional numbers but without density stratification. As we can see, the differential rotation is still cylindrical and its amplitude is similar to the Boussinesq case. This property is confirmed at other density ratios (not shown). The main difference in the AM distribution involves the local maximum of the rotation rate along the tangent cylinder, which is present in the Boussinesq case but absent in the anelastic case. The differences in the meridional circulation are more striking because, in the density stratified case, the circulation is no longer characterised by a strong vertical jet towards the inner equator.
Fig. 9. Meridional cuts showing the coloured contours of the steady differential rotation normalised to the top value in the left panel, and the norm of the total meridional velocity field with its associated streamlines (in black) in the right panel, in the Taylor–Proudman regime. As previously, outside the tangent cylinder the solid lines describe a clockwise circulation, and the dashed lines an anticlockwise circulation; inside the tangent cylinder the solid lines represent an upward meridional flow and the dashed lines a downward flow. The parameters of the simulation are E = 10^{−4}, P_{r}(N_{0}/Ω_{0})^{2} = 10^{−4}, Re_{c} = 10^{−2}, and ρ_{i}/ρ_{0} = 20.9 (run 1.1 in Table 6). 
We first interpret the fact that the amplitude of the differential rotation appears to be independent of the density contrast between the two spheres. As in the Boussinesq case, we expect the AM conservation to drive a meridional circulation that balances the effect of the imposed contraction field. In Appendix C a linear and inviscid solution of such meridional flow is given. As in the Boussinesq case, this circulation does not satisfy the boundary conditions at the outer sphere, and an Ekman boundary layer takes charge of it. This induces a differential rotation across the Ekman layer, which is determined by the relation Eq. (30) between the azimuthal velocity jump and the interior radial velocity field. As the meridional velocity field near the outer sphere should not depend on the density contrast, the jump in rotation rate is not affected either. This differential rotation is then communicated to the rest of the flow through the Taylor–Proudman constraint ∂U_{ϕ}/∂z = 0. This reasoning explains why the amplitude of the differential rotation is not affected by the density stratification and why its value is similar to the Boussinesq case.
We further show in Appendix C that the analytical solution found in the Boussinesq case for the differential rotation is also a solution in the anelastic case. Away from the tangent cylinder this solution indeed provides a good approximation to the numerical results. However, the lack of a local maximum along the tangent cylinder observed in the anelastic numerical simulations makes it less relevant than in the Boussinesq case.
The marked difference between the meridional flows in the Boussinesq and anelastic cases is first due to the dependence of the imposed contraction field on the density stratification, with radial velocities decreasing inwards for high enough density ratios. The linear and inviscid solution of the induced meridional circulation (see Appendix C) shows the same behaviour, thus explaining the global distribution of the total velocity observed in the right panel of Fig. 9. We also find that the equatorial boundary layer is practically suppressed at high densitystratification, which in turn explains the lack of a strong vertical jet near the equator. The behaviour close to the inner sphere of three different fields (the radial gradient of angular velocity ∂Ω(r, θ)/∂r, and the cylindrical and vertical velocity fields) is displayed in Fig. 10. In the Boussinesq case we observe that an equatorial boundary layer is present and, as in the classical spherical Couette flow, it is associated with a strong vertical jet along the tangent cylinder. However, when the density contrast increases, the jumps of the three fields near the inner sphere are practically suppressed together with the vertical jet.
Fig. 10. Radial partial derivative of the ratio of the azimuthal velocity field U_{ϕ} to the radius r (left panel), cylindrical radial velocity field U_{s} = U_{r} sin θ + U_{θ} cos θ (middle panel), and vertical velocity field U_{z} = U_{r} cos θ − U_{θ} sin θ (right panel) as a function of radius. In each panel the top row corresponds to the latitude θ = π/4. In the left and middle panels the bottom row corresponds to the equator, but since the vertical velocity field U_{z} is zero at the equator, the bottom right panel is obtained close to this location, at about 85degree angle with the vertical axis. In each case E = 10^{−5}, P_{r}(N_{0}/Ω_{0})^{2} = 10^{−4} and Re_{c} = 10^{−2}. In all panels the plain blue curve corresponds to the Boussinesq case (run 1.4 in Table 2), while the other curves are obtained by varying the density contrast between the inner and outer spheres. Shown in the left panel are ρ_{i}/ρ_{0} = 10.8 (dashed lines), 20.9 (dashdotted lines), and 60.2 (dotted lines), corresponding to runs 2.2 to 2.4 in Table 6; in the middle and right panels the two black curves are respectively obtained for ρ_{i}/ρ_{0} = 1.75 (plain curve) and ρ_{i}/ρ_{0} = 20.9 (dashed lines), corresponding to runs 2.1 and 2.3 in Table 6. 
The absence of the equatorial boundary layer could be explained by the fact that the inviscid interior solution connects in a smoother way to the boundary conditions at the inner sphere. This is indeed the case for the cylindrical radial velocity because the mismatch at the inner equator falls from ΔU_{s} ≈ 11.1 for the Boussinesq inviscid solution to ΔU_{s} ≈ 0.18 when ρ_{i}/ρ_{0} = 60.2.
7.2. Anelastic viscous case
In this section we study, in the viscous regime, the modification to the differential rotation and the meridional flow induced by the introduction of a varying density across the shell. Simulations are performed for the various density contrasts and contraction Reynolds numbers Re_{c} listed in Table 7.
Parameters of the simulations in the anelastic viscous regime.
Figure 11 shows the coloured contours of the rotation rate Ω and the streamlines of the meridional circulation in the steady state for two anelastic simulations with ρ_{i}/ρ_{0} = 20.9. The other parameters are the same as those used in the Boussinesq cases (see Fig. 3). Similar to the Boussinesq case, the differential rotation profile is still radial and the amplitude of differential rotation still increases with Re_{c} (first and third panels of Fig. 11). However, the amplitude of the differential rotation is decreased compared to the Boussinesq case. For ρ_{i}/ρ_{0} = 20.9 the maximum value of the rotation contrast between the inner and the outer spheres is nearly four times smaller than in the Boussinesq case when Re_{c} = 10^{−2}, and almost twice as small as when Re_{c} = 10. The total meridional circulation, visible in the second and fourth panels of Fig. 11, is dominated by the radial contraction field. As in the Boussinesq case, a circulation confined in a thin layer close to the outer boundary can be seen when the forced radial flow is subtracted. However, the circulation at the inner sphere boundary is virtually suppressed by the density stratification.
Fig. 11. Stationary differential rotation δΩ(r, θ)/Ω_{0} (first and third panels) and norm of the total meridional velocity field (second and fourth panels) with its associated streamlines in black, in the viscous anelastic regime. The parameters of the simulation shown in the two left panels are E = 10^{−5}, P_{r}(N_{0}/Ω_{0})^{2} = 10^{ 4}, Re_{c} = 10^{−2}, and ρ_{i}/ρ_{0} = 20.9 (run 5.2 in Table 7), and in the two right panelsE = 10^{−4}, P_{r}(N_{0}/Ω_{0})^{2} = 10^{ 4}, Re_{c} = 10, and ρ_{i}/ρ_{0} = 20.9 (run 4.3 in Table 7). In the first and third panels, the coloured contours are normalised to the top value Ω_{0}. 
We now derive an analytical expression of the differential rotation in the linear regime ΔΩ/Ω_{0} ≪ 1 that clearly shows how the density stratification lowers the level of differential rotation compared to the Boussinesq case. In the anelastic approximation the AM equation reads
where D^{2} is the vector Laplace operator in the azimuthal direction D^{2} = (∇^{2}−1/r^{2}sin^{2}θ). When the Coriolis and the nonlinear advection terms are neglected in Eq. (54), and when the differential rotation is sufficiently weak (i.e. ΔΩ/Ω_{0} ≪ 1, and thus the first part of the contraction term can be neglected), the steady balance is the linear equation
whose solution is
It simply corresponds to the Boussinesq linear viscous solution Eq. (49) weighted by the inverse of the background density profile.
In Fig. 12 we show the agreement between the analytical expression (56) and the numerical simulations. In all panels of this figure the profiles of the differential rotation found in the numerical calculations for various values of the parameters are plotted at a particular latitude (here θ = π/8) as a function of radius. We also overplot in dashed lines the analytical profile found by solving Eq. (55). In the left panel the linear regime (with Re_{c} = 10^{−2}) is represented at different density contrasts. In the other panels the density contrast between the inner and outer spheres is fixed to ρ_{i}/ρ_{0} = 20.9, but Re_{c} is increased. It is quite clear from the left panel that the agreement between the numerical and analytical solutions is perfect for the low Reynolds number cases. When the Reynolds number is increased, a departure from the analytical solution appears, as observed in the middle panel of Fig. 12 where the profile of Eq. (56), scaled with Re_{c}, is shown in dashed lines. In this case the level of differential rotation becomes comparable to Ω_{0} and the first part of the contraction term in Eq. (54) is no longer negligible, leading to a correction to Eq. (56). When this correction is applied, the agreement between the numerical solutions at higher Re_{c} and the new analytical expression (not given here) is recovered, as shown in the last panel of Fig. 12.
Fig. 12. Differential rotation normalised to the rotation rate of the outer sphere δΩ/Ω_{0}, rescaled with Re_{c} and plotted at the latitude θ = π/8 as a function of radius. In each case E = 10^{−4} and P_{r}(N_{0}/Ω_{0})^{2} = 10^{ 4}. Left panel: analytical solution Eq. (56) is represented by dashed lines and numerical solutions (runs 1.1–1.4 in Table 7) as solid coloured lines for the linear regime Re_{c} = 10^{−2}. The colours green, blue, purple, and red correspond respectively to the values of the density contrasts between the inner and outer spheres ρ_{i}/ρ_{0} = 1.75, 10.8, 20.9, and 60.2. Middle panel: numerical solutions are now represented for Re_{c} = 10^{−2}, 10^{−1}, 1, and 10 (red, purple, blue, and green, respectively) for ρ_{i}/ρ_{0} = 20.9 (runs 1.3, 2.3, 3.3, and 4.3 in Table 7). Right panel: same, but when Re_{c} ≥ 1 the solution in dashed lines is numerically estimated by taking the full contraction term in Eq. (55) into account. 
To summarise, we find that in the anelastic viscous regime the presence of a density stratification leads to a weaker differential rotation between the inner and outer spheres. In the linear regime, ΔΩ/Ω_{0} ≪ 1, the differential rotation is given by Eq. (56), i.e. the analytical solution previously derived in the Boussinesq case but weighted by the inverse of the background density.
7.3. Anelastic Eddington–Sweet regime
The role of the density stratification is now studied in the regime where the transport of AM is dominated by an Eddington–Sweet type meridional circulation. The density ratios ρ_{i}/ρ_{0} considered together with the other parameters in the simulations are listed in Table 8. The parameters E, P_{r}(N_{0}/Ω_{0})^{2}, and Re_{c} being defined in the same way and covering the same range as in the Boussinesq case (see Table 4), we are able to simply compare the anelastic simulations to the Boussinesq ones.
Parameters of the simulations in the Eddington–Sweet regime under the anelastic approximation.
The meridional cut displayed in the left panel of Fig. 13 shows, as in the Boussinesq case, that the differential rotation profile is neither radial nor cylindrical. However, a comparison of Fig. 13 to its Boussinesq counterpart in Fig. 5 reveals clear differences between the two simulations that only differ in their density stratification. Firstly, the amplitude of the rotation contrast taken between the inner and outer spheres is smaller in the anelastic case. Secondly, the distribution of the differential rotation within the spherical shell is smoother overall and predominantly radial; in particular, the region of high rotation rates localised at the inner sphere close to the equator that was clearly visible in the Boussinesq case is now absent. Quantitatively, along the inner sphere, the rotation rates of the equator and the pole only differ by 11%, while this contrast was 67% in the Boussinesq case. From the right panel of Fig. 13 it is obvious that the meridional flow in the anelastic case is also smoother and not focused towards the equator at the inner sphere. This comes as no surprise since the meridional flow is similar in the Eddington–Sweet and Taylor–Proudman regimes, and we have already found in the anelastic Taylor–Proudman case that the circulation is weighted by the inverse of the background density profile and is no longer affected by an equatorial boundary layer.
Fig. 13. Steady differential rotation (left panel) and norm of the total meridional velocity field with the associated streamlines in black (right panel), in the Eddington–Sweet regime obtained for E = 10^{−6}, P_{r}(N_{0}/Ω_{0})^{2} = 10^{−1}, Re_{c} = 10^{−2}, and ρ_{i}/ρ_{0} = 20.9 (run 10.3 in Table 8). In the left panel the coloured contours are normalised to the rotation rate of the outer sphere Ω_{0}. In the right panel the streamlines are as follows: outside the tangent cylinder the dashed lines describe an anticlockwise circulation and the solid lines a clockwise one; inside the tangent cylinder the dashed lines correspond to a downward meridional flow and the solid lines to an inward flow. 
We have computed the overall differential rotation between the inner and the outer spheres for different density ratios, and we propose a scaling relation to describe how it decreases with higher density ratios. As in Sect. 6.3 we estimate the dominant terms in the governing equations, but we pay attention to the effects of the density stratification. In particular, in determining the contraction timescale we take into account that the contraction velocity field is weighted by the inverse of the density. This allows us a more accurate determination of the contraction timescale:
Here the superscript ‘A’ stands for anelastic, and differentiates this contraction timescale from τ_{c} = r_{0}/V_{0} used in the Boussinesq case. The second expression is obtained using the definition of V_{f} (r) and assuming that r ∼ r_{0}.
Another specificity of the anelastic equations is the presence of an additional term in the thermal balance,
namely the first term in the brackets. However, in our simulations we found that this term is of the same order as the second term, so we can use only this second term to derive a scaling relation. This scaling relation takes the form
and provides a reasonable approximation for the differential rotation between the inner and outer spheres obtained in the anelastic Eddington–Sweet regime. However, as for the Boussinesq case, we found that the agreement is better when the contribution of the outer Ekman layer is taken away by subtracting the numerical Taylor–Proudman solution. This is done in Fig. 14 which shows that the global differential rotation scales approximately with the inverse of the integral of the density between the two spheres. It also shows that the rotation contrast is practically independent of the latitude except for the lowest density ratio ρ_{i}/ρ_{0} = 1.75.
Fig. 14. Differential rotation between the inner and outer spheres ΔΩ = Ω(r_{i})−Ω(r_{0}) normalised to the top value Ω_{0} obtained after subtracting the numerical Taylor–Proudman contribution, plotted as a function of the inverse of the integral of the background density profile normalised to the outer sphere value ρ_{0} for different latitudes: θ = π/2, π/4, π/8, π/16, and π/32 shown respectively in green, cyan, light blue, purple, and blue. The parameters are E = 10^{−5}, P_{r}(N_{0}/Ω_{0})^{2} = 10^{−1}, and Re_{c} = 10^{−2}. The symbols, circles, diamonds, squares, and triangles, respectively correspond to the different density contrasts between the inner and outer spheres ρ_{i}/ρ_{0} = 1.75, 10.8, 20.9, and 60.2 (runs 1.1–1.4 in Table 8). 
The latitudinal variation of the inner sphere rotation rate found at ρ_{i}/ρ_{0} = 1.75 is also found in the Boussinesq simulations and is a manifestation of the local rotation rate maximum found near the equator in these simulations. This local maximum seems to be related to the strong vertical jet which advects AM towards the equator, this jet being in turn associated with the presence of an equatorial boundary layer. Figure 15 shows the radial profiles of the gradient of the angular velocity together with the cylindrical radial and vertical velocity fields (respectively left, middle, and right panels) near the inner sphere at θ = π/4 and at the equator (respectively top and bottom rows), both for the Boussinesq (blue curve) and anelastic cases (black curves). As in the Taylor–Proudman regime, we observe that the strong equatorial boundary layer and the vertical jet that are present in the Boussinesq case nearly disappear in the anelastic simulations. Thus, contrary to the Boussinesq case, the meridional flow is no longer focused towards the equator. The different AM advections between the Boussinesq and anelastic cases thus appear to explain the smoother, more radial differential rotation observed in the anelastic case.
Fig. 15. Same as Fig. 10, but for the Eddington–Sweet regime with parameters E = 10^{−5}, P_{r}(N_{0}/Ω_{0})^{2} = 10^{−1}, and Re_{c} = 10^{−2}. The plain blue curve corresponds to the Boussinesq case (run 1.3 in Table 4), while the other curves are obtained by varying the density contrast between the inner and outer spheres in the same way as previously, but for runs 1.2–1.4 in Table 8 in the left panel, and for cases 1.1 and 1.3 in Table 8 in the middle and right panels. 
8. Summary and conclusions
In this paper we have presented the results of our extensive parametric study of an axisymmetric rotating spherical layer undergoing contraction, mimicking a stellar radiative zone during a rapid contraction phase. Our goal was to understand the motions induced by the contraction, in particular the amplitude and the spatial distribution of the differential rotation and of the meridional circulation. The contraction, modelled through an imposed massconserving radial velocity field V_{f}, produces an inward transport of AM that is balanced, in a stationary state, by an outward viscous or advective transport of AM. The full compressible gas dynamics has been approximated using either the Boussinesq or the anelastic equations.
We show that the parameter P_{r}(N_{0}/Ω_{0})^{2} controls the amplitude and the distribution of the differential rotation. When P_{r}(N_{0}/Ω_{0})^{2} ≫ 1, the buoyancy force inhibits the circulation and the viscosity dominates the AM transport. For a weak differential rotation ΔΩ/Ω ≪ 1, the amplitude of the differential rotation is proportional to the ratio between the viscous and the contraction timescales and its purely radial profile can be derived analytically. This regime is relevant in the core of contracting subgiants. If P_{r}(N_{0}/Ω_{0})^{2} ≪ 1, an Eddington–Sweettype circulation dominates the AM transport. Then, for a weak differential rotation and in the absence of density stratification, the amplitude of the differential rotation is approximately the ratio of the Eddington–Sweet timescale to the contraction timescale. In the anelastic simulations this amplitude decreases approximately as the inverse of the integral of the background density across the layer. The rotation is neither cylindrical nor radial, but when the density stratification is taken into account it tends to be predominantly radial and smoothly distributed. This regime holds for PMS stars and outside the degenerate core of subgiants.
Although less relevant for stars, a third regime of weak stratification was studied to understand the effects of the small though finite Ekman numbers in our simulations. It enabled us to understand the effects of the boundary layers formed at the inner and outer limits of the numerical domain. In particular, we found that in the Boussinesq simulations the presence of an equatorial boundary layer and of an associated vertical jet along the tangent cylinder produces a localised region of high rotation rates near the equator. This equatorial layer and the vertical jet are absent from the more realistic density stratified anelastic simulations leading to a smoother distribution of the differential rotation. This is at odds with the conclusions of Hypolite & Rieutord (2014) where such jets are invoked as a way to connect the core and the convective envelope of the contracting subgiants.
In this study we have been mostly concerned with the weak differential rotation regime, although the limits of this linear regime ΔΩ/Ω_{0} ∼ 1 have been exceeded in various cases. It would be interesting to further investigate the nonlinear regime, although this may be numerically challenging even for axisymmetric studies. In our setup the rotation of the frame of reference has been kept constant, while we expect it to increase as the star contracts. This effect could be introduced by adding a Euler term that describes the time evolution of rotating frame of reference (e.g., Hypolite & Rieutord 2014). However, in density stratified cases we expect that this temporal variation is slow compared to the contraction time in most parts of the domain.
We also want to note two important points about the anelastic approximation. Firstly, our formulation of the momentum equation Eq. (2) neglects the term; this is the LBR approximation (Lantz 1992; Braginsky & Roberts 1995; Lantz & Fan 1999). Physically, it means that the density scale height is smaller than the potential temperature scale height, which is approximately satisfied in the Earth’s atmosphere (Vallis 2017). This hypothesis is ensured in our simulations by keeping a low deviation from isentropy (ϵ_{s} ≪ 1). The legitimate question is, what is the validity of this approximation for a stellar interior? For a convective zone this assumption is obviously correct since . For a radiative interior a model of subgiant computed with the stellar evolution code MESA shows that outside the degenerate core the ratio of the potential temperature to density scale heights (h_{θ}/h_{ρ}) is only four during the evolution of the subgiant up to the red giant branch (RGB). Inside the core, the condition h_{ρ} ≪ h_{θ} can be satisfied at the terminal age main sequence (TAMS) and at the base of the RGB, but never during the subgiant phase where h_{θ} and h_{ρ} are of the same order. Nevertheless, using a soundproof approach other than the LBR approximation would lead to problems associated with the nonconservation of energy (Brown et al. 2012).
Secondly, in the thermal balance equation, we used the entropy diffusion term instead of the temperature diffusion because temperature diffusion leads to numerical difficulties with the code MagIC. In principle, these numerical difficulties could be relaxed by ensuring a strong coupling between the pressure and velocity fields, either by using staggered grids (Harlow & Welch 1965) or through an interpolation method (Rhie & Chow 1983).
In this work we present the Eddington–Sweet regime as a possible axisymmetric steady state in a contracting stellar radiative zone. This might seem at odds with previous studies that questioned the existence of an Eddington–Sweet meridional circulation (Busse 1981, 1982), but as we discuss now there is actually no contradiction between these studies and the present one. Busse (1981) showed that for an inviscid baroclinic stably stratified rotating fluid the only possible steady state is a purely zonal flow and that any meridional motions should be only transients. A stationary meridional flow would indeed produce an uncompensated advection in the AM equation. Consequently, if an Eddington–Sweet circulation initially existed, it would be modified on a rotation timescale (Busse 1982) because of the lack of AM conservation. Such a system would then evolve thermally so that the divergence of the heat flux balances the temperature advection leading, after an Eddington–Sweet timescale, to a steady state where the meridional circulation and the divergence of the heat flux are both zero. Our work does not disagree with these results since we introduce a contraction term in the AM equation which compensates for the advection by a steady meridional flow, thus ensuring AM conservation. We then get a meridional flow whose amplitude is of the order of the forcing amplitude and whose timescale corresponds to the Eddington–Sweet timescale.
The existence of the axisymmetric states described in this paper can also be debated against potential hydrodynamical instabilities. Axisymmetric instabilities like the centrifugal instability (Zahn 1974) or the GoldreichSchubertFricke instability (Goldreich & Schubert 1967) are unlikely because they should have been observed in our numerical simulations. However, the differentially rotating solutions might still be subject to nonaxisymmetric instabilities either barotropic instabilities driven by the radial (Lignières et al. 1999) or the horizontal shear (Deloncle et al. 2007), or to baroclinic instabilities that occur because of the noncoincidence of the isobaric and isentropic surfaces (Spruit & Knobloch 1984).
In addition to the hydrodynamical stability of the axisymmetric steady states described in this work, contracting radiative zones could also potentially trigger magnetohydrodynamical instabilities (Spruit 1999). The AM transport associated with the development of such instabilities tends to decrease the level of differential rotation imposed by the contraction such as required by subgiant seismic data (Deheuvels et al. 2012, 2014). For example, Fuller et al. (2019) propose a revised mechanism for the saturation of the Tayler instability that reproduces fairly well the core rotation rates along the red giant branch during the red clump and up to the white dwarf phase. The magnetorotational instability (MRI; Rüdiger et al. 2015; Jouve et al. 2020) could also be at play in contracting radiative zones.
In intermediatemass PMS stars, such instabilities have also been invoked as a way to explain the observed dichotomy between the magnetism of Ap/Bp stars and the weak complex magnetic fields discovered on bright A and Am stars (Lignieres et al. 2014; Gaurat et al. 2015).
Outside the tangent cylinder (i.e. for s > r_{i}) this analytical solution is similar to the expression derived by Hypolite & Rieutord (2014) up to a factor related to the aspect ratio.
Acknowledgments
The authors acknowledge the developers of MagIC (https://github.com/magicsph/magic) for the opensource code thanks to which the simulations were performed. This work was granted access to the HPC resources of CALMIP supercomputing center under the allocation P1118. The authors also wish to thank Sébastien Deheuvels for providing us with stellar parameters from the stellar evolution code MESA. L. J. acknowledges funding by the Institut Universitaire de France.
References
 Aerts, C., Mathis, S., & Rogers, T. M. 2019, ARA&A, 57, 35 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Alecian, E., Wade, G., Catala, C., et al. 2013a, MNRAS, 429, 1001 [NASA ADS] [CrossRef] [Google Scholar]
 Alecian, E., Wade, G., Catala, C., et al. 2013b, MNRAS, 429, 1027 [NASA ADS] [CrossRef] [Google Scholar]
 Barcilon, V., & Pedlosky, J. 1967, J. Fluid Mech., 29, 1 [Google Scholar]
 Beck, P. G., Montalban, J., Kallinger, T., et al. 2012, Nature, 481, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Böhm, T., & Catala, C. 1995, A&A, 301, 155 [Google Scholar]
 Bouvier, J., Bertout, C., Benz, W., & Mayor, M. 1986, A&A, 165, 110 [NASA ADS] [Google Scholar]
 Braginsky, S. I., & Roberts, P. H. 1995, Geophys. Astrophys. Fluid Dyn., 79, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, B. P., Vasil, G. M., & Zweibel, E. G. 2012, ApJ, 756, 109 [Google Scholar]
 Busse, F. 1981, Geophys. Astrophys. Fluid Dyn., 17, 215 [NASA ADS] [CrossRef] [Google Scholar]
 Busse, F. 1982, ApJ, 259, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Ceillier, T., Eggenberger, P., García, R., & Mathis, S. 2013, A&A, 555, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Charbonnel, C., & Talon, S. 2005, Science, 309, 2189 [Google Scholar]
 Clune, T. C., Elliott, J., Miesch, M., Toomre, J., & Glatzmaier, G. A. 1999, Parallel Comput., 25, 361 [CrossRef] [Google Scholar]
 Deheuvels, S., García, R., Chaplin, W., et al. 2012, ApJ, 756, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Deheuvels, S., Doğan, G., Goupil, M., et al. 2014, A&A, 564, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Deloncle, A., Chomaz, J.M., & Billant, P. 2007, J. Fluid Mech., 570, 297 [Google Scholar]
 den Hartogh, J., Eggenberger, P., & Deheuvels, S. 2020, A&A, 634, L16 [CrossRef] [EDP Sciences] [Google Scholar]
 Dietrich, W., & Wicht, J. 2018, Front. Earth Sci., 6, 189 [NASA ADS] [CrossRef] [Google Scholar]
 Eddington, A. 1925, Observatory, 48, 73 [Google Scholar]
 Eggenberger, P., Montalbán, J., & Miglio, A. 2012, A&A, 544, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Eggenberger, P., den Hartogh, J., Buldgen, G., et al. 2019, A&A, 631, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fuller, J., Piro, A. L., & Jermyn, A. S. 2019, MNRAS, 485, 3661 [NASA ADS] [Google Scholar]
 Gallet, F., & Bouvier, J. 2013, A&A, 556, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Garaud, P., Medrano, M., Brown, J., Mankovich, C., & Moore, K. 2015, ApJ, 808, 89 [Google Scholar]
 Gastine, T., & Wicht, J. 2012, Icarus, 219, 428 [NASA ADS] [CrossRef] [Google Scholar]
 Gaurat, M., Jouve, L., Lignières, F., & Gastine, T. 2015, A&A, 580, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Goldreich, P., & Schubert, G. 1967, ApJ, 150, 571 [Google Scholar]
 Harlow, F. H., & Welch, J. E. 1965, Phys. Fluids, 8, 2182 [Google Scholar]
 Hypolite, D., & Rieutord, M. 2014, A&A, 572, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jones, C., Boronski, P., Brun, A., et al. 2011, Icarus, 216, 120 [NASA ADS] [CrossRef] [Google Scholar]
 Jouve, L., Lignières, F., & Gaurat, M. 2020, A&A, 641, A13 [EDP Sciences] [Google Scholar]
 Lantz, S. R. 1992, PhD thesis, Dissertation Abstracts International, 5212, 6460, Cornell University, USA [Google Scholar]
 Lantz, S., & Fan, Y. 1999, ApJS, 121, 247 [Google Scholar]
 Lara, F. E., & Rieutord, M. 2013, A&A, 552, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lignières, F., Califano, F., & Mangeney, A. 1999, ArXiv eprints [arXiv:astroph/9908184] [Google Scholar]
 Lignieres, F., Petit, P., Auriere, M., Wade, G. A., & Böhm, T. 2014, Proc. Int. Astron. Union, 9, 338 [Google Scholar]
 Maeder, A. 2008, Physics Formation and Evolution of Rotating Stars (Springer Science& Business Media) [Google Scholar]
 Maeder, A., & Meynet, G. 2001, A&A, 373, 555 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marques, J., Goupil, M., Lebreton, Y., et al. 2013, A&A, 549, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pinçon, C., Belkacem, K., Goupil, M., & Marques, J. 2017, A&A, 605, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rhie, C., & Chow, W. L. 1983, AIAA J., 21, 1525 [NASA ADS] [CrossRef] [Google Scholar]
 Rieutord, M. 2006, A&A, 451, 1025 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rieutord, M., & Beth, A. 2014, A&A, 570, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rüdiger, G., Gellert, M., Spada, F., & Tereshin, I. 2015, A&A, 573, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schou, J., Antia, H., Basu, S., et al. 1998, ApJ, 505, 390 [NASA ADS] [CrossRef] [Google Scholar]
 Spruit, H. 1999, Arxiv eprints [arxiv:astroph/9907138] [Google Scholar]
 Spruit, H. 2002, A&A, 381, 923 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Spruit, H. C., & Knobloch, E. 1984, A&A, 132, 89 [NASA ADS] [Google Scholar]
 Sweet, P. 1950, MNRAS, 110, 548 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Talon, S., & Charbonnel, C. 2005, A&A, 440, 981 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Talon, S., & Charbonnel, C. 2008, A&A, 482, 597 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Talon, S., Kumar, P., & Zahn, J.P. 2002, ApJ, 574, L175 [Google Scholar]
 Vallis, G. K. 2017, Atmospheric and Oceanic Fluid Synamics (Cambridge University Press) [Google Scholar]
 Von Zeipel, H. 1924, MNRAS, 84, 665 [NASA ADS] [CrossRef] [Google Scholar]
 Wade, G. A., Drouin, D., Bagnulo, S., et al. 2005, A&A, 442, L31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wicht, J. 2002, Phys. Earth Planet. Inter., 132, 281 [NASA ADS] [CrossRef] [Google Scholar]
 Wicht, J., French, M., Stellmach, S., et al. 2018, Magnetic Fields in the Solar System, (Springer), 7 [Google Scholar]
 Zahn, J. P. 1974, SymposiumInternational Astronomical Union (Cambridge University Press), 59, 185 [NASA ADS] [CrossRef] [Google Scholar]
 Zahn, J.P. 1992, A&A, 265, 115 [NASA ADS] [Google Scholar]
Appendix A: Hydrostatic state forced by the contraction
In order to derive the different timescales of AM transport in Sect. 3, we subtract the hydrostatic state forced by the contraction velocity field from our equations.
We first decompose the temperature fluctuations as the sum of a spherical contribution and a deviation to the spherical symmetry,
where
and
We now rewrite Eq. (8) by using the decomposition of Eq. (A.1), where T′(r) is adimensionalised in the same way as and δΘ(r, θ) is nondimensionalised as Θ′(r, θ):
If R_{o} ≪ (N_{0}/Ω_{0})^{2}, the steady case of the foregoing equation reads
Using both Ostrogradsky’s theorem and the continuity equation, we can write
which then allows us to deduce
Furthermore,
According to Eq. (A.3) the first integral in brackets is zero. Since the last integral is also zero, we conclude the following:
Consequently, the latitudinal mean of Eq. (A.5) gives
The radial derivative of is given by solving and leads to
Finally, is derived by solving
using the boundary conditions . The analytical solution is readily
We thus define
and likewise
After subtracting the spherically symmetric temperature field induced by the contraction velocity field, the evolution equation of temperature fluctuations finally reads
while the momentum equation Eq. (7) becomes
in which we have subtracted the hydrostatic state induced by the radial contraction
with , and where is the contractioninduced hydrostatic pressure nondimensionalised with the hydrostatic pressure scale r_{0} ρ_{0} g_{0}.
Appendix B: Analytical solution for the differential rotation in the Taylor–Proudman regime (Boussinesq case)
In this appendix we seek an analytical solution for the differential rotation in the Taylor–Proudman regime under the Boussinesq approximation.
The linear steady inviscid AM equation Eq. (21) provides a balance between the Coriolis and contraction terms. This balance allows us to derive an analytical solution for the cylindrical radial velocity field
where and . In cylindrical coordinates the steady axisymmetric velocity field reads
where is the dimensionless stream function. By using Eq. (B.1) we can then derive an analytical solution for ,
while the use of Eq. (B.2) allows us to determine the vertical velocity field:
By using the classical relations between the cylindrical and spherical coordinates, we deduce
The noslip condition at the inner sphere leads to
By symmetry, at the equator, which means that when , . This allows us to find the analytical solutions of the meridional velocity fields:
By integrating Eq. (B.7), we deduce the stream function (up to a constant):
Since the stream function is zero at the poles, we deduce that cte = 0.
The mass conservation leads to the relationship
which gives in spherical coordinates:
From the analysis of the Ekman layer located at the outer sphere, we obtain the nondimensional form of the Ekman condition:
We now rewrite Eq. (B.13) in cylindrical coordinates using the Taylor–Proudman constraint and by noting that , which leads to at the outer sphere:
In addition,
thus, with the following substitutions and dz = cos θ dr, the relationship Eq. (B.12) becomes
The choice of the closed surface is a cylinder of radius s = r sin θ between and (see Fig. B.1).
By noting the dimensionless outward mass flux crossing this cylinder, we have
Likewise, with the inward dimensionless mass flux crossing the cylinder, we obtain
Finally, from Eq. (B.16) we get
from which we conclude
which can then be rewritten under dimensional form to lead to the solution Eq. (45).
Fig. B.1. Sketch of a mass–flux balance inside a cylinder of radius s = r sin θ between and . 
Appendix C: Analytical solution for the differential rotation in the Taylor–Proudman regime (anelastic case)
We now show that the analytical solution previously derived, namely Eq. (B.20), is still valid independently of the amplitude of the density contrast applied between the inner and outer spheres.
In the anelastic approximation, the dimensionless balance between the Coriolis and contraction terms now provides
The poloidaltoroidal decomposition now relates the stream function to the meridional velocity fields as
which is equivalent to
From these relationships, it is obvious that all the analytical solutions derived for the meridional fields will be identical to the Boussinesq case except that they will now be weighted by the inverse of the background density profile:
Instead, the previous analytical solution derived for the stream function is unmodified:
In the anelastic approximation, Eq. (B.11) is rewritten as follows:
Since the Ekman condition Eq. (B.13) is given in , and the background density profile at this location is such that , we obtain the following relationship:
In other words, this is the same dimensional analytical solution for the differential rotation as in the Boussinesq case, namely Eq. (B.20) (or under dimensional form Eq. (45)).
Appendix D: Solution in the Boussinesq Eddington–Sweet regime
In this appendix we show, in the Eddington–Sweet regime, that a solution for the differential rotation can be expressed as the sum of an Eddington–Sweettype solution plus a Taylor–Proudman solution. The second term comes from the Ekman boundary layer at the outer sphere and cannot be entirely neglected near this boundary even at very small Ekman numbers.
The dimensionless linear axisymmetric steady equations of momentum, temperature, and continuity read
For small Ekman numbers, solutions are sought as the sum of an interior inviscid solution and a boundary layer solution, that is where by construction the boundary layer solution must vanish in the interior. For an E^{1/2} wide boundary layer at the outer sphere, a dimensionless stretched coordinate is introduced and the previous condition reads .
For the interior solution we assume the expansion
where E^{−1/2} ≪ P_{r}(N_{0}/Ω_{0})^{2}/E. Substitution of this expansion in Sys. (D.1) yields at the dominant order the following system:
From the expression of , the radial velocity field can be determined using mass conservation and the boundary condition . Then we assume that the system of equations for the remaining unknowns and admits a solution for the full boundary conditions , in , and .
However, the meridional flow does not satisfy the noslip condition at the outer sphere. An Ekman layer is therefore needed to correct the 𝒪(1) jump on this field and this will induce the term. Inserting the expansion Sys. (D.2) in Sys. (D.1), we see at 𝒪(E^{−1/2}) order that such a term must satisfy
Meanwhile, within the Ekman layer, the boundary layer solution must be of the following order to compensate the 𝒪(1) meridional velocity jump:
By rewriting Sys. (D.1) with the stretched coordinate ξ and by inserting the previous expressions, the boundary layer equations read as follows:
The boundary condition leads to at the order E^{−1/2}. Together with the other boundary condition and the evanescent condition when , the first two equations are solved to give
Then, integrating the last equation of Sys. (D.6) and using both the evanescent condition and the boundary condition leads to the Ekman pumping condition Eq. (B.13) in its dimensionless form. From this condition on the interior flow, one can show, as in Appendix B, that is the Taylor–Proudman solution. In the interior, its contribution to the total differential rotation decreases as increases. However, near the outer sphere its contribution can never be neglected because the Eddington–Sweet term vanishes there .
All Tables
Density contrasts between the inner and outer spheres investigated in the anelastic approximation, and their associated parameters.
Parameters of the simulations in the Eddington–Sweet regime under the anelastic approximation.
All Figures
Fig. 1. Steady axisymmetric flow in the Taylor–Proudman regime. In the left panel, the coloured contours represent the differential rotation δΩ(r, θ) = Ω(r, θ)−Ω_{0} normalised to the rotation rate of the outer sphere Ω_{0} while in the right panel, we show the norm of the total meridional velocity field . The parameters of the simulation are E = 10^{−4}, P_{r}(N_{0}/Ω_{0})^{2} = 10^{−4} and Re_{c} = 10^{−2} (simulation 1.3 in Table 2). The black lines show the streamlines of the total meridional circulation i.e. taking the effect of the contraction into account. Outside the tangent cylinder, the dashed lines correspond to an anticlockwise circulation and the solid lines to a clockwise one while inside the tangent cylinder, the dashed lines represent a downward circulation and the full ones, an upward circulation. 

In the text 
Fig. 2. Axisymmetric steady solutions of the differential rotation as a function of radius at the latitude θ = π/4 in the Taylor–Proudman regime. Left panel: analytical solution Eq. (45) is plotted in black dashed lines and compared with the numerical simulations obtained for fixed P_{r}(N_{0}/Ω_{0})^{2} parameter (10^{−4}) and contraction Reynolds number Re_{c} (10^{−2}) with various Ekman numbers: 10^{−2} (in grey), 10^{−3} (in cyan), 10^{−4} (in light blue), 10^{−5} (in dark blue), and 10^{−6} (in purple) respectively for runs 1.1–1.5 in Table 2. For E = 10^{−5} and E = 10^{−6}, additional contraction Reynolds numbers are studied, namely 10^{−1} (runs 2.1 and 3.1 in Table 2), 1 (runs 2.2 and 3.2 in Table 2), and 10 (runs 2.3 and 3.3 in Table 2). All curves are rescaled by . Right panel: differential rotation between the inner and outer spheres ΔΩ = Ω(r_{i})−Ω(r_{0}) at the same latitude and normalised with the rotation rate taken at the outer sphere Ω_{0}, now plotted as a function of . The previous numerical solutions obtained for E = 10^{−5} (runs 1.4 and 2.1–2.3 in Table 2) are plotted in dark blue, while the simulations performed for E = 10^{−6} (runs 1.5 and 3.1–3.3 in Table 2) are plotted in purple. 

In the text 
Fig. 3. Differential rotation normalised to the rotation rate of the outer sphere δΩ(r, θ)/Ω_{0} in the viscous regime (first and third panels) and norm of the total meridional velocity field (second and fourth panels) with the streamlines of the associated total meridional circulation (in black). In the two left panels the parameters of the simulation are P_{r}(N_{0}/Ω_{0})^{2} = 10^{ 4}, Re_{c} = 10^{−2}, and E = 10^{−5} (run 3.4 in Table 3), and in the two right panels the parameters are P_{r}(N_{0}/Ω_{0})^{2} = 10^{4}, Re_{c} = 10, and E = 10^{−4} (run 4.3 in Table 3). In the second and fourth panels the black solid lines represent the streamlines of an ascending and inwardly directed circulation, while the black dashed lines correspond to a descending and inwardly directed one. 

In the text 
Fig. 4. Differential rotation in the viscous regime scaled by Re_{c} as a function of radius at the latitude θ = π/8. The analytical solution Eq. (48) (black dashed lines) is compared with the numerical results. Low Reynolds number runs (1.1–4.1 in Table 3, with Re_{c} = 10^{−2} − 10^{−1}) are plotted with violet dotted lines and are all very close to the low Reynolds number version of the analytical solution, Eq. (49). Simulations at Re_{c} = 1 (runs 4.2 and 5.1 in Table 3) and Re_{c} = 10 (run 4.3 Table 3) are also in agreement with the analytical solution. 

In the text 
Fig. 5. Steady differential rotation (left panel) normalised to the top value, and norm and streamlines of the total meridional velocity field (right panel). The parameters of the simulation (run 1.4 in Table 4) are E = 10^{−6}, P_{r}(N_{0}/Ω_{0})^{2} = 10^{−1}, and Re_{c} = 10^{−2}. Outside the tangent cylinder the black dashed lines represent the streamlines of an anticlockwise circulation, and the solid lines those of a clockwise circulation. Inside the tangent cylinder the dashed lines simply correspond to the streamlines of a downward circulation, while the solid lines are for an upward circulation. 

In the text 
Fig. 6. Differential rotation between the inner equator and the outer sphere normalised to the outer value of the rotation rate Ω_{0} in the Eddington–Sweet regime, plotted as a function of P_{r}(N_{0}/Ω_{0})^{2}R_{o}/E. The Rossby number R_{o} is varied to encompass the linear and nonlinear regimes. The Ekman number is fixed to E = 10^{−6}. The green curve corresponds to simulations with P_{r}(N_{0}/Ω_{0})^{2} = 10^{−1} (runs 1.4, 3.1, 4.1, 5.1, and 6.1 in Table 4) and the blue curve to simulations with P_{r}(N_{0}/Ω_{0})^{2} = 10^{−2} (runs 2.2, 3.2, 4.2, and 5.2 in the same table). 

In the text 
Fig. 7. Differential rotation in the Eddington–Sweet linear regime for P_{r}(N_{0}/Ω_{0})^{2} = 10^{−1} (top panel) and 10^{−2} (bottom panel), the other two parameters E = 10^{−6} and Re_{c} = 10^{−2} being identical (runs 1.4 and 2.2 in Table 4). The coloured contours are normalised to the rotation rate of the outer sphere Ω_{0}. Left column: differential rotation δΩ(r, θ)/Ω_{0} rescaled by E (P_{r}(N_{0}/Ω_{0})^{2}R_{o})^{−1}. Middle column: differential rotation obtained after subtracting the unstratified Taylor–Proudman solution given by Eq. (45) and rescaling by E (P_{r}(N_{0}/Ω_{0})^{2}R_{o})^{−1} (see Eq. (53)). Right column: ratio . 

In the text 
Fig. 8. Contraction radial velocity field V_{f}, Eq. (19), normalised to the top value V_{0} and plotted as a function of radius r. The lines (from light to dark) show the weakest to the strongest density contrasts between the inner and outer spheres, which are respectively ρ_{i}/ρ_{0} = 1.75, 10.8, 20.9, and 60.2. The contraction radial velocity field in the Boussinesq case given by Eq. (20) (i.e. corresponding to a density ratio of 1) is plotted as a dashed line as a reference. 

In the text 
Fig. 9. Meridional cuts showing the coloured contours of the steady differential rotation normalised to the top value in the left panel, and the norm of the total meridional velocity field with its associated streamlines (in black) in the right panel, in the Taylor–Proudman regime. As previously, outside the tangent cylinder the solid lines describe a clockwise circulation, and the dashed lines an anticlockwise circulation; inside the tangent cylinder the solid lines represent an upward meridional flow and the dashed lines a downward flow. The parameters of the simulation are E = 10^{−4}, P_{r}(N_{0}/Ω_{0})^{2} = 10^{−4}, Re_{c} = 10^{−2}, and ρ_{i}/ρ_{0} = 20.9 (run 1.1 in Table 6). 

In the text 
Fig. 10. Radial partial derivative of the ratio of the azimuthal velocity field U_{ϕ} to the radius r (left panel), cylindrical radial velocity field U_{s} = U_{r} sin θ + U_{θ} cos θ (middle panel), and vertical velocity field U_{z} = U_{r} cos θ − U_{θ} sin θ (right panel) as a function of radius. In each panel the top row corresponds to the latitude θ = π/4. In the left and middle panels the bottom row corresponds to the equator, but since the vertical velocity field U_{z} is zero at the equator, the bottom right panel is obtained close to this location, at about 85degree angle with the vertical axis. In each case E = 10^{−5}, P_{r}(N_{0}/Ω_{0})^{2} = 10^{−4} and Re_{c} = 10^{−2}. In all panels the plain blue curve corresponds to the Boussinesq case (run 1.4 in Table 2), while the other curves are obtained by varying the density contrast between the inner and outer spheres. Shown in the left panel are ρ_{i}/ρ_{0} = 10.8 (dashed lines), 20.9 (dashdotted lines), and 60.2 (dotted lines), corresponding to runs 2.2 to 2.4 in Table 6; in the middle and right panels the two black curves are respectively obtained for ρ_{i}/ρ_{0} = 1.75 (plain curve) and ρ_{i}/ρ_{0} = 20.9 (dashed lines), corresponding to runs 2.1 and 2.3 in Table 6. 

In the text 
Fig. 11. Stationary differential rotation δΩ(r, θ)/Ω_{0} (first and third panels) and norm of the total meridional velocity field (second and fourth panels) with its associated streamlines in black, in the viscous anelastic regime. The parameters of the simulation shown in the two left panels are E = 10^{−5}, P_{r}(N_{0}/Ω_{0})^{2} = 10^{ 4}, Re_{c} = 10^{−2}, and ρ_{i}/ρ_{0} = 20.9 (run 5.2 in Table 7), and in the two right panelsE = 10^{−4}, P_{r}(N_{0}/Ω_{0})^{2} = 10^{ 4}, Re_{c} = 10, and ρ_{i}/ρ_{0} = 20.9 (run 4.3 in Table 7). In the first and third panels, the coloured contours are normalised to the top value Ω_{0}. 

In the text 
Fig. 12. Differential rotation normalised to the rotation rate of the outer sphere δΩ/Ω_{0}, rescaled with Re_{c} and plotted at the latitude θ = π/8 as a function of radius. In each case E = 10^{−4} and P_{r}(N_{0}/Ω_{0})^{2} = 10^{ 4}. Left panel: analytical solution Eq. (56) is represented by dashed lines and numerical solutions (runs 1.1–1.4 in Table 7) as solid coloured lines for the linear regime Re_{c} = 10^{−2}. The colours green, blue, purple, and red correspond respectively to the values of the density contrasts between the inner and outer spheres ρ_{i}/ρ_{0} = 1.75, 10.8, 20.9, and 60.2. Middle panel: numerical solutions are now represented for Re_{c} = 10^{−2}, 10^{−1}, 1, and 10 (red, purple, blue, and green, respectively) for ρ_{i}/ρ_{0} = 20.9 (runs 1.3, 2.3, 3.3, and 4.3 in Table 7). Right panel: same, but when Re_{c} ≥ 1 the solution in dashed lines is numerically estimated by taking the full contraction term in Eq. (55) into account. 

In the text 
Fig. 13. Steady differential rotation (left panel) and norm of the total meridional velocity field with the associated streamlines in black (right panel), in the Eddington–Sweet regime obtained for E = 10^{−6}, P_{r}(N_{0}/Ω_{0})^{2} = 10^{−1}, Re_{c} = 10^{−2}, and ρ_{i}/ρ_{0} = 20.9 (run 10.3 in Table 8). In the left panel the coloured contours are normalised to the rotation rate of the outer sphere Ω_{0}. In the right panel the streamlines are as follows: outside the tangent cylinder the dashed lines describe an anticlockwise circulation and the solid lines a clockwise one; inside the tangent cylinder the dashed lines correspond to a downward meridional flow and the solid lines to an inward flow. 

In the text 
Fig. 14. Differential rotation between the inner and outer spheres ΔΩ = Ω(r_{i})−Ω(r_{0}) normalised to the top value Ω_{0} obtained after subtracting the numerical Taylor–Proudman contribution, plotted as a function of the inverse of the integral of the background density profile normalised to the outer sphere value ρ_{0} for different latitudes: θ = π/2, π/4, π/8, π/16, and π/32 shown respectively in green, cyan, light blue, purple, and blue. The parameters are E = 10^{−5}, P_{r}(N_{0}/Ω_{0})^{2} = 10^{−1}, and Re_{c} = 10^{−2}. The symbols, circles, diamonds, squares, and triangles, respectively correspond to the different density contrasts between the inner and outer spheres ρ_{i}/ρ_{0} = 1.75, 10.8, 20.9, and 60.2 (runs 1.1–1.4 in Table 8). 

In the text 
Fig. 15. Same as Fig. 10, but for the Eddington–Sweet regime with parameters E = 10^{−5}, P_{r}(N_{0}/Ω_{0})^{2} = 10^{−1}, and Re_{c} = 10^{−2}. The plain blue curve corresponds to the Boussinesq case (run 1.3 in Table 4), while the other curves are obtained by varying the density contrast between the inner and outer spheres in the same way as previously, but for runs 1.2–1.4 in Table 8 in the left panel, and for cases 1.1 and 1.3 in Table 8 in the middle and right panels. 

In the text 
Fig. B.1. Sketch of a mass–flux balance inside a cylinder of radius s = r sin θ between and . 

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.