© B. Gouhier et al. 2022

1. Introduction

Rotation is ubiquitous at every stage of stellar evolution. Yet, it is still often considered as a second order effect in stellar evolution models. A full description necessarily entails a thorough study of differential rotation and meridional circulation induced by rotation, how they interact with other physical processes (e.g. magnetic field or internal gravity waves), and the ensuing potential instabilities. In his seminal work, Zahn (1992) proposed a model for the transport of chemical elements and angular momentum (AM) (without a magnetic field and internal gravity waves) that was later implemented in state-of-the-art stellar evolution codes. One of the strong assumptions is that the differential rotation in radiative zones is close to shellular (i.e. constant on an isobar) because of the anisotropic turbulence induced by shear instabilities in stably stratified conditions. This formalism has been successful at explaining a large number of stellar observations (see Maeder 2008 for a review). However, a growing body of evidence shows that additional AM transport mechanisms are still needed. This is particularly true for contracting stars. During the pre-main sequence (PMS) or post-MS evolution, a spin-up of stellar cores is naturally produced. However, observations of surface rotation of stars in young clusters (Gallet & Bouvier 2013) as well as asteroseismic studies of red giant stars (Eggenberger et al. 2012; Ceillier et al. 2013; Marques et al. 2013) tend to show that rotation rates are in fact strongly overestimated in models. In order to improve this formalism, the first thing to ask is what kind of flows are really expected in those contracting stars.

In a previous work (Gouhier et al. 2021), we investigated the differential rotation and meridional circulation produced in a modelled contracting stellar radiative zone. This axisymmetric study included the effects of stable stratification and was conducted both under the Boussinesq and the anelastic approximations, but the effects of a magnetic field were ignored. We showed that a radial differential rotation should be expected only in strongly stably stratified radiative zones such as the degenerate cores of subgiants. Indeed, any meridional circulation is inhibited by the strong buoyancy force, and the characteristic amplitude of the differential rotation in the linear regime is found to be proportional to the ratio of the viscous to contraction timescales ΔΩ/Ω₀ ∝ τ_ν/τ_c. In conditions relevant for the outside of the degenerate cores of subgiants and for PMS stars though, thermal diffusion weakens the stable stratification and allows a meridional circulation to exist, with a typical amplitude of the order of the contraction speed. The differential rotation profile then exhibits both a dependence in latitude and radius and its characteristic amplitude in the linear regime is found to be ΔΩ/Ω₀ ∝ τ_ED/τ_c, where τ_ED is the Eddington–Sweet timescale associated with the AM transport by the meridional flow. Both estimates, assuming a contraction timescale comparable to the Kelvin-Helmholtz timescale, tend to indicate that the amplitude of the differential rotation can be quite strong and that another process of AM transport should be invoked to reproduce the rotation rates of pre-MS and post-MS stars.

The presence of a large-scale magnetic field in such contracting radiative zones can drastically modify this picture. It is commonly agreed that radiative zones can host such fields (see Braithwaite & Spruit 2017 for a review). Magnetic fields with surface intensities of a few hundred Gauss have indeed been detected in a fraction of intermediate-mass PMS stars (Alecian et al. 2013). These stars are most probably the progenitors of the MS chemically peculiar Ap/Bp stars that host large-scale mostly dipolar fields with intensities ranging from 300 G to 30 kG (Donati & Landstreet 2009). Meanwhile, a distinct population of MS intermediate-mass stars, including the A-type star Vega and the Am-type stars Sirius, β Ursae Majoris and θ Leonis, exhibits much weaker (∼1 G) multi-polar magnetic fields (Lignières et al. 2009; Petit et al. 2010, 2011; Blazère et al. 2016). This magnetic dichotomy could be explained if, during the PMS, contraction forces a differential rotation that destroys pre-existing large-scale weak magnetic fields through magnetohydrodynamic (MHD) instabilities (Aurière et al. 2007; Lignières et al. 2014; Jouve et al. 2015, 2020).

Unlike PMS stars, the radiative zone of post-MS stars is overlaid by a large convective envelope preventing any direct measurement of a magnetic field in the convectively stable regions. Recently, asteroseismology revealed a class of red giants exhibiting dipolar oscillation modes with a very weak visibility (Mosser et al. 2012). This phenomenon has been attributed to the presence of an internal magnetic field modifying the angular structure of the dipole waves, thus leading to their trapping in the radiative core (Fuller et al. 2015). This so-called greenhouse effect, supported by other authors (Stello et al. 2016a,b; Cantiello et al. 2016) is, however, seriously questioned because these modes preserve their mixed character at odds with Fuller’s scenario (Mosser et al. 2017). In parallel, the possibility to detect magnetic fields through their effect on the oscillation frequencies is under investigation (see e.g. Bugnet et al. 2021). Awaiting observational evidence, theoretical work strongly supports the presence of magnetic fields in stars that possess a convective core (MS stars with M ≳ 1.2 M_⊙). The numerical simulations of core-dynamos of MS A- and B-type stars (Brun et al. 2005; Augustson et al. 2016) indicate that magnetic fields with intensities ranging from 0.1 − 1.0 MegaGauss can indeed be generated. Such a magnetic field could then relax into a large-scale stable configuration in the radiative interior of the post-MS stars where it would be buried for the rest of its evolution (Braithwaite & Spruit 2004).

Large-scale magnetic fields are able to impose a quasi-solid rotation on very short timescales, even for very weak intensities (Ferraro 1937; Mestel & Weiss 1987). Besides, enforcing a quasi-solid rotation during ∼1 Gyr after the end of the MS enabled Spada et al. (2016) to reproduce the rotation rates of subgiants, as measured by asteroseismology (Deheuvels et al. 2014). This result was also obtained by Eggenberger et al. (2019) and an observational support was then provided by Deheuvels et al. (2020) who measured a near solid rotation in two young subgiants. Interestingly, the efficiency of the AM transport then seems to decrease up to the tip of the red giant branch (RGB) before increasing again (Spada et al. 2016; Eggenberger et al. 2019; Deheuvels et al. 2020). Those works suggest a scenario broken down into several key phases. First, a quasi-solid rotation is maintained for some time through an efficient AM transport mechanism (possibly magnetic tension imposed by a large-scale field). Then, this mechanism would become inefficient and differential rotation would build up again before another AM transport mechanism, such as turbulent transport induced by MHD instabilities (Spruit 2002; Rüdiger et al. 2015; Fuller et al. 2019; Jouve et al. 2020) takes over.

In this work, we intend to study the flows induced by a contracting radiative zone in the presence of a large-scale magnetic field, through axisymmetric MHD simulations. In particular, we focus on the structure of the steady-states differential rotation and on the ability of the magnetic field to transport AM. The paper is organised as follows: in Sect. 2 we present the mathematical model, then in Sect. 3, the different timescales involved in our problem. The initial and boundary conditions as well as the numerical method are described in Sects. 4 and 5 respectively. In Sect. 6 we provide the reader with the relevant timescales in the stellar context and the consequences for our numerical study. The results of the simulations in the viscous and Eddington–Sweet regimes are finally given in Sect. 7, and are then summarised in Sect. 8 where the astrophysical implications are also discussed.

2. Mathematical formulation

In our previous work (Gouhier et al. 2021), we investigated the differential rotation and meridional flows produced in a contracting stellar radiative zone. In this follow-up work, we add the effect of an initial large-scale magnetic field. To do so, we numerically solve the Boussinesq or anelastic magnetohydrodynamical (MHD) equations in a spherical shell filled with a stably-stratified fluid subject to a radial contraction and embedded in a magnetic field. In this section we present the governing equations that we numerically solve in the two aforementioned approximations.

In this study, the fluid contraction is modelled using a mass-conserving contraction velocity field defined by

$$\begin{array}{c}\hfill {\displaystyle {\mathit{V}}_f=V_f(r){\mathit{e}}_r=-\frac{V_0{\rho }_{0}r{}_0^2}{\overline{\rho }{r}^2}{\mathit{e}}_r,}\end{array}$$(1)

where r is the radius and $\overline{\rho }$ the background density, r₀ and ρ₀ their respective values at the outer sphere and V₀ is the amplitude of the contraction velocity at the outer sphere. Using the Lantz-Braginsky-Roberts (LBR) approximation (Lantz 1992; Braginsky & Roberts 1995), assuming a uniform kinematic viscosity ν, thermal diffusion κ, magnetic diffusion η, and neglecting the centrifugal effects and local sources of heat, the dimensionless axisymmetric anelastic equations of a magnetised fluid (Jones et al. 2011) undergoing contraction read

$$\begin{array}{c}\hfill \begin{array}{c}\mathbf{\nabla }·\left[\stackrel{\sim }{\overline{\rho }}(\stackrel{\mathbf{\sim }}{\mathit{U}}+{\stackrel{\mathbf{\sim }}{\mathit{V}}}_f)\right]=0\mathrm{and}\mathbf{\nabla }·\stackrel{\mathbf{\sim }}{\mathit{B}}=0,\hfill \end{array}\end{array}$$(2)

$$\begin{array}{cc}& R_o\left[{\displaystyle \frac{\partial \stackrel{\mathbf{\sim }}{\mathit{U}}}{\partial t}+((\stackrel{\mathbf{\sim }}{\mathit{U}}+{\stackrel{\mathbf{\sim }}{\mathit{V}}}_f)·\mathbf{\nabla })(\stackrel{\mathbf{\sim }}{\mathit{U}}+{\stackrel{\mathbf{\sim }}{\mathit{V}}}_f)}\right]+2{\mathit{e}}_z\times (\stackrel{\mathbf{\sim }}{\mathit{U}}+{\stackrel{\mathbf{\sim }}{\mathit{V}}}_f)\hfill \\ \hfill & {\displaystyle =-\mathbf{\nabla }\left({\displaystyle \frac{{\stackrel{\sim }{\mathrm{\Pi }}}^{\prime }}{\stackrel{\sim }{\overline{\rho }}}}\right)+\frac{{\stackrel{\sim }{S}}^{\prime }}{{\stackrel{\sim }{r}}^2}{\mathit{e}}_r+\frac{E}{\stackrel{\sim }{\overline{\rho }}}\mathbf{\nabla }·\stackrel{\mathbf{\sim }}{\mathit{\sigma }}+{\left({\displaystyle \frac{L_{u}E}{P_m}}\right)}^2\frac{R_o^{-1}}{\stackrel{\sim }{\overline{\rho }}}[(\mathbf{\nabla }\times \stackrel{\mathbf{\sim }}{\mathit{B}})\times \stackrel{\mathbf{\sim }}{\mathit{B}}],}\hfill \end{array}$$(3)

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial \stackrel{\mathbf{\sim }}{\mathit{B}}}{\partial t}=\mathbf{\nabla }\times [(\stackrel{\mathbf{\sim }}{\mathit{U}}+{\stackrel{\mathbf{\sim }}{\mathit{V}}}_f)\times \stackrel{\mathbf{\sim }}{\mathit{B}}]+\frac{1}{R_m}{\mathbf{\nabla }}^2\stackrel{\mathbf{\sim }}{\mathit{B}},}\end{array}$$(4)

$$\begin{array}{cc}& \stackrel{\sim }{\overline{\rho }}\stackrel{\sim }{\overline{T}}[P_r{R}_o({\displaystyle \frac{\partial {\stackrel{\sim }{S}}^{\prime }}{\partial t}}+((\stackrel{\mathbf{\sim }}{\mathit{U}}+{\stackrel{\mathbf{\sim }}{\mathit{V}}}_f)·\mathbf{\nabla }){\stackrel{\sim }{S}}^{\prime })+P_r{\left({\displaystyle \frac{N_0}{{\mathrm{\Omega }}_0}}\right)}^2({\stackrel{\sim }{U}}_r+{\stackrel{\sim }{V}}_f){\displaystyle \frac{\mathrm{d}\stackrel{\sim }{\overline{S}}}{\mathrm{d}r}}]\hfill \\ \hfill & =E\mathbf{\nabla }·\left(\stackrel{\sim }{\overline{\rho }}\stackrel{\sim }{\overline{T}}\mathbf{\nabla }{\stackrel{\sim }{S}}^{\prime }\right)+D_{i}P{e}_c{E}^2[\stackrel{\sim }{\mathcal{Q}}{}_{\nu }+P_m^{-1}{\left({\displaystyle \frac{L_{u}E}{P_m{R}_o}}\right)}^2{(\mathbf{\nabla }\times \stackrel{\mathbf{\sim }}{\mathit{B}})}^2],\hfill \end{array}$$(5)

where the diffusion of entropy is introduced in the energy equation instead of the diffusion of temperature (Braginsky & Roberts 1995; Clune et al. 1999). The non-dimensional form (identified by the tilde variables) of these equations is obtained using the radius of the outer sphere r₀ as the reference lengthscale, the value of the contraction velocity field at the outer sphere V₀ as a characteristic velocity, and τ_c = r₀/V₀ as the reference timescale of the contraction. The frame rotates at Ω₀, the rotation rate of the outer sphere and all the thermodynamics variables are expanded as a background value plus fluctuations, respectively denoted with an overbar and 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 $\mathrm{\Delta }\overline{T}$ and $\mathrm{\Delta }\overline{S}$ 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. Finally, we use the value of the surface poloidal field at the poles B₀ as the reference scale for the magnetic field.

In Eqs. (2)–(5), U is the velocity field, ${\stackrel{\sim }{\sigma }}_{\mathit{ij}}$ is the dimensionless stress tensor, ${\stackrel{\sim }{\mathcal{Q}}}_{\nu }$ is the dimensionless viscous heating and the gravity profile is ∝r⁻². The reference state is non-adiabatic and a uniform positive entropy gradient is used to produce a stable stratification. It is related to the Brunt-Väisälä frequency defined by:

$$\begin{array}{c}\hfill N_0=\sqrt{{\displaystyle \frac{g_0}{C_{\mathrm{p}}}\frac{\mathrm{\Delta }\overline{S}}{r_0}}},\end{array}$$(6)

and which controls the amplitude of this stable stratification. The magnitude of the deviation to the isentropic state is controlled by the parameter ${ϵ}_s=\mathrm{\Delta }\overline{S}/C_{\mathrm{p}}$, chosen sufficiently small to ensure the validity of the anelastic approximation. With the dissipation number D_i = g₀r₀/T₀C_p, it sets the background temperature and density profiles (see Gouhier et al. 2021).

These anelastic equations involve six independent dimensionless numbers, a Rossby number based on the amplitude of the contraction velocity R_o = V₀/Ω₀r₀, the Ekman number $E=\nu /{\mathrm{\Omega }}_0{r}_0^2$, the Prandtl number P_r = ν/κ, the ratio between the reference Brunt-Väisälä frequency and the rotation rate of the outer sphere N₀/Ω₀, the magnetic Prandtl number P_m = ν/η, and the Lundquist number $L_u=B_0{r}_0/\sqrt{{\mu }_0{\rho }_0}\eta$ where μ₀ is the vacuum permeability. From these dimensionless numbers, three additional parameters can then be defined: a contraction Reynolds number Re_c = R_o/E, a Péclet number Pe_c = P_rRe_c, and a magnetic Reynolds number R_m = P_mRe_c.

When the compressibility effects are neglected, except in the buoyancy term of the momentum equation, the Boussinesq approximation is recovered. In that case, using Ω₀V₀T₀/g₀ as the scale of the temperature deviations Θ′, the dimensionless governing equations read

$$\begin{array}{c}\hfill \mathbf{\nabla }·(\stackrel{\mathbf{\sim }}{\mathit{U}}+{\stackrel{\mathbf{\sim }}{\mathit{V}}}_f)=0\mathrm{and}\mathbf{\nabla }·\stackrel{\mathbf{\sim }}{\mathit{B}}=0,\end{array}$$(7)

$$\begin{array}{cc}& R_o\left[{\displaystyle \frac{\partial \stackrel{\mathbf{\sim }}{\mathit{U}}}{\partial t}+((\stackrel{\mathbf{\sim }}{\mathit{U}}+{\stackrel{\mathbf{\sim }}{\mathit{V}}}_f)·\mathbf{\nabla })(\stackrel{\mathbf{\sim }}{\mathit{U}}+{\stackrel{\mathbf{\sim }}{\mathit{V}}}_f)}\right]+2{\mathit{e}}_z\times (\stackrel{\mathbf{\sim }}{\mathit{U}}+{\stackrel{\mathbf{\sim }}{\mathit{V}}}_f)\hfill \\ \hfill & =-\mathbf{\nabla }{\stackrel{\sim }{\mathrm{\Pi }}}^{\prime }+{\stackrel{\sim }{\mathrm{\Theta }}}^{\prime }\stackrel{\mathbf{\sim }}{\mathit{r}}+E{\mathbf{\nabla }}^2(\stackrel{\mathbf{\sim }}{\mathit{U}}+{\stackrel{\mathbf{\sim }}{\mathit{V}}}_f)+{\left({\displaystyle \frac{L_{u}E}{P_m}}\right)}^2{R}_o^{-1}[(\mathbf{\nabla }\times \stackrel{\mathbf{\sim }}{\mathit{B}})\times \stackrel{\mathbf{\sim }}{\mathit{B}}],\hfill \end{array}$$(8)

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial \stackrel{\mathbf{\sim }}{\mathit{B}}}{\partial t}=\mathbf{\nabla }\times [(\stackrel{\mathbf{\sim }}{\mathit{U}}+{\stackrel{\mathbf{\sim }}{\mathit{V}}}_f)\times \stackrel{\mathbf{\sim }}{\mathit{B}}]+\frac{1}{R_m}{\mathbf{\nabla }}^2\stackrel{\mathbf{\sim }}{\mathit{B}},}\end{array}$$(9)

