Axisymmetric investigation of differential rotation in contracting stellar radiative zones

Context. Stars experience rapid contraction or expansion at di ﬀ erent 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 di ﬀ erential rotation and meridional circulation. Methods. We consider a rotating spherical layer crossed by an imposed radial velocity ﬁeld that mimics the contraction, and numeri- cally 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 stratiﬁcation, and density stratiﬁcation that are relevant for stars. Results. The di ﬀ erential rotation and the meridional circulation result from a competition between the contraction-driven inward transport of angular momentum and an outward transport dominated by either viscosity or an Eddington–Sweet-type 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 stratiﬁcation into account is important to study more realistic radial contraction ﬁelds, and also because the resulting ﬂow is less a ﬀ ected by unwanted e ﬀ ects of the boundary conditions. In these di ﬀ erent regimes and for a weak di ﬀ erential rotation we derive scaling laws that relate the amplitude of the di ﬀ erential rotation to the contraction timescale.


Introduction
The transport of angular momentum (AM) and chemical elements in stars strongly affects their evolution, from pre-main 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 self-consistent meridional flow and models of the turbulent transport driven by hydrodynamical instabilities. The importance of additional processes such as internal waves (e.g.,  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 blue-to-red 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 inferred through helioseismology (Schou et al. 1998), as well as the lithium dip in solar-like stars . Despite these early encouraging results many stellar observations remain unexplained, especially the internal rotation rates revealed by asteroseismology in evolved stars and in intermediate-mass 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 large-scale 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.
Solar-like 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 break-up 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 break-up 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 break-up 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 solar-type 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 break-up velocity on the ZAMS while a 3 M star should rotate at around 65% of its break-up velocity and a 5 M at around 52% break-up. Again, observations indicate lower surface rotation rates on the ZAMS, corresponding to 30% break-up 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 wind-driven 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 contraction-driven 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, up-to-date 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 large-scale 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 contraction-driven 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.

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 sound-proof 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.

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 non-dimensional 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 non-dimensionalised 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 ∆T and ∆S between the two spheres. The pressure fluctuations are non-dimensionalised 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 mass-conserving contraction velocity field whose expression is given in Sect. 2.3. The dimensionless momentum equation 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σ i j is the dimensionless stress tensor andQ ν 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 E = ν/Ω 0 r 2 0 , 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 N 2 0 /Ω 2 0 . Instead, only the last four are necessary in the Boussinesq approximation.

Background state
In the anelastic approximation we use a non-adiabatic 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 dS /dr 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 dS /dr, 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 T is the solution of the conduction equation for prescribed temperatures at the spherical boundaries: Here ∆T = T (r = r 0 ) − T (r = r i ) > 0, again ensuring stable stratification. The background temperature gradient dT /dr 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

Contraction velocity field
The specificity of our model is to introduce a contraction process that will force large-scale 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 ∇ · ρ V f = 0. 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.

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 non-linear 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 τ L c = (r 0 /V 0 )(∆Ω/Ω 0 ), 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 well-known types of circulation, the Eddington-Sweet and the Ekman circulations, that will be relevant to analyse our simulation results. In this paper the steady-state 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, large-scale 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-Sweet-type 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 T m (r) = T (r) + T (r) 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 large-scale 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 no-slip 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 well-known 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 U r ∼ √ EL∆Ω J 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 spin-up 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 large-scale 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 non-negligible 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, √ E 1, 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 non-linear regime, ∆Ω/Ω 0 1.
A non-linear 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 set-up with r i = 0.3 r 0 , the maximum value of the differential rotation is which is clearly in the non-linear 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.

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.

Pre-main sequence stars
Among PMS stars, low-mass stars are characterised by a bimodal distribution of their rotational period (P rot ≈ 8 days for those experiencing a disk-locking phase and P rot ≈ 1 day for the others; Maeder 2008), while intermediate-and high-mass 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 N 2 T = 10 −7.5 s −2 at the beginning of the PMS to N 2 T = 10 −5.8 s −2 at the ZAMS. The Prandtl number is very low in non-degenerate 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.
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 τ κ = r 2 0 /κ. In this case Re c ≈ P −1 r , 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 solar-mass star, the radius ranges from 2.6 R during the Table 1. Estimates of P r (N 0 /Ω 0 ) 2 in PMS stars.

Status
Period 1 day 8 days N 2 T = 10 −7.5 s −2 5.98 · 10 −6 3.83 · 10 −4 (Beginning of the PMS) N 2 T = 10 −5.8 s −2 3.00 · 10 −4 1.91 · 10 −2 (End of the PMS) 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 non-linear Eddington-Sweet regime defined by

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 N 2 T = 10 −4.5 s −2 at the beginning of the subgiant phase to N 2 T = 10 −3 s −2 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 non-linear viscous regime is the relevant one inside the degenerate core of the subgiants, whereas outside the degenerate core the non-linear Eddington-Sweet regime dominates:

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 pseudo-spectral 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 1 https://github.com/magic-sph/magic 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 stress-free 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 no-slip 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.

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 non-linear 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.

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, A109, page 6 of 22 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 U tot = U r + V f (r) e r +U θ e θ . 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. Simulation E P r (N 0 /Ω 0 ) 2 Re c 1.1 10 −2 10 −4 10 −2 1.2 10 −3 10 −4 10 −2 1.3 10 −4 10 −4 10 −2 1.4 10 −5 10 −4 10 −2 1.5 10 −6 10 −4 10 −2 2.1 10 −5 10 −4 10 −1 2.2 10 −5 10 −4 1 2.3 10 −5 10 −4 10 3.1 10 −6 10 −4 10 −1 3.2 10 −6 10 −4 1 3.3 10 −6 10 −4 10 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 U tot = U r + V f (r) e r + U θ e θ . 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.
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: 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 R o / √ E 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 Finally, the existence of the Taylor-Proudman regime in a stably stratified atmosphere where P r (N 0 /Ω 0 ) 2 √ E can  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 √ E/R o . 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 R o / √ E. 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.
be explained as follows. From Eq. (44) we see that U r ≈ V 0 . Thus, according to the thermal balance Eq. (25) The ratio of these two quantities is U φ /δΘ ≈ √ E/P r (N 0 /Ω 0 ) 2 . As a consequence, in the regime P r (N 0 /Ω 0 ) 2 √ E, 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.

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.
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.
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.
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.

Eddington-Sweet regime
When the stratification is not too strong, √ E P r (N 0 /Ω 0 ) 2 1, 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.
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 .
This figure first shows that a linear regime exists as long as ∆Ω/Ω 0 is low enough. The differential rotation is clearly lin- 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 non-linear 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). early dependent on R o until ∆Ω/Ω 0 10 −1 , where it starts to deviate. The differential rotation then falls below the linear relation because the non-linear 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 ∆Ω/Ω 0 ∝ P r (N 0 /Ω 0 ) 2 R o /E. 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 by-product of the Ekman layer induced by the mismatch between the inviscid contraction-driven 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 O 1/ √ E jump in differential rotation across the O √ E 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 √ E/P r (N 0 /Ω 0 ) 2 , 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).
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 ∝1/ √ E 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 √ E/P r (N 0 /Ω 0 ) 2 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.
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 √ E/P r (N 0 /Ω 0 ) 2 (such as runs 1.3 and 2.3 in Table 4) have similar scaled differential rotation δΩ ED /Ω 0 E P r (N 0 /Ω 0 ) 2 . While the Eddington-Sweet scaling and its O √ E/P r (N 0 /Ω 0 ) 2 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 O 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 non-dimensional Table 5. Density contrasts between the inner and outer spheres investigated in the anelastic approximation, and their associated parameters.  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. numbers might be tricky and is beyond the scope of the present study.

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)). 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.

