© B. Gouhier et al. 2021

1. 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 large-scale 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., 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 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 Charbonnel & Talon 2005) inferred through helioseismology (Schou et al. 1998), as well as the lithium dip in solar-like 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 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⁵³ g⋅cm² 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² 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.

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

2.1. Governing equations

To model a stellar radiative zone we consider a spherical shell of inner radius r_i and outer radius r₀ filled with a self-gravitating 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₀ as the reference length, the value of the contraction velocity field at the outer sphere V₀ as a characteristic velocity, and the contraction timescale τ_c = r₀/V₀ as the reference timescale. The frame rotates at Ω₀, 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 ρ₀ and temperature T₀, 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 non-dimensionalised by ρ₀r₀Ω₀V₀ and the entropy fluctuations by C_pΩ₀V₀/g₀, where C_p is the heat capacity and g₀ 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)

(1)

(2)

(3)

where V_f(r) is the mass-conserving contraction velocity field whose expression is given in Sect. 2.3. The dimensionless momentum equation Eq. (2) was derived using the Lantz-Braginsky-Roberts (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:

(4)

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₀/Ω₀r₀, E is the Ekman number , P_r is the Prandtl number P_r = ν/κ, and N₀ 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:

(5)

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 Ω₀V₀T₀/g₀ as the temperature scale. The dimensionless governing equations for continuity, momentum evolution, and temperature evolution become

(6)

(7)

(8)

In Eq. (7) we use a gravity profile corresponding to a uniform fluid density ρ_c:

(9)

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 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 is related to the Brunt–Väisälä frequency through

(10)

from which we define a reference Brunt–Väisälä frequency:

(11)

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:

(12)

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

(13)

as well as the background density profile

(14)

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

(15)

for prescribed temperatures at the spherical boundaries:

(16)

Here , again ensuring stable stratification. The background temperature gradient is then related to the Brunt–Väisälä frequency through the relationship

(17)

from which the reference Brunt–Väisälä frequency is defined as

(18)

2.3. 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 . In the anelastic case the contraction velocity field thus simply reads:

(19)

Under the Boussinesq approximation, since the background density is uniform, this velocity field simplifies into

(20)

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

(21)

where D² = (∇²−1/r²sin²θ) 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₀/V₀, when ΔΩ/Ω₀ ≳ 1, or its linear version , when ΔΩ/Ω₀ ≪ 1 so that rU_ϕ ≪ r² sin θ Ω₀.

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:

(22)

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

(23)

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:

(24)

We then use the steady linear dimensional form of the thermal balance Eq. (A.16)

(25)

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

(26)

(27)

from which we get the following relationship between the radial velocity field and the level of differential rotation:

(28)

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:

(29)

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,

(30)

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

(31)

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,

(32)

are controlled by two dimensionless numbers, the Ekman number E and the P_r(N₀/Ω₀)² parameter. As we are interested in cases with low Ekman numbers, , we can thus define three regimes of interest, namely

(33)

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:

(34)

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, ΔΩ/Ω₀ ≪ 1, and the beginning of the non-linear regime, ΔΩ/Ω₀ ≳ 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

(35)

which can be solved to obtain

(36)

Thus, for our set-up with r_i = 0.3 r₀, the maximum value of the differential rotation is

(37)

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 ΔΩ/Ω₀ ≪ 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₀/Ω₀)² parameter, and the contraction Reynolds number Re_c in stars undergoing contraction.

4.1. 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 at the beginning of the PMS to at the ZAMS. The Prandtl number is very low in non-degenerate stellar interiors, typically P_r ≈ 10⁻⁶ (Garaud et al. 2015). We thus obtain the orders of magnitude of P_r(N₀/Ω₀)² 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 . In this case , thus τ_ED/τ_c = Re_c P_r(N₀/Ω₀)² ≈ (N₀/Ω₀)² 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¹ cm² ⋅ s⁻¹ in the core to ν ≈ 5 ⋅ 10² cm² ⋅ s⁻¹ 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³ cm² ⋅ s⁻¹). 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 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⁻¹⁵ to E = 10⁻¹⁷. We thus conclude that the dominant parameter regime for PMS stars is the non-linear Eddington–Sweet regime defined by

(38)

4.2. Subgiants

In the degenerate cores of subgiants, the Prandtl number increases to ∼10⁻³ 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⁻⁶ rad ⋅ s⁻¹, we find that the degenerate core is characterised by very large values of the parameter P_r(N₀/Ω₀)² up to ∼10⁵.