$$\begin{array}{cc}& {\displaystyle P_r{R}_o\left[{\displaystyle \frac{\partial {\stackrel{\sim }{\mathrm{\Theta }}}^{\prime }}{\partial t}+((\stackrel{\mathbf{\sim }}{\mathit{U}}+{\stackrel{\mathbf{\sim }}{\mathit{V}}}_f)·\mathbf{\nabla }){\stackrel{\sim }{\mathrm{\Theta }}}^{\prime }}\right]+P_r{\left({\displaystyle \frac{N_0}{{\mathrm{\Omega }}_0}}\right)}^2({\stackrel{\sim }{U}}_r+{\stackrel{\sim }{V}}_f)\frac{\mathrm{d}\stackrel{\sim }{\overline{T}}}{\mathrm{d}r}}\hfill \\ \hfill & =E{\mathbf{\nabla }}^2{\stackrel{\sim }{\mathrm{\Theta }}}^{\prime },\hfill \end{array}$$(10)

where the gravity profile is now ∝r, the reference Brunt-Väisälä frequency is defined using the correspondence $\mathrm{\Delta }\overline{S}/C_{\mathrm{p}}=\mathrm{\Delta }\overline{T}/r_0$ in Eq. (6), and the contraction velocity field Eq. (1) is simplified using $\overline{\rho }={\rho }_0$.

To conclude this section, we note that when the anelastic approximation is used the parameter space is defined by eight dimensionless numbers: ϵ_s, D_i, Re_c, E, P_r, P_m, L_u and $N_0^2/{\mathrm{\Omega }}_0^2$. Instead, only the last six are necessary in the Boussinesq approximation. In this study, only Re_c, L_u and the product P_r(N₀/Ω₀)² will be varied.

3. Timescales of physical processes

In this section, we describe the various timescales involved in the transport of AM in our problem. We start by briefly recalling the relevant hydrodynamical timescales (see Gouhier et al. 2021), then we introduce two timescales associated with the presence of a magnetic field.

In this work, the dimensional form of the AM equation under the Boussinesq approximation reads

$$\begin{array}{cc}& {\displaystyle \frac{\partial U_{\varphi }}{\partial t}+\mathrm{NL}+\underset{\mathrm{Coriolis}\mathrm{term}}{\underbrace{2{\mathrm{\Omega }}_0{U}_s}}-\underset{\mathrm{Lorentz}\mathrm{force}}{\underbrace{{\displaystyle \frac{1}{{\mu }_0{\rho }_0}({\mathit{B}}_p·\mathbf{\nabla })B_{\varphi }}}}-\underset{\mathrm{Viscous}\mathrm{term}}{\underbrace{\nu D^2{U}_{\varphi }}}}\hfill \\ \hfill & =\underset{\mathrm{Contraction}}{\underbrace{{\displaystyle \frac{V_{0}r{}_0^2}{r^3}\frac{\partial }{\partial r}(r{U}_{\varphi }+r^{2}sin\theta {\mathrm{\Omega }}_0)}}},\hfill \end{array}$$(11)

where D² is the azimuthal component of the vector Laplacian operator, U_s = cos θU_θ + sin θU_r is the component of the velocity field, perpendicular to the rotation axis, and NL denotes the non-linear advection term. By balancing the partial time derivative with the contraction term in Eq. (11), we recover the contraction timescale used to non-dimensionalise the governing equations in Sect. 2

$$\begin{array}{c}\hfill {\displaystyle {\tau }_{\mathrm{c}}=\frac{r_0}{V_0},}\end{array}$$(12)

as well as its linear form

$$\begin{array}{c}\hfill {\displaystyle {\tau }_{\mathrm{c}}^{\mathrm{L}}=\frac{r_0}{V_0}\frac{\mathrm{\Delta }{\mathrm{\Omega }}_0}{{\mathrm{\Omega }}_0},}\end{array}$$(13)

when ΔΩ/Ω₀ ≪ 1. In the anelastic approximation, the contraction term in Eq. (11) is multiplied by ${\rho }_0/\overline{\rho }$ and the resulting timescale is then weighted by the background density profile

$$\begin{array}{c}\hfill {\tau }_{\mathrm{c}}^{\mathrm{A}}=\left({\displaystyle {\int }_{r_i/r_0}^1\frac{\overline{\rho }}{{\rho }_0}\mathrm{d}(r/r_0)}\right){\tau }_{\mathrm{c}}\mathrm{or}{\tau }_{\mathrm{c}}^{{\mathrm{L}}^{\mathrm{A}}}=\left({\displaystyle {\int }_{r_i/r_0}^1\frac{\overline{\rho }}{{\rho }_0}\mathrm{d}(r/r_0)}\right){\tau }_{\mathrm{c}}^{\mathrm{L}},\end{array}$$(14)

where ${\tau }_{\mathrm{c}}^{{\mathrm{L}}^{\mathrm{A}}}$ denotes the linear version. The AM transport by contraction can be balanced either by the viscous processes on a viscous timescale ${\tau }_{\nu }=r_0^2/\nu$, or by a meridional circulation of Eddington–Sweet type, in which case it redistributes the AM on the following timescale

$$\begin{array}{c}\hfill {\displaystyle {\tau }_{\mathrm{ED}}=\frac{r_0^2}{\kappa }{\left({\displaystyle \frac{N_0}{{\mathrm{\Omega }}_0}}\right)}^2.}\end{array}$$(15)

Ekman layers tend to develop at the spherical boundaries to accommodate the interior flow to the boundary conditions. In unstratified flows, they drive a global circulation. In stars, the stable stratification efficiently opposes this global flow although it can still exist in numerical simulations because the Ekman numbers can not reach stellar values. It then transports the AM on a spin-up timescale defined by

$$\begin{array}{c}\hfill {\tau }_{\mathrm{E}}=\sqrt{{\displaystyle \frac{r_0^2}{{\mathrm{\Omega }}_0\nu }}}.\end{array}$$(16)

The relative importance of these AM redistribution processes is given by the ratio of the above timescales, namely:

$$\begin{array}{c}\hfill {\displaystyle \frac{{\tau }_{\mathrm{ED}}}{{\tau }_{\nu }}=P_r{\left({\displaystyle \frac{N_0}{{\mathrm{\Omega }}_0}}\right)}^2;\frac{{\tau }_{\mathrm{E}}}{{\tau }_{\nu }}=\sqrt{E};\frac{{\tau }_{\mathrm{E}}}{{\tau }_{\mathrm{ED}}}=\frac{\sqrt{E}}{P_r{\left({\displaystyle \frac{N_0}{{\mathrm{\Omega }}_0}}\right)}^2}.}\end{array}$$(17)

Two main dimensionless numbers thus appear, the Ekman number E and the P_r(N₀/Ω₀)² parameter. As first noticed by Garaud (2002), the latter is of prime interest since it controls the flow dynamics. Thus, at fixed P_r and N₀, the author shows that depending on the rotation rate value this parameters defines two rotation regimes (slow or fast). More in line with our work, Garaud & Brummell (2008), Garaud & Garaud (2008), Garaud & Acevedo-Arreguin (2009), Wood & Brummell (2012, 2018), Acevedo-Arreguin et al. (2013) have shown that it also controls the efficiency of the burrowing of the meridional circulation in radiative layers adjacent to convective regions. In particular, when P_r(N₀/Ω₀)² ≫ 1, this circulation is suppressed by the stable stratification, a situation similar to the one that we encountered in the viscous regime described in Gouhier et al. (2021). In our case, E and P_r(N₀/Ω₀)² allow us to distinguish three regimes of interest:

$$\begin{array}{cc}& P_r{\left({\displaystyle \frac{N_0}{{\mathrm{\Omega }}_0}}\right)}^2\ll \sqrt{E}\ll 1;\sqrt{E}\ll P_r{\left({\displaystyle \frac{N_0}{{\mathrm{\Omega }}_0}}\right)}^2\ll 1;\hfill \\ \hfill & \sqrt{E}\ll 1\ll P_r{\left({\displaystyle \frac{N_0}{{\mathrm{\Omega }}_0}}\right)}^2.\hfill \end{array}$$(18)

As discussed in Gouhier et al. (2021), the two last regimes, namely the Eddington–Sweet regime τ_E ≪ τ_ED ≪ τ_ν and the viscous regime τ_E ≪ τ_ν ≪ τ_ED, are the most relevant for stars. We thus focus on them for the magnetic study.

The magnetic field introduces two new timescales: the magnetic diffusion timescale

$$\begin{array}{c}\hfill {\displaystyle {\tau }_{\eta }=\frac{r_0^2}{\eta },}\end{array}$$(19)

and the Alfvén timescale

$$\begin{array}{c}\hfill {\displaystyle {\tau }_{{\mathrm{A}}_{\mathrm{p}}}=\frac{r_0\sqrt{{\mu }_0{\rho }_0}}{B_0}.}\end{array}$$(20)

When the density contrast is taken into account (anelastic approximation), a new Alfvén timescale can be defined as

$$\begin{array}{c}\hfill {\tau }_{{\mathrm{A}}_{\mathrm{p}}}^{\mathrm{A}}=\left({\displaystyle {\int }_{r_i/r_0}^1\sqrt{{\displaystyle \frac{\overline{\rho }}{{\rho }_0}}}\mathrm{d}(r/r_0)}\right){\tau }_{{\mathrm{A}}_{\mathrm{p}}}.\end{array}$$(21)

4. Initial and boundary conditions

Initially, we impose a dipolar or a quadrupolar poloidal magnetic field. In both cases, the radial distribution of the magnetic field is such that the azimuthal current density does not depend on r, (∂j_ϕ/∂r = 0), to avoid possible numerical instabilities resulting from strong current sheets at the boundaries. For the dipole topology (left panel in Fig. 1), the initial field reads

$$\begin{array}{cc}& {\displaystyle \mathit{B}(r,\theta ,t=0)=\frac{3r{B}_0}{r_0(1-{(r_i/r_0)}^4)}cos\theta \left({\displaystyle 1+\frac{r_i^4}{3r^4}-\frac{4r_0}{3r}}\right){\mathit{e}}_r}\hfill \\ \hfill & {\displaystyle -\frac{3r{B}_0}{2r_0(1-{(r_i/r_0)}^4)}sin\theta \left({\displaystyle 3-\frac{r_i^4}{3r^4}-\frac{8r_0}{3r}}\right){\mathit{e}}_{\theta }.}\hfill \end{array}$$(22)

For the quadrupole topology (right panel in Fig. 1), it takes the following form

$$\begin{array}{cc}& {\displaystyle \mathit{B}(r,\theta ,t=0)=\frac{r{B}_0}{2r_0(1-{(r_i/r_0)}^5)}·(2{cos}^2\theta -{sin}^2\theta )}\hfill \\ \hfill & {\displaystyle \left({\displaystyle \frac{r_i^5}{r^5}-1-5ln\left({\displaystyle \frac{r_0}{r}}\right)}\right){\mathit{e}}_r-\frac{r{B}_0}{r_0(1-{(r_i/r_0)}^5)}sin\theta cos\theta }\hfill \\ \hfill & \left({\displaystyle 1-\frac{r_i^5}{r^5}-\frac{15}{2}ln\left({\displaystyle \frac{r_0}{r}}\right)}\right){\mathit{e}}_{\theta }.\hfill \end{array}$$(23)

Figure 1 shows that some field lines connect to the inner and outer boundaries while others are only connected to the outer boundary, or even loop-back on themselves inside the domain. In both configurations, the norm of the magnetic field at the poles and at the outer sphere is B₀.

Insulated boundary conditions are imposed at the inner and outer spheres. For our axisymmetric setup, these conditions translate into

$$\begin{array}{c}\hfill {\mathit{B}}_{\mathit{p}}=\mathbf{\nabla }\mathrm{\Phi }\mathrm{and}{B}_{\varphi }=0\mathrm{at}r=r_i,r_0,\end{array}$$(24)

where Φ is a potential field.

The rotation rate is chosen to be initially uniform Ω(r, θ, t = 0) = Ω₀ in the Boussinesq approximation. In the anelastic case, the initial profile is

$$\begin{array}{c}\hfill \mathrm{\Omega }(r,\theta ,t=0)=\overline{\rho }(r){\mathrm{\Omega }}_{0}exp\left({\displaystyle \frac{-(r-r_0)}{\sigma }}\right),\end{array}$$(25)

where σ controls the amplitude of the differential rotation.

For all simulations, we impose stress-free conditions at the inner sphere for the latitudinal and azimuthal components of the velocity field, and impermeability condition for the radial component of the velocity field:

$$\begin{array}{c}\hfill {\displaystyle U_r=\frac{\partial }{\partial r}\left({\displaystyle \frac{U_{\varphi }}{r}}\right)=\frac{\partial }{\partial r}\left({\displaystyle \frac{U_{\theta }}{r}}\right)=0\mathrm{at}r=r_i.}\end{array}$$(26)

At the outer sphere, we impose an impermeability condition on the radial component of velocity field, and no-slip conditions on the latitudinal and azimuthal components of the velocity field:

$$\begin{array}{c}\hfill U_r=U_{\theta }=U_{\varphi }=0\mathrm{at}r=r_0,\end{array}$$(27)

the rotation rate of the outer sphere being thus fixed to Ω₀. The boundary layers induced by these conditions in the absence of magnetic field were analysed in Gouhier et al. (2021).

Finally, in the Boussinesq approximation the temperature is prescribed at the inner and outer spheres, and the initial temperature field is the purely radial solution of the conduction equation. In the anelastic case, the entropy is also fixed at the boundaries. The initial stably stratified background density and temperature profiles are displayed in Gouhier et al. (2021).

5. Numerical method

The numerical study was carried out using the fully documented, publicly available code MagIC¹ to solve the set of axisymmetric magneto-hydrodynamical equations in a spherical shell under the anelastic approximation (Gastine & Wicht 2012) (Eqs. (3)–(5)), or under the Boussinesq approximation (Wicht 2002) (Eqs. (8)–(10)). The solenoidal condition of Eqs. (2) and (7) is ensured by a poloidal–toroidal decomposition for the mass flux and the magnetic field. Then, the different fields are expanded on the basis of the spherical harmonics for the horizontal direction, and on the set of the Chebyshev polynomials for the radial direction. In particular, the Chebyshev discretisation guarantees a better resolution near the boundaries. The extent of the spherical shell can be reduced to a two-dimensional domain such as 𝒟 = {r_i = 0.3 ≤ r ≤ r₀ = 1.0;0 ≤ θ ≤ π}.

In the viscous regime, most of the simulations were performed with N_r × N_θ = 127 × 256, while for the most resolved cases N_r × N_θ = 193 × 512. In the Eddington–Sweet regime, a higher resolution was needed and for most of the simulations N_r × N_θ = 193 × 512, where N_r could be increased when the boundary layers needed to be carefully resolved.

6. Space of parameters

In this section, we estimate the relevant timescales in stars in order to constrain the regime of parameters of our numerical study.

6.1. The stellar context

The magnetic fields considered in this paper are large-scale fossil fields. They can be remnants of the proto-stellar phase, or the product of a core-dynamo buried in the radiative zone. At large scales, the ohmic diffusion timescale is around one Gyr (Braithwaite & Spruit 2017), that is longer than the MS lifetime of the intermediate mass-stars. Dynamic processes such as the Alfvén waves propagation occur over much shorter timescales, and we will always have

$$\begin{array}{c}\hfill {\tau }_{{\mathrm{A}}_{\mathrm{p}}}\ll {\tau }_{\eta }.\end{array}$$(28)

Typical magnetic Prandtl numbers in stellar plasmas are ∼5 × 10⁻³ − 10⁻² (Rüdiger et al. 2016) so that τ_η ≪ τ_ν. There is an exception though in the core of subgiants where the electrons are partially or fully degenerate. In that case, the Prandtl and magnetic Prandtl numbers increase because the thermal and magnetic diffusivities as well as the kinematic viscosity are dominated by the electron conduction (Garaud et al. 2015). The magnetic Prandtl number then ranges from 0.1 to 10 (Cantiello & Braithwaite 2011; Garaud et al. 2015; Rüdiger et al. 2015) and consequently, τ_η ≤ τ_ν or τ_ν ≤ τ_η.