Anelastic Taylor-Proudman regime
We first study the effects of the density stratification in the regime P r (N 0 /Ω 0 ) 2 √ E 1 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 Table 6. Parameters of the simulations in the anelastic Taylor-Proudman regime.
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. 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 non-dimensional 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.
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 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 85-degree 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 (dash-dotted 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. 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.

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. 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 Fig. 11. Stationary differential rotation δΩ(r, θ)/Ω 0 (first and third panels) and norm of the total meridional velocity field U tot = U r + V f (r) e r + U θ e θ (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 panels E = 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 .  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.
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.
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 non-linear 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 V 0 r 2 0 ρ 0 /ρr 3 ∂(rU φ )/∂r 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.
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.

Anelastic Eddington-Sweet regime
The role of the density stratification is now studied in the regime √ E P r (N 0 /Ω 0 ) 2 1 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.
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.
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. 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.

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 mass-conserving 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-Sweet-type 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 P r (N 0 /Ω 0 ) 2 ≤ √ E was studied to understand the effects A109, page 16 of 22 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.
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 non-linear regime, although this may be numerically challenging even for axisymmetric studies. In our set-up 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 P /ρ C p dS /dr 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 dS /dr ≈ 0. 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 sound-proof 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(Busse , 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 Goldreich-Schubert-Fricke 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 non-coincidence 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 magneto-hydrodynamical 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(Deheuvels et al. , 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 magneto-rotational instability (MRI; Rüdiger et al. 2015;Jouve et al. 2020) could also be at play in contracting radiative zones.
In intermediate-mass 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).

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, Θ (r, θ) = T (r) + δΘ(r, θ), where T (r) is adimensionalised in the same way as T (r), and δΘ(r, θ) is non-dimensionalised as Θ (r, θ): gives (A.10) The radial derivative ofT is given by solving ∇ 2T = 0 and leads to dT dr = −r ĩ r 2 (r i − 1) . (A.11) Finally,T (r) is derived by solving d 2T dr 2 + 1 r Pe c r + 2 dT dr = Pe crĩ r 4 (r i − 1) (A.12) using the boundary conditionsT (1) =T (r i ) = 0. The analytical solution is readilỹ We thus definẽ with δΠ =Π −P c , and whereP c is the contraction-induced hydrostatic pressure non-dimensionalised with the hydrostatic pressure scale r 0 ρ 0 g 0 .