Outside the degenerate core, typical P_r(N₀/Ω₀)² are 10⁻³, 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,

(39)

whereas outside the degenerate core the non-linear Eddington–Sweet regime dominates:

(40)

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 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¹ 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₀ = 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:

(41)

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:

(42)

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) = Ω₀ in the Boussinesq case and a fixed level of differential rotation

(43)

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⁻⁷), 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 non-negligible. 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⁻⁷ case to 12 grid cells in the E = 10⁻⁵ 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₀/Ω₀)² 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.

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⁻⁴, P_r(N₀/Ω₀)² = 10⁻⁴, Re_c = 10⁻²). 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,θ)−Ω₀)/Ω₀ 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.

In order to characterise this regime we performed several simulations for a range of Ekman and contraction Reynolds numbers, E = 10⁻² − 10⁻⁶, Re_c = 10⁻² − 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:

(44)

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²:

(45)

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

(46)

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₀. Thus, according to the thermal balance Eq. (25), δΘ ≈ P_r(N₀/Ω₀)² E⁻¹ V₀, while according to Eq. (30), U_ϕ ≈ E^−1/2 V₀. 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₀/Ω₀)². To investigate this regime, numerical simulations have been performed for values of P_r(N₀/Ω₀)² varying between 10² and 10⁴. The different runs characterised by the three parameters, P_r(N₀/Ω₀)², 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₀/Ω₀)² = 10⁴ 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⁻² for Re_c = 10⁻² 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.

The results of the numerical simulations are now compared with the analytical solution that can be derived by neglecting the Coriolis and the non-linear advection terms in the AM conservation Eq. (21). The steady balance between the contraction and viscous terms indeed reads

(47)

whose solution is given by the following expression:

(48)

If the differential rotation δΩ_ν(r)/Ω₀ is less than one, the forcing term in Eq. (47) is reduced to , and the analytical solution for the differential rotation simplifies into

(49)

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₀/Ω₀)² ≥ 10² considered. We also find that the linear relation ΔΩ/Ω₀ ∝ 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₀/Ω₀)² ≳ 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₀/Ω₀)², and Re_c were varied; the values of these parameters are listed in Table 4.

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:

(50)

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

(51)

which is equivalent to

(52)

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₀/Ω₀)² = 10⁻¹ and 10⁻², and a fixed low Ekman number E = 10⁻⁶.

This figure first shows that a linear regime exists as long as ΔΩ/Ω₀ is low enough. The differential rotation is clearly linearly dependent on R_o until ΔΩ/Ω₀ ≳ 10⁻¹, 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₀/Ω₀)² 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,

(53)

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 jump in differential rotation across the outer Ekman layer. The second term scales with P_r(N₀/Ω₀)²/E and is of the Eddington–Sweet type. It connects smoothly to the outer sphere boundary condition, that is δΩ_ED(r₀, θ) = 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₀/Ω₀)² = 10⁻¹ and 10⁻², and the same Ekman number E = 10⁻⁶, 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 non-negligible 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⁻² 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 (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₀/Ω₀)², although determining their exact scaling with these non-dimensional 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/ρ₀, 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/ρ₀ > (r₀/r_i)² ≈ 11.11) the absolute value of the contraction velocity field even decreases with depth.

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₀/Ω₀)² parameter, and the contraction Reynolds number are fixed to 10⁻⁵, 10⁻⁴, and 10⁻², 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⁻⁴ and E = 10⁻⁶, 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/ρ₀ = 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 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 density-stratification, 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.

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/ρ₀ = 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.

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/ρ₀ = 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/ρ₀ = 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⁻², 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 ΔΩ/Ω₀ ≪ 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

(54)

where D² is the vector Laplace operator in the azimuthal direction D² = (∇²−1/r²sin²θ). When the Coriolis and the non-linear advection terms are neglected in Eq. (54), and when the differential rotation is sufficiently weak (i.e. ΔΩ/Ω₀ ≪ 1, and thus the first part of the contraction term can be neglected), the steady balance is the linear equation

(55)

whose solution is

(56)

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⁻²) is represented at different density contrasts. In the other panels the density contrast between the inner and outer spheres is fixed to ρ_i/ρ₀ = 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 Ω₀ 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.

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, ΔΩ/Ω₀ ≪ 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/ρ₀ considered together with the other parameters in the simulations are listed in Table 8. The parameters E, P_r(N₀/Ω₀)², 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.

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:

(57)