In Gouhier et al. (2021), we already found that contracting stars (PMS or subgiant stars) always lie in the regime

$$\begin{array}{c}\hfill {\tau }_{\mathrm{c}}\ll {\tau }_{\nu },{\tau }_{\mathrm{ED}}.\end{array}$$(29)

In addition, we showed that the Eddington–Sweet regime is relevant for PMS stars and outside the degenerate core of subgiants. We thus have P_r(N₀/Ω₀)² ≪ P_m ≪ 1, or τ_ED ≪ τ_η ≪ τ_ν, in these cases. By contrast, the degenerate cores of subgiants experience a viscous regime with higher P_m, such that 1 ∼ P_m ≪ P_r(N₀/Ω₀)², that is τ_ν ∼ τ_η ≪ τ_ED.

We can now wonder how the Alfvén time compares to the contraction time in stars. We lack precise information about the magnetic field intensities within stars, but we can get some insight from the spectropolarimetric data of Herbig stars, or from the asteroseismology of red giant stars combined to numerical simulations. On the one hand, high-resolution spectropolarimetric surveys show that a small fraction of HAeBes hosts large-scale dipolar fields stronger than a hundred Gauss (Wade et al. 2005; Alecian et al. 2013). For a typical Herbig star of 3 M_⊙ and 3 R_⊙ hosting a magnetic field of 300 G, the Alfvén poloidal timescale, computed using the mean density of the star, would be of the order of a few tens of years. For these PMS stars, the mass is between 2 M_⊙ and 5 M_⊙ and the Kelvin-Helmholtz timescale τ_KH typically ranges from 23 to 1.2 Myr (Maeder 2008). Then, assuming τ_c ≈ τ_KH, implies that τ_Ap ≪ τ_c.

On the other hand, the recent discovery of depressed dipole oscillation modes in red giants (Mosser et al. 2012) has been assigned to a greenhouse effect resulting from a strong magnetic field ∼1 MG trapping the gravity waves in the radiative core (Fuller et al. 2015). Although this scenario is controversial (see e.g. Mosser et al. 2017), the three-dimensional MHD simulations of convective core dynamos of Brun et al. (2005), Augustson et al. (2016) also point towards magnetic field intensities of the order of 10⁵ − 10⁶ G. Such intensities again lead to an Alfvén timescale much smaller than the contraction timescale. Even for a 1 G magnetic field in a subgiant such as KIC 5955122, which is a 1.1 M_⊙ star of 2 R_⊙ with a radiative interior extending to 0.74 R_⋆, we get an Alfvén timescale of 8 × 10³ yr. According to Deheuvels et al. (2020), the instantaneous contraction time defined as Ω_core/(dΩ_core/dt), where Ω_core is the mean rotation rate of the core, varies between 100 Myr and 3 Gyr in the subgiant phase. Its average value during this phase is ∼1 Gyr, thus much higher than the Alfvén timescale corresponding to a 1 G field.

6.2. The numerical study

We intend to perform numerical simulations in the timescale regimes thought to exist in stars. However, the regime τ_c ≪ τ_ν, τ_ED is strongly non-linear and too challenging numerically (Gouhier et al. 2021). The ratio τ_ν/τ_c = Re_c will instead vary in the range 0.1 − 5 which will allow us to study the viscous regime in the linear and non-linear regimes, while the Eddington–Sweet regime will be studied in the linear and weakly non-linear regimes as τ_ED/τ_c = P_r(N₀/Ω₀)²Re_c varies between 10⁻³ and 0.5.

The fact that τ_ν/τ_c ∼ 1 also constrains the magnetic Prandtl number. Indeed, to avoid a significant dissipation of the initial poloidal field during the simulation, the diffusion time τ_η must exceed the timescale τ_c ∼ τ_ν for the establishment of the stationary flow, which implies that P_m has to be larger than one. Such magnetic Prandtl numbers are expected in the degenerate cores of subgiants but are not realistic in the radiative envelope of subgiants or PMS stars. As in Charbonneau & MacGregor (1993), to prevent the diffusion of the poloidal field, an alternative option would have been to fix it. This is however not suited for the present problem where at large-scale, the field topology is modified by the contraction.

Finally, simulations were run for τ_Ω/τ_c = Re_cE, τ_Ω/τ_ν = E, and τ_Ω/τ_ED = E/P_r(N₀/Ω₀)² far from typical stellar values since realistic Ekman numbers are numerically unreachable. However, as shown in Gouhier et al. (2021), the flow dynamics do not critically depend on these ratios and we thus expect the model associated with our numerical simulations to remain valid for stars. Indeed, the first important parameter is τ_ED/τ_ν = P_r(N₀/Ω₀)² which determines if we are in the Eddington–Sweet or viscous regime. Realistic values can be used for this parameter. Another important ratio is τ_ν/τ_c = Re_c which governs the level of differential rotation. Realistic values for this last parameter are more difficult to reach numerically but we come back on the implications of this discrepancy in Sect. 8.

To conclude, all the simulations performed in the viscous regime fulfil the following conditions:

$$\begin{array}{c}\hfill {\tau }_{\mathrm{\Omega }}\ll {\tau }_{{\mathrm{A}}_{\mathrm{p}}}\ll {\tau }_{\mathrm{c}}\le {\tau }_{\nu }\ll {\tau }_{\eta }\ll {\tau }_{\mathrm{ED}},\end{array}$$(30)

or in terms of dimensionless numbers

$$\begin{array}{c}\hfill E\ll {\left({\displaystyle \frac{L_u}{P_m}}\right)}^{-1}\ll R{e}_c^{-1}\le 1\ll P_m\ll P_r{\left({\displaystyle \frac{N_0}{{\mathrm{\Omega }}_0}}\right)}^2.\end{array}$$(31)

For the Eddington–Sweet regime, we shall have

$$\begin{array}{c}\hfill {\tau }_{\mathrm{\Omega }}\ll {\tau }_{{\mathrm{A}}_{\mathrm{p}}}\ll {\tau }_{\mathrm{ED}}<{\tau }_{\mathrm{c}}\ll {\tau }_{\nu }\ll {\tau }_{\eta },\end{array}$$(32)

or equivalently

$$\begin{array}{c}\hfill E\ll {\left({\displaystyle \frac{L_u}{P_m}}\right)}^{-1}\ll P_r{\left({\displaystyle \frac{N_0}{{\mathrm{\Omega }}_0}}\right)}^2<R{e}_c^{-1}\ll 1\ll P_m.\end{array}$$(33)

7. Numerical results

We are mostly interested in the steady state differential rotation produced by the contracting flow in the stably stratified magnetised layer. We shall investigate separately the viscous regime $\sqrt{E}\ll 1\ll P_r{(N_0/{\mathrm{\Omega }}_0)}^2$ and the Eddington–Sweet regime $\sqrt{E}\ll P_r{(N_0/{\mathrm{\Omega }}_0)}^2\ll 1$. For the initial magnetic field we consider either a dipole or a quadrupole, and we include, or not, the effect of the density stratification. For each configurations we vary the Lundquist and the contraction Reynolds numbers to study the effect of the amplitudes of the initial poloidal field and of the contraction. To help understand physically our numerical results we also performed simulations where the contraction term is artificially removed from the induction equation, thus preventing the advection of the magnetic field by the contraction velocity field. All the relevant simulations and their associated parameters are listed in Table 1.

7.1. Unsteady evolution

Before we focus on the steady states, we start with a brief description of the unsteady phase. Figure 2 illustrates for a particular run, the typical evolution of the toroidal field B_ϕc (middle panel) and of the normalised differential rotation δΩ_c = (Ω_c−Ω₀)/Ω₀ (right panel) at some fixed locations in the domain (black points in the left panel). The run corresponds to a dipolar field in a viscous regime without advection of the field lines (run D5 in Table 1). We observe that as soon as the contraction produces a differential rotation, the initially zero toroidal field linearly grows by the Ω-effect. After ∼1τ_Ap, it saturates because the amplification of the toroidal field produces a Lorentz force that back-reacts on the differential rotation, thus counteracting further Ω-effect. This leads to the propagation of Alfvén waves along the poloidal field lines, illustrated by the oscillations of the control points in Fig. 2. As these waves oscillate independently from each other, they quickly get out of phase. This builds gradients of the toroidal magnetic field and of the differential rotation on sufficiently small scales so that they can be efficiently dissipated by the diffusion processes: this is the so-called phase mixing mechanism (see Heyvaerts & Priest 1983 or Cally 1991). In the absence of a process forcing the differential rotation, this would lead to a uniformly rotating steady state.

7.2. Steady state in the viscous regime with a dipolar field

In this section we investigate, in the viscous regime and for an initial dipolar field, the steady states obtained after the Alfvén waves have dissipated and the contraction has imposed some level of differential rotation. These states are actually quasi-steady because the initial poloidal field continues to slowly decrease through magnetic diffusion. Three parameters are fixed P_r(N₀/Ω₀)² = 10⁴, E = 10⁻⁴, and P_m = 10², while the contraction Reynolds number Re_c varies between 10⁻¹ to 5, and the Lundquist number between 10³ and 5 × 10⁴ (see details in Table 1). The anelastic simulations were performed with a density contrast ρ_i/ρ₀ = 20.85.

Representative results are displayed in Fig. 3, for a case without advection of the field lines at Re_c = 1 and L_u = 5 × 10⁴ (first two panels, run D5 of Table 1), and for a case with advection at Re_c = 1 and L_u = 10⁴ (last two panels, run D13 of Table 1). From this figure we clearly distinguish two magnetically decoupled regions with different levels of differential rotation. In the first region, the poloidal field lines are connected to the outer sphere and the flow is in quasi-solid rotation. In the second region, either the poloidal field lines loop-back on themselves inside the spherical shell (case without advection of the field lines) or they loop-back on the inner sphere (case with advection of the field lines). As in Charbonneau & MacGregor (1993), we shall refer to these particular regions as the ‘dead zone’ (DZ). The level of differential rotation is always significant in the DZs. With contraction of the field lines, the maximum differential rotation max(δΩ(r,θ)/Ω₀) = max((Ω(r,θ)−Ω₀)/Ω₀) is 30%, while it is reduced to ∼1% in the other case. The snapshots of the toroidal magnetic field (second and fourth panels in Fig. 3), show the presence of two anti-symmetric lobes of toroidal field in both hemispheres due to the Ω-effect acting on the dipolar poloidal field. Interestingly, the DZs are surrounded by a region of strong toroidal field. This corresponds to a Shercliff boundary layer (Shercliff 1956) that develops between two magnetic regions that are forced to rotate at different rates. Indeed, in our simulations the first poloidal line delimiting the DZ and the neighbouring lines connected to the outer sphere are forced to rotate differentially and a layer involving strong toroidal fields accommodates this jump in rotation rate. We verified that, as expected (Roberts 1967), the thickness of this Shercliff boundary layer scales with the inverse of the square root of the Hartmann number defined by $H_a=L_u/\sqrt{P_m}$. When we prevent the advection of the field lines (first panel in Fig. 3), a second Shercliff layer is visible at the separation between the region where the field lines are connected to both the outer and inner spheres, and the region where they are connected to the outer sphere only.

To give a more detailed description of the dynamics outside and inside the DZ, we compare in Fig. 4 the relative amplitudes of the different terms of the AM balance Eq. (11). From the first two panels of this figure, we see that outside the DZ, the quasi-steady configuration is characterised by a balance between the contraction term and the Lorentz force, that is:

$$\begin{array}{c}\hfill {\displaystyle -2sin\theta {\mathrm{\Omega }}_0\frac{V_0{r}_0^2}{r^2}=\frac{1}{{\mu }_0{\rho }_0}\left[{\displaystyle \frac{B_{\theta }}{rsin\theta }\frac{\partial }{\partial \theta }(sin\theta B_{\varphi })+\frac{B_r}{r}\frac{\partial }{\partial r}\left(r{B}_{\varphi }\right)}\right].}\end{array}$$(34)

On the contrary, the last two panels of Fig. 4 show that inside the DZ the viscous term balances the contraction term, thus leading to

$$\begin{array}{cc}& \nu \left[{\displaystyle r\frac{{\partial }^2}{\partial r^2}\left({\displaystyle \frac{\delta \mathrm{\Omega }}{{\mathrm{\Omega }}_0}}\right)+4\frac{\partial }{\partial r}\left({\displaystyle \frac{\delta \mathrm{\Omega }}{{\mathrm{\Omega }}_0}}\right)+\frac{1}{r}\frac{{\partial }^2}{\partial {\theta }^2}\left({\displaystyle \frac{\delta \mathrm{\Omega }}{{\mathrm{\Omega }}_0}}\right)+\frac{3cot\theta }{r}\frac{\partial }{\partial \theta }\left({\displaystyle \frac{\delta \mathrm{\Omega }}{{\mathrm{\Omega }}_0}}\right)}\right]\hfill \\ \hfill & {\displaystyle =\frac{-2V_0{r}_0^2}{r^2},}\hfill \end{array}$$(35)

where only the linear part of the contraction term has been retained, an approximation only valid if δΩ/Ω₀ ≪ 1. In the following two sub-sections, the differential rotation resulting from these two different balances is analysed.

7.2.1. Region outside the dead zone

Although at first glance the flow seems to be in solid rotation outside the DZ, there is a rotation rate jump across a boundary layer at the outer sphere as well as a residual differential rotation along the poloidal field lines.

By rescaling the colour range of Fig. 3, the top panel of Fig. 5 indeed shows that outside the DZ, the differential rotation between the interior flow and the outer sphere (Ω(r,θ)−Ω₀)/Ω₀) is ∼6 × 10⁻⁴. In order to understand that value, we first estimate the toroidal field amplitude outside the DZ using Eq. (34). Assuming B_r ∼ B_θ ∼ B₀ and r ∼ r₀, we get

$$\begin{array}{c}\hfill {\displaystyle \frac{B_{\varphi }}{B_0}\approx \frac{r_0{\mu }_0{\rho }_0{V}_0{\mathrm{\Omega }}_0}{B_0^2}={\left({\displaystyle \frac{P_m}{L_u}}\right)}^2\frac{R{e}_c}{E}.}\end{array}$$(36)

This estimate is confirmed in Fig. 6 where B_ϕ/B₀ determined at a particular location, θ = π/6, r = 0.65r₀, is compared with the right-hand side of Eq. (36) for the runs D1–D7.

This toroidal field is however unable to naturally match the vacuum condition imposed at the outer sphere (B_ϕ = 0 at r = r₀). This is done across a $H_a^{-1}$ thickness magnetic boundary layer known as an Hartmann boundary layer. This layer is analysed in Appendix B and the conclusion is that the $\mathcal{O}\left(P_m^{2}R{e}_c/E{L}_u^2\right)$ jump on B_ϕ/B₀, induces a $\mathcal{O}(\sqrt{P_m}R{e}_c/L_u)$ jump on the differential rotation ΔΩ/Ω₀ across the layer. This last scaling is confirmed in Fig. 7 by computing ΔΩ/Ω₀ at a particular location for the runs D1–D7.

Figure 5b shows a zoom on the zone delimited in black lines in Fig. 5a. It illustrates the residual differential rotation that exists along each poloidal field lines. We note it as ΔΩ_pol/Ω₀. Such a contrast of differential rotation is at odds with Ferraro’s law of isorotation (Ferraro 1937) requiring a constant angular velocity along each poloidal field line. Ferraro’s isorotation state is indeed found in our cases with no contraction in the induction equation. However, taking field advection into account, the Ω-effect term can be balanced by the radial advection term in the induction equation, that is:

$$\begin{array}{c}\hfill {\displaystyle V_f(r)\frac{\partial }{\partial r}\left({\displaystyle \frac{B_{\varphi }}{r}}\right)=sin\theta ({\mathit{B}}_p·\mathbf{\nabla })\delta \mathrm{\Omega }.}\end{array}$$(37)

Using the order of magnitude of the toroidal field derived in Eq. (36) together with Eq. (37), we end up with an estimate of this differential rotation

$$\begin{array}{c}\hfill {\displaystyle \frac{\mathrm{\Delta }{\mathrm{\Omega }}_{\mathrm{pol}}}{{\mathrm{\Omega }}_0}\approx {\left({\displaystyle \frac{R{e}_c{P}_m}{L_u}}\right)}^2={\left({\displaystyle \frac{{\tau }_{{\mathrm{A}}_{\mathrm{p}}}}{{\tau }_{\mathrm{c}}}}\right)}^2.}\end{array}$$(38)

Figure 8 shows the differential rotation taken between two ends of a poloidal field line as a function of the right-hand side of Eq. (38), for various Re_c and L_u respectively ranging from 10⁻¹ to 2 and from 5 × 10³ to 10⁴. We can see that the level of differential rotation along a field line is indeed consistent with our estimate.

Another consequence of the force balance Eq. (34) is that the toroidal field amplitude slowly increases over time. Figure 9 illustrates this by showing that in order to counteract the magnetic diffusion of the poloidal field (blue curve) and maintain this balance (black curve), the toroidal field amplitude is forced to increase (green curve). A similar situation was reported by Charbonneau & MacGregor (1993) where the magnetic stresses were in their case balanced by a wind torque.

7.2.2. Dead zone

We now turn to investigate the dynamics of the DZ. We already found that in this region, the dominant balance is between the contraction and the viscous effects, and is given by Eq. (35) in the linear regime ΔΩ/Ω₀  ≪  1. This implies that the order of magnitude of the differential rotation in the DZ, denoted ΔΩ_DZ/Ω₀, is ∼L_DZV_f(r)/ν, where L_DZ is a typical lengthscale of the DZ. While the magnetic field does not explicitly enter in this expression, it is of prime importance because its interaction with the contraction determines the shape and the location of the DZ.

The differences observed in Fig. 3 concerning the level of differential rotation between the simulations with and without contraction of the field lines can thus be explained by the location and size of the DZ. Indeed, in the case without field line contraction, the DZ is smaller and located closer to the outer sphere, that is in a area where the contraction velocity is weaker, which leads to a smaller level of differential rotation. Another difference between these two cases is the time taken to reach the stationary solution. As this time is the viscous time based on the size of the DZ, it is not surprising that the stationary state is reached much more rapidly when the DZ occupies a small fraction of the spherical shell. We note that in the case with contraction of the field lines, the viscous time is of the same order as the contraction time τ_c because we are in the regime ΔΩ/Ω₀ ∼ 1.

It is also possible to obtain a more precise determination of the DZ differential rotation, by deriving an approximate analytical solution of Eq. (35). To this end, the DZ is first assimilated to a conical domain, r ∈ [r_iDZ,r₀], θ ∈ [ − θ₀, θ₀], outside of which the flow is in solid rotation. In the case where the domain connects to the inner sphere, we adopt a stress-free condition on the azimuthal component of the velocity field at r = r_i. We then neglect the last term of the left-hand side in Eq. (35) as we found in our simulations that its contribution is small, particularly at low latitude. The homogeneous problem, which is separable in two radial and latitudinal eigenvalue sub-problems, allows us to construct a basis of orthogonal eigenfunctions that satisfy the boundary conditions on the conical domain. The details of the method are provided in Appendices C.1 and C.2.

In the case where the contraction does not advect the field lines, the approximate analytical solution takes the following form

$$\begin{array}{cc}& {\displaystyle \frac{\delta \mathrm{\Omega }(r,\theta )}{{\mathrm{\Omega }}_0}=R{e}_c{\displaystyle \sum _{n=1}^{\infty }\sum _{k=1}^{\infty }{A}_{\mathit{nk}}{\left({\displaystyle \frac{r_0}{r}}\right)}^{3/2}sin(\left({\displaystyle \frac{n\pi }{ln\left({\displaystyle \frac{r_{i_{\mathrm{DZ}}}}{r_0}}\right)}}\right)ln\left({\displaystyle \frac{r}{r_0}}\right))}}\hfill \\ \hfill & cos\left(\left({\displaystyle \frac{(2k-1)\pi }{2{\theta }_0}}\right)\theta \right),\hfill \end{array}$$(39)

wherein the expression of the coefficient A_nk can be found in Appendix C.1. This solution, computed for a chosen conical domain, is compared in Fig. 10 with the results of the numerical simulations performed at various Re_c and L_u respectively ranging from 10⁻¹ to 5 and from 5 × 10³ to 5 × 10⁴. A good agreement is found between the numerical and analytical solutions. The differential rotation scales as Re_c and this was expected because the condition δΩ/Ω₀  ≪  1 for the linear approximation of the contraction term is fulfilled in the simulations. We also find that the differential rotation is almost independent of the initial amplitude of the magnetic field: an increase of L_u of one order of magnitude causes only a small change on ΔΩ_DZ/Ω₀. This is consistent with the fact that, in the regime considered here, the shape and location of the DZ do not depend on the magnetic field. We rather observe that a higher magnetic field affects the rotation rate outside the DZ, or equivalently the differential rotation outside the DZ. Again this is expected as the jump in rotation rate across the outer sphere Hartmann layer decreases with the field amplitude.

When the contraction term is introduced in the induction equation, the solution of Eq. (35) becomes

$$\begin{array}{cc}& {\displaystyle \frac{\delta \mathrm{\Omega }(r,\theta )}{{\mathrm{\Omega }}_0}=R{e}_c{\displaystyle \sum _{n=1}^{\infty }\sum _{k=1}^{\infty }{A}_{\mathit{nk}}{\left({\displaystyle \frac{r_0}{r}}\right)}^{3/2}sin\left({\displaystyle \frac{1}{2}\sqrt{|9+4{\mu }_n|}ln\left({\displaystyle \frac{r}{r_0}}\right)}\right)}}\hfill \\ \hfill & cos\left(\left({\displaystyle \frac{(2k-1)\pi }{2{\theta }_0}}\right)\theta \right),\hfill \end{array}$$(40)

where the expressions of the coefficients A_nk and μ_n are given in Appendix C.2. This solution, computed for a chosen conical domain, is compared in Fig. 11 with the numerical results obtained at various Re_c and L_u. As previously, the initial amplitude of the magnetic field causes only small changes on the level of differential rotation. However, some discrepancies are visible. First of all, the numerical curves obtained at different Re_c do not overlap, which indicates that for these higher levels of differential rotation (δΩ/Ω₀ ∼ 30%) the contraction term can no longer be linearised, a necessary condition to derive the analytical solution. A second limitation is the assumed conical shape of the DZ. This shape is indeed more representative of the actual DZ for run D14 performed at Re_c = 2, where we can observe that the analytical solution tends to reproduce the expected differential rotation, than for simulations D10–D13 obtained at Re_c = 0.5 and 1.

7.2.3. Effect of the density stratification

We have performed simulations in the anelastic approximation for a fixed density contrast ρ_i/ρ₀ = 20.85. Two typical results are presented in Fig. 12 at Re_c = 1, L_u = 10⁵ and Re_c = 5, L_u = 5 × 10⁴ (left and right panels respectively). For these runs where the contraction advects the field lines, we find again two separate regions, with a differential rotation still mostly located in the DZ. This zone is now more extended in latitude compared to the Boussinesq case and the main difference between the calculations at two different Re_c lies in the level of differential rotation in the DZ. Indeed, when this level in the left panel is compared to its Boussinesq counterpart in the third panel of Fig. 3, we notice that it has been divided by a factor ∼3. This is explained by the density stratification effect on the contraction velocity field. Indeed, as we already found in Gouhier et al. (2021), the normalised differential rotation ΔΩ_DZ/Ω₀ resulting from the balance between the viscous and contraction effects is weighted by the inverse of the background density profile:

$$\begin{array}{c}\hfill {\displaystyle \frac{\mathrm{\Delta }{\mathrm{\Omega }}_{\mathrm{DZ}}^{\mathrm{A}}}{{\mathrm{\Omega }}_0}\approx \left({\displaystyle {\int }_1^{r_i/r_0}\frac{{\rho }_0}{\overline{\rho }}\mathrm{d}(r/r_0)}\right)\frac{\mathrm{\Delta }{\mathrm{\Omega }}_{\mathrm{DZ}}^{\mathrm{B}}}{{\mathrm{\Omega }}_0},}\end{array}$$(41)

where the index A stands for ‘anelastic’ and B for ‘Boussinesq’. For ρ_i/ρ₀ = 20.85 and r ∈ [r_i = 0.3;r₀ = 1], using the characteristic amplitude of the differential rotation given in the third panel of Fig. 3, we obtain $\mathrm{\Delta }{\mathrm{\Omega }}_{\mathrm{DZ}}^{\mathrm{A}}/{\mathrm{\Omega }}_0=0.095$, in agreement with the rotation rate shown in the left panel of Fig. 12. A good estimate of the rotation rate displayed in the right panel is then readily obtained by multiplying this result by Re_c.

7.3. Quadrupolar field in the viscous regime

We now consider an initial quadrupolar magnetic field as given by Eq. (23) and displayed in the right panel of Fig. 1. The outcome will be completely different as the differential rotation resulting from the contraction is subject to an instability that will be discussed in detail below.

Figure 13 displays typical states obtained after a contraction timescale. The two left panels show a case without contraction of the field lines while the two right panels correspond to simulations where contraction is included in the induction equation. At this stage, the similarities with the dipolar case are striking. First, we observe the quasi-solid rotation region outside the DZs. When the advection of the field lines is prevented, these DZs are located near the outer sphere while they connect to the inner sphere when contraction is present. In addition, for a given set of parameters, when the dipolar and quadrupolar cases are compared with each other, the same levels of differential rotation are found. For the runs presented here, the maximum amplitude of the normalised differential rotation inside the DZs is again ∼30% when the contraction acts on the field lines, and falls to ∼1% otherwise, similar to the cases presented in Fig. 3. By observing the second and fourth panels of Fig. 13, we note, again, the presence of magnetic boundary layers namely, the Shercliff layers separating the poloidal field lines constrained to rotate at different rates, and the Hartmann layer at the outer sphere. However, compared to the dipolar case, after ∼1τ_c the quadrupolar configuration is the seat of an axisymmetric instability.

The evolution of the maximum differential rotation is shown in Fig. 14 where the different steps that will be described thereafter are highlighted. In particular, the differential rotation first builds up before an instability starts to kick in at t ∼ 1τ_c (red dashed lines). Between ∼1τ_c and ∼3.5τ_c (blue dashed lines), the instability grows, saturates and strongly modifies the flow and field as we shall see later. Finally, after 3.5τ_c, the configuration evolves more smoothly until a final steady state is reached at ∼4.7τ_c (purple dashed lines).

7.3.1. Description of the instability

Figure 15 shows, in colour, the structure of the unstable modes when the contraction does not act on the field lines (first panel) and when it acts on them (second panel). In these meridional cuts are also represented, in black, the contours of the latitudinal shear ∂lnΩ/∂θ. Dashed (solid) lines represent negative (positive) values of this gradient. Contrary to the dipolar case, we now have a region of negative shear in the northern hemisphere and positive shear in the south. As can be observed, this is precisely at these locations that the perturbations grow. We already note that this instability is not of the centrifugal type because the shear in these regions is not strong enough, that is |∂lnΩ/∂lns| < 2. Figure 15 also shows that the unstable modes are characterised by small radial length scales and large horizontal length scales, which implies that the meridional motions are predominantly horizontal. The last panel of Fig. 15 displays the local evolution of the square of the latitudinal component of the velocity field at the fixed points marked on the second panel (in light blue, blue and black). We can see that the kinetic energy of the perturbations grows exponentially for about ∼20τ_Ap before saturation occurs. This time interval corresponds to the linear phase of the instability. The oscillating behaviour that adds to the exponential growth is due to the fact that the perturbations propagate.

Growth rates have been determined from these plots. Their values, normalised by the product of the surface-averaged shear parameter ⟨q⟩ = 1/S∬(∂lnΩ/∂θ)dS to the mean local rotation rate, are listed in Table 2. We observe that when the contraction Reynolds number is multiplied by two, the shear rate and the growth rate are both doubled. This shows that the growth rate of the instability seems to be proportional to the shear rate. Moreover, for the simulation Q5 performed at higher L_u = 5 × 10⁴, the instability is not triggered, clearly indicating that a strong enough poloidal field has a stabilising effect. Finally, while runs Q4, Q10 and Q11 carried out for L_u = 10⁴ are unstable, runs Q1 and Q8 performed for the same L_u, but at lower Re_c, are stable. This shows that, at a lower shear rate, a lower poloidal field is required to stabilise the flow. These findings, namely the requirement that the rotation decreases away from the rotation axis, and the facts that the growth rates are proportional to the shear and that stabilisation occurs above a certain magnetic tension, are all in agreement with an MRI-type instability (Balbus & Hawley 1991).

Our numerical results also exhibit more subtle effects that are not accounted for in the local WKB approach of the standard MRI (SMRI) (Balbus & Hawley 1991, 1994; Menou et al. 2004). First, the growth rates determined in Table 2 are significantly lower than the maximum growth rate σ_max = |q|Ω/2, where q = dlnΩ/dlns, predicted by these studies (e.g. Balbus & Hawley 1994). Second, a comparison between runs Q9 and Q10 shows that the growth rate increases with L_u whereas the predicted σ_max does not depend on the poloidal field. Finally, as noted earlier, the perturbations propagate in the domain while the phase velocity of the modes is zero in the SMRI.

We now argue that these differences are due to the effect of the stable stratification and to the presence of a toroidal field. In Balbus & Hawley (1994), the maximum growth rate is determined by assuming a zero latitudinal wavenumber. This avoids any stabilising effect of the stratification because the buoyancy force has no effect on purely horizontal motions. It follows that the most unstable radial scale is inversely proportional to the poloidal field amplitude and does not depend on the stable stratification. Even if the unstable motions found in our simulations are predominantly horizontal, assuming a zero latitudinal wavenumber is too extreme because our background flow is not uniform in latitude, so that the perturbation must have, and indeed, has a finite wavelength in this direction. As a consequence, the stable stratification is expected to play a role in determining the most unstable radial lengthscale and the associated maximum growth rate. We indeed found that the radial wavenumbers of the unstable modes are always larger the theoretical value from Balbus & Hawley (1994) and that it is very little dependent on the amplitude of the poloidal field. The effect of the stable stratification may thus potentially explain why the growth rates found in our simulations are significantly smaller than σ_max. We note that the thermal diffusion must be taken into account to consider the effect of the stable stratification. Indeed, we estimated that the thermal diffusion time scale ${\tau }_{\kappa }^{\ast }={\kappa }^{-1}/(k_r^2+k_{\theta }^2)$ associated with the observed unstable modes is about one order of magnitude smaller than the buoyancy time scale ${\tau }_B^{\ast }=(\sqrt{k_r^2+k_{\theta }^2}/k_{\theta })N^{-1}$, which implies that thermal diffusion will play an important role in determining the amplitude of the buoyancy force (e.g. Lignières et al. 1999).

To interpret the increase of the growth rate with L_u (runs Q9 and Q10), we first recall that in the context of the SMRI the toroidal field is assumed to be zero. This is not the case in our simulations where, according to Eq. (34), the amplitude of the stationary toroidal field decreases when the initial poloidal field (and thus L_u) increases. As reviewed in Rüdiger et al. (2018), introducing a toroidal field leads to the so-called helical MRI (HMRI) (Hollerbach & Rüdiger 2005). Following the dispersion relation derived by Kirillov et al. (2014), the toroidal field can be either stabilising or destabilising depending on the sign of $V_{A_{\varphi }}/s=B_{\varphi }/(rsin\theta \sqrt{{\mu }_0{\rho }_0})$. Moreover, the HMRI differs from the SMRI through the non-vanishing phase velocity of the unstable modes and a different phase shift between the perturbed fields. This phenomenon has been mentioned for the first time by Knobloch (1992), then found experimentally by Stefani et al. (2006) and in numerical simulations by Petitdemange et al. (2013).

In our simulations, the quantity V_Aϕ/s changes sign in the region where the instability is triggered so that the WKB-type analysis of Kirillov et al. (2014) is not directly conclusive. Our results nevertheless would imply that the toroidal field has a stabilising effect since the growth rate increases when the intensity of B_ϕ decreases. As noted earlier, we did observe that the perturbations propagate in our simulations. Their phase velocity V_phase, normalised by the contraction velocity field, are listed in Table 3. As expected, comparing the runs Q9 and Q10, we observe that the phase velocity is higher when the toroidal field is higher. Moreover, contrary to a SMRI, we observe neither exact phase quadrature between the latitudinal velocity perturbation and the latitudinal and azimuthal magnetic perturbations, nor exact opposition phase between the perturbed latitudinal and toroidal magnetic fields (Petitdemange et al. 2013). This is visible in Fig. 16 where the different perturbed quantities, the latitudinal velocity (purple) as well as the latitudinal (green) and azimuthal (blue) magnetic fields, are plotted as a function of the radius.

We conclude that the instability triggered in the DZ of the quadrupole is of the MRI-type and that its full description would require a detailed modelling of the effects of the stable stratification and of the toroidal field.

7.3.2. Non-linear evolution

After the exponential growth phase, the evolution becomes non-linear and the instability saturates. Figure 17 displays the flow structure obtained at ∼3τ_c for the run Q10 (Re_c = 1, L_u = 10⁴), through the rotation rate (colour) and the meridional circulation (black) in the left panel, as well as the norm of the poloidal field (colour) and its associated field lines (black) in the right panel. We observe that the instability proceeds via a multi-cellular meridional circulation, radially confined by the stable stratification and latitudinally extended in both hemispheres. From Fig. 17, we see that the poloidal field lines (right panel) are dragged around by this meridional circulation everywhere it is present (left panel), and then warped. The poloidal field thus behaves like a passive scalar advected and mixed by the multi-cellular circulation. This process creates small scales on which the magnetic diffusion can efficiently act to dissipate the poloidal field.