Here the superscript ‘A’ stands for anelastic, and differentiates this contraction timescale from τ_c = r₀/V₀ used in the Boussinesq case. The second expression is obtained using the definition of V_f (r) and assuming that r ∼ r₀.

Another specificity of the anelastic equations is the presence of an additional term in the thermal balance,

(58)

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

(59)

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/ρ₀ = 1.75.

The latitudinal variation of the inner sphere rotation rate found at ρ_i/ρ₀ = 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.

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 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₀/Ω₀)² controls the amplitude and the distribution of the differential rotation. When P_r(N₀/Ω₀)² ≫ 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₀/Ω₀)² ≪ 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 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 ΔΩ/Ω₀ ∼ 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 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 sound-proof approach other than the LBR approximation would lead to problems associated with the non-conservation 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 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 non-axisymmetric 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, 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).

---

Acknowledgments

The authors acknowledge the developers of MagIC (https://github.com/magic-sph/magic) for the open-source 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

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,

(A.1)

where

(A.2)

and

(A.3)

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 non-dimensionalised as Θ′(r, θ):

(A.4)

If R_o ≪ (N₀/Ω₀)², the steady case of the foregoing equation reads

(A.5)

Using both Ostrogradsky’s theorem and the continuity equation, we can write

(A.6)

which then allows us to deduce

(A.7)

Furthermore,

(A.8)

According to Eq. (A.3) the first integral in brackets is zero. Since the last integral is also zero, we conclude the following:

(A.9)

Consequently, the latitudinal mean of Eq. (A.5) gives

(A.10)

The radial derivative of is given by solving and leads to

(A.11)

Finally, is derived by solving

(A.12)

using the boundary conditions . The analytical solution is readily

(A.13)

We thus define

(A.14)

and likewise

(A.15)

After subtracting the spherically symmetric temperature field induced by the contraction velocity field, the evolution equation of temperature fluctuations finally reads

(A.16)

while the momentum equation Eq. (7) becomes

(A.17)

in which we have subtracted the hydrostatic state induced by the radial contraction

(A.18)

with , and where is the contraction-induced hydrostatic pressure non-dimensionalised with the hydrostatic pressure scale r₀ ρ₀ g₀.

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

(B.1)

where and . In cylindrical coordinates the steady axisymmetric velocity field reads

(B.2)

where is the dimensionless stream function. By using Eq. (B.1) we can then derive an analytical solution for ,

(B.3)

while the use of Eq. (B.2) allows us to determine the vertical velocity field:

(B.4)

By using the classical relations between the cylindrical and spherical coordinates, we deduce

(B.5)

(B.6)

The no-slip condition at the inner sphere leads to

(B.7)

By symmetry, at the equator, which means that when , . This allows us to find the analytical solutions of the meridional velocity fields:

(B.8)

(B.9)

By integrating Eq. (B.7), we deduce the stream function (up to a constant):

(B.10)

Since the stream function is zero at the poles, we deduce that cte = 0.

The mass conservation leads to the relationship

(B.11)

which gives in spherical coordinates:

(B.12)

From the analysis of the Ekman layer located at the outer sphere, we obtain the non-dimensional form of the Ekman condition:

(B.13)

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:

(B.14)

In addition,

(B.15)

thus, with the following substitutions and dz = cos θ dr, the relationship Eq. (B.12) becomes

(B.16)

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

(B.17)

Likewise, with the inward dimensionless mass flux crossing the cylinder, we obtain

(B.18)

Finally, from Eq. (B.16) we get

(B.19)

from which we conclude

(B.20)

which can then be rewritten under dimensional form to lead to the solution Eq. (45).

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

(C.1)

The poloidal-toroidal decomposition now relates the stream function to the meridional velocity fields as

(C.2)

which is equivalent to

(C.3)

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:

(C.4)

(C.5)

(C.6)

Instead, the previous analytical solution derived for the stream function is unmodified:

(C.7)

In the anelastic approximation, Eq. (B.11) is rewritten as follows:

(C.8)

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:

(C.9)

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–Sweet-type 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

(D.1)

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

(D.2)

where E^−1/2 ≪ P_r(N₀/Ω₀)²/E. Substitution of this expansion in Sys. (D.1) yields at the dominant order the following system:

(D.3)

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 no-slip 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

(D.4)

Meanwhile, within the Ekman layer, the boundary layer solution must be of the following order to compensate the 𝒪(1) meridional velocity jump:

(D.5)

By rewriting Sys. (D.1) with the stretched coordinate ξ and by inserting the previous expressions, the boundary layer equations read as follows:

(D.6)

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

(D.7)

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

All Figures