In Fig. 18, we compare the evolution of the norm of the poloidal field in an unstable case (in black) and in a stable case (in blue). The field norm is determined either through a volume average (solid line) or at a fixed point (displayed in black in Fig. 17) in the middle of the unstable region (dashed lines). The slope of these curves enables us to estimate a diffusion rate, and so a diffusion lengthscale, of the poloidal field. For the unstable configuration, this rate is determined during the saw-tooth evolution ranging from ∼2.25τ_c to ∼3.5τ_c. By denoting ${\omega }_{\text{stab}}=\eta /L_{\text{stab}}^2$ and ${\omega }_{\text{unst}}=\eta /L_{\text{unst}}^2$ the diffusion rates of the stable and unstable runs respectively, we obtain ω_unst/ω_stab ≈ 42, hence a diffusion of the poloidal field 42 times faster in the unstable case. In other words, the motions driven by the instability induce an effective diffusion at a lengthscale 6.5 smaller than the diffusion acting in the stable case (L_stab/L_unst = 6.5).

Although we shall see below that the differential rotation increases as a result of the instability, the toroidal field does not grow through the Ω-effect because the poloidal field is too weak in this phase. Then, the toroidal field experiences a diffusive-like decay, similar to the poloidal field.

7.3.3. Post-instability description

By destroying the poloidal field, the instability allowed a reconfiguration of the flow structure. This is shown at t = 3.5τ_c in the first two panels of Fig. 19, then at t = 4.66τ_c in the last two. From the first panel, we observe that the maximum level of differential rotation is now three times higher than before the instability. The reason is as follows: from a DZ to another, the large-scale structure of the poloidal field has been destroyed by the instability. As a result, even if a significant level of toroidal field still exists in this region, as displayed in the second panel of Fig. 19, the Lorentz force remains weak between the two DZs. The domain within which the contraction is balanced by the viscous effects thus becomes larger and, as expected, the differential rotation increases. By contrast, the poloidal field amplitude is still significant near the rotation axis and the Lorentz force imposes a very weak differential rotation 𝒪((Re_cP_m/L_u)²) (see Sect. 7.2.1) in this region.

Coming back to the numerical results presented in Fig. 19. The first two panels show that the magnetic topology has also completely changed after the development of the instability. A comparison between the third panel of Fig. 13 and the first panel of Fig. 19 shows that the field lines which looped back on themselves before the instability have now been moved towards the poles. In addition, the toroidal field is now very weak close to the outer sphere. An Hartmann layer is thus no longer needed to connect to the vacuum condition at the outer sphere. Likewise, the Shercliff layers have been removed with the dissipation of the poloidal field. Interestingly, we note that the new magnetic configuration, characterised by its positive lobe of toroidal field located near the rotation axis in both hemispheres, is from now on likely to be unstable to a non-axisymmetric instability of Tayler-type (see e.g. Spruit 1999).

After ∼4.66τ_c, the third panel of Fig. 19 shows that the differential rotation is mostly radial and occupies the whole shell. As a consequence, its amplitude further increased. In Fig. 20, we plotted in black dashed line the analytical solution corresponding to the balance on the full sphere between the viscous and contraction terms of Eq. (11). This solution, derived in Gouhier et al. (2021), perfectly matches the numerical solution in blue, thus showing that the hydrodynamic steady state is recovered. In conclusion, the magnetic field now has a negligible effect on the flow dynamics as supported by the fourth panel of Fig. 19 where we can see that the amplitude of the toroidal field has been divided by more than 10.

7.4. Steady state in the Eddington–Sweet regime

We now focus on the Eddington–Sweet regime ($\sqrt{E}\ll P_r{(N_0/{\mathrm{\Omega }}_0)}^2\ll 1$) considering a dipolar or a quadrupolar field as the pre-existing field. All simulations were performed in the anelastic approximation and included the contraction term in the induction equation. The contrast of density between the inner and the outer spheres was fixed to 20.85, and the Ekman and magnetic Prandtl numbers were respectively equal to 10⁻⁵ and 10². Our parametric study in this regime study consists in varying Re_c from 10⁻¹ to 5 and L_u from 5 × 10³ to 10⁵ for P_r(N₀/Ω₀)² = 10⁻² and 10⁻¹ (see details in Table 1).

We basically found three types of steady states. The first one is characterised, as in the viscous case, by two magnetically decoupled regions, one of them including a DZ where the contraction enforces a certain level of differential rotation. This state is the most relevant in a stellar context because it is obtained for the lowest values of τ_Ap/τ_ED and τ_Ap/τ_c, and these ratios are a priori very small in magnetic contracting stars as discussed in Sect. 6.1. As we increase these ratios, we find another type of solution where the differential rotation and the meridional circulation are no longer confined within the DZ while the field topology is unchanged (for higher τ_Ap/τ_ED), and finally a state where the advection of the poloidal field destroys the DZ and significantly reconfigures the magnetic field and the rotation profile (for higher τ_Ap/τ_c). In the following, these last two solutions will be described briefly as they are thought to be less relevant for our purpose, although physically interesting.

7.4.1. Meridional circulation and differential rotation confined to the dead zone

Figure 21 displays the typical structure of the quasi-steady flows and fields obtained for a dipolar (top row) or a quadrupolar (bottom row) initial field. These simulations were performed at Re_c = 1, L_u = 10⁵ and P_r(N₀/Ω₀)² = 10⁻¹ (runs D28 and Q13 of Table 1), and thus satisfy τ_Ap/τ_ED = 10⁻² and τ_Ap/τ_c = 10⁻³. There are many similarities with the viscous case. From the panels of the first column, we again observe two regions that are magnetically decoupled. One occupies the major part of the spherical shell and is in quasi-solid rotation while the other, the DZ, exhibits a certain level of differential rotation. The amplitude of this differential rotation is similar for the dipolar and quadrupolar cases. We also note the presence of magnetic boundary layers: the Hartmann layer at the outer sphere, and the Shercliff layers wherever adjacent poloidal field lines are forced to rotate differently, namely along the tangent cylinder and around the DZ. As in the viscous case, the toroidal field is characterised by a strong amplitude at these locations, as indicated by the panels of the second column.

We also see major differences with the viscous regime. First, although the contraction of the field lines is allowed, the DZ is confined near the outer sphere, whether a dipolar or quadrupolar field is initially imposed. This can be compared to the third panels of Figs. 3 and 13 in the viscous regime where the DZ was clearly advected towards the inner sphere. This difference is attributed to the effect of the contraction-induced meridional flow which now plays a significant role in the DZ advection. This flow is illustrated in the fourth column of Fig. 21 displaying the norm (colour) and the streamlines (black) of the poloidal velocity field U_p = U_re_r + U_θe_θ. For an initial dipolar field, this circulation is characterised by the presence of one cell of anticlockwise (clockwise) circulation in the northern (southern) hemisphere. This contraction-induced flow contributes to the total meridional circulation ${\mathit{U}}_p^{\mathrm{tot}}=(U_r+V_f){\mathit{e}}_r+U_{\theta }{\mathit{e}}_{\theta }$ displayed in the third column. As seen in the two top right panels, the induced flow inside the DZ tends to oppose contraction and the resulting total circulation becomes very weak, thus preventing the inward advection of the DZ. Outside the DZ, the total meridional flow is approximately parallel to the poloidal field lines close to the outer sphere, where the contraction velocity is maximum. In the deeper regions close to the inner sphere, the advection of poloidal field lines by the weaker contraction field is balanced by magnetic diffusion.

For an initial quadrupolar field, we observe a strong circulation around the DZ while inside the DZ the flow is predominantly vertical, downwards (upwards) in the northern (southern) hemisphere. Again, away from the DZ, the contraction-induced meridional flow has only a negligible contribution to the total meridional circulation. Finally, in contrast to the viscous regime, for the parameters numerically reachable in this study, the quadrupolar configurations are stable with respect to MRI because the shear built in the DZs is not strong enough to counteract the stabilising effect of the poloidal field. By comparison, the contrast of differential rotation in run Q5 of the viscous regime is ∼130 times larger and is not even unstable despite a weaker L_u of 5 × 10⁴.

In order to understand the flow dynamics inside and outside the DZ, we now examine the force balance in the AM equation Eq. (11), as was done in the viscous regime. The force amplitudes are analysed in the 2D maps of Fig. 22 where we display the ratio of the Lorentz force (left panel) and of the Coriolis force (right panel) to the contraction. We can observe that inside the DZ, the contraction is now balanced by the Coriolis force because the toroidal component of the magnetic field tends to zero and the Lorentz force becomes negligible accordingly (see right panel of Fig. 22). This implies that

$$\begin{array}{c}\hfill {\displaystyle U_s=sin\theta \frac{V_0{\rho }_0{r}_0^2}{\overline{\rho }{r}^2}.}\end{array}$$(42)

Here, contrary to the viscous case, the thermal diffusion weakens the stable stratification and enables a contraction-driven meridional circulation to exist. Inside the DZ, we also find that the thermal balance

$$\begin{array}{c}\hfill {\displaystyle U_r\frac{\mathrm{d}{S}_m}{\mathrm{d}r}=\kappa \left[{\displaystyle \left({\displaystyle \frac{\mathrm{d}ln\overline{\rho }}{\mathrm{d}r}+\frac{\mathrm{d}ln\overline{T}}{\mathrm{d}r}}\right)\frac{\partial \delta S}{\partial r}+{\mathbf{\nabla }}^2\delta S}\right],}\end{array}$$(43)

and the thermal wind balance

$$\begin{array}{c}\hfill {\displaystyle 2r{\mathrm{\Omega }}_0\frac{\partial U_{\varphi }}{\partial z}=\frac{g_0{r}_0^2}{C_{\mathrm{p}}{r}^2}\frac{\partial \delta S}{\partial \theta },}\end{array}$$(44)

are both satisfied. In these expressions, $S_m(r)=\overline{S}(r)+S_c(r)$ where ${\displaystyle S_c(r)=1/2{\int }_0^{\pi }{S}^{\prime }(r,\theta )sin\theta \mathrm{d}\theta }$ is the spherically symmetric entropy field induced by the contraction, and δS(r, θ) = S′(r, θ)−S_c(r) (see Gouhier et al. 2021 for details). According to Eq. (42), contraction then drives a meridional circulation U_s ∼ 𝒪(V_f) which redistributes AM on an Eddington–Sweet timescale. We note that the circulation timescale τ_c can be quite different from the Eddington–Sweet timescale. Indeed, as stated in Sect. 3, the ratio τ_ED/τ_c is measured by the dimensionless quantity Re_cP_r(N₀/Ω₀)². In this numerical study of the Eddington–Sweet regime, τ_ED ≪ τ_c because P_r(N₀/Ω₀)² ≪ 1 and because large contraction Reynolds numbers are too difficult to reach numerically.

Outside the DZ, the timescale of AM transport by the Alfvén waves is much shorter than the Eddington–Sweet timescale, and the Alfvén waves impose their dynamics. The left panel of Fig. 22 thus shows that the Lorentz force balances the contraction and Eq. (34) holds, as in the viscous regime. In this case, a quasi-isorotation state along the field lines is obtained, verifying:

$$\begin{array}{c}\hfill V_f(r)\left[{\displaystyle \frac{\partial }{\partial r}\left({\displaystyle \frac{B_{\varphi }}{r}}\right)-\frac{B_{\varphi }}{r}\frac{\mathrm{d}ln\overline{\rho }}{\mathrm{d}r}}\right]=sin\theta ({\mathit{B}}_p\times \mathbf{\nabla })\delta \mathrm{\Omega }.\end{array}$$(45)

As a result, the estimate of the characteristic amplitude of the differential rotation along the field lines Eq. (38) still holds, except that it must be weighted by ${\left({\displaystyle {\int }_{r_i/r_0}^1\frac{\overline{\rho }}{{\rho }_0}\mathrm{d}(r/r_0)}\right)}^{-1}$ accounting for the effect of the density stratification in the domain.

In the hydrodynamical case (Gouhier et al. 2021), we showed that the characteristic amplitude of the steady differential rotation resulting from the balance between the inward AM transport by the contraction and the AM redistribution by the Eddington–Sweet circulation should be $\mathcal{O}\left(P_r{(N_0/{\mathrm{\Omega }}_0)}^{2}R{e}_c{\left({\int }_{r_i/r_0}^1(\overline{\rho }/{\rho }_0)\mathrm{d}\stackrel{\sim }{r}\right)}^{-1}\right)$. This global analysis does not apply directly in the present situation where the DZ is reduced to a small fraction of the spherical shell, confined near the outer sphere. As in Sect. 7.2.2, and following Oglethorpe & Garaud (2013), to account for the DZ size and its effect on the differential rotation induced by the Eddington–Sweet circulation, we introduce the lengthscale L_DZ = 0.1r₀. Because of the density stratification, the contraction velocity is not very different between the outer and inner spheres. Thus, after using Eq. (42) and the continuity equation, we have U_r ≈ V₀. Then, from Eq. (43) we get $\delta S\approx \left(V_0{L}_{\text{DZ}}^2/\kappa \right)\left(\text{d}{S}_m/\text{d}r\right)$. Injecting this estimate in Eq. (44) yields finally to $\mathrm{\Delta }{\mathrm{\Omega }}_{\text{DZ}}/{\mathrm{\Omega }}_0\approx \left(g(r)/C_{\text{p}}\right)\left(\text{d}{S}_m/\text{d}r\right)\left(V_0{r}_0/\kappa {\mathrm{\Omega }}_0^2\right){\left(L_{\text{DZ}}/r_0\right)}^2\approx P_r{\left(N_0/{\mathrm{\Omega }}_0\right)}^{2}R{e}_c{\left(L_{\text{DZ}}/r_0\right)}^2$, thus enabling us to recover the level of differential rotation inside the DZ up to a factor two.

In Fig. 23, we plotted the maximum amplitude of the differential rotation inside the DZ as a function of P_r(N₀/Ω₀)²Re_c, for the runs D21–D24 and D27–D30 of Table 1 performed with Re_c ranging from 10⁻¹ to 2 (identified with symbols) and L_u from 5 × 10⁴ (light blue) to 10⁵ (blue). The other parameters are fixed to P_r(N₀/Ω₀)² = 10⁻¹, E = 10⁻⁵, P_m = 10², and ρ_i/ρ₀ = 20.85. As expected, the maximum contrast of differential rotation follows a linear relation with Re_c. Moreover, we also observe that this level is almost independent of the Lundquist number, consistent with the balance Eq. (42) inside the DZ. However, for the highest contraction Reynolds number Re_c = 2, we observe a clear discrepancy between runs performed at L_u = 5 × 10⁴ and L_u = 10⁵. This deviation can be attributed to the fact that, in the lower L_u case, the magnetic tension no longer prevents the advection of the magnetic field by the meridional flows. As shown in Fig. 24, this produces a significant deformation of the DZ geometry and a related expulsion of the magnetic flux outside the DZ. This phenomenon is discussed in more details below.

7.4.2. Meridional circulation and differential rotation not confined to a dead zone

Simulations carried out at a smaller P_r(N₀/Ω₀)² parameter (runs D32-D39 of Table 1) or at higher Re_c (runs D38, Q14, Q16 and Q17 of Table 1), that is at higher values of τ_Ap/τ_ED and τ_Ap/τ_c, exhibit different features. In the former case (smaller P_r(N₀/Ω₀)²), the differential rotation is no longer confined to the DZ as an AM redistribution by the Eddington–Sweet circulation occurs outside the DZ. In the second case (higher Re_c), the amplitude of the meridional circulation is strong enough to warp the DZ and expel the magnetic flux there. These two phenomena are now discussed.

The left panel of Fig. 25 displays a meridional cut of the quasi-steady differential rotation (in colour) obtained for Re_c = 1, L_u = 5 × 10⁴, and P_r(N₀/Ω₀)² = 10⁻² (run D36 of Table 1), on which are represented the vector lines of the total meridional velocity field as arrows and the poloidal field lines (in black). Compared to the simulation shown in Fig. 21, the ratio τ_Ap/τ_ED has been increased by a factor 20 (from 10⁻² to 2  ×  10⁻¹). Actually, if local values of this ratio are considered, a value of order 1 can be reached, in particular in the vicinity of the DZ. This is shown in the right panel of Fig. 25 that displays the distribution of the ratio of the Eddington–Sweet time to the local Alfvén time. In this regime, the differential rotation and the meridional circulation have spread out away from the DZ while the poloidal field lines, and thus the DZ, have not been affected by the circulation flow.

As presented in Fig. 26, another regime is encountered at sufficiently high Re_c. The left panel of this figure displays the differential rotation in colour with the streamlines of the meridional circulation in black. The right panel shows, in colour, the norm of the poloidal field with the poloidal field lines in black. We observe that this meridional circulation significantly warps the DZ, thus leading to a reconfiguration of the magnetic field and the differential rotation. Compared to the simulation shown in the top panel of Fig. 21, the ratio τ_Ap/τ_c has been increased by a factor of 6 (from 10⁻³ to 6 × 10⁻³). This is apparently sufficient for the Lorentz force not to be able to confine the circulation in the DZ. The magnetic field is then advected and partly dissipated in the vicinity of the original DZ. The dissipation process is reminiscent of the phenomenon of magnetic flux expulsion studied by Weiss (1966), whereby an eddy advects the magnetic field to such small scales that magnetic diffusion is efficient. We indeed observe that the magnetic flux ends up being expelled from the regions where the meridional circulation exists, and the poloidal field is found concentrated in free layers separating the quasi-solid rotation region from the one in differential rotation.

8. Summary and conclusions

In this work, we investigated the dynamics of a contracting radiative spherical layer embedded in a large-scale magnetic field. The aim was to determine the differential rotation that results from the combined effects of contraction and magnetic fields. The contraction was modelled through an imposed radial velocity field V_f and the gas dynamics was modelled using either the Boussinesq or the anelastic approximations. The parametric study was guided by the results obtained without magnetic field (Gouhier et al. 2021) highlighting two hydrodynamical regimes, namely the viscous regime in the strongly stratified cases and the Eddington–Sweet regime in the weakly stratified cases.

We find that the contracting layer first evolves towards a quasi-steady configuration characterised by two magnetically decoupled regions. In the first region, all the poloidal field lines connect to the outer sphere. The rotation is quasi-uniform in this region because the contraction only allows very small deviations from Ferraro’s isorotation law along the field lines and the outer sphere rotates uniformly. The second region, called the DZ, is decoupled from the first one as the poloidal field lines loop-back on themselves or connect to the inner sphere. In addition, the poloidal field amplitude vanishes at some point within the DZ. A significant level of differential rotation can be produced in these DZs, the inward AM transport by the contraction being balanced either by a viscous transport or by an Eddington–Sweet circulation. The exact amplitude of the differential rotation also depends on the size, the shape and the location of the DZ.

In a second step, after a time of the order of the contraction time, the shear built in the DZ can trigger a powerful axisymmetric instability that profoundly modifies the subsequent evolution of the flow. Indeed, for an initial quadrupolar field in the viscous regime, we observe that if the field strength is low enough, an MRI-type instability grows and produces a multi-cellular meridional circulation organised at small scales in the radial direction. This flow advects and eventually enables to efficiently dissipate the magnetic energy. The new field configuration is strongly modified, and the differential rotation which is no longer constrained to the DZ spreads to most of the spherical layer while its amplitude increases. This instability has not been observed for the quadrupolar field in the Eddington–Sweet regime because numerical limitations did not allow us to reach significant levels of differential rotation. However, we anticipate that for realistic contraction Reynolds numbers and Lundquist numbers, the differential rotation in the DZ of the quadrupolar field will also trigger an instability in the Eddington–Sweet regime. By contrast, the dipolar field configuration does not lead to an instability. Indeed, in this case, the DZ is symmetric with respect to the equator and the contraction produces maximum rotation rates along the equator. The latitudinal differential rotation thus increases away from the rotation axis which implies stability with respect to the MRI. We note that the same configuration in an expanding layer would lead to minimum rotation rates along the equator and thus to differential rotations potentially unstable to MRI.

If we intent to extrapolate to a more complex geometry of the initial poloidal field, the dipolar topology with a single equatorially symmetric DZ appears exceptional. Thus, we expect that generically negative latitudinal gradients of the rotation rate, potentially unstable to the MRI, are present in DZs. Rather than the topology of the poloidal field, what can prevent the MRI to develop is its intensity. Indeed, according to Balbus & Hawley (1998), the magnetic tension stabilises the flow if the perturbation length scales λ_r are smaller than $(B_r/\mathrm{\Omega })·(\sqrt{2}\pi /\sqrt{{\mu }_0\overline{\rho }\left|q\right|})$. Applying this criteria to the degenerate core of a typical subgiant of 1.1 M_⊙ and 2 R_⊙, we find that, assuming a 𝒪(1) shear |q|, a rotation rate Ω = 3.1 × 10⁻⁶ rad s⁻¹ and a mean core density ${\overline{\rho }}_{\mathrm{c}}=2.1\times {10}^5$ kg m⁻³, fields higher than 3 × 10⁵ G will stabilise all the perturbations smaller than the degenerate core size of 0.06 R_⊙. In practice, the radial wavelength of the unstable modes is constrained by the stable stratification rather than by the core size and thus even lower field intensities will be stabilising. For example, in our simulations λ_r is ∼44 times smaller than the outer radius of the spherical layer. The critical field in our numerical simulations is reached for a Lorentz number $L_o=B_0/\sqrt{{\mu }_0{\rho }_0}{r}_0{\mathrm{\Omega }}_0$ equal to ∼10⁻². For the subgiant core rotation and density given above, this corresponds to a ∼10⁴ G critical field strength. As P_r(N₀/Ω₀)² and P_m of the simulations are not too far from realistic values in subgiant cores, and the shear should remain limited to 𝒪(1) even for more realistic contraction Reynolds numbers, this critical field extrapolated from the simulations might be of the right order of magnitude. To be more precise, a closer look at the MRI driven by a negative rotation latitudinal gradient in a radiative zone will be necessary.

Our numerical study thus points towards the following scenario: during a first period of the order of the contraction timescale, a contracting radiative layer embedded in a large-scale poloidal field tends to rotate uniformly, except in localised DZs where the contraction induces a significant differential rotation. If the field is weak enough and not purely dipolar, the development of a powerful axisymmetric MRI reconfigures the field and diminishes its intensity. The magnetic coupling then becomes inefficient in the major part of the radiative layer and the contraction can force the differential rotation there.

Such a scenario could potentially explain the evolution of the rotation of the subgiants between the end of the MS and the tip of the RGB. As mentioned in the introduction, the asteroseismic data can be reproduced by assuming a uniform rotation during a first period after the end of the MS followed by a second period where the contraction is left free to enforce differential rotation (Spada et al. 2016). This is consistent with the two young subgiants in near solid-body found by Deheuvels et al. (2020). At their age, the post-MS contraction should have increased their core rotation by a factor of four which means that the period of uniform rotation should last at least a contraction timescale. This timescale is compatible with our scenario.

Our simulations are nevertheless far to describe the full complexity of a magnetic and contracting subgiant. In particular, an expanding layer and boundary conditions mimicking the effect of a convective envelope should be added. The role of non-axisymmetric instabilities should also be considered, especially in the magnetic configuration that results from the axisymmetric instability. Non-axisymmetric MRI or Tayler instability might indeed be present and take part to the AM transport particularly along the giant branch as already invoked (Cantiello et al. 2014; Fuller et al. 2019).

As far as intermediate-mass stars are concerned, the occurrence of a powerful contraction-driven instability could help explain the dichotomy between Ap/Bp and Vega-like magnetisms. Indeed, strong Ap/Bp-like magnetic fields are expected to survive the instability during the PMS while below a certain field intensity, the axisymmetric MRI would change the pre-existing large-scale field into a small-scale field of smaller amplitude, thus leading to Vega-like magnetism. This is in line with the scenario proposed by Aurière et al. (2007), except that the instability invoked in this paper was a non-axisymmetric instability produced during the initial winding-up of the poloidal field by the differential rotation. Numerical investigations of this process confirmed the presence of such instabilities but not their ability to profoundly modify the pre-existing poloidal field (Jouve et al. 2015, 2020). By contrast, the axisymmetric MRI found in the present paper has a very strong impact on the initial poloidal field, destroying its large-scale structure and even diminishing its amplitude. To test this scenario, the threshold field strength that separates MRI stable and from MRI unstable configurations is crucial because it should be compatible with the observed 300 G lower bound of Ap/Bp surface magnetic fields. Again, this calls for further numerical and theoretical investigations of the critical field of the MRI driven by rotation latitudinal gradients in radiative zones.

Part of the above discussion is based on extrapolations of our numerical results to stellar conditions. Our simulations are indeed a simplified version of a contracting star. Among the simplifications, the ratio between the contraction time and the rotation time is larger in stars than in our simulations (τ_c/τ_Ω ∼ 3.4 × 10⁸ − 1.1 × 10¹¹ in stars while this ratio is comprised between 10³ and 5 × 10⁵ in our simulations). However, the physical model derived from our simulation does not depend critically on this ratio. Indeed, by running our simulations for 5 − 6 contraction times, we observed that, after a contraction time, a powerful axisymmetric MHD instability develops. This leads to a complete reconfiguration of the initial magnetic field and to the subsequent development of differential rotation in most part of the spherical shell. This process should not be affected by increasing the ratio τ_c/τ_Ω to stellar values. In the same spirit, the ratio τ_Ω/τ_ν, the Ekman number, is much lower in stars than in numerical simulations. But as shown in Gouhier et al. (2021), the hydrodynamical AM transport is not affected when this ratio is decreased by various orders of magnitude. Thus, despite the simplifications inherent to numerical simulations, the physical model derived from these simulations seems robust enough to apply to stars, especially in the viscous regime where the MRI has been observed. A question that remains to be addressed in future works concerns the occurrence of the MRI in more realistic Eddington–Sweet regimes which will require to explore the strongly non-linear regime corresponding to very large ratio τ_c/τ_ν.

---

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. L.J. acknowledges funding by the Institut Universitaire de France. The authors wish to thank Sébastien Deheuvels for very fruitful discussions.

References

Appendix A: Electrical-current function

In this appendix, we determine the stream function that is constant along the streamlines of the poloidal component of the current density j_p. The toroidal component of the magnetic field is related to j_p through

$$\begin{array}{c}\hfill {\displaystyle {\mathit{j}}_{\mathit{p}}=\mathbf{\nabla }\times \left(B_{\varphi }{\mathit{e}}_{\mathit{\varphi }}\right)=\frac{1}{rsin\theta }\frac{\partial }{\partial \theta }(sin\theta B_{\varphi }){\mathit{e}}_{\mathit{r}}-\frac{1}{r}\frac{\partial }{\partial r}\left(r{B}_{\varphi }\right){\mathit{e}}_{\mathit{\theta }}.}\end{array}$$(A.1)

As the divergence of the curl is zero and the problem is axisymmetric, we can define a vector potential χ as

$$\begin{array}{c}\hfill {\displaystyle j_r=\frac{1}{rsin\theta }\frac{\partial }{\partial \theta }(sin\theta \chi );j_{\theta }=\frac{-1}{r}\frac{\partial }{\partial r}\left(r\chi \right),}\end{array}$$(A.2)

from which we define the electrical-current stream function 𝒥_p = r sin θχ:

$$\begin{array}{c}\hfill {\displaystyle j_r=\frac{1}{r^{2}sin\theta }\frac{\partial {\mathcal{J}}_p}{\partial \theta };j_{\theta }=\frac{-1}{rsin\theta }\frac{\partial {\mathcal{J}}_p}{\partial r}.}\end{array}$$(A.3)

We thus obtain a relationship between the toroidal field and the electrical-current stream function:

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{rsin\theta }\frac{\partial }{\partial \theta }(sin\theta B_{\varphi })=\frac{1}{r^{2}sin\theta }\frac{\partial {\mathcal{J}}_p}{\partial \theta };-\frac{1}{r}\frac{\partial }{\partial r}\left(r{B}_{\varphi }\right)=\frac{-1}{rsin\theta }\frac{\partial {\mathcal{J}}_p}{\partial r}.}\end{array}$$(A.4)

By integrating the first equation we have

$$\begin{array}{c}\hfill {\displaystyle {\mathcal{J}}_p=r{\displaystyle \int }\frac{\partial }{\partial \theta }(sin\theta B_{\varphi })\mathrm{d}\theta =rsin\theta B_{\varphi }+f(r),}\end{array}$$(A.5)

and substituting in the second equation we deduce

$$\begin{array}{c}\hfill {\displaystyle \frac{f^{\prime }(r)}{-rsin\theta }=0⟹f(r)=\mathrm{cte}.}\end{array}$$(A.6)

As the electrical-current stream function is zero at the poles, we conclude that cte = 0 hence:

$$\begin{array}{c}\hfill {\mathcal{J}}_p=rsin\theta B_{\varphi }.\end{array}$$(A.7)

Appendix B: Hartmann boundary layer equations

In the present appendix, the Hartmman boundary layer equations are derived with the aim of relating the jump of the toroidal field across the layer with the differential rotation jump. The Hartmann layer occurs at the electrically insulating outer boundary where the azimuthal velocity is fixed and a magnetic field perpendicular to the flow is present. When the contraction term is removed from the induction equation, the governing equations for the azimuthal flow and the toroidal field component in the steady linear limit read:

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle 2{\mathrm{\Omega }}_0\left({\displaystyle U_s-sin\theta \frac{V_0{r}_0^2}{r^2}}\right)=\frac{1}{{\mu }_0{\rho }_0}\left[{\displaystyle \frac{B_{\theta }}{r}\frac{\partial B_{\varphi }}{\partial \theta }+\frac{cot\theta }{r}{B}_{\theta }{B}_{\varphi }+B_r\frac{\partial B_{\varphi }}{\partial r}+\frac{B_r{B}_{\varphi }}{r}}\right]+}\hfill \\ \hfill \\ \nu \left[{\displaystyle \frac{{\partial }^2{U}_{\varphi }}{\partial r^2}+\frac{2}{r}\frac{\partial U_{\varphi }}{\partial r}+\frac{1}{r^2}\frac{{\partial }^2{U}_{\varphi }}{\partial {\theta }^2}+\frac{cot\theta }{r^2}\frac{\partial U_{\varphi }}{\partial \theta }-\frac{U_{\varphi }}{r^2{sin}^2\theta }}\right],\hfill \end{array}\end{array}$$(B.1)

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle U_r\frac{\partial B_{\varphi }}{\partial r}+\frac{U_{\theta }}{r}\frac{\partial B_{\varphi }}{\partial \theta }-B_r\frac{\partial U_{\varphi }}{\partial r}-\frac{B_{\theta }}{r}\frac{\partial U_{\varphi }}{\partial \theta }+U_{\varphi }\left({\displaystyle \frac{B_r}{r}+\frac{cot\theta B_{\theta }}{r}}\right)-B_{\varphi }\left({\displaystyle \frac{U_r}{r}+\frac{cot\theta U_{\theta }}{r}}\right)=}\hfill \\ \hfill \\ \eta \left[{\displaystyle \frac{{\partial }^2{B}_{\varphi }}{\partial r^2}+\frac{2}{r}\frac{\partial B_{\varphi }}{\partial r}+\frac{1}{r^2}\frac{{\partial }^2{B}_{\varphi }}{\partial {\theta }^2}+\frac{cot\theta }{r^2}\frac{\partial B_{\varphi }}{\partial \theta }-\frac{B_{\varphi }}{r^2{sin}^2\theta }}\right].\hfill \end{array}\end{array}$$(B.2)

Since in our simulations the Elsasser number $\mathrm{\Lambda }=B_0^2/{\mu }_0{\rho }_0{\mathrm{\Omega }}_0\eta$ is ≫1, the Coriolis term becomes negligible with respect to the Lorentz force (Acheson & Hide 1973), then, assuming U_r ≪ U_ϕ and conserving the highest radial derivatives in each term, we get the Hartmann boundary layer equations (see also reviews by Roberts (1967), Acheson & Hide (1973), Dormy et al. (1998)):

$$\begin{array}{c}\hfill {\displaystyle \frac{B_r}{{\mu }_0{\rho }_0}\frac{\partial B_{\varphi }^{\mathrm{H}}}{\partial r}=-\nu \frac{{\partial }^2{U}_{\varphi }^{\mathrm{H}}}{\partial r^2},B_r\frac{\partial U_{\varphi }^{\mathrm{H}}}{\partial r}=-\eta \frac{{\partial }^2{B}_{\varphi }^{\mathrm{H}}}{\partial r^2},}\end{array}$$(B.3)

where the index H denotes a boundary layer flow. As B_r(r₀, θ) = − B₀ cos θ (see Eq. (22)), then after introducing the stretched coordinate $\xi =(B_0{r}_0/\sqrt{{\mu }_0{\rho }_0\eta \nu })(r_0-r)$, Sys. (B.3) is rewritten as:

$$\begin{array}{c}\hfill {\displaystyle \frac{cos\theta }{\sqrt{{\mu }_0{\rho }_0}}\frac{\partial B_{\varphi }^{\mathrm{H}}}{\partial \xi }=-r_0\sqrt{{\displaystyle \frac{\nu }{\eta }}}\frac{{\partial }^2{U}_{\varphi }^{\mathrm{H}}}{\partial {\xi }^2},\sqrt{{\mu }_0{\rho }_0}cos\theta \frac{\partial U_{\varphi }^{\mathrm{H}}}{\partial \xi }=-r_0\sqrt{{\displaystyle \frac{\eta }{\nu }}}\frac{{\partial }^2{B}_{\varphi }^{\mathrm{H}}}{\partial {\xi }^2}.}\end{array}$$(B.4)

By deriving the second of these two equations, $U_{\varphi }^{\mathrm{H}}$ can be eliminated to yield:

$$\begin{array}{c}\hfill {\displaystyle \frac{{\partial }^3{B}_{\varphi }^{\mathrm{H}}}{\partial {\xi }^3}-\frac{{cos}^2\theta }{r_0^2}\frac{\partial B_{\varphi }^{\mathrm{H}}}{\partial \xi }=0,}\end{array}$$(B.5)

whose the solution is

$$\begin{array}{c}\hfill {\displaystyle B_{\varphi }^{\mathrm{H}}(\xi ,\theta )=C_1+\frac{C_2{r}_0}{cos\theta }exp\left({\displaystyle -cos\theta \frac{\xi }{r_0}}\right)+\frac{C_3{r}_0}{cos\theta }exp\left({\displaystyle cos\theta \frac{\xi }{r_0}}\right).}\end{array}$$(B.6)

From the evanescent condition when ξ → ∞, C₁ = C₃ = 0, and from the vacuum condition at the outer sphere $B_{\varphi }^{\mathrm{H}}(0,\theta )+B_{\varphi }^{\mathrm{I}}(r_0,\theta )=0$, $C_2=-(cos\theta /r_0)B_{\varphi }^{\mathrm{I}}(r_0,\theta )$, hence:

$$\begin{array}{c}\hfill B_{\varphi }^{\mathrm{H}}(\xi ,\theta )=-B_{\varphi }^{\mathrm{I}}(r_0,\theta )exp\left({\displaystyle -cos\theta \frac{\xi }{r_0}}\right).\end{array}$$(B.7)

This time, the index I stands for interior flow. We now integrate the second equation of Sys. (B.4):

$$\begin{array}{c}\hfill {\displaystyle U_{\varphi }^{\mathrm{H}}(\xi ,\theta )=-\frac{1}{\sqrt{{\mu }_0{\rho }_0}}\sqrt{{\displaystyle \frac{\eta }{\nu }}}{B}_{\varphi }^{\mathrm{I}}(r_0,\theta )exp\left({\displaystyle -cos\theta \frac{\xi }{r_0}}\right)+C_4.}\end{array}$$(B.8)

Again, from the evanescent equation when ξ → ∞, we readily get C₄ = 0. Then, as $U_{\varphi }^{\mathrm{H}}(\xi ,\theta )+U_{\varphi }^{\mathrm{I}}(r_0,\theta )=0$ at the outer sphere, we obtain the following relationship between the interior flows:

$$\begin{array}{c}\hfill {\displaystyle U_{\varphi }^{\mathrm{I}}(r_0,\theta )=\frac{1}{\sqrt{{\mu }_0{\rho }_0}}\sqrt{{\displaystyle \frac{\eta }{\nu }}}{B}_{\varphi }^{\mathrm{I}}(r_0,\theta ),}\end{array}$$(B.9)

and,

$$\begin{array}{c}\hfill U_{\varphi }^{\mathrm{H}}(\xi ,\theta )=-U_{\varphi }^{\mathrm{I}}(r_0,\theta )exp\left({\displaystyle -cos\theta \frac{\xi }{r_0}}\right).\end{array}$$(B.10)

Thus, using the estimate Eq. (36), we deduce from Eq. (B.9) the order of magnitude of the jump on the differential rotation across the Hartmann layer

$$\begin{array}{c}\hfill {\displaystyle U_{\varphi }^{\mathrm{I}}(r_0,\theta )\approx \mathcal{O}\left({\displaystyle \frac{r_0{\mu }_0{\rho }_0{V}_0{\mathrm{\Omega }}_0}{\sqrt{{\mu }_0{\rho }_0}{B}_0}\sqrt{{\displaystyle \frac{\eta }{\nu }}}}\right)\Rightarrow \frac{\mathrm{\Delta }{\mathrm{\Omega }}^{\mathrm{I}}(r_0,\theta )}{{\mathrm{\Omega }}_0}\approx \mathcal{O}\left({\displaystyle \frac{\sqrt{{\mu }_0{\rho }_0}{V}_0}{B_0}\sqrt{{\displaystyle \frac{\eta }{\nu }}}}\right)=\mathcal{O}\left({\displaystyle \frac{\sqrt{P_m}R{e}_c}{L_u}}\right).}\end{array}$$(B.11)

In Fig. B.1, we plotted the normalised toroidal field B_ϕ/B₀ (left panel) and the normalised differential rotation δΩ/Ω₀ (right panel) as a function of the stretched coordinate ξ at the particular location θ = π/4. This was done for runs obtained at various L_u ranging from 5 ⋅ 10³ to 5 ⋅ 10⁴, both for Re_c = 1 and 5 (runs D3 to D7 of Table 1). These quantities have been rescaled with their characteristic amplitude as given by Eqs. (36) and (B.11). We can observe that the different curves overlap and are of 𝒪(1) after rescaling. This enables us to conclude that the 𝒪(r₀μ₀ρ₀V₀Ω₀/B₀) jump on the toroidal magnetic field at the outer sphere induces a $\mathcal{O}(r_0\sqrt{{\mu }_0{\rho }_0}{V}_0{\mathrm{\Omega }}_0\sqrt{\eta }/B_0\sqrt{\nu })$ jump on the azimuthal component of the velocity field across the Hartmann layer or equivalently, a $\mathcal{O}(\sqrt{{\mu }_0{\rho }_0}{V}_0\sqrt{\eta }/B_0\sqrt{\nu })$ jump on the normalised differential rotation. As a result, the quasi-solid rotation region is in differential rotation with the outer sphere and the characteristic amplitude of this differential rotation is given by Eq. (B.11).

Appendix C: Approximate solutions of differential rotation in the dead zone

C.1. Case 1 - No effect of the contraction on the field lines

This appendix is intended to solve in the DZ the viscous balance Eq. (35), whose dimensionless form reads

$$\begin{array}{c}\hfill {\displaystyle \stackrel{\sim }{r}\frac{{\partial }^2\delta \stackrel{\sim }{\mathrm{\Omega }}}{\partial r^2}+4\frac{\partial \delta \stackrel{\sim }{\mathrm{\Omega }}}{\partial r}+\frac{1}{\stackrel{\sim }{r}}\frac{{\partial }^2\delta \stackrel{\sim }{\mathrm{\Omega }}}{\partial {\theta }^2}+\frac{3cot\theta }{\stackrel{\sim }{r}}\frac{\partial \delta \stackrel{\sim }{\mathrm{\Omega }}}{\partial \theta }=\frac{-2}{E{\stackrel{\sim }{r}}^2}.}\end{array}$$(C.1)

To do so, the DZ is assimilated to a conical domain as sketched in Fig. C.1. This domain, denoted 𝒟, is such that:

$$\begin{array}{c}\hfill \mathcal{D}=\{r_{i_{\mathrm{DZ}}}\le r\le r_0;-{\theta }_0\le \theta \le {\theta }_0\},\end{array}$$(C.2)

where r_iDZ = 0.77r₀ and the latitude θ₀ is ≈π/14. In addition, we assume that the term $(3cot\theta /\stackrel{\sim }{r})\partial \delta \stackrel{\sim }{\mathrm{\Omega }}/\partial \theta$ is negligible as compared to the others. This assumption has been verified a posteriori. Thus, Eq. (C.1) is rewritten as follows:

$$\begin{array}{c}\hfill {\displaystyle {\stackrel{\sim }{r}}^3\frac{{\partial }^2\delta \stackrel{\sim }{\mathrm{\Omega }}}{\partial r^2}+4{\stackrel{\sim }{r}}^2\frac{\partial \delta \stackrel{\sim }{\mathrm{\Omega }}}{\partial r}+\stackrel{\sim }{r}\frac{{\partial }^2\delta \stackrel{\sim }{\mathrm{\Omega }}}{\partial {\theta }^2}=\frac{-2}{E}.}\end{array}$$(C.3)

Outside the domain 𝒟, the flow is assumed to be in solid rotation and symmetrical with respect to the equatorial plane so that we adopt the following homogeneous boundary conditions:

$$\begin{array}{c}\hfill \delta \mathrm{\Omega }(r_0,\theta )=\delta \mathrm{\Omega }(r_{i_{\mathrm{DZ}}},\theta )=\delta \mathrm{\Omega }(r,{\theta }_0)={{\displaystyle \frac{\partial \delta \mathrm{\Omega }(r,\theta )}{\partial \theta }}|}_{\theta =0}=0.\end{array}$$(C.4)

The method for solving Eq. (C.3) was excerpted from a lecture of the Ira A. Fulton College of Engineering and Technology at Brigham Young University ². The broad lines are now given. We first build a set of basis functions that will be used to express the solution. These are obtained by looking for separable solutions of the eigenvalue problem

$$\begin{array}{c}\hfill {\displaystyle {\stackrel{\sim }{r}}^3\frac{{\partial }^2\delta {\stackrel{\sim }{\mathrm{\Omega }}}_{\mathrm{h}}}{\partial r^2}+4{\stackrel{\sim }{r}}^2\frac{\partial \delta {\stackrel{\sim }{\mathrm{\Omega }}}_{\mathrm{h}}}{\partial r}+\stackrel{\sim }{r}\frac{{\partial }^2\delta {\stackrel{\sim }{\mathrm{\Omega }}}_{\mathrm{h}}}{\partial {\theta }^2}=\stackrel{\sim }{\lambda }\stackrel{\sim }{r}\delta {\stackrel{\sim }{\mathrm{\Omega }}}_{\mathrm{h}},}\end{array}$$(C.5)

satisfying the boundary conditions Eq. (C.4). We finally determine a general solution using the orthogonality properties of the basis functions.

Separable solutions, $\delta {\stackrel{\sim }{\mathrm{\Omega }}}_{\mathrm{h}}=\stackrel{\sim }{g}(\theta )\stackrel{\sim }{f}(r)$, of Eq. (C.5) must verify:

$$\begin{array}{c}\hfill {\displaystyle {\stackrel{\sim }{r}}^3{\stackrel{\sim }{f}}^{″}(r)+4{\stackrel{\sim }{r}}^2{\stackrel{\sim }{f}}^{\prime }(r)+\stackrel{\sim }{r}{\stackrel{\sim }{g}}^{″}(\theta )=\stackrel{\sim }{\lambda }\stackrel{\sim }{r}\stackrel{\sim }{g}(\theta )\stackrel{\sim }{f}(r)⟹\frac{{\stackrel{\sim }{r}}^2{\stackrel{\sim }{f}}^{″}(r)+4\stackrel{\sim }{r}{\stackrel{\sim }{f}}^{\prime }(\stackrel{\sim }{r})}{\stackrel{\sim }{f}(r)}+\frac{{\stackrel{\sim }{g}}^{″}(\theta )}{\stackrel{\sim }{g}(\theta )}=\stackrel{\sim }{\lambda },}\end{array}$$(C.6)

with the boundary conditions

$$\begin{array}{c}\hfill g({\theta }_0)=g^{\prime }(0)=0\mathrm{and},f(r_0)=f(r_{i_{\mathrm{DZ}}})=0.\end{array}$$(C.7)

The problem thus reduces to the solving the two following sub-eigenvalue problems:

$$\begin{array}{c}\hfill {\stackrel{\sim }{g}}^{″}(\theta )-\stackrel{\sim }{\nu }\stackrel{\sim }{g}(\theta )=0\mathrm{and},{\stackrel{\sim }{r}}^2{\stackrel{\sim }{f}}^{″}(r)+4\stackrel{\sim }{r}{\stackrel{\sim }{f}}^{\prime }(r)-\stackrel{\sim }{\mu }\stackrel{\sim }{f}(r)=0\mathrm{with}\stackrel{\sim }{\nu }+\stackrel{\sim }{\mu }=\stackrel{\sim }{\lambda }.\end{array}$$(C.8)

The differential equation in θ is a classical Sturm-Liouville problem while the differential equation in r is known as the Euler’s problem. We first deal with the Sturm-Liouville problem.

In order to avoid unphysical solutions, the eigenvalues must be strictly negative. The solution then reads

$$\begin{array}{c}\hfill {\stackrel{\sim }{g}}_k(\theta )={\stackrel{\sim }{A}}_{k}cos\left(\sqrt{\left|{\stackrel{\sim }{\nu }}_k\right|}\theta \right)+{\stackrel{\sim }{B}}_{k}sin\left(\sqrt{\left|{\stackrel{\sim }{\nu }}_k\right|}\theta \right).\end{array}$$(C.9)

Its latitudinal derivative is readily

$$\begin{array}{c}\hfill {\stackrel{\sim }{g}}_k^{\prime }(\theta )=-{\stackrel{\sim }{A}}_k\sqrt{\left|\stackrel{\sim }{{\nu }_k}\right|}sin\left(\sqrt{\left|\stackrel{\sim }{{\nu }_k}\right|}\theta \right)+{\stackrel{\sim }{B}}_k\sqrt{\left|\stackrel{\sim }{{\nu }_k}\right|}cos\left(\sqrt{\left|\stackrel{\sim }{{\nu }_k}\right|}\theta \right).\end{array}$$(C.10)

From the condition of symmetry, ${\stackrel{\sim }{B}}_k=0$ and ${\stackrel{\sim }{g}}_k(\theta )={\stackrel{\sim }{A}}_{k}cos\left(\sqrt{\left|{\stackrel{\sim }{\nu }}_k\right|}\theta \right)$. Besides, ${\stackrel{\sim }{g}}_k({\theta }_0)={\stackrel{\sim }{A}}_{k}cos\left(\sqrt{\left|{\stackrel{\sim }{\nu }}_k\right|}{\theta }_0\right)=0$. Since we are looking for a non-trivial solution (i.e. ${\stackrel{\sim }{A}}_k\ne 0$) we obtain

$$\begin{array}{c}\hfill {\displaystyle \sqrt{\left|{\stackrel{\sim }{\nu }}_k\right|}{\theta }_0=k\pi -\frac{\pi }{2}\mathrm{with}k\in \mathbb{Z},}\end{array}$$(C.11)

thus concluding that

$$\begin{array}{c}\hfill {\stackrel{\sim }{\nu }}_k=-{\left({\displaystyle \frac{(2k-1)\pi }{2{\theta }_0}}\right)}^2\mathrm{with}k\in \mathbb{Z},\end{array}$$(C.12)

are the sought eigenvalues and are associated with the eigenfunctions

$$\begin{array}{c}\hfill {\stackrel{\sim }{g}}_k(\theta )={\stackrel{\sim }{A}}_{k}cos\left(\left({\displaystyle \frac{(2k-1)\pi }{2{\theta }_0}}\right)\theta \right)\mathrm{with}k\in \mathbb{Z}.\end{array}$$(C.13)

We are now going to solve the Euler’s problem. Here the eigenvalues must be <  − 9/4 to avoid unphysical solutions. In that case, the solution reads:

$$\begin{array}{c}\hfill {\stackrel{\sim }{f}}_n(r)={\stackrel{\sim }{C}}_n{\stackrel{\sim }{r}}^{{\displaystyle -(3+i\sqrt{|9+4{\stackrel{\sim }{\mu }}_n|})/2}}+{\stackrel{\sim }{D}}_n{\stackrel{\sim }{r}}^{{\displaystyle -(3-i\sqrt{|9+4{\stackrel{\sim }{\mu }}_n|})/2}},\end{array}$$(C.14)

which can be rewritten as

$$\begin{array}{c}\hfill {\displaystyle {\stackrel{\sim }{f}}_n(r)=\frac{{\stackrel{\sim }{C}}_n}{{\stackrel{\sim }{r}}^{3/2}}cos\left({\displaystyle \frac{1}{2}\sqrt{|9+4{\stackrel{\sim }{\mu }}_n|}ln\left(\stackrel{\sim }{r}\right)}\right)+\frac{{\stackrel{\sim }{D}}_n}{{\stackrel{\sim }{r}}^{3/2}}sin\left({\displaystyle \frac{1}{2}\sqrt{|9+4{\stackrel{\sim }{\mu }}_n|}ln\left(\stackrel{\sim }{r}\right)}\right).}\end{array}$$(C.15)

From the condition ${\stackrel{\sim }{f}}_n(1)=0$ we have ${\stackrel{\sim }{C}}_n=0$, and

$$\begin{array}{c}\hfill {\displaystyle {\stackrel{\sim }{f}}_n(r)=\frac{{\stackrel{\sim }{D}}_n}{{\stackrel{\sim }{r}}^{3/2}}sin\left({\displaystyle \frac{1}{2}\sqrt{|9+4{\stackrel{\sim }{\mu }}_n|}ln(\stackrel{\sim }{r})}\right).}\end{array}$$(C.16)

The second boundary condition leads to

$$\begin{array}{c}\hfill {\displaystyle \frac{{\stackrel{\sim }{D}}_n}{{\stackrel{\sim }{r}}_{i_{\mathrm{DZ}}}^{3/2}}sin\left({\displaystyle \frac{1}{2}\sqrt{|9+4{\stackrel{\sim }{\mu }}_n|}ln\left({\stackrel{\sim }{r}}_{i_{\mathrm{DZ}}}\right)}\right)=0.}\end{array}$$(C.17)

Excluding the non-trivial solution ${\stackrel{\sim }{D}}_n=0$ we obtain

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{2}\sqrt{|9+4{\stackrel{\sim }{\mu }}_n|}ln\left({\stackrel{\sim }{r}}_{i_{\mathrm{DZ}}}\right)=n\pi \mathrm{with}n\in \mathbb{Z}.}\end{array}$$(C.18)

The sought eigenvalues are thus

$$\begin{array}{c}\hfill {\stackrel{\sim }{\mu }}_n=-{\left({\displaystyle \frac{3}{2}}\right)}^2-{\left({\displaystyle \frac{n\pi }{ln\left({\stackrel{\sim }{r}}_{i_{\mathrm{DZ}}}\right)}}\right)}^2\mathrm{with}n\in \mathbb{Z},\end{array}$$(C.19)

and are associated with the eigenfunctions

$$\begin{array}{c}\hfill {\displaystyle {\stackrel{\sim }{f}}_n(r)=\frac{{\stackrel{\sim }{D}}_n}{{\stackrel{\sim }{r}}^{3/2}}sin\left({\displaystyle \frac{n\pi ln\left(\stackrel{\sim }{r}\right)}{ln\left({\stackrel{\sim }{r}}_{i_{\mathrm{DZ}}}\right)}}\right)\mathrm{with}n\in \mathbb{Z}.}\end{array}$$(C.20)

The general solution δΩ(r, θ) is now expanded over the basis of the eigenfunctions satisfying the boundary conditions Eq. (C.7):

$$\begin{array}{c}\hfill {\displaystyle \delta \stackrel{\sim }{\mathrm{\Omega }}(r,\theta )={\displaystyle \sum _{n=1}^{\infty }\sum _{k=1}^{\infty }\frac{{\stackrel{\sim }{A}}_{\mathit{nk}}}{{\stackrel{\sim }{r}}^{3/2}}sin(\left({\displaystyle \frac{n\pi }{ln\left({\stackrel{\sim }{r}}_{i_{\mathrm{DZ}}}\right)}}\right)ln\left(\stackrel{\sim }{r}\right))cos\left(\left({\displaystyle \frac{(2k-1)\pi }{2{\theta }_0}}\right)\theta \right)\mathrm{with}n,k\in \mathbb{Z}.}}\end{array}$$(C.21)

Thus, this solution verifies the boundary conditions Eq. (C.4). We now need to find the A_nk coefficients. These are determined using the orthogonality properties of the eigenfunctions. After determining the different partial derivatives of the function δΩ(r, θ) and after substituting their expressions in Eq. (C.3), we obtain the following relationship:

$$\begin{array}{c}\hfill {\displaystyle \sum _{n=1}^{\infty }\sum _{k=1}^{\infty }\frac{{\stackrel{\sim }{A}}_{\mathit{nk}}}{\sqrt{\stackrel{\sim }{r}}}[-{\left({\displaystyle \frac{3}{2}}\right)}^2-{\left({\displaystyle \frac{n\pi }{ln\left({\stackrel{\sim }{r}}_{i_{\mathrm{DZ}}}\right)}}\right)}^2-{\left({\displaystyle \frac{(2k-1)\pi }{2{\theta }_0}}\right)}^2]sin(\left({\displaystyle \frac{n\pi }{ln\left({\stackrel{\sim }{r}}_{i_{\mathrm{DZ}}}\right)}}\right)ln\left(\stackrel{\sim }{r}\right))cos\left(\left({\displaystyle \frac{(2k-1)\pi }{2{\theta }_0}}\right)\theta \right)=\frac{-2}{E}.}\end{array}$$(C.22)

As

$$\begin{array}{c}\hfill {\stackrel{\sim }{\lambda }}_{\mathit{nk}}=\left[\underset{\left|{\stackrel{\sim }{\mu }}_n\right|}{\underbrace{-{\left({\displaystyle \frac{3}{2}}\right)}^2-{\left({\displaystyle \frac{n\pi }{ln\left({\stackrel{\sim }{r}}_{i_{\mathrm{DZ}}}\right)}}\right)}^2}}\underset{\left|{\stackrel{\sim }{\nu }}_k\right|}{\underbrace{-{\left({\displaystyle \frac{(2k-1)\pi }{2{\theta }_0}}\right)}^2}}\right]\mathrm{with}n,k\in \mathbb{Z},\end{array}$$(C.23)

Eq. (C.22) is rewritten as follows

$$\begin{array}{c}\hfill {\displaystyle \sum _{n=1}^{\infty }\sum _{k=1}^{\infty }\frac{{\stackrel{\sim }{A}}_{\mathit{nk}}{\stackrel{\sim }{\lambda }}_{\mathit{nk}}}{{\stackrel{\sim }{r}}^{3/2}}sin(\left({\displaystyle \frac{n\pi }{ln\left({\stackrel{\sim }{r}}_{i_{\mathrm{DZ}}}\right)}}\right)ln\left(\stackrel{\sim }{r}\right))cos\left(\left({\displaystyle \frac{(2k-1)\pi }{2{\theta }_0}}\right)\theta \right)=\frac{-2}{\stackrel{\sim }{r}E}.}\end{array}$$(C.24)

Multiplying each members of this relationship by ${\stackrel{\sim }{r}}^{-1/2}sin(\left({\displaystyle \frac{n\pi }{ln\left({\stackrel{\sim }{r}}_{i_{\mathrm{DZ}}}\right)}}\right)ln\left(\stackrel{\sim }{r}\right))cos\left(\left({\displaystyle \frac{(2k-1)\pi }{2{\theta }_0}}\right)\theta \right)$, then integrating over the domain 𝒟 the resulting relationship, we obtain the coefficients A_nk:

$$\begin{array}{c}\hfill {\stackrel{\sim }{A}}_{\mathit{nk}}=\left({\displaystyle \frac{4{\left(n\pi \right)}^2+{(ln\left({\stackrel{\sim }{r}}_{i_{\mathrm{DZ}}}\right))}^2}{E{\stackrel{\sim }{\lambda }}_{\mathit{nk}}({\stackrel{\sim }{r}}_{i_{\mathrm{DZ}}}-1){\left(n\pi \right)}^2{\theta }_0}}\right)\left({\displaystyle \frac{-4{\theta }_0}{(2k-1)\pi }{(-1)}^k}\right)\left({\displaystyle \frac{4n\pi ln\left({\stackrel{\sim }{r}}_{i_{\mathrm{DZ}}}\right)}{4{\left(n\pi \right)}^2+9{(ln\left({\stackrel{\sim }{r}}_{i_{\mathrm{DZ}}}\right))}^2}[{(-1)}^n{\stackrel{\sim }{r}}_{i_{\mathrm{DZ}}}^{3/2}-1]}\right)\mathrm{with}n,k\in \mathbb{Z},\end{array}$$(C.25)

and the final solution, which reads under dimensional form

$$\begin{array}{cc}\hfill {\displaystyle \frac{\delta \mathrm{\Omega }(r,\theta )}{{\mathrm{\Omega }}_0}}& =R{e}_c{\displaystyle \sum _{n=1}^{\infty }\sum _{k=1}^{\infty }[\left({\displaystyle \frac{{\left(2n\pi \right)}^2+{(ln\left({\displaystyle \frac{r_{i_{\text{DZ}}}}{r_0}}\right))}^2}{{\lambda }_{\mathit{nk}}\left({\displaystyle \frac{r_{i_{\text{DZ}}}}{r_0}-1}\right){\left(n\pi \right)}^2{\theta }_0}}\right)\left({\displaystyle \frac{-4{\theta }_0}{(2k-1)\pi }{(-1)}^k}\right)\left({\displaystyle \frac{4n\pi ln\left({\displaystyle \frac{r_{i_{\text{DZ}}}}{r_0}}\right)}{{\left(2n\pi \right)}^2+{(3ln\left({\displaystyle \frac{r_{i_{\text{DZ}}}}{r_0}}\right))}^2}({(-1)}^n{\left({\displaystyle \frac{r_{i_{\text{DZ}}}}{r_0}}\right)}^{3/2}-1)}\right)}\hfill \\ \hfill & {\left({\displaystyle \frac{r_0}{r}}\right)}^{3/2}sin(\left({\displaystyle \frac{n\pi }{ln\left({\displaystyle \frac{r_{i_{\text{DZ}}}}{r_0}}\right)}}\right)ln\left({\displaystyle \frac{r}{r_0}}\right))cos\left(\left({\displaystyle \frac{(2k-1)\pi }{2{\theta }_0}}\right)\theta \right)]\mathrm{with}n,k\in \mathbb{Z}.\hfill \end{array}$$(C.26)

C.2. Case 2 - Advection of the field lines by the contraction

When the contraction term is introduced in the induction equation, the poloidal field lines are advected and the DZ connects to the inner sphere (see Fig. C.2). This modifies the boundary conditions Eq. (C.4) since at the inner sphere, the azimuthal component of the velocity field now satisfies a stress-free condition hence:

$$\begin{array}{c}\hfill \delta \mathrm{\Omega }(r_0,\theta )={{\displaystyle \frac{\partial \delta \mathrm{\Omega }(r,\theta )}{\partial r}}|}_{r=r_i}=\delta \mathrm{\Omega }(r,{\theta }_0)={{\displaystyle \frac{\partial \delta \mathrm{\Omega }(r,\theta )}{\partial \theta }}|}_{\theta =0}=0.\end{array}$$(C.27)

Thus, the solution of the Sturm-Liouville problem is unchanged and the sought eigenvalues are again Eq. (C.12) and are associated with the eigenfunctions Eq. (C.13). For the Euler’s problem, the eigenvalues must still be < 9/4 to avoid unphysical solutions. By taking Eq. (C.15), and after using the condition ${\stackrel{\sim }{f}}_n(1)=0$, we have

$$\begin{array}{c}\hfill {\displaystyle {\stackrel{\sim }{f}}_n(r)=\frac{{\stackrel{\sim }{D}}_n}{{\stackrel{\sim }{r}}^{3/2}}sin\left({\displaystyle \frac{1}{2}\sqrt{|9+4{\stackrel{\sim }{\mu }}_n|}ln\left(\stackrel{\sim }{r}\right)}\right).}\end{array}$$(C.28)

By deriving this expression to apply the stress-free condition to it, we get the transcendental equation

$$\begin{array}{c}\hfill 3tan\left({\displaystyle \frac{1}{2}\sqrt{|9+4{\stackrel{\sim }{\mu }}_n|}ln({\stackrel{\sim }{r}}_i)}\right)-\sqrt{|9+4{\stackrel{\sim }{\mu }}_n|}=0.\end{array}$$(C.29)

After numerically solving this equation, we obtain the eigenvalues μ_n. The general solution δΩ(r, θ) of Eq. (C.3) is then expanded on the basis of the eigenfunctions satisfying the boundary conditions Eq. (C.27) on the domain 𝒟

$$\begin{array}{c}\hfill {\displaystyle \delta \stackrel{\sim }{\mathrm{\Omega }}(r,\theta )={\displaystyle \sum _{n=1}^{\infty }\sum _{k=1}^{\infty }\frac{{\stackrel{\sim }{A}}_{\mathit{nk}}}{{\stackrel{\sim }{r}}^{3/2}}sin\left({\displaystyle \frac{1}{2}\sqrt{|9+4{\stackrel{\sim }{\mu }}_n|}ln\left(\stackrel{\sim }{r}\right)}\right)cos\left(\left({\displaystyle \frac{(2k-1)\pi }{2{\theta }_0}}\right)\theta \right)\mathrm{with}n,k\in \mathbb{Z}.}}\end{array}$$(C.30)

As previously, after determining the different partial derivatives of the function δΩ(r, θ), and after substituting them in Eq. (C.3), we get:

$$\begin{array}{c}\hfill {\displaystyle \sum _{n=1}^{\infty }\sum _{k=1}^{\infty }\frac{{\stackrel{\sim }{A}}_{\mathit{nk}}}{\sqrt{\stackrel{\sim }{r}}}\left[\underset{{\stackrel{\sim }{\lambda }}_{\mathit{nk}}}{\underbrace{{\stackrel{\sim }{\mu }}_n-{\left({\displaystyle \frac{(2k-1)\pi }{2{\theta }_0}}\right)}^2}}\right]sin\left({\displaystyle \frac{1}{2}\sqrt{|9+4{\stackrel{\sim }{\mu }}_n|}ln\left(\stackrel{\sim }{r}\right)}\right)cos\left(\left({\displaystyle \frac{(2k-1)\pi }{2{\theta }_0}}\right)\theta \right)=\frac{-2}{E}.}\end{array}$$(C.31)

Again, using the orthogonality properties enables us to determine the A_nk coefficients

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle {\stackrel{\sim }{A}}_{\mathit{nk}}=\frac{[2\sqrt{|9+4{\stackrel{\sim }{\mu }}_n|}({\stackrel{\sim }{r}}_i^{3/2}cos\left({\displaystyle \frac{1}{2}\sqrt{|9+4{\stackrel{\sim }{\mu }}_n|}ln({\stackrel{\sim }{r}}_i)}\right)-1)-6{\stackrel{\sim }{r}}_i^{3/2}sin\left({\displaystyle \frac{1}{2}\sqrt{|9+4{\stackrel{\sim }{\mu }}_n|}ln({\stackrel{\sim }{r}}_i)}\right)]}{[|9+4{\stackrel{\sim }{\mu }}_n|(1-{\stackrel{\sim }{r}}_i)+{\stackrel{\sim }{r}}_i(cos(\sqrt{|9+4{\stackrel{\sim }{\mu }}_n|}ln({\stackrel{\sim }{r}}_i))+\sqrt{|9+4{\stackrel{\sim }{\mu }}_n|}sin(\sqrt{|9+4{\stackrel{\sim }{\mu }}_n|}ln({\stackrel{\sim }{r}}_i))-1)]}·}\hfill \\ \hfill \\ {\displaystyle \frac{16(1+|9+4{\stackrel{\sim }{\mu }}_n|){(-1)}^k}{E{\stackrel{\sim }{\lambda }}_{\mathit{nk}}(2k-1)(9+|9+4{\stackrel{\sim }{\mu }}_n|)\pi },}\hfill \end{array}\end{array}$$(C.32)

and so the final solution which, under dimensional form, reads as follows:

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle \frac{\delta \mathrm{\Omega }(r,\theta )}{{\mathrm{\Omega }}_0}=R{e}_c{\displaystyle {\sum }_{n=1}^{\infty }{\sum }_{k=1}^{\infty }\frac{[2\sqrt{|9+4{\mu }_n|}({\left({\displaystyle \frac{r_i}{r_0}}\right)}^{3/2}cos\left({\displaystyle \frac{1}{2}\sqrt{|9+4{\mu }_n|}ln\left({\displaystyle \frac{r_i}{r_0}}\right)}\right)-1)-6{\left({\displaystyle \frac{r_i}{r_0}}\right)}^{3/2}sin\left({\displaystyle \frac{1}{2}\sqrt{|9+4{\mu }_n|}ln\left({\displaystyle \frac{r_i}{r_0}}\right)}\right)]}{\left[{\displaystyle |9+4{\mu }_n|\left({\displaystyle 1-\frac{r_i}{r_0}}\right)+\frac{r_i}{r_0}(cos(\sqrt{|9+4{\mu }_n|}ln\left({\displaystyle \frac{r_i}{r_0}}\right))+\sqrt{|9+4{\mu }_n|}sin(\sqrt{|9+4{\mu }_n|}ln\left({\displaystyle \frac{r_i}{r_0}}\right))-1)}\right]}}}\hfill \\ \hfill \\ {\displaystyle \frac{16(1+|9+4{\mu }_n|){(-1)}^k}{{\lambda }_{\mathit{nk}}(2k-1)(9+|9+4{\mu }_n|)\pi }{\left(\frac{r_0}{r}\right)}^{3/2}sin\left({\displaystyle \frac{1}{2}\sqrt{|9+4{\mu }_n|}ln\left({\displaystyle \frac{r}{r_0}}\right)}\right)cos\left(\left({\displaystyle \frac{(2k-1)\pi }{2{\theta }_0}}\right)\theta \right),}\hfill \end{array}\end{array}$$(C.33)

where n, k ∈ ℤ.

All Tables

All Figures