© L. Jouve et al. 2020

1. Introduction

Considerable progress has recently been made in the understanding of magnetic fields at the surface of stars, mostly thanks to the ground-based instruments NARVAL at the Pic du Midi observatory in France and ESPaDOnS at the Mauna Kea Observatory in Hawaï. It has been known for more than a century now that the Sun harbors a strong magnetic field which manifests itself as spots popping-up at the solar surface (Hale 1908). Currently, there is also a general consensus on the fact that this magnetic field is produced by dynamo action inside the convective envelope of the Sun and that such a process should be quite general for all solar-like stars (Parker 1955; Moffatt 1978). The magnetism of intermediate-mass and massive stars has also been thoroughly investigated. It is however expected to differ strongly from that of cool stars because of the presence of the outer radiative zone. Indeed, if a convective dynamo is at play in the core of hot stars, it might be more difficult for the magnetic field created in the convective core to travel all the way to the surface so that observers from Earth could see it. In intermediate-mass and massive stars, the magnetism is indeed quite different from what is observed on cool stars: 5–10% of these stars do exhibit a strong surface magnetic field above 300 G and these stars are also the ones which show chemical peculiarities in their spectra (Ap/Bp stars). Thanks to recent spectropolarimetric observations, detections of a much smaller amplitude field (at the sub-Gauss level) have been obtained on stars like Vega, Sirius A, Alhena, β-Uma or θ-Leo (Lignières et al. 2009; Blazère et al. 2016a,b), leading to the idea that two classes of magnetism could exist in intermediate-mass and massive stars: the strong dipolar field of Ap/Bp stars and the ultra-weak Vega-like magnetic field. A sound explanation for the existence of these two types of magnetism and the absence of stars possessing fields with amplitudes between approximately 1 G and 300 G is still lacking. A possible scenario was proposed by Aurière et al. (2007), relying on the existence of a magnetic instability which could develop only for weak enough dipolar fields and which would lead to the disruption of the axisymmetric magnetic configuration. In that scenario, a crucial role is given to the differential rotation which acts on the dipolar magnetic field to produce a configuration dominated by the toroidal component. This toroidal field is very likely to become unstable to a magnetohydrodynamical (MHD) instability. If such an instability existed, not only would it possibly explain the minimum field of Ap/Bp stars, but it could also be at the origin of dynamo action in the radiative zones of Vega-like stars. Various studies have indeed recently focused on the appealing idea that dynamo action would not require convective motions but, rather, non-axisymmetric hydro or MHD instabilities which, in conjunction with the differential rotation, would produce the electromotive force needed to close the dynamo loop (Spruit 2002; Braithwaite 2006; Zahn et al. 2007; Guervilly & Cardin 2010; Marcotte & Gissinger 2016). For now, it is still debated whether such a radiative zone dynamo could exist in stars.

The interplay between differential rotation and magnetic fields, which is at the heart of the Aurière et al. (2007) explanation of the magnetism of hot stars, is also invoked to interpret the recent asteroseismic observations of more than 300 red giants provided by the Kepler satellite in the last decade. Indeed, in those stars, the radiative zone contracts below the H-burning shell and expands above, naturally leading to a spin-up of the innermost regions and a braking of the layers above. This is indeed what is observed, a differential rotation is established between the inner and outer shells in these stars because of the contraction or expansion phenomena (Deheuvels et al. 2012, 2014). However, simple models assuming conservation of angular momentum considerably overestimate the level of differential rotation produced. More puzzling is the fact that even sophisticated stellar evolution models including the rotationally-induced transport of angular momentum fail at reproducing the observations (Eggenberger et al. 2012a,b; Ceillier et al. 2013; Marques et al. 2013). In this case, a more efficient transport of angular momentum seems to be at play in those stellar radiative zones and magnetic fields are seriously considered as possible candidates for this role. In particular, the transport by traveling Alfvén waves could strongly modify the level of differential rotation, through, for example, the phase-mixing mechanism (Ionson 1978; Spruit 1999). Moreover, the development of magnetohydrodynamical (MHD) instabilities could lead to a turbulent transport which would efficiently redistribute the angular momentum. This possibility has been studied recently (Cantiello et al. 2014; Fuller et al. 2019; Eggenberger et al. 2019) with the conclusions unclear so far.

Instabilities of a differentially rotating stellar radiative zone with or without the presence of a magnetic field have also been widely investigated theoretically, experimentally and numerically. In hydrodynamical situations, differential rotation can be unstable to various types of instabilities, such as centrifugal or shear instabilities. Centrifugal (or inertial) instabilities require strong enough gradient while weak shear instabilities tend to be stabilized by the Coriolis force (e.g., Knobloch & Spruit 1982). In the MHD case, a shear flow which is hydrodynamically stable can become unstable because of the presence of a large-scale magnetic field. This has been studied in various configurations and in particular when the differential rotation is forced through the boundaries. This is the case of the Taylor Couette flow in cylindrical geometry (or the equivalent spherical Couette flow in spherical geometry). A detailed review of the various MHD instabilities which can arise in Taylor-Couette flows for different rotation rates of the inner and outer cylinders was published recently by Rüdiger et al. (2018). The main instabilities described in that review are the current-driven Tayler instability (Tayler 1973; Markey & Tayler 1973), which is purely magnetic, and the magnetorotational instability (Velikhov 1959; Chandrasekhar 1960; Acheson 1978; Balbus & Hawley 1992) which necessitates a gradient of rotation and is thus shear-driven. As described in Rüdiger et al. (2018), the MRI exists for various large-scale magnetic field geometries: the standard MRI is found for purely axial fields, the so-called azimuthal-MRI for purely azimuthal fields and the so-called helical-MRI for a mixed axial or azimuthal configuration. It could be argued that the Tayler instability is the most relevant for stellar interiors since it only requires a magnetic configuration sufficiently dominated by its toroidal or its poloidal component and a rather weak rotation or differential rotation (Spruit 1999). Detailed studies have been conducted using linear stability analysis for purely toroidal fields with various latitudinal dependences in rotating or differentially rotating radiative zones (Kitchatinov & Rüdiger 2008; Rüdiger & Kitchatinov 2010; Rüdiger et al. 2016). These analyses were local in radius but global in the horizontal directions and took into account the effects of stratification, focusing particularly on a realistic stellar regime where the heat conductivity is high. The m = 1 Tayler instability was found to develop even for a large rotation rate compared to the toroidal Alfvén frequency but with very weak growth rates.

In this work, we do not focus on the instability of a purely toroidal field but, rather, we wish to study the global 3D evolution of an initally poloidal field wound-up into a toroidal field by an initial differential rotation. The system containing all the physical ingredients of a stellar radiative zone (i.e., stratification, axisymmetric meridional flow, shear, global rotation, a mixed poloidal and toroidal magnetic field configuration, heat conductivity, viscosity, and magnetic diffusivity) is then let free to evolve into potentially unstable equilibria. Our recent numerical studies (Jouve et al. 2015; Meduri et al. 2019) show that it is in fact the MRI which is favored in these specific conditions. In these calculations, the initial poloidal field is wound-up by the differential rotation imposed initially for Jouve et al. (2015) and forced trough the boundaries in Meduri et al. (2019) until the Maxwell stresses feed back on the flow. In these situations, the toroidal Alfvén frequency always remains small compared to the rotation frequency and the dynamics associated with the rotation and the shear dominate. The growth rate of the Tayler instability is thus probably strongly reduced by the rotation, as shown by Pitts & Tayler (1985) or Kitchatinov & Rüdiger (2008) but the conditions for the development of the MRI are gathered so that the instability grows on a rotation time-scale. However, the important effects of stable stratification are omitted in the 3D numerical calculations cited above. Only a few recent 3D global numerical studies have focused on the effect of stable stratification on MHD instabilities in specific cases, like for example Philidet et al. (2020) for spherical Couette flows, Guerrero et al. (2019) for the Tayler instability in a non-rotating spherical shell or Szklarski & Arlt (2013) for the Tayler instability of a toroidal field produced by the winding-up of an initial poloidal field. In this last study, very similarly to what is presented in this paper, the wound-up magnetic field is found to be unstable only if the feedback on the differential rotation is inhibited until the ratio of toroidal Alfvén frequency to rotation frequency becomes sufficiently large that the Tayler instability sets in. The MRI has thus probably been stabilized by the stable stratification in these particular calculations. In this paper, we investigate the possibility that high heat conductivities could let the MRI develop again in the same type of numerical setup.

In fact, in most studies dedicated to instabilities of MHD flows with differential rotation, the effect of the stable stratification is often neglected. For the application to stellar interiors, this is a crucial ingredient which may suppress a large number of instabilities, in particular the MRI (Spruit 1999). Indeed, in order to avoid doing work against the stable stratification, the unstable displacements must be nearly horizontal and thus the vertical wavenumber must be high, at which point the diffusive effects will act to make the perturbations decay away. However, in the hydrodynamical case, it has been shown that the largest growth rates of the instability of a horizonthal shear flow would be mostly unaffected by the presence of a large Brunt-Väisälä frequency N Deloncle et al. (2007) for the inflectional instability, Kloosterziel & Carnevale (2008) for the inertial instability). Moreover, non-adiabatic effects should also be considered: if the thermal diffusivity is large, which is the case for stellar radiative zones, the effect of the stable stratification can be strongly reduced and some instabilities may survive for higher values of N (see Townsend (1958), Zahn (1992) for the case of a vertical shear in a stably stratified atmosphere). The possible effects of a high thermal diffusion on MHD instabilities have been discussed theoretically for example by Acheson (1978) or Spruit (1999) but very few global numerical simulations exist where MHD states containing mixed poloidal/toroidal fields and meridional flows and differential rotation subject to the Lorentz force feedback in a stably stratified environment with a varying thermal diffusivity have been analysed in detail. This is what we present in this article. This work is a follow-up on Jouve et al. (2015), where the following intial value problem was considered: an initially imposed large-scale poloidal field is wound-up by an initially imposed differential rotation to produce an axisymmetric toroidal field. After approximately an Alfvén time-scale, the magnetic field back-reacts on the differential rotation and the dynamics is dominated by Alfvén waves which progressively damp the differential rotation. We focused in this last work on the possible development of non-axisymmetric instabilities during this whole process. We now study the effects of the stable stratification with various values of the Brunt-Väisälä frequency, when the thermal diffusivity is also allowed to vary. In particular, we wish to determine the characteristics of the new axisymmetric MHD states and whether the MRI found in Jouve et al. (2015) can survive in a stably stratified environment.

The paper is organized as follows: in Sect. 2 we present the model and the numerical code used to solve the MHD equations. Section 3 then discusses the axisymmetric joint evolution of the magnetic field and the flow. We then investigate the stability of this axisymmetric configuration in Sect. 4 and we present our conclusions in Sect. 5.

2. Numerical model

We wish to explore the interplay between magnetic fields and differential rotation in a 3D spherical shell with stable stratification to mimic the physical processes at play in a differentially rotating stellar radiative zone. To do so, we choose to focus on an initial value problem where a magnetic field and differential rotation will be initially prescribed and then leave it free to evolve with time, according to the MHD equations in the Boussinesq approximation. Indeed, for now, we neglect the effects of a varying density. This will be considered in future works. The details of the equations are given in Sect. 2.1, the initial and boundary conditions are then discussed in Sects. 2.2 and 2.3, and the numerical method is finally briefly described in Sect. 2.4.

2.1. Governing equations

Assuming uniform dynamic viscosity μ, magnetic diffusivity η, thermal conductivity χ, and neglecting the local sources of heat and the centrifugal force, the governing equations under the Boussinesq approximation of a magnetized flow are

$$\begin{array}{cc}& \mathbf{\nabla }·\mathit{v}=0,\hfill \end{array}$$(1)

$$\begin{array}{cc}& \frac{D\mathit{v}}{\mathit{Dt}}=-2{\mathbf{\Omega }}_{\mathbf{0}}\times \mathit{v}-\alpha T_1\mathit{g}-\frac{1}{\rho }\mathbf{\nabla }(p_1+\frac{B^2}{8\pi })\hfill \\ \hfill & +\frac{1}{4\pi \rho }(\mathit{B}·\mathbf{\nabla })\mathit{B}+\nu \mathrm{\Delta }\mathit{v},\hfill \end{array}$$(2)

$$\begin{array}{cc}& \frac{D{T}_1}{\mathit{Dt}}+\mathit{v}·\mathbf{\nabla }\overline{T}=\kappa \mathrm{\Delta }{T}_1,\hfill \end{array}$$(3)

$$\begin{array}{cc}& \frac{\partial \mathit{B}}{\partial t}=\mathbf{\nabla }\times (\mathit{v}\times \mathit{B})+\eta \mathrm{\Delta }\mathit{B},\hfill \end{array}$$(4)

where v is the velocity field, B is the magnetic field, Ω₀ is the rotation rate at the rotation axis, $T(r,\theta ,t)=\overline{T}(r)+T_1(r,\theta ,t)$ is the temperature field with $\overline{T}(r)$ the temperature of the reference state and T₁ its fluctuation, ρ is the uniform density of the reference state, p₁ is the pressure fluctuation, gravity is proportional to 1/r², α is the coefficient of thermal expansion, ν = μ/ρ is the kinematic viscosity and κ = χ/(ρ c_p) is the thermal diffusivity where c_p is the heat capacity at constant pressure.

These equations are then non-dimensionalised using d = r_o − r_i (where r_i and r_o are, respectively, the inner and outer radii of the spherical shell) the thickness of the spherical domain, as the length unit, the poloidal Alfvén time $t_{\mathrm{Ap}}=d\sqrt{4\pi \rho }/B_0$ as the time unit where the surface radial magnetic field at the poles B₀ is the poloidal magnetic field unit, $d{\mathrm{\Omega }}_0\sqrt{4\pi \rho }$ as the toroidal magnetic field unit, V_Ap = d/t_Ap as the meridional circulation unit, dΩ₀ as the azimuthal velocity flow unit, ΔT = T_o − T_i as the temperature unit where T_o and T_i are respectively the temperature at the outer and at the inner radius of the spherical shell and $d^2{\mathrm{\Omega }}_0^2\rho$ as the pressure unit. The full set of governing equations of the problem is given in Appendix A, namely the equations for the three components of the velocity field, for the three components of the magnetic field, and for the temperature field.

Five dimensionless numbers appear in the set of equations:

$$\begin{array}{cc}& \mathrm{Lo}=\frac{t_{\mathrm{\Omega }}}{t_{\mathrm{Ap}}}=\frac{B_0}{d{\mathrm{\Omega }}_0\sqrt{4\pi \rho }},\hfill \end{array}$$(5)

$$\begin{array}{cc}& \frac{N}{{\mathrm{\Omega }}_0}=\frac{1}{{\mathrm{\Omega }}_0}\sqrt{\frac{\alpha g\mathrm{\Delta }T}{d}},\hfill \end{array}$$(6)

$$\begin{array}{cc}& \mathrm{Lu}=\frac{t_{\eta }}{t_{\mathrm{Ap}}}=\frac{d{B}_0}{\eta \sqrt{4\pi \rho }},\hfill \end{array}$$(7)

$$\begin{array}{cc}& \mathrm{Pr}=\frac{t_{\kappa }}{t_{\nu }}=\frac{\nu }{\kappa },\hfill \end{array}$$(8)

$$\begin{array}{cc}& \mathrm{Pm}=\frac{t_{\eta }}{t_{\nu }}=\frac{\nu }{\eta }·\hfill \end{array}$$(9)

The Lorentz number Lo measures the ratio between the rotation time-scale t_Ω and the Alfvén time-scale based on the poloidal field t_Ap, N/Ω₀ is the ratio between the Brunt-Väisälä frequency and the rotation frequency, the Lundquist number Lu measures the ratio between the poloidal Alfvén time-scale and the magnetic diffusion time t_η = d²/η and finally the Prandtl numbers quantify the ratio of diffusivities or the ratio of diffusive time-scales where t_ν = d²/ν is the viscous time-scale and t_κ = d²/κ is the thermal diffusive time-scale.

We can add to these numbers, the definition of the Ekman number, which will be mentioned in the text, and measures the ratio of rotation to viscous time-scales: E_k = ν/Ωd² = Lo/(PmLu). We chose in this study to fix the values of two dimensionless numbers, namely the Lundquist number Lu and the magnetic Prandtl number Pm. We then focus on the effects of the three other parameters: the Lorentz number Lo, the ratio N/Ω₀ and the Prandtl number Pr. We shall see that in the axisymmetric case, the number of relevant dimensionless parameters can in fact be, in some limit cases, reduced to only one.

2.2. Initial and boundary conditions

In this work, we focused on initial conditions which would produce a large-scale magnetic field, likely to be unstable to MHD instabilities under certain circumstances. To do so, we started from a poloidal field which was acted upon by an initial differential rotation. The winding-up of the initial poloidal field by the differential rotation naturally produced a toroidal magnetic field. We propose to focus on the conditions for stability of such a magnetic configuration embedded in a stably stratified atmosphere.

Initially, we thus choose the magnetic field to be axisymmetric, purely poloidal with a constant current density. The detailed expression of the initial magnetic field then reads:

$$\begin{array}{cc}\hfill \mathit{B}(r,\theta ,t=0)& ={\mathit{B}}_{\mathrm{p}}(r,\theta ,t=0)\hfill \\ \hfill & =\frac{3rcos\theta B_0}{r_o}\frac{1-4r_o/3r+r_i^4/3r^4}{1-{(r_i/r_o)}^4}{\mathit{e}}_{\mathit{r}}\hfill \\ \hfill & -\frac{3rsin\theta B_0}{2r_o}\frac{3-8r_o/3r-r_i^4/3r^4}{1-{(r_i/r_o)}^4}{\mathit{e}}_{\mathit{\theta }}.\hfill \end{array}$$(10)

With this choice of normalization, B₀ is the value of the radial field on the axis of rotation at the outer shell r = r_o. For the boundary conditions, we impose that the magnetic field matches continuously to a potential field at both inner and outer boundaries.

The velocity field is also initially axisymmetric but purely azimuthal,

$$\begin{array}{c}\hfill \mathit{v}(r,\theta ,t=0)=v_{\phi }(r,\theta ,t=0){\mathit{e}}_{\mathit{\phi }}=rsin\theta \mathrm{\Omega }(r,\theta ,t=0){\mathit{e}}_{\mathit{\phi }},\end{array}$$(11)

and two different initial rotation profiles will be used. They are discussed in the following section. The boundary conditions for the velocity field are chosen to be impenetrable and stress-free at both inner and outer shells.

The initial temperature field is a purely radial solution satisfying the thermal equilibrium ${\mathbf{\nabla }}^2\overline{T}=0$. Fixed values are imposed for the temperature at both boundaries:

$$\begin{array}{c}\hfill T(r,t=0)=\overline{T}(r)=T_i+\mathrm{\Delta }T\frac{1-r_i/r}{1-r_i/r_o}·\end{array}$$(12)

Finally, to ensure that the flow is stable with respect to convection, the Brunt-Väisälä frequency N must be real which means that ΔT >  0.

2.3. Radial vs. cylindrical differential rotation

The evolution of the toroidal field originating from the winding-up of an initial poloidal field by the differential rotation is expected to be strongly dependent on the differential rotation profile and magnetic configuration. Indeed, the term producing the toroidal field, known as the Ω-effect is proportional to B_p ⋅ ∇Ω and the angle between the poloial field lines and the isocontours of Ω will thus determine the amount of toroidal field created. The efficiency of this Ω-effect is quite important for our study since the ratio between toroidal and poloidal fields is known to be crucial for the stability of the magnetic configurations. That is why we chose to study two different profiles for the initial differential rotation, namely one dependent on the cylindrical radius only and the other one dependent on the spherical radius only. The expressions of both rotation profiles are given below:

$$\begin{array}{cc}& \mathrm{\Omega }(rsin\theta ,t=0)={\mathrm{\Omega }}_0\sqrt{\frac{2}{1+{(rsin\theta /r_o)}^4}},\hfill \end{array}$$(13)

$$\begin{array}{cc}& \mathrm{\Omega }(r,t=0)={\mathrm{\Omega }}_0\frac{1-c_1{(r-r_i)}^2{r}_o/r^3-c_2{(r-r_i)}^2/(r{r}_o)}{1-(c_1+c_2){(1-r_i/r_o)}^2},\hfill \end{array}$$(14)

where Ω₀ is the rotation rate at the equator at r = r_o and where c₁ = 0.980 and c₂ = 0.214 are chosen such that the contrast in the rotation rate between the inner and outer shells is approximately the same for both profiles, namely (Ω_i − Ω₀)/Ω₀ ≈ 1. We thus choose to initially impose a strong differential rotation which is then let free to evolve without any forcing. The transport of angular momentum resulting from the system dynamics will then naturally modify this initial profile.

Figure 1 illustrates these two different profiles of differential rotation and enables us to envision the interaction with the initial poloidal field since it is overplotted in black dashed lines. For example, we can tell that the radial differential rotation profile is likely to produce a strong toroidal field in the bulk of our domain since this is where the poloidal magnetic field lines are almost perpendicular to the isocontours of Ω. On the contrary, in the cylindrical case, the orthogonality is more confined to a region close to the top boundary at mid-latitudes and the toroidal field will thus be mostly produced in this region. Another consequence of our initial setup is that the resulting toroidal field will be antisymmetric with respect to the equator, positive in the Northern hemisphere, negative in the Southern hemisphere, and vanishing at the equator.

Fig. 1. Initial configurations for the cylindrical rotation profile (left) and radial one (right). Ω is scaled with Ω0 and superimposed are the poloidal magnetic field lines.

2.4. Numerical method

The numerical simulations were computed with the numerical code MagIC (Wicht 2002; Gastine & Wicht 2012). MagIC is a fully documented, publicly available code¹ which solves the MHD equations in a spherical shell using a poloidal toroidal decomposition for the mass flux and the magnetic fields:

$$\begin{array}{cc}& \mathit{\rho }\mathit{u}=\mathrm{\nabla }\times \mathrm{\nabla }\times (W{\mathit{e}}_{\mathit{r}})+\mathrm{\nabla }\times (Z{\mathit{e}}_{\mathit{r}}),\hfill \end{array}$$(15)

$$\begin{array}{cc}& \mathit{B}=\mathrm{\nabla }\times \mathrm{\nabla }\times (C{\mathit{e}}_{\mathit{r}})+\mathrm{\nabla }\times (D{\mathit{e}}_{\mathit{r}}),\hfill \end{array}$$(16)

where W (C) and Z (D) are the poloidal and toroidal potentials. The scalar potentials W, Z, C, D and the pressure p are further expanded in spherical harmonic functions up to degree l_max in colatitude θ and longitude φ and in Chebyshev polynomials up to degree N_r in the radial direction. An exhaustive description of the complete numerical technique can be found in Gilman & Glatzmaier (1981). We also make use of the spherical harmonic transforms contained in the SHTns library (Schaeffer 2013) which greatly decreases the computational time for our calculations. Typical numerical resolutions employed in this study range from (N_r = 65, l_max = 170) for the more diffusive cases to (N_r = 129, l_max = 341) for the less diffusive ones. The considered spherical shell extends from longitude φ = 0 to φ = 2π, from colatitude θ = 0 to θ = π and from radius r = r_i = η/(1 − η) to r = r_o = 1/(1 − η) where η = r_i/r_o is the aspect ratio of the shell, equal to 0.3 in all our calculations.

3. Axisymmetric evolution

In this section, we consider the axisymmetric configurations resulting from the evolution of a poloidal field in a differentially rotating stably stratified spherical shell. We performed a parametric study to get an overview of the different configurations that can be reached. In analyzing the simulation results, we shall benefit from previous studies by Gaurat et al. (2015) and Jouve et al. (2015) where the same problem has been considered although with simplifying assumptions. In Gaurat et al. (2015) meridional motions were neglected altogether, while they were taken into account in Jouve et al. (2015) but without the effects of a stable stratification.

As shown on Fig. 1, the initial differential rotation is either cylindrical or radial, but has the same maximum contrast ΔΩ/Ω. We vary the non-dimensional numbers Lo = t_Ω/t_Ap, N/Ω₀ and Pr, while the two others Lu and Pm are fixed to Lu = 50 and Pm = 1. We also restrict ourselves to initial poloidal fields that are weak enough (that is low Lo value) so that the Ω-effect produces magnetic configurations dominated by the toroidal component of the magnetic field. Indeed, the toroidal field will grow linearly with time until the magnetic field back-reacts on the differential rotation. The winding time-scale for the toroidal field to get to the same amplitude as the poloidal field is defined as t_w = 1/qΩ where q = r||∇Ω||/Ω is the shear parameter, of order 1 in our case, such that for our situation, t_w ≈ t_Ω. For the Maxwell stresses to back-react on the flow, a time-scale equal to t_Ap is needed. If t_Ω ≪ t_Ap and thus Lo ≪ 1, the toroidal field will then have time to grow significantly above the poloidal field value before the differential rotation profile is affected by the magnetic field.

In the following, we first discuss the range of parameters that is relevant to stellar radiative zones and then specify the parameters of the present simulations (Sect. 3.1). The typical magnetic configurations obtained from radial vs cylindrical initial differential rotation are described in Sect. 3.2. The effects of the stable stratification on these configurations are analyzed in Sect. 3.3.

3.1. Physical parameters in stellar radiative zones and the parameter range of our numerical simulations

Parameters such as the Brunt-Väisälä frequency and the Prandtl number come from stellar evolution models, whereas rotation rates are obtained from observations. For main sequence massive and intermediate-mass stars, the ratio N/Ω₀ is typically much larger than 1, except for stars rotating near break-up velocity for which N/Ω₀ ∼ 1, meanwhile the Prandtl number Pr is always much smaller than 1. For example, stellar structure models of a 3 M_⊙ star indicate that, during the main sequence, N ∼ 1 − 2 × 10⁻³ s⁻¹ away from the convective core (Talon & Charbonnel 2008) while Pr ≤ 4 × 10⁻⁶ (Garaud et al. 2015). For rotation periods between 1 and 2.7 days, N/Ω₀ is then comprised between ∼10 and 75. For comparison, this ratio is much higher, N/Ω₀ ∼ 300, in the radiative zone of the Sun. As we shall see below, the product Pr(N/Ω₀)² is also a relevant parameter and its value is typically much smaller than one in main sequence massive and intermediate-mass stars. Taking a rotation period of 2.7 days, we find that Pr(N/Ω₀)² ≤ 0.02 for a 3 M_⊙ main sequence star. The situation is different in the solar radiative zone where Pr(N/Ω₀)² is rather of the order of 1 (Garaud & Acevedo Arreguin 2009). Another important parameter is the Ekman number, that compares the rotation and viscous time scales. Its very small value, E_k = ν/(R²Ω)∼10⁻¹⁴, is out of reach in direct numerical simulations. We nevertheless intend to consider small enough Ekman numbers to respect the order of the characteristic times involved, if not the actual time scale ratio.

The typical conditions at large length scales d in radiative zones of intermediate-mass and massive stars thus read E_k ≪ Pr <  Pr(N/Ω₀)² ≪ 1. It corresponds to the following time scale order: t_ν ≫ t_es >  t_κ ≫ t_Ω >  t_B where t_ν = d²/ν, t_es = (d²/κ)(N/Ω₀)², t_κ = d²/κ, t_Ω = 1/Ω₀ and t_B = 1/N. As shown in Table 1, the simulations performed respect that ordering.

The timescale associated with the initial poloidal field is $t_{\mathrm{Ap}}=d\sqrt{4\pi \rho }/B_0$, the poloidal Alfvén time. We have no direct constrain on the field intensity within stars, but we can use spectropolarimetric observations to get surface values. In particular, the lower limit of the dipolar field of Ap/Bp magnetic stars is close to 300 Gauss and, for a rotation period of 5 days, this field corresponds to a Lorentz number Lo = t_Ω/t_Ap close to 1 (Aurière et al. 2007). This number is expected to decrease strongly towards the stellar interior as the variation of the Alfvén speed $v_{\mathrm{Ap}}=B_{\mathrm{p}}/\sqrt{4\pi \rho }$ is dominated by the density increase. For example, at a radius r = R_⋆/3, the Lorentz number would be 2.7 × 10⁻³ assuming a density ratio of 10⁸ and a dipolar-like radial increase B_p ∝ 1/r³. Even lower Lorentz numbers are expected in Vega-like magnetic stars with 1 Gauss surface field and 1-day rotation period. In the radiative interior of intermediate-mass and massive stars, the magnetic field could also result from a convective core dynamo. Numerical simulations of A and B-type star convective cores (Brun et al. 2005; Augustson et al. 2016) indicate that, in the low Rossby number regime characterizing these convective motions, the generated fields have low Lorentz numbers. Indeed, in Brun et al. (2005) simulation of a 7-days rotating A star, a ratio $B_{\text{rms}}^2/(4\pi \rho r_{\text{c}}^2{\mathrm{\Omega }}^2)~2\times {10}^{-4}$ is found at mid-depth of the convective core. Finally, the Lundquist number measures the ratio of the magnetic diffusion time scale to the poloidal Alfvén time and is expected to be very large, much larger than the value attainable in numerical simulations.

Table 1 lists the parameters used in our simulations. The cylindrical and radial cases are respectively labelled C and R. The effect of varying the profile of differential rotation is first studied, keeping the other parameters fixed (cases C1 and R1). Both for the cylindrical and radial cases, we vary Pr (cases 2–4), then Pr and N/Ω₀ while keeping $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2$ fixed (cases 2, 5, 6 and R9 for the radial case). We also decrease the value of $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2$ (cases 7 and 8) and finally consider lower Lo (case C9 and R1, to be compared with C2 and R2). In all cases the Lundquist number is maintained equal to 50 and the magnetic Prandtl number to 1.

3.2. Influence of a radial vs cylindrical initial differential rotation

Figure 2 displays the evolution of the ratio of the total (integrated over the whole spherical shell) dimensioned azimuthal magnetic energy E_mφ to the total dimensioned poloidal magnetic energy E_mp in two axisymmetric simulations in which the cylindrical differential rotation defined by Eq. (13) (top panel) and the radial differential rotation defined by Eq. (14) (bottom panel) are used. The other parameters of these simulations, denoted C1 and R1 in Table 1, are identical. In both simulations, the ratio E_mφ/E_mp initially increases quadratically before reaching a maximal value in a fraction of t_Ap, namely around 0.70t_Ap for the cylindrical case and 0.35t_Ap in the radial case. The maximal value of the quantity E_mφ/E_mp is equal to 1600 in the cylindrical case and to 620 in the radial case, showing that the magnetic configurations are dominated by the toroidal component.

Fig. 2. Evolution of the ratio between the toroidal and poloidal magnetic energies for cases C1 (top) and R1 (bottom). Eφ seems much stronger in the cylindrical case but it’s because Bφ is more extended and located closer to the surface where Bp is much smaller. If we look locally, we also find that Bφ/Bp is stronger in the cylindrical case but mostly because of the location of Bφ (where Bp is small).

This evolution of the magnetic field, namely a near-linear growth of B_φ followed by a maximum reached at t ∼ t_Ap, is the same as observed in simulations where only the coupled equations for the azimuthal magnetic and velocity fields were solved (e.g. Charbonneau & MacGregor 1992 or Gaurat et al. 2015 for exactly identical initial conditions). This means that the avdection by meridional flows and the diffusive decay of poloidal fields only have a weak effect in this case. The physical explanation is rather simple and well-known in this context of purely azimuthal dynamics: the toroidal magnetic field B_φ, initially set to zero, increases through the winding-up of the initial poloidal magnetic field by the initial differential rotation. This growth is linear as long as the differential rotation and the poloidal field are not modified by the back-reaction of the Lorentz force on the flow. When the Maxwell stress associated with the magnetic field becomes sufficiently strong to change the differential rotation, the Ω-effect is modified accordingly and the growth of the toroidal field stops. This results in a maximum in the evolution of B_φ and of the ratio B_φ/B_p, that occurs after a time of the order of the poloidal Alfvén time t_Ap. Locally, the maximum ratio should then be of the order of (B_φ/B_p)_max ∼ ΔΩt_Ap where ΔΩ is the initial differential rotation. The exact value depends on the shear along the poloidal field line during the linear growth phase (see a more detailed model in Gaurat et al. 2015) and this explains why the maxima reached for the cylindrical and radial differential rotation are different. On longer timescales after the maximum (not shown in Fig. 2), damped oscillations of the global magnetic energy are observed and are due to torsional Alfvén waves whose damping through the so-called phase-mixing mechanism (see e.g. Ionson 1978; Spruit 1999) finally leads to uniform rotation as nothing enforces differential rotation in our set-up. We note that we are not interested in the final uniformly rotating state of the axisymmetric evolution as we focus on MHD instabilities that are likely to be triggered before it is reached. The simplified problem considered in Gaurat et al. (2015) is still relevant to describe the global evolution of the present simulations even though we included meridional motions and stable stratification. However the detailed distributions of the differential rotation and the magnetic field, that are crucial for the occurrence and the nature of possible MHD instabilities, will depend on these processes.

Figure 3 shows the distribution of the toroidal field (together with the contours (in black) of the poloidal field) during the linear growth of magnetic field, that is, at 0.25t_ap for the cylindrical case (left panel) and at 0.125t_ap for the radial case (right panel). The value of B_φ is normalized by $d{\mathrm{\Omega }}_0\sqrt{4\pi \rho }$. We clearly see that the distribution of B_φ is quite different in both cases, the field being mostly confined close to the upper boundary in the cylindrical case and to the bottom boundary in the radial case. This is due to the Ω-effect which acts differently because the angle between the isocontours of Ω and the poloidal field lines is maximal at very different locations in the two cases.

Fig. 3. Contours of toroidal magnetic field (colors) for cases C1 (left) and R1 (right). Time is taken in the middle of the linear growth of Bφ, that is at 0.125tap for case R1 and 0.25tap for case C1. Superimposed are the poloidal magnetic field lines. Units are in Loφ = ωAφ/Ω. Here, the strength of Bφ compared to the global rotation is similar, we always have Loφ <  1.

Beside these spatial distributions, the ratio of the azimuthal Alfvén frequency, ${\omega }_{A_{\phi }}=B_{\phi }/r\sqrt{4\pi \rho }$ to the rotation rate Ω, denoted Lo_φ = ω_Aφ/Ω, is another relevant quantity as it allows to distinguish between the two types of instabilities likely to be triggered: the Tayler instability (TI) or the azimuthal-magnetorotational instability (A-MRI) (Jouve et al. 2015). For sufficiently large values of Lo_φ (i.e. when the toroidal field dominates over rotation), the TI should be favoured because the A-MRI is suppressed when the magnetic field becomes too strong. On the contrary, the A-MRI is favoured for small values of Lo_φ (when the rotation is fast compared to the toroidal Alfvén time) since the growth rate of the TI is strongly reduced by a fast rotation (Pitts & Tayler 1985).

With the chosen normalization, the azimuthal field presented in Fig. 3 thus provides the values of the Lorentz number Lo_φ, showing that it is everywhere lower than one. More generally, in our simulations, the maximum value of Lo_φ reached in the whole computational domain is always of the order of 0.2 − 0.3. Indeed, arbitrarily strong values of Lo_φ cannot be reached because only a certain amount of toroidal field can be built before the Lorentz force back-reacts on the differential rotation. As a consequence, we can already expect that the A-MRI will be the favored instability likely to develop in our simulations. This point will be addressed in Sect. 4.

3.3. Influence of the stable stratification

In the previous subsection, the initial growth of the toroidal field and the subsequent regime of damped oscillations have been explained by the winding-up of poloidal field induced by the initial differential rotation followed by the back-reaction of the magnetic field through Alfvén waves. While stable stratification does not seem to play a role in this process, we know from simulations performed in uniform density background (Jouve et al. 2015) that its presence is in fact crucial when the initial differential rotation is radial. Indeed, in the absence of stable stratification, radial differential rotation drives a fast meridional circulation (of time scale 1/ΔΩ) that strongly redistributes the initial poloidal field and angular momentum, before any coherent and therefore efficient winding-up can happen. As a consequence, the build-up of a magnetic configuration dominated by a toroidal component does not happen. With stable stratification, particularly in a N/Ω₀ ≫ 1 regime, such a fast adiabiatic meridional circulation is efficiently suppressed. Instead, on a larger thermal diffusion time scale, the latitudinal temperature perturbations generated by the differential rotation (through the so-called thermal wind equation) drive an Eddington–Sweet type circulation of time scale t_es = (d²/κ)(N/Ω₀)² (see e.g. Spiegel & Zahn 1992)². As long as t_Ap is of the same order or smaller than t_es, the Eddington–Sweet circulation should not prevent the winding-up process to occur. This is indeed what occurs in our numerical simulations as the ratio $t_{\text{Ap}}/t_{\text{es}}=\text{Pm}/(\text{LuPr}{N}^2/{\mathrm{\Omega }}_0^2)$ varies between 8 × 10⁻⁴ and 1 (according to Table 1). In this regime, the detailed distribution of the angular momentum and of the azimuthal field can nevertheless be affected by the meridional circulation, as we shall see below.

The effect of the stable stratification on the axisymmetric evolution has been analyzed by varying the parameters Pr and N/Ω₀. A striking feature is that, after a transient phase, most of our axisymmetric solutions are controlled by the product $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2$, rather than by the two parameters Pr and N/Ω₀ independently. Indeed, as seen in Fig. 4, the evolution of the kinetic and magnetic energies is very similar after about t = 0.2t_Ap for three simulations having the same $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2=0.25$ but different Pr and N/Ω₀, namely Pr = 10⁻², N/Ω₀ = 5 (blue), Pr = 10⁻³, N/Ω₀ = 15.8 (green) and Pr = 10⁻⁴, N/Ω₀ = 50 (red). The two cases Pr = 10⁻³ and Pr = 10⁻⁴ are even indistinguishable on this plot after about t = 10⁻³t_Ap. Figure 4 also displays the energy evolution of two simulations at $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2=0.04$, the Pr = 10⁻³ (cyan) and Pr = 10⁻⁴ (magenta) cases, again showing a very similar behaviour. All the runs of Fig. 4 (R2, R5, R6, R7, R8) have been performed with the initial radial differential rotation but using the initial cylindrical profile (runs C2, C5 and C6) leads to the same conclusion regarding the dependence on $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2$.

Fig. 4. Temporal evolution of the kinetic energy (on a long timescale, left and zoomed in, mid-panel) and magnetic energy ratio (right) for five different cases: three cases with $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2=0.25$ (blue: R2, green: R5, red: R6 and black: R9) and two cases with $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2=0.04$ (cyan: R7 and magenta: R8). For the magnetic energy plot, the cyan and magenta curves are almost superimposed, as well as the red and green curves. The long term evolution is similar for the cases with the same $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2$, as long as Pr is small enough, but different for different values of $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2$. In particular, the amount of toroidal field produced is much less in the case where $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2=0.04$.

On a short timescale however, flows having the same $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2$ can evolve differently. This is illustrated in the mid-panel of Fig. 4 where a zoom is made between t = 0 and t = 0.2t_Ap. Indeed, the initial conditions generate gravity waves that propagate and oscillate until they are damped by thermal diffusion. The oscillations are clearly visible in the Pr = 10⁻² case (run R2). The thermal damping is so efficient in the case Pr = 10⁻⁴ (for which thermal diffusivity is 100 times larger than for the case Pr = 10⁻²) that the oscillations of the kinetic energy are not even visible. The invariance of the solution with $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2$ is thus relevant only after this initial transient phase. In addition, we observe that run R9 deviates also at late time from the other runs with $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2=0.25$.

Then, not only the kinetic and magnetic energies evolution are close for identical value of $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2$ but in fact the whole solutions are very similar. This is illustrated in Fig. 5 where the spatial structures of the flow and the magnetic field are shown at t = 0.1t_Ap at $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2=0.04$ (two left panels) and $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2=0.25$ (two right panels), in both cases for 2 different values of Pr, namely 1.6 × 10⁻³ and 1.6 × 10⁻⁵ for the left panels (cases R7 and R8) and 10⁻² and 10⁻⁴ for the right panels (cases R2 and R6). The top panels show the rotation rate in color and the meridional flow contours while the bottom panels present the toroidal magnetic field in color and the poloidal field lines in dashed lines.

Fig. 5. Structure of the flow (top panels: rotation rate in color, meridional flow contours in black lines) and of the magnetic field (bottom panels: toroidal field in color and poloidal field lines in dashed lines) at time t = 0.1tap, for R7 and R8 (two left panels) and for R2 and R6 (two right panels).

In Appendix B, we perform a scale analysis of the Boussinesq MHD equations in the parameter regime of the simulations. It shows that for a time scale ordering, t_ν ≫ t_Ap ≫ t_κ ≫ t_Ω ≫ t_B, the evolution of the system only depends on $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2$, Lu and Pm if time is scaled by t_Ap and B_φ by $d\mathrm{\Delta }\mathrm{\Omega }\sqrt{4\pi \rho }$. The Lorentz number Lo only appears as a scaling factor of the ratio $B_{\phi }/B_{\mathrm{p}}={\mathrm{Lo}}^{-1}f(t/t_{\mathrm{Ap}},\mathit{r}/d,\mathrm{Pr}{N}^2/{\mathrm{\Omega }}_0^2,\mathrm{Lu},\mathrm{Pm})$. The simulations that verify the required time ordering do show this $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2$ dependence. The deviations at small time due to the initially excited gravity waves are expected because gravity waves are filtered out by the scale analysis. The strong deviation observed at late time for run R9 is also expected since $\frac{t_{\kappa }}{t_{\mathrm{Ap}}}=\frac{\mathrm{Lu}\mathrm{Pr}}{\mathrm{Pm}}=3.125$ is larger than 1 in this case (in Appendix B, the regime t_κ ≫ t_Ap ≫ t_Ω ≫ t_B is shown to be dominated by waves with negligible effect of the meridional circulation). Besides the $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2$ dependence, the expression of B_φ/B_p is fully compatible with the maximum toroidal to poloidal magnetic energy ratio found in simulations performed for the same $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2=0.25$ but two different Lorentz numbers (runs R1 and R2). The ratio indeed increased by a factor ≈2.5², as Lo was reduced from 2.5 × 10⁻³ (right panel of Fig. 2) to Lo = 10⁻³ (right panel of Fig. 4). In addition to a simplification in the physical interpretation, the scale analysis allows us to conduct the parametric study of the flow by varying only one non-dimensional number $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2$ instead of the three Lo, N/Ω₀ and Pr.

We thus consider how the rotation and magnetic configurations depend on the values of $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2$. As shown on the right panel of Fig. 4, the maximum ratio of toroidal to poloidal magnetic energy is 1500 for $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2=0.04$ while it is close to 4000 for $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2=0.25$ (and is reached at an earlier time). The stable stratification is thus favorable to the creation of a magnetic field configuration more strongly dominated by its toroidal component. According to Fig. 5, in the most stably stratified case (right panels), the differential rotation remains close to its initial profile and then mostly dependent on radius, contrary to the less stably stratified case (left panels), for which the differential rotation is reduced and tends to become cylindrical. This tendency is expected because when stable stratification is less efficient the system can evolve more freely towards a flow satisfying the Taylor–Proudman constraint, valid for unstratified systems: $\frac{\partial \mathrm{\Omega }}{\partial z}=0$, with the z-direction parallel to the rotation axis. The reduced level of differential rotation can also be explained by an efficient meridional transport of angular momentum in the less stratified case. The ratio Ω_i/Ω_o indeed decreases from 2 at t = 0 to 1.6 at t = 0.1t_Ap in the less stratified case whereas it remains close its initial value in the more stratified case.

This difference in the level of differential rotation then naturally explains why a weaker toroidal magnetic field is produced in the less stratified case. This is visible on the bottom left panels of Fig. 5 where the maximum value of Lo_φ only reaches 0.15 compared to the other cases where it is already close to 0.2. In addition, we observe that the poloidal field configuration has been significantly altered in the less stratified cases compared to the initial condition. The poloidal field tends to align on the cylindrical isocontours of Ω at mid-latitudes, again preventing a strong Ω-effect to be at play. This significant change of the poloidal field is due to its advection by the meridional circulation which is more efficient in the less stratified case.

From the axisymmetric numerical simulations performed for different stable stratifications, we conclude that in the regime considered, the effect of the stable stratification is controlled by the product $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2$ and that stable stratification favors the creation of magnetic field configurations more strongly dominated by their toroidal component.

4. Stability of the magnetic configurations

We now turn to investigate the stability of the axisymmetric magnetic configurations determined in the first part of this work. We perturb the magnetic field by adding a random noise on the axisymmetric poloidal field and then follow the temporal evolution of the various non-axisymmetric modes m ≠ 0, in the same way as was done in Jouve et al. (2015). We first consider the stability of the system with radial and cylindrical differential rotations and a fixed $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2=0.25$ (cases R2 and C2) and argue that the observed instability is of the MRI type, this is presented in Sect. 4.1. In Sect. 4.2, the effect of varying the thermal diffusivity on the instability is then studied (cases C3, C4 and R3, R4). To help us understand the characteristics of the instabilities, we compare our results with a local stability analysis in Sect. 4.3. While the Lorentz number Lo has been fixed to a small value (namely Lo = 5 × 10⁻³ for the cylindrical case and Lo = 10⁻³ for the radial case) to maximize the possibility for a non-axisymmetric instability to fully develop (see Jouve et al. 2015), we investigate in Sect. 4.4 the effects of increasing Lo (cases C9 and R1).

4.1. Radial vs. cylindrical differential rotation

We first investigate the typical evolution of an unstable situation and compare the behaviour of the simulations initialized with the cylindrical and radial differential rotation profiles. Figure 6 shows the evolution of the poloidal magnetic energy contained in the first 11 azimuthal wavenumbers, including the axisymmetric m = 0 mode, which is approximately steady during the time considered. In both cases, a non-axisymmetric instability grows exponentially in a fraction of Alfvén time to quickly reach the level of the axisymmetric energy at about 0.6t_ap. The right panels of Fig. 6 enable us to visualize the location and structure of the unstable modes, by showing the amplitude of the fluctuations of the radial component of the magnetic field. In both cases, the instability develops preferentially where the azimuthal magnetic field is maximum (see Fig. 3 for the axisymmetric configuration which was perturbed) and where and when a significant amount of differential rotation exists. In both cases, the growth rate of the most unstable mode is a fraction of the rotation rate. It is approximately equal to 5 × 10⁻²Ω for the m = 4 mode in the cylindrical case and 2.5 × 10⁻²Ω for the m = 1 mode in the radial case.

Fig. 6. Instability in case C2 (top) and R2 (bottom). Shown are the temporal evolution of the poloidal magnetic energy in the first 11 azimuthal wavenumbers (averaged in r and θ) (left), the fluctuations of the radial component of the magnetic field at a particular longitude (mid) and the rotation profile (right) when the instability starts to develop, that is at t = 0.15tap for C2 and t = 0.35tap for R2. The color bar applies to the rotation rate.

We now emphasize the differences between the two cases. First, the time at which the instability starts to grow is quite different. Indeed, in the cylindrical case, the axisymmetric equilibrium which is perturbed is already unstable as soon as the perturbation is introduced, leading to the exponential growth of the non-axisymmetric modes from approximately t = 0.1t_ap. At this stage, the axisymmetric evolution is still in its linear growth of azimuthal magnetic field, as shown in Fig. 2. On the contrary, in the radial case, the instability develops only later, at about t = 0.3t_ap, approximately when the maximum of axisymmetric B_φ is reached and thus when a strong back-reaction of the Lorentz force on the differential rotation profile has acted. To illustrate this, the right panels of Fig. 6 show the profiles of differential rotation at the time when the unstable non-axisymmetric modes start to grow. It is clear that the cylindrical differential rotation is still mostly identical to its initial condition whereas the radial case has been significantly modified by the back-reaction of the Lorentz force. In particular, a latitudinal differential rotation appears here, which was not present initially since the rotation rate was dependent on radius only. It is exactly at the location where the latitudinal shear is the strongest that the unstable modes are confined.

Another major difference between cases R2 and C2 lies in the structure of the unstable modes. From the middle panels of Fig. 6, it is clear that the displacement of the perturbations is not in the same direction in both cases. Let us express the perturbation as proportional to:

$$\begin{array}{c}\hfill exp[\mathrm{i}(k_{r}r+k_{\theta }\theta +m\phi -\sigma t)],\end{array}$$

with k_r, k_θ and m the radial, latitudinal and azimutal wavenumbers, and σ the complex growth rate. Then in the cylindrical case, the latitudinal wavenumber k_θ is large compared to k_r and the displacement is thus mainly in the radial direction, parallel to gravity. We argue that the radial extent of the unstable mode is in fact mostly due to the structure of the axisymmetric background and not due to the effect of stable stratification. Indeed, in our previous study where the effects of stable stratification were not included (Jouve et al. 2015), the structure of the unstable modes and the growth rates in the equivalent of case C2 where very close to the ones found here.

On the other hand, in case R2, k_r is now dominant compared to k_θ, so that the displacement is mainly in the latitudinal direction, perpendicular to gravity. We thus anticipate that in this simulation, the stable stratification, which is much less effective if the displacement is horizontal, only affects the geometry of the unstable mode and not its growth rate. This is investigated in the next section, where the effect of the stable stratification is increased for the two initial differential rotations.

In both cases, we argue that the instability found here is of MRI type. First, we checked that the flow is hydrodynamically stable by perturbing the flow when the magnetic field is set to 0 at the time where the instability develops in the MHD case. The instability could nevertheless be a current-driven instability of the Tayler type since the magnetic configuration contains current and is strongly dominated by the toroidal component, as we can clearly see on Fig. 7. This figure shows the magnetic field lines of the background axisymmetric magnetic field traced around the location of the instability. This 3D view enables to clearly see the dominance of the toroidal component of the field and also allows us to see that the maximum amplitude of the unstable modes (shown on Fig. 7 by the isosurfaces of the fluctuating axial component of the field) is mainly located where the toroidal field is maximum. The location of the maximum B_φ naturally corresponds to the region where the shear is also maximum since the shear is responsible for the generation of toroidal field through the Ω-effect. We thus find here that the instability develops mainly where the shear is concentrated. This would be different if the instability was of the Tayler type, because then the location of the unstable modes would be correlated with the gradients of toroidal field, where the currents are maximal. Moreover, the most unstable mode in the cylindrical case is not the m = 1 as expected for the Tayler instability. It is however the m = 1 mode which is the most unstable in the radial case as seen on the figure but this is not incompatible with an MRI instability in the fast thermal diffusion case. In Acheson (1978), a detailed theoretical description is made of all the various MHD instabilities which can arise in stellar radiative zones. In this seminal paper, the MRI is not explicitly quoted but an instability associated with a shear and which necessitates the presence of a magnetic field is studied, when the thermal diffusivity is high and in the limit where ${({\omega }_{A_{\phi }}/\mathrm{\Omega })}^2{E}_k^{-1}\text{Pm}$ is also high. In this situation, he argues that the most favored unstable mode is precisely the m = 1 mode. The values in our simulations of the parameter ${({\omega }_{A_{\phi }}/\mathrm{\Omega })}^2{E}_k^{-1}\text{Pm}$ when the instability develops in the radial case is of the order of 2 × 10³ so that the limit studied by Acheson (1978) does apply here. In both cases started with a cylindrical or a radial differential rotation, we thus observe the presence of a MRI which is driven mostly by the initial radial differential rotation in the cylindrical case and driven by the latitudinal shear that is produced by the back-reaction of the Lorentz force in the radial case. The instability is allowed to exist in both cases with a relatively high thermal diffusivity (Prandtl number of 10⁻²). We now wish to investigate the effect of varying the thermal diffusion on the instability.

Fig. 7. 3D views of the instability in the two cases: Case C2 on top and case R2 below. The magnetic field lines of the background axisymmetric magnetic field are plotted around the location of the instability and colored with the values of the toroidal magnetic field (in the Y-direction in the Cartesian frame shown at the bottom left) and isosurfaces of the axial component of the fluctuating magnetic field are overplotted.

4.2. Effect of the thermal diffusivity

The stable stratification has the tendency to strongly reduce the development of non-axisymmetric instabilities, as shown for example in Spruit (1999). In particular, as the stable stratification limits radial displacements, it will strongly affect instabilities that require them to develop. By damping temperature deviations, thermal diffusion diminishes the amplitude of the restoring buoyancy force and thus the effect of the stable stratification. In order to study these effects, we therefore decrease the thermal diffusivity, and thus increase the Prandtl number Pr, keeping the same value for N/Ω₀. These correspond to cases C3, C4 and R3, R4. Figure 8 shows the temporal evolution of the poloidal magnetic energy decomposed into the first 11 azimuthal wavenumbers in cases where Pr = 0.1 and Pr = 1. The left panels correspond to cases C3 and C4 and the right panels cases R3 and R4. Compared to the magnetic energy evolution of Fig. 6 where we had Pr = 10⁻², it is clear that for the cylindrical case, the instability is largely suppressed by the increase of the stable stratification effect. In particular, for Pr = 1, the axisymmetric solution becomes completely stable to any non-axisymmetric perturbation. In other words, the preferentially radial displacements that were unstable at Pr = 10⁻² are inhibited at Pr = 1. The transition between Pr = 10⁻² and Pr = 1 can be linked to the value of the critical lengthscale above which the effects of the stable stratification are not diminished by thermal diffusion. This critical lengthscale l_c is determined by equating the buoyancy and the thermal diffusion time scales:

$$\begin{array}{c}\hfill \frac{l_{\mathrm{c}}^2}{\kappa }=\frac{1}{N}\mathrm{a}ndthus{l}_{\mathrm{c}}=\sqrt{\frac{\kappa }{N}}·\end{array}$$

Fig. 8. Temporal evolution of the poloidal magnetic energy in the first 11 azimuthal wavenumbers for cases C3 (top left), C4 (bottom left), R3 (top right), and R4 (bottom right).

With the dimensionless parameters used in our calculations, the critical lengthscale reads:

$$\begin{array}{c}\hfill \frac{l_{\mathrm{c}}}{d}=\sqrt{\frac{E_k}{\mathrm{Pr}}\frac{\mathrm{\Omega }}{N}}·\end{array}$$

The computation of this quantity gives a value ranging from 4% to 0.4% when Pr goes from 10⁻² to Pr = 1, for these cases where $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2=0.25$. The unstable radial lengthscale seen in Fig. 6 being of the order of a few percent of the computational domain, we argue that this case is only marginally affected by the stratification. This is consistent with the fact that, as mentioned above, a very similar unstable mode was found in the corresponding unstratified simulation by Jouve et al. (2015). However, the reduction of the critical lengthscale causes the instability to disappear in the Pr = 1 case. Such a behaviour where increasing the stable stratification removes the instability is reminiscent of the vertical shear instability in a vertically stratified medium (Dudis 1974; Lignières et al. 1999).

The situation is quite different in the radial case (right panels of Fig. 8). Now, the instability survives even with the increase of the effects of stable stratification, and grows on time scales similar as in the Pr = 10⁻² case. As we show below, this comes with the fact that the unstable displacements become more and more horizontal. In Fig. 9, we show the structure of the unstable mode for cases R3 and R4 at two different longitudes and we plot the profile of the magnetic field and the rotation rate, averaged in longitude. In both cases, the background flow and field are quite similar even if the value of $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2$ differs. This is also true for the cylindrical case (not shown here), which also confirms that the absence of an instability in cases C3 and C4 is mostly due to the effect of stable stratification on the characteristics of the instability (namely the lengthscale) and not on the background flow and field. We recover the fact that the displacement is indeed mostly horizontal, with a latitudinal lengthscale extremely dominant in comparison to the radial scale in case R4 where $\text{Pr}{N}^2/{\Omega }_0^2=25$. The location of the instability is still mostly where the latitudinal gradient of Ω lie, as seen on the right panels. It is thus clear here that the effect of the strong stratification is to force the unstable modes to become more horizontal and since their origin is the latitudinal gradient of Ω, the instability survives even when the degree of stratification is increased. We note that the most unstable azimuthal wavenumber is still m = 1 so that the strong stratification does not seem to significantly affect the azimuthal scale.

Fig. 9. Fluctuations of the radial component of the magnetic field at a particular longitude, axisymmetric magnetic field (toroidal field in color and poloidal field lines) and rotation rate at approximately time t = 0.5tap (during the linear phase of the instability) for cases R3 (Pr = 0.1, top) and R4 (Pr = 1, bottom). We clearly see that the perturbation direction is orthogonal to gravity, especially for the most stratified case R4. We note that the background flow and field are very similar in both cases.

It is quite striking here that the growth rates of the unstable modes do not seem to be strongly affected by the stratification. Indeed, the growth rates of the cases R2, R3 and R4 are similar. Meanwhile, the ratio of the radial to the latitudinal wavenumbers increases with the increased stratification. This can be understood by the fact that the instability here is driven by the latitudinal (or horizontal) gradient of Ω. Indeed, this behaviour is reminiscent of previously studied hydrodynamical instabilities driven by an horizontal shear in a vertically stratified medium. In the case of the centrifugal (or inertial) instability studied by Kloosterziel & Carnevale (2008), the dispersion relation of the unstable modes clearly shows that the growth rate of a mode with given latitudinal and azimuthal wavenumbers can be made invariant to a stratification increase by adapting (that is increasing) the vertical wavenumber accordingly. The possibility to adapt the vertical lengthscale to get the same growth rate also exists when the shear instability is of the inflectional type (Deloncle et al. 2007). The effect of the stable stratification on a vertical shear instability is very different. In an inviscid and adiabatic case there is simply no instability when the Richardson number exceeds 1/4, while a high thermal diffusivity can potentially destabilize predominantly horizontal perturbations but then the growth rates are vanishingly small (Lignières et al. 1999). A simple physical interpretation is that the most unstable modes of a vertical shear necessarily involve vertical motions, such as for example in Kelvin–Helmholtz billows. Thus, by opposing vertical motions, stable stratification either kills the instability or reduces it strongly. On the contrary, for an horizontal shear, the stable stratification may affect the preferred vertical wavelength of the perturbation but this does not prevent the unaffected horizontal motions to efficiently draw energy from the horizontal shear.

While these purely hydrodynamical cases help interpret the effect of the stratification, both the centrifugal and the inflectional instabilities are absent from our simulations since the differential rotation does not fulfill the inviscid and unstratified criteria for these instabilities. We thus expect that the observed instability is a magnetorotational instability due to the latitudinal shear and supported by the magnetic field. In the next section, we check the consistency of our interpretation using a local linear stability analysis in the MHD case.

4.3. Comparison to the Acheson dispersion relation

We now wish to analyze our numerical results at the light of a local linear instability analysis, strongly inspired by the work of Acheson (1978) where various types of MHD instabilities in different regimes were investigated, as already quoted at the end of Sect. 4.1. We are particularly interested in the impact of stable stratification on our instabilities and on the differences found between the cylindrical and radial cases. We recall here the various steps of the establishment of the Acheson dispersion relation of interest in our case (Eq. (3.20) in Acheson 1978) without indicating all the details, which can be found in Appendix C.

First, the MHD equations governing the system with thermal, viscous and magnetic diffusion are linearised around the background axisymmetric state in cylindrical geometry (which is assumed to be purely toroidal both for the magnetic and the velocity fields) and by considering small amplitude harmonic perturbations in space and time of the form

$$\begin{array}{c}\hfill exp[\mathrm{i}(k_{s}s+k_{z}z+m\phi -\sigma t)].\end{array}$$(17)

Here k_s = 2π/λ_s (k_z = 2π/λ_z) is the radial (axial) wavenumber of the instability and m its azimuthal order which is an O(1) integer. When the imaginary part of σ is positive, the applied perturbation is unstable and grows exponentially at a rate γ = ℑ(σ). Then, we assume here that the thermal diffusivity κ is much higher than the magnetic diffusivity η, which is the case in our setup where Pm = 1 and Pr ≪ 1. In this situation, the dispersion relation of Acheson is reduced to a simpler expression: a polynomial equation of degree 4 in the dimensionless frequency $\stackrel{\sim }{\omega }=\omega /{\mathrm{\Omega }}_0$.

We solve numerically that polynomial Eq. (C.4) by choosing as background axisymmetric profiles our numerical solutions Ω(r, θ) and B_φ(r, θ) at the time where the instability develops in the simulations. The various parameters defined in Appendix B and which play a role in the calculation of the instability growth rate are the ratio of poloidal wavenumbers β, the azimuthal wavenumber m, the shear parameter q, a parameter b quantifying the gradient of B_φ, the azimuthal Lorentz number Lo_φ, the stratification parameter $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2$ and the Reynolds numbers Re, Rm and Rt. For all these parameters, we take the values estimated or calculated from the simulations. With this procedure, we obtain a 2D map of the theoretical growth rate σ(r, θ) at the time where the instability starts to grow in the simulation. The aim is then to compare the location and the value of the maximum theoretical growth rate in the 2D domain with the location of the unstable mode and the growth rate estimated from the simulation.

An example of such a map is given in Fig. 10, where the azimuthal wavenumber was chosen to be m = 1 and the ratio of poloidal wavenumbers such that k_θ ≪ k_r. This case corresponds to case R2 where the instability was clearly present in the numerical simulation and mostly on the m = 1. On this map of σ/Ω₀, we superimpose the isocontours of the fluctuating component of the radial magnetic field coming from the 3D simulation. We find that the location of the unstable mode coincides well with the theoretical location of the maximum growth rate. The value of the maximum growth rate reaches σ/Ω₀ ≈ 0.1, compared to 2.5 × 10⁻² in the simulation. We do not expect to recover exactly the same growth rates because of the various assumptions underlying the derivation of the dispersion relation, which might not be entirely fulfilled in our simulations. In particular, the analysis of Acheson (1978) is local and we are comparing it here with global numerical simulations, with possible effects of the boundary conditions, especially for the cylindrical case where the instability develops very close to the top boundary of our computational domain. Then, as already pointed out in a previous work (Meduri et al. 2019), the use of the short wavelength approximation (meridional perturbation wavelength much smaller than the typical scale of variation of the background) could also be questioned here, in particular for the radial direction. Anyhow, we do not try here to understand in detail the discrepancy in the values of the growth rate obtained in the local analysis and in the numerical simulations, we just aim at gaining some insight from the local analysis on the possible causes for instability observed in the simulations.

Fig. 10. Cases R2 (left) and R3 (right): contours of the growth rate σ/Ω0 of the m = 1 mode obtained through the Acheson dispersion relation (colors) and superimposed in black lines are the contours of the fluctuating radial magnetic field coming from the 3D simulation. The agreement for the location of the instability is quite satisfactory in both cases.

In the radial cases, the local dispersion relation helps us to determine that it is mostly the latitudinal gradient of Ω which is responsible for the instability and that the presence of the background magnetic field is needed. Indeed, the local analysis predicts that the background flow is hydrodynamically stable (the growth rate is negative when the magnetic field is set to 0). Moreover, when the latitudinal gradient of Ω is set to 0, the growth rate drops dramatically while it stays around the same maximum value of σ/Ω₀ ≈ 0.1 when the radial gradient is set to 0. We thus confirm here the argument developed in the previous section: the instability is here driven by the gradient of rotation in the θ-direction, that is, orthogonal to the stable stratification. As discussed above, this is also probably the reason for the persistence of the instability when the stability of the stratification is increased. The right panel of Fig. 10 shows the location of the instability when the effect of the stable stratification is increased, namely, with Pr = 10⁻¹ instead of Pr = 10⁻². The local analysis still predicts a significant growth rate, again located around the maximum latitudinal gradient of rotation. The value of the growth rate itself is reduced by about 20% but the instability still exists and indeed also observed in the 3D simulation at approximately the same location for the m = 1 mode. In fact, in this case where we chose the poloidal wavenumber β such that the displacement is nearly horizontal (k_θ ≪ k_r), it is expected that the local dispersion relation predicts a small effect of the stable stratification on the growth rate. Indeed, if we look at the coefficients of Eq. (C.4), we see that all the terms involving $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2$ are multiplied by the quantity sin θ − β cos θ. In the limit case where k_θ is 0 and k_r is large (but finite), β reduces to tanθ and the term multiplying $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2$ vanishes. Of course, we are not really at this limit here but there can be a factor of at least ten between the poloidal wavenumbers such that the effect of the stable stratification becomes very weak on the value of the growth rates.

In Fig. 11, we illustrate the fact that the linear analysis also predicts that the geometry of the unstable mode in the radial case should change as the stable stratification increases. The figure shows the maximum growth rate reached in the (r, θ) plane for the background flow and field of case R2 as a function of the poloidal wavenumber ratio when the value of N/Ω₀ is varied from 5 to 100. We clearly see that when the level of stratification is increased, the most unstable mode adapts its radial to horizontal wavenumber ratio: the most unstable mode becomes more and more horizontal when the stratification is increased, as also seen in the 3D simulation and as observed in the hydrodynamical studies discussed in the previous section. We also note that the maximum growth rate always tends to the same value as the level of stratification is increased, so that, theoretically, all modes with a sufficiently large radial to horizontal wavenumber ratio should be equally unstable. The unstable mode seen in the 3D simulation of course possesses a finite k_r/k_θ, probably chosen to minimize the stable stratification effects while fitting in the extension of the background field.

Fig. 11. Spatial maximum of the instability growth rate for the m = 1 mode calculated from Eq. (C.4), as a function of the ratio of poloidal wavenumbers, for different values of N/Ω0, keeping Pr = 10−2 and for the background azimuthal velocity and magnetic fields of case R2.

The situation is different in the cylindrical cases C2 and C3 where we chose, on the contrary, to calculate the growth rate as a function of r and θ but using a ratio β such that k_θ ≫ k_r, as seen in the 3D simulation. And we now choose to focus on the mode m = 4 which is one of the most unstable ones. The results would be similar for the equally unstable m = 2 and m = 3 modes. Figure 12 shows the map of the theoretical growth rate obtained from the dispersion relation for the m = 4 mode, together with the contours of the radial component of the magnetic fluctuations coming from the 3D simulation. Again, the location of the maximum growth rate coincides quite well with the position where the instability is observed in the simulation but the expected growth rate is larger (2.7 × 10⁻¹ compared to 5 × 10⁻² in the simulation). In this case, our procedure enables us to attest that it is now the radial gradient of Ω which is responsible for the instability found here, the growth rates (value and location) being very similar when the $\frac{\partial \mathrm{\Omega }}{\partial \theta }$ is set to 0. When the Prandtl number is now increased to Pr = 0.1, the instability completely vanishes, showing the very strong effect of the stable stratification on the instability in this cylindrical case, as observed in the simulation.

Fig. 12. Same as Fig. 10 but for Case C2 and for the m = 4 mode. Again, the location of the instability in the simulation corresponds quite well with the position of the expected maximum growth rate from the local analysis. Case C3 is not shown since the local analysis does not predict any instability in this case, in agreement with the 3D simulation.

4.4. Effect of Lo

In this last section, we investigate the effect of varying the Lorentz number, which measures the ratio between the dynamical timescales of interest in this study: the rotation time scale to the poloidal Alfvén time scale. From our previous study (Jouve et al. 2015), we know that this parameter is crucial to the full development of the instability. Indeed, since we identify our instability here of the magneto-rotational type with a typical growth rate of the order of the rotation frequency, the Lorentz number quantifies the time it takes for the instability to grow compared to the typical lifetime of the background toroidal magnetic field on which it grows. The optimal case for the full development of the instability is consequently when the Lorentz number is small. To test this argument with the simulations performed in this work, we increased the Lorentz number both in the radial and the cylindrical cases. The growth in time of the magnetic energy contained in the first 11 azimuthal modes for an increased Lorentz number is shown in Fig. 13, both for the cylindrical (left panel) and the radial case (right panel).

Fig. 13. Temporal evolution of the poloidal magnetic energy in the first 11 azimuthal wavenumbers for case C9 (left) and R1 (right): compared to cases C2 and R2, the value of Lo was increased to Lo = 10−2 for the cylindrical case and Lo = 2.5 × 10−3 for the radial case.

As expected, the main effect of increasing Lo in both cases is to suppress the instability in the cylindrical case and drastically decrease its impact on the axisymmetric field in the radial case. To be more precise, the Lorentz number was increased here by decreasing the rotation rate and thus increasing the rotation time. The evolution of the axisymmetric magnetic field is then similar to what is shown in Fig. 2 but with a smaller value for the ratio between toroidal and poloidal magnetic energies. In particular, the maximum toroidal field will still peak at approximately t = 0.7t_Ap for the cylindrical case and at t = 0.3t_Ap for the radial case, but the growth rates are divided by approximately 2 since the Lorentz number was doubled in case C9 and multiplied by 2.5 in case R1 compared to C2 and R2 respectively. As a consequence, the instability does not have time to sufficiently develop to reach the level of the axisymmetric field. The magnetic field after a few Alfvén times will remain mostly unaffected by the presence of non-axisymmetric components.

The conclusion here is similar to the unstratified case (Jouve et al. 2015) and is still valid for the radial case. This is not surprising since the instability is also of MRI type and thus the ratio between the instability growth time and the background magnetic field lifetime will still control the ability of the non-axisymmetric unstable modes to reach the energy of the axisymmetric field. We then anticipate that for the mean axisymmetric field to be significantly modified by the development of the instability, the system must be at low Lorentz number, that is, a relatively weak poloidal magnetic field embedded in a fastly rotating environment. To be more quantitative, in the cylindrical case, Lo must be weaker than 5 × 10⁻³ while in the radial case, it must be even weaker, of the order of 10⁻³. Both values are compatible with the values expected to be found in stellar interiors (see Sect. 3.1), especially for Vega-like stars which possess a weak surface field and a rapid rotation. The difference between the radial and cylindrical cases can be understood by the fact that for the instability to be triggered in the radial case, we first need to wait for the magnetic field to back-react on the flow to produce the latitudinal shear. The instability starts to develop when the toroidal field has already reached its maximum value and begins its decay. The instability thus needs to grow quite fast so that the toroidal field keeps approximately its maximum value during the whole development of the instability. The cylindrical case is different since the instability is able to grow right away on the existing radial shear and consequently while the toroidal field is building up. In the cylindrical case, the instability is thus allowed more time to grow and the range of Lorentz numbers allowing the instability to fully develop is thus extended.

5. Conclusion

In this work, we studied the effects of the stable stratification in the non-adiabatic case on instabilities which can develop when an initial poloidal field is wound up by an initial differential rotation. Two different profiles for the differential rotation were considered, both likely to exist in stellar radiative zones: one, cylindrical, which satisfies the Taylor–Proudman constrain and the other, shellular, which corresponds to what could be expected in a strongly stably stratified layer.

The axisymmetric solutions of this initial value problem were first investigated. We showed that, for fixed Lu and Pm, the axisymmetric evolution depends only on one dimensionless parameter $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2$ measuring the level of stratification, instead of the three independent parameters Lo, Pr and N/Ω₀. This result is found to be consistent with a scale analysis of the Boussinesq MHD equations performed for a time ordering t_ν ≫ t_Ap ≫ t_κ ≫ t_Ω ≫ t_B. In this simplified form, the gravity waves are filtered out, and the system evolves through Alfvén waves and an Eddington Sweet circulation prescribed by a magneto-thermal wind equilibrium and a thermal equilibrium. The parameter Lo = t_Ω/t_Ap only controls the ratio between toroidal and poloidal field. An interesting feature of the strongly stably stratified cases (with relatively high values of $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2$) is that the toroidal to poloidal field ratio becomes higher since the transport of angular momentum through meridional flows is inhibited. Indeed, in this situation, the initial differential rotation is not modified before the Lorentz force starts to back-react on the flow and the Ω-effect is more efficient at producing a toroidal field component.

In stars, the ordering in time-scales given above may not apply since the Alfvén time-scale may be small compared to the thermal diffusion time-scale. Then, the Eddington–Sweet circulation becomes negligible and the system evolution is dominated by Alfvén waves as in Gaurat et al. (2015) where only the coupled evolution of v_φ and B_φ was analysed. In our calculations, we considered a situation where t_Ap ∼ t_κ and found that the gravity waves existing in the transient phase do persist during the whole winding-up process and significantly perturb the flow and field. However, this initial gravity wave transient is a direct consequence of our initial condition which is far from an equilibrium and such a transient is not likely to be present in stars.

When axisymmetric solutions strongly dominated by their toroidal component exist, there are expected to be unstable. This is indeed what was found in Jouve et al. (2015) in the non-stratified case. We tested here the effects of the stable stratification on the instability. It turns out that the situations involving two different initial differential rotation profiles respond quite differently to perturbations. When non-adiabatic effects are important, that is when a large thermal diffusivity is considered, both cases are unstable to a magneto-rotational instability. However, when the thermal diffusivity is reduced and thus when the effects of the stable stratification are increased, the instability disappears in the cylindrical case while the unstable displacements become more and more horizontal in the radial case, with similar growth rates. We argue that this is due to the fact that the radial shear is responsible for the instability in the cylindrical case while it is driven by the latitudinal shear in the other. This latitudinal shear does not exist initially, it is produced by the back-reaction of the magnetic field on the flow. The situation may appear quite specific since it is here the magnetic field itself which creates the conditions for its own instability. However, such phenomena could occur in stellar radiative zones where angular momentum is permanently redistributed by meridonal flows or Aflvén waves. In our case, the level of latitudinal shear which produces the instability does not need to be very high (ΔΩ/Ω <  1) and can be localized in space. If such a gradient appears in a stellar radiative zone and persists for a few hundreds of rotation periods, we predict that an instability could develop and strongly modify the axisymmetric magnetic field despite the stable stratification.

As far as Ap and Bp stars are concerned, we predict here that the instability could appear in stars for which the Lorentz number is less than 10⁻³, meaning that the Alfvén frequency should be 1000 times smaller than the rotation frequency. As argued in Sect. 3.1, small Lo are indeed expected in stellar interiors especially for Vega-like stars which rotate rapidly and exhibit a small surface magnetic field. The Lorentz number is even smaller when deep layers of the stars are considered, where latitudinal shears could be locally generated and likely to be unstable. The existence of an instability for low Lo stars would then possibly explain why strong fields are observed only for about 10% of intermediate-mass and massive stars, these stars having potentially sufficiently high magnetic frequency compared to their rotation frequency so that the instability does not reach the level of the axisymmetric field. The present study also potentially applies to the angular momentum transport in evolved stars. Although the turbulent transport associated with the MRI is not quantified here, various studies (Rüdiger et al. 2014, 2015; Jouve et al. 2015) have shown that the MRI of a toroidal field as seen here could produce a significant transport of angular momentum, which could possibly help to reconcile models and observations of the differential rotation of sub-giant and red giant stars observed with Kepler.

Acknowledgments

The authors acknowledge financial support from the Agence Nationale de la Recherche (ANR) through the project IMAGINE (Investigating MAGnetism of INtErmediate-mass and massive stars). This work was granted access to the HPC resources of CALMIP supercomputing center under the allocation P1118. LJ acknowledges funding by the Institut Universitaire de France.

References

Appendix A: Full set of non-axisymmetric MHD Boussinesq equations

We give in this appendix the full set of non-axisymmetric Boussinesq equations solved in this work, using the adimensionalisation detailed in the main text in Sect. 2.

For the three components of the velocity field, separating the toroidal and poloidal dynamics, the equations read:

$$\begin{array}{cc}\hfill \frac{\partial {\stackrel{\sim }{v}}_r}{\partial \stackrel{\sim }{t}}& +{\stackrel{\sim }{\mathit{v}}}_{\mathit{m}}·\stackrel{\sim }{\mathbf{\nabla }}{\stackrel{\sim }{v}}_r+\frac{1}{\mathrm{Lo}}\frac{{\stackrel{\sim }{v}}_{\phi }}{\stackrel{\sim }{r}sin\theta }\frac{\partial {\stackrel{\sim }{v}}_r}{\partial \phi }-\frac{{\stackrel{\sim }{v}}_{\theta }^2}{\stackrel{\sim }{r}}-\frac{1}{{\mathrm{Lo}}^2}\frac{{\stackrel{\sim }{v}}_{\phi }^2}{\stackrel{\sim }{r}}=\frac{2}{{\mathrm{Lo}}^2}sin\theta {\stackrel{\sim }{v}}_{\phi }\hfill \\ \hfill & +\frac{1}{{\mathrm{Lo}}^2}{\left(\frac{N}{{\mathrm{\Omega }}_0}\right)}^2{\stackrel{\sim }{T}}_1-\frac{1}{{\mathrm{Lo}}^2}\frac{\partial {\stackrel{\sim }{p}}_1}{\partial \stackrel{\sim }{r}}+\frac{{\stackrel{\sim }{B}}_{\theta }}{\stackrel{\sim }{r}}(\frac{\partial {\stackrel{\sim }{B}}_r}{\partial \theta }-\frac{\partial (\stackrel{\sim }{r}{\stackrel{\sim }{B}}_{\theta })}{\partial \stackrel{\sim }{r}})\hfill \\ \hfill & -\frac{1}{{\mathrm{Lo}}^2}\frac{{\stackrel{\sim }{B}}_{\phi }}{\stackrel{\sim }{r}}\frac{\partial (\stackrel{\sim }{r}{\stackrel{\sim }{B}}_{\phi })}{\partial \stackrel{\sim }{r}}+\frac{1}{\mathrm{Lo}}\frac{{\stackrel{\sim }{B}}_{\phi }}{\stackrel{\sim }{r}sin\theta }\frac{\partial {\stackrel{\sim }{B}}_r}{\partial \phi }+\frac{\mathrm{Pm}}{\mathrm{Lu}}{\stackrel{\sim }{\mathrm{\Delta }}\stackrel{\sim }{\mathit{v}}|}_r,\hfill \end{array}$$(A.1)

$$\begin{array}{cc}\hfill \frac{\partial {\stackrel{\sim }{v}}_{\theta }}{\partial \stackrel{\sim }{t}}& +{\stackrel{\sim }{\mathit{v}}}_{\mathit{m}}·\stackrel{\sim }{\mathbf{\nabla }}{\stackrel{\sim }{v}}_{\theta }+\frac{1}{\mathrm{Lo}}\frac{{\stackrel{\sim }{v}}_{\phi }}{\stackrel{\sim }{r}sin\theta }\frac{\partial {\stackrel{\sim }{v}}_{\theta }}{\partial \phi }+\frac{{\stackrel{\sim }{v}}_{\theta }{\stackrel{\sim }{v}}_r}{\stackrel{\sim }{r}}-\frac{1}{{\mathrm{Lo}}^2}cot\theta \frac{{\stackrel{\sim }{v}}_{\phi }^2}{\stackrel{\sim }{r}}=\frac{2}{{\mathrm{Lo}}^2}cos\theta {\stackrel{\sim }{v}}_{\phi }\hfill \\ \hfill & -\frac{1}{{\mathrm{Lo}}^2}\frac{1}{\stackrel{\sim }{r}}\frac{\partial {\stackrel{\sim }{p}}_1}{\partial \theta }+\frac{{\stackrel{\sim }{B}}_r}{\stackrel{\sim }{r}}(\frac{\partial (\stackrel{\sim }{r}{\stackrel{\sim }{B}}_{\theta })}{\partial \stackrel{\sim }{r}}-\frac{\partial {\stackrel{\sim }{B}}_r}{\partial \theta })\hfill \\ \hfill & -\frac{1}{{\mathrm{Lo}}^2}\frac{{\stackrel{\sim }{B}}_{\phi }}{\stackrel{\sim }{r}sin\theta }\frac{\partial (sin\theta {\stackrel{\sim }{B}}_{\phi })}{\partial \theta }+\frac{1}{\mathrm{Lo}}\frac{{\stackrel{\sim }{B}}_{\phi }}{\stackrel{\sim }{r}sin\theta }\frac{\partial {\stackrel{\sim }{B}}_{\theta }}{\partial \phi }+\frac{\mathrm{Pm}}{\mathrm{Lu}}{\stackrel{\sim }{\mathrm{\Delta }}\stackrel{\sim }{\mathit{v}}|}_{\theta },\hfill \end{array}$$(A.2)

$$\begin{array}{cc}\hfill \frac{\partial {\stackrel{\sim }{v}}_{\phi }}{\partial \stackrel{\sim }{t}}& +\frac{1}{\stackrel{\sim }{r}sin\theta }({\stackrel{\sim }{\mathit{v}}}_{\mathit{m}}·\stackrel{\sim }{\mathbf{\nabla }})(\stackrel{\sim }{r}sin\theta {\stackrel{\sim }{v}}_{\phi })+\frac{1}{{\mathrm{Lo}}^2}\frac{{\stackrel{\sim }{v}}_{\phi }}{\stackrel{\sim }{r}sin\theta }\frac{\partial {\stackrel{\sim }{v}}_{\phi }}{\partial \phi }=-2(sin\theta {\stackrel{\sim }{v}}_r+cos\theta {\stackrel{\sim }{v}}_{\theta })-\frac{1}{\stackrel{\sim }{r}sin\theta }\frac{\partial {\stackrel{\sim }{p}}_1}{\partial \phi }\hfill \\ \hfill & +\frac{1}{\stackrel{\sim }{r}sin\theta }({\stackrel{\sim }{\mathit{B}}}_{\mathit{p}}·\stackrel{\sim }{\mathbf{\nabla }})(\stackrel{\sim }{r}sin\theta {\stackrel{\sim }{B}}_{\phi })-\mathrm{Lo}\frac{{\stackrel{\sim }{B}}_{\theta }}{\stackrel{\sim }{r}sin\theta }\frac{\partial {\stackrel{\sim }{B}}_{\theta }}{\partial \phi }-\mathrm{Lo}\frac{{\stackrel{\sim }{B}}_r}{\stackrel{\sim }{r}sin\theta }\frac{\partial {\stackrel{\sim }{B}}_r}{\partial \phi }+\frac{\mathrm{Pm}}{\mathrm{Lu}}{\stackrel{\sim }{\mathrm{\Delta }}\stackrel{\sim }{\mathit{v}}|}_{\phi }.\hfill \end{array}$$(A.3)

The equations for the 3 components of the magnetic field then read:

$$\begin{array}{cc}\hfill \frac{\partial {\stackrel{\sim }{B}}_r}{\partial \stackrel{\sim }{t}}=& \frac{1}{\stackrel{\sim }{r}sin\theta }\frac{\partial }{\partial \theta }(sin\theta {\stackrel{\sim }{v}}_r{\stackrel{\sim }{B}}_{\theta }-sin\theta {\stackrel{\sim }{v}}_{\theta }{\stackrel{\sim }{B}}_r)\hfill \\ \hfill & -\frac{1}{\mathrm{Lo}}\frac{1}{\stackrel{\sim }{r}sin\theta }\frac{\partial }{\partial \phi }({\stackrel{\sim }{v}}_{\phi }{\stackrel{\sim }{B}}_r-{\stackrel{\sim }{v}}_r{\stackrel{\sim }{B}}_{\phi })+\frac{1}{\mathrm{Lu}}{\stackrel{\sim }{\mathrm{\Delta }}\stackrel{\sim }{\mathit{B}}|}_r,\hfill \end{array}$$(A.4)

$$\begin{array}{cc}\hfill \frac{\partial {\stackrel{\sim }{B}}_{\theta }}{\partial \stackrel{\sim }{t}}=& -\frac{1}{\stackrel{\sim }{r}}\frac{\partial }{\partial \theta }(\stackrel{\sim }{r}{\stackrel{\sim }{v}}_r{\stackrel{\sim }{B}}_{\theta }-\stackrel{\sim }{r}{\stackrel{\sim }{v}}_{\theta }{\stackrel{\sim }{B}}_r)\hfill \\ \hfill & +\frac{1}{\mathrm{Lo}}\frac{1}{\stackrel{\sim }{r}sin\theta }\frac{\partial }{\partial \phi }({\stackrel{\sim }{v}}_{\theta }{\stackrel{\sim }{B}}_{\phi }-{\stackrel{\sim }{v}}_{\phi }{\stackrel{\sim }{B}}_{\theta })+\frac{1}{\mathrm{Lu}}{\stackrel{\sim }{\mathrm{\Delta }}\stackrel{\sim }{\mathit{B}}|}_{\theta },\hfill \end{array}$$(A.5)

$$\begin{array}{cc}\hfill \frac{\partial {\stackrel{\sim }{B}}_{\phi }}{\partial \stackrel{\sim }{t}}=& \frac{1}{\stackrel{\sim }{r}}\frac{\partial }{\partial \stackrel{\sim }{r}}(\stackrel{\sim }{r}{\stackrel{\sim }{v}}_{\phi }{\stackrel{\sim }{B}}_r-\stackrel{\sim }{r}{\stackrel{\sim }{v}}_r{\stackrel{\sim }{B}}_{\phi })\hfill \\ \hfill & -\frac{1}{\stackrel{\sim }{r}}\frac{\partial }{\partial \theta }({\stackrel{\sim }{v}}_{\theta }{\stackrel{\sim }{B}}_{\phi }-{\stackrel{\sim }{v}}_{\phi }{\stackrel{\sim }{B}}_{\theta })+\frac{1}{\mathrm{Lu}}{\stackrel{\sim }{\mathrm{\Delta }}\stackrel{\sim }{\mathit{B}}|}_{\phi }.\hfill \end{array}$$(A.6)

And finally the temperature equation reads:

$$\begin{array}{c}\hfill \frac{\partial {\stackrel{\sim }{T}}_1}{\partial \stackrel{\sim }{t}}+{\stackrel{\sim }{\mathit{v}}}_{\mathit{m}}·\stackrel{\sim }{\mathbf{\nabla }}{\stackrel{\sim }{T}}_1+\frac{1}{\mathrm{Lo}}\frac{{\stackrel{\sim }{v}}_{\phi }}{\stackrel{\sim }{r}sin\theta }\frac{\partial \stackrel{\sim }{T_1}}{\partial \phi }+\stackrel{\sim }{v_r}\frac{\partial \overline{T}}{\partial \stackrel{\sim }{r}}=\frac{\mathrm{Pm}}{\mathrm{Lu}\mathrm{Pr}}\stackrel{\sim }{\mathrm{\Delta }}{\stackrel{\sim }{T}}_1,\end{array}$$(A.7)

where v_m = v_re_r + v_θe_θ is the meridional velocity field and B_p = B_re_r + B_θe_θ is the poloidal magnetic field. The tildes indicate the dimensionless quantities. We note that the choice of reference scales in this appendix is slightly different from the one chosen in the next appendix where a scaling analysis of the axisymmetric version of the equations is performed. The variables with a tilde in this appendix are thus different from the tilde-variables of Appendix B.

Appendix B: Scaling analysis of the axisymmetric MHD Boussinesq equations

In the following, we present a scale analysis of the axisymmetric MHD equations with the aim of finding a simplified form of these equations that approximates the evolution of our system. We note that the choice of reference scales to make the axisymmetric equations dimensionless will be slightly different in this appendix than the choice given in Sect. 2 which enabled to produce the full non-axisymmetric set of equations of Appendix A.

The initial conditions provide the characteristic magnitude of some variables: the poloidal field B₀, the rotation rate Ω₀, the stable stratification $N=\sqrt{\alpha g_0\frac{\mathrm{\Delta }T}{d}}$, the azimuthal velocity $U_{\phi }^{*}=d\mathrm{\Delta }\mathrm{\Omega }$, the domain size and also the lengthscale of the initial gradients d = r_o − r_i. We restrict our analysis to the regime $t_{\mathrm{Ap}}=\frac{d\sqrt{4\pi \rho }}{B_0}\gg t_{\mathrm{\Omega }}=\frac{1}{{\mathrm{\Omega }}_0}$. The toroidal field has no initially prescribed amplitude and there is no physical reason to choose B₀. We anticipate instead that a characteristic amplitude is $B_{\phi }^{\ast }=d\mathrm{\Delta }\mathrm{\Omega }\sqrt{4\pi \rho }$, the magnetic field resulting from the winding-up of the initial poloidal field by the differential rotation ΔΩ over an Alfvén time t_Ap. We also need to choose a typical amplitude for the meridional motions U_m. Due to the strong stable stratification N ≥ Ω₀, we argue that U_m should be small because radial motions are efficiently limited and the mass conservation ensures that latitudinal velocities are of the same order as radial velocities, v_θ ∼ v_r. In practice assuming U_m ≪ dΔΩ ≲ dΩ₀ allows us to simplify the system of equation and to obtain U_m as a result of the scale analysis. The consistency of the assumption U_m ≪ dΔΩ ≲ dΩ₀ is verified afterwards. As demonstrated below, such small meridional velocities lead to a thermal-wind balance which in turn determines a typical amplitude for the temperature fluctuations, $T^{\ast }=\frac{{\mathrm{\Omega }}_0^2}{N^2}\frac{\mathrm{\Delta }\mathrm{\Omega }}{{\mathrm{\Omega }}_0}\mathrm{\Delta }T$, and the pressure fluctuations P^* = ρd²Ω₀ΔΩ. Finally, as we are interested in the evolution of the angular momentum, the characteristic time scale is chosen from the equation governing this evolution:

$$\begin{array}{c}\hfill \frac{\partial M}{\partial t}+({\mathit{v}}_m·\mathbf{\nabla })M=\frac{1}{4\pi \rho }({\mathit{B}}_{\mathit{p}}·\mathbf{\nabla })(rsin\theta B_{\phi })+rsin\theta \nu (\mathrm{\Delta }-\frac{1}{r^2{sin}^2\theta })v_{\phi },\end{array}$$(B.1)

where M = r² sin ²θΩ₀ + r sin θv_φ is the specific angular momentum. The meridional velocity, v_m = v_re_r + v_θe_θ, advects the angular momentum on a time scale t_c = (ΔΩ/Ω)(d/U_m), where the factor $\mathrm{Ro}=\frac{\mathrm{\Delta }\mathrm{\Omega }}{\mathrm{\Omega }}$ accounts for the effect of the Coriolis force that speed-up the transport when Ro <  1. In our simulations, the initial differential rotation is such that Ro ∼ 1 while in the following we consider more generally Ro ≤ 1 regimes. The other time scale that controls the angular momentum evolution is the poloidal Alfvén time t_Ap as the time over which the toroidal field produced by the Ω-effect back reacts onto the rotation. The third time scale is the viscous time t_ν and it is supposed to be always larger than t_Ap. Consequently, the relevant time scale to study the angular momentum evolution should be either t_c or t_Ap. We don’t have to choose between these two times yet. But as we already assumed that t_c ≫ t_Ω (as a consequence of U_m ≪ dRoΩ₀) and t_Ap ≫ t_Ω, we can safely assume that the characteristic time of the angular momentum evolution, denoted t_*, verifies t_* ≫ t_Ω.

With these choices, the scaled version of the radial and latitudinal components of the MHD Boussinesq equations read:

$$\begin{array}{cc}& \frac{t_{\mathrm{\Omega }}^2}{t_{\ast }{t}_{\mathrm{c}}}\frac{\partial {\stackrel{\sim }{v}}_r}{\partial \stackrel{\sim }{t}}+\mathrm{Ro}{\left(\frac{t_{\mathrm{\Omega }}}{t_{\mathrm{c}}}\right)}^2({\stackrel{\sim }{\mathit{v}}}_{\mathit{m}}·\stackrel{\sim }{\mathbf{\nabla }}{\stackrel{\sim }{v}}_r-\frac{{\stackrel{\sim }{v}}_{\theta }^2}{\stackrel{\sim }{r}})-(2sin\theta {\stackrel{\sim }{v}}_{\phi }+\mathrm{Ro}\frac{{\stackrel{\sim }{v}}_{\phi }^2}{\stackrel{\sim }{r}})={\stackrel{\sim }{T}}_1-\frac{\partial {\stackrel{\sim }{p}}_1}{\partial \stackrel{\sim }{r}}\hfill \\ \hfill & +\frac{1}{\mathrm{Ro}}{\left(\frac{t_{\mathrm{\Omega }}}{t_{\mathrm{Ap}}}\right)}^2\left[\frac{{\stackrel{\sim }{B}}_{\theta }}{\stackrel{\sim }{r}}(\frac{\partial {\stackrel{\sim }{B}}_r}{\partial \theta }-\frac{\partial (\stackrel{\sim }{r}{\stackrel{\sim }{B}}_{\theta })}{\partial \stackrel{\sim }{r}})\right]-\mathrm{Ro}\frac{{\stackrel{\sim }{B}}_{\phi }}{\stackrel{\sim }{r}}\frac{\partial (\stackrel{\sim }{r}{\stackrel{\sim }{B}}_{\phi })}{\partial \stackrel{\sim }{r}}+\frac{t_{\mathrm{\Omega }}^2}{t_{\nu }{t}_{\mathrm{c}}}{\stackrel{\sim }{\mathrm{\Delta }}\stackrel{\sim }{\mathit{v}}|}_r,\hfill \end{array}$$(B.2)

$$\begin{array}{cc}& \frac{t_{\mathrm{\Omega }}^2}{t_{\ast }{t}_{\mathrm{c}}}\frac{\partial {\stackrel{\sim }{v}}_{\theta }}{\partial \stackrel{\sim }{t}}+\mathrm{Ro}{\left(\frac{t_{\mathrm{\Omega }}}{t_{\mathrm{c}}}\right)}^2({\stackrel{\sim }{\mathit{v}}}_{\mathit{m}}·\stackrel{\sim }{\mathbf{\nabla }}{\stackrel{\sim }{v}}_{\theta }+\frac{{\stackrel{\sim }{v}}_{\theta }{\stackrel{\sim }{v}}_r}{\stackrel{\sim }{r}})-(2cos\theta {\stackrel{\sim }{v}}_{\phi }+\mathrm{Ro}cot\theta \frac{{\stackrel{\sim }{v}}_{\phi }^2}{\stackrel{\sim }{r}})=-\frac{1}{r}\frac{\partial {\stackrel{\sim }{p}}_1}{\partial \theta }\hfill \\ \hfill & +\frac{1}{\mathrm{Ro}}{\left(\frac{t_{\mathrm{\Omega }}}{t_{\mathrm{Ap}}}\right)}^2\left[\frac{{\stackrel{\sim }{B}}_r}{\stackrel{\sim }{r}}(\frac{\partial (\stackrel{\sim }{r}{\stackrel{\sim }{B}}_{\theta })}{\partial \stackrel{\sim }{r}}-\frac{\partial {\stackrel{\sim }{B}}_r}{\partial \theta })\right]-\mathrm{Ro}\frac{{\stackrel{\sim }{B}}_{\phi }}{\stackrel{\sim }{r}sin\theta }\frac{\partial (sin\theta {\stackrel{\sim }{B}}_{\phi })}{\partial \theta }+\frac{t_{\mathrm{\Omega }}^2}{t_{\nu }{t}_{\mathrm{c}}}{\stackrel{\sim }{\mathrm{\Delta }}\stackrel{\sim }{\mathit{v}}|}_{\theta }.\hfill \end{array}$$(B.3)

From these expressions, the inertial terms that do not involve the azimuthal velocity can be neglected because t_c ≫ t_Ω and t_* ≫ t_Ω. Moreover, the viscous terms is negligible if t_ν ≫ t_Ω, and, as long as Ro is finite and non-zero, the term of the Lorentz force that contains the poloidal field is very small because t_Ap ≫ t_Ω. We thus simplify Eqs. (B.2), (B.3) into:

$$\begin{array}{cc}& -2sin\theta {\stackrel{\sim }{v}}_{\phi }-\frac{\mathrm{Ro}}{\stackrel{\sim }{r}}({\stackrel{\sim }{v}}_{\phi }^2-{\stackrel{\sim }{B}}_{\phi }^2)={\stackrel{\sim }{T}}_1-\frac{\partial }{\partial \stackrel{\sim }{r}}({\stackrel{\sim }{p}}_1+\frac{\mathrm{Ro}}{2}{\stackrel{\sim }{B}}_{\phi }^2),\hfill \end{array}$$(B.4)

$$\begin{array}{cc}& -2cos\theta {\stackrel{\sim }{v}}_{\phi }-\mathrm{Ro}\frac{cot\theta }{\stackrel{\sim }{r}}({\stackrel{\sim }{v}}_{\phi }^2-{\stackrel{\sim }{B}}_{\phi }^2)=-\frac{1}{\stackrel{\sim }{r}}\frac{\partial }{\partial \theta }({\stackrel{\sim }{p}}_1+\frac{\mathrm{Ro}}{2}{\stackrel{\sim }{B}}_{\phi }^2).\hfill \end{array}$$(B.5)

The pressure terms, including the magnetic pressure, can be eliminated to get a magneto-thermal wind equation that relates the temperature fluctuations to the differential rotation and the azimuthal field. This relation has been anticipated to determine the characteristic temperature fluctuation $T^{\ast }=\frac{{\mathrm{\Omega }}_0^2}{N^2}\frac{\mathrm{\Delta }\mathrm{\Omega }}{{\mathrm{\Omega }}_0}\mathrm{\Delta }T$ associated with the differential rotation. We now turn to the thermal energy equation that relates temperature fluctuations and meridional velocities:

$$\begin{array}{c}\hfill \frac{t_{\kappa }}{t_{\ast }}\frac{\partial {\stackrel{\sim }{T}}_1}{\partial \stackrel{\sim }{t}}+\mathrm{Ro}\frac{t_{\kappa }}{t_{\mathrm{c}}}{\stackrel{\sim }{\mathit{v}}}_{\mathit{m}}·\stackrel{\sim }{\mathbf{\nabla }}{\stackrel{\sim }{T}}_1+\frac{N^2}{{\mathrm{\Omega }}^2}\frac{t_{\kappa }}{t_{\mathrm{c}}}\stackrel{\sim }{v_r}\frac{d\overline{T}}{d\stackrel{\sim }{r}}=\stackrel{\sim }{\mathrm{\Delta }}{\stackrel{\sim }{T}}_1.\end{array}$$(B.6)

The advection of the temperature has been split into the advection of temperature fluctuations by meridional motions and the radial advection against the background stratification. This last term is expected to dominate the advection if $\mathrm{Ro}\frac{{\mathrm{\Omega }}^2}{N^2}\ll 1$. Then, depending on the ratio t_κ/t_*, it can be balanced either by the time derivative of temperature fluctuations or by the thermal diffusion term. The two cases are now considered separately

B.1. Alfvén waves and Eddington–Sweet circulation

If t_* ≫ t_κ the thermal diffusion term dominates over the temperature time variation in Eq. (B.6). Thus the balance between the thermal diffusion transport and the radial advection against the background stratification determines the circulation time $t_{\mathrm{c}}=t_{\kappa }\frac{N^2}{{\mathrm{\Omega }}^2}$ and the characteristic meridional velocity $U_m=\frac{\kappa }{d}\frac{{\mathrm{\Omega }}^2}{N^2}\frac{\mathrm{\Delta }\mathrm{\Omega }}{\mathrm{\Omega }}$. The scaled thermal energy equation is then:

$$\begin{array}{c}\hfill \stackrel{\sim }{v_r}\frac{d\overline{T}}{d\stackrel{\sim }{r}}=\stackrel{\sim }{\mathrm{\Delta }}{\stackrel{\sim }{T}}_1,\\ \hfill \end{array}$$(B.7)

where the circulation appears driven by the thermal diffusion of the temperature deviations, that were produced by the differential rotation. It is an Eddington–Sweet type circulation of time scale $t_{\mathrm{c}}=t_{\mathrm{es}}=\frac{d^2}{\kappa }\frac{N^2}{{\mathrm{\Omega }}^2}$. We can now verify that the meridional circulation satisfies the condition t_c ≫ t_Ω necessary to simplify Eqs. (B.4), (B.5) if t_es ≫ t_Ω. This is satisfied in stars because $t_{\kappa }=\frac{d^2}{\kappa }\gg t_{\mathrm{\Omega }}$ and N ≥ Ω. The system of equation is completed by three prognostic equations for v_φ, B_φ and the potential A, defined by B_p = ∇ × Ae_φ.

Their scaled form is:

$$\begin{array}{cc}& \frac{\partial {\stackrel{\sim }{v}}_{\phi }}{\partial \stackrel{\sim }{t}}+\mathrm{Ro}\frac{t_{\mathrm{Ap}}}{t_{\mathrm{es}}}\frac{1}{\stackrel{\sim }{r}sin\theta }({\stackrel{\sim }{\mathit{v}}}_{\mathit{m}}·\stackrel{\sim }{\mathbf{\nabla }})(\stackrel{\sim }{r}sin\theta {\stackrel{\sim }{v}}_{\phi })+2\frac{t_{\mathrm{Ap}}}{t_{\mathrm{es}}}(sin\theta {\stackrel{\sim }{v}}_r+cos\theta {\stackrel{\sim }{v}}_{\theta })=\frac{1}{\stackrel{\sim }{r}sin\theta }({\stackrel{\sim }{\mathit{B}}}_{\mathit{p}}·\stackrel{\sim }{\mathbf{\nabla }})(\stackrel{\sim }{r}sin\theta {\stackrel{\sim }{B}}_{\phi })+\frac{t_{\mathrm{Ap}}}{t_{\nu }}(\stackrel{\sim }{\mathrm{\Delta }}-\frac{1}{{\stackrel{\sim }{r}}^2{sin}^2\theta }){\stackrel{\sim }{v}}_{\phi },\hfill \end{array}$$(B.8)

$$\begin{array}{cc}& \frac{\partial {\stackrel{\sim }{B}}_{\phi }}{\partial \stackrel{\sim }{t}}+\mathrm{Ro}\frac{t_{\mathrm{Ap}}}{t_{\mathrm{es}}}\stackrel{\sim }{r}sin\theta ({\stackrel{\sim }{\mathit{v}}}_{\mathit{m}}·\stackrel{\sim }{\mathbf{\nabla }})\left(\frac{{\stackrel{\sim }{B}}_{\phi }}{\stackrel{\sim }{r}sin\theta }\right)=\stackrel{\sim }{r}sin\theta ({\stackrel{\sim }{\mathit{B}}}_{\mathit{p}}·\stackrel{\sim }{\mathbf{\nabla }})\left(\frac{{\stackrel{\sim }{v}}_{\phi }}{\stackrel{\sim }{r}sin\theta }\right)+\frac{t_{\mathrm{Ap}}}{t_{\eta }}(\stackrel{\sim }{\mathrm{\Delta }}-\frac{1}{{\stackrel{\sim }{r}}^2{sin}^2\theta }){\stackrel{\sim }{B}}_{\phi },\hfill \end{array}$$(B.9)

$$\begin{array}{cc}& \frac{\partial \stackrel{\sim }{A}}{\partial \stackrel{\sim }{t}}+\mathrm{Ro}\frac{t_{\mathrm{Ap}}}{t_{\mathrm{es}}}\frac{1}{\stackrel{\sim }{r}sin\theta }({\stackrel{\sim }{\mathit{v}}}_{\mathit{m}}·\stackrel{\sim }{\mathbf{\nabla }})(\stackrel{\sim }{r}sin\theta \stackrel{\sim }{A})=\frac{t_{\mathrm{Ap}}}{t_{\eta }}(\stackrel{\sim }{\mathrm{\Delta }}-\frac{1}{{\stackrel{\sim }{r}}^2{sin}^2\theta })\stackrel{\sim }{A},\hfill \end{array}$$(B.10)

where we used t_* = t_Ap (but could also have used t_* = t_es).

The scale analysis thus led to a simplified system formed by Eqs. (B.4), (B.5) and (B.7)–(B.10), plus the mass conservation equation, $\stackrel{\sim }{\mathbf{\nabla }}·{\stackrel{\sim }{\mathit{v}}}_{\mathit{m}}=0$. To be consistent the approximations requires $t_{\mathrm{Ap}}\gg t_{\kappa }\gg t_{\mathrm{\Omega }}\gg {\mathrm{Ro}}^{\frac{1}{2}}{t}_B$ together with t_ν ≫ t_Ω and ${\mathrm{Ro}}^{\frac{1}{2}}{t}_{\mathrm{Ap}}\gg t_{\mathrm{\Omega }}$. It intends to describe axisymmetric motions for time scale of the order of t_* = t_Ap. In particular, it should fail when solid body rotation is reached because the Lorentz force term involving the poloidal field component in Eqs. (B.4) and (B.5) will no longer be negligible. Also, short time dynamics like gravity waves have been filtered out by the approximation of the scaling analysis.

In this simplified form, the system is fully determined by the azimuthal velocity and the two components of the magnetic field. The meridional velocity components and the temperature fluctuations are intermediate variables determined by the magneto-thermal wind equilibrium and the thermal equilibrium. Flows, where such equilibrium equation reduces the number of independent variables, are said to have balanced dynamics (e.g., Vallis 2006). Physically, the flow evolves through Alfvén wave dynamics and an Eddington Sweet circulation prescribed by the instantaneous angular momentum and azimuthal field distributions.

As compared to the full MHD problem that depends on five non-dimensional numbers (plus Ro = ΔΩ/Ω₀), this simplified system has the advantage to depend only on three non-dimensional numbers $\frac{t_{\mathrm{es}}}{t_{\mathrm{Ap}}}=\frac{\mathrm{Lu}}{\mathrm{Pm}}\mathrm{Pr}{N}^2/{\mathrm{\Omega }}_0^2$, $\frac{t_{\eta }}{t_{\mathrm{Ap}}}=\mathrm{Lu}$ and $\frac{t_{\nu }}{t_{\mathrm{Ap}}}=\frac{\mathrm{Lu}}{\mathrm{Pm}}$, or equivalently on $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2$, Lu and Pm. Consequently, for given initial conditions and thus a given Ro, solutions can be expressed in the general form ${\stackrel{\sim }{v}}_{\phi }=\frac{v_{\phi }}{d\mathrm{\Delta }\mathrm{\Omega }}=f_0(t/t_{\mathrm{Ap}},\mathit{r}/d,\mathrm{Pr}{N}^2/{\mathrm{\Omega }}_0^2,\mathrm{Lu},\mathrm{Pm})$, ${\stackrel{\sim }{B}}_{\phi }=\frac{B_{\phi }}{d\mathrm{\Delta }\mathrm{\Omega }\sqrt{4\pi \rho }}=f_1(t/t_{\mathrm{Ap}},\mathit{r}/d,\mathrm{Pr}{N}^2/{\mathrm{\Omega }}_0^2,\mathrm{Lu},\mathrm{Pm})$, ${\stackrel{\sim }{\mathit{B}}}_{\mathrm{p}}=\frac{{\mathit{B}}_{\mathrm{p}}}{B_0}=f_2(t/t_{\mathrm{Ap}},\mathit{r}/d,\mathrm{Pr}{N}^2/{\mathrm{\Omega }}_0^2,\mathrm{Lu},\mathrm{Pm})$ from which we deduce $B_{\phi }/B_{\mathrm{p}}={\mathrm{Lo}}^{-1}f(t/t_{\mathrm{Ap}},\mathit{r}/d,\mathrm{Pr}{N}^2/{\mathrm{\Omega }}_0^2,\mathrm{Lu},\mathrm{Pm})$, that is the expression given in Sect. 3.3.

Most of the numerical simulations listed in Table 1 meet the requirement of the scaling analysis as they verify t_ν ≫ t_Ap ≫ t_κ ≫ t_Ω ≫ t_B together with Ro = (Ω_i − Ω₀)/Ω₀ ≈ 1. Except for the transient period during which initially excited gravity waves are dissipated, their dependence on $\text{Pr}{N}^2/{\mathrm{\Omega }}_0^2$ and Lo indicate that they are indeed governed by the simplified equations derived from the present scale analysis.

Below, we consider the case t_Ap ≤ t_κ. It holds in particular for the run R9 of Table 1 for which $\frac{t_{\kappa }}{t_{\mathrm{Ap}}}=\frac{\mathrm{Lu}\mathrm{Pr}}{\mathrm{Pm}}=3.125$.

B.2. Alfvén waves

If t_* ≪ t_κ, the balanced thermal energy equation is:

$$\begin{array}{c}\hfill \frac{\partial {\stackrel{\sim }{T}}_1}{\partial \stackrel{\sim }{t}}+\stackrel{\sim }{v_r}\frac{\mathrm{d}\overline{T}}{\mathrm{d}\stackrel{\sim }{r}}=0,\end{array}$$(B.11)

with $t_{\mathrm{c}}=t_{\ast }\frac{N^2}{{\mathrm{\Omega }}^2}$ or equivalently $U_m=\frac{{\mathrm{\Omega }}^2}{N^2}\frac{\mathrm{\Delta }\mathrm{\Omega }}{\mathrm{\Omega }}\frac{d}{t_{\ast }}$. Then, the conditions t_* ≫ t_Ω and t_c ≫ t_Ω necessary to simplify Eqs. (B.4), (B.5), now read $t_{\ast }\gg \frac{t_B^2}{t_{\mathrm{\Omega }}}$. As N ≥ Ω, we have t_c >  t_* which implies that t_* = t_Ap is the more relevant choice for the time scale characterizing the angular momentum evolution. The condition t_* = t_Ap ≫ t_Ω is met because t_Ap ≫ t_Ω and t_Ω ≥ t_B. The amplitude of the meridional motion is now $U_m=\frac{{\mathrm{\Omega }}^2}{N^2}\frac{\mathrm{\Delta }\mathrm{\Omega }}{\mathrm{\Omega }}{v}_{\mathrm{Ap}}$. The scaled version of the three prognostic equations for v_φ, B_φ and B_p = ∇ × Ae_φ simplifies into:

$$\begin{array}{cc}& \frac{\partial {\stackrel{\sim }{v}}_{\phi }}{\partial \stackrel{\sim }{t}}+2\frac{{\mathrm{\Omega }}^2}{N^2}(sin\theta {\stackrel{\sim }{v}}_r+cos\theta {\stackrel{\sim }{v}}_{\theta })=\frac{1}{\stackrel{\sim }{r}sin\theta }({\stackrel{\sim }{\mathit{B}}}_{\mathit{p}}·\stackrel{\sim }{\mathbf{\nabla }})(\stackrel{\sim }{r}sin\theta {\stackrel{\sim }{B}}_{\phi })+\frac{t_{\mathrm{Ap}}}{t_{\nu }}(\stackrel{\sim }{\mathrm{\Delta }}-\frac{1}{{\stackrel{\sim }{r}}^2{sin}^2\theta }){\stackrel{\sim }{v}}_{\phi },\hfill \end{array}$$(B.12)

$$\begin{array}{cc}& \frac{\partial {\stackrel{\sim }{B}}_{\phi }}{\partial \stackrel{\sim }{t}}=\stackrel{\sim }{r}sin\theta ({\stackrel{\sim }{\mathit{B}}}_{\mathit{p}}·\stackrel{\sim }{\mathbf{\nabla }})\left(\frac{{\stackrel{\sim }{v}}_{\phi }}{\stackrel{\sim }{r}sin\theta }\right)+\frac{t_{\mathrm{Ap}}}{t_{\eta }}(\stackrel{\sim }{\mathrm{\Delta }}-\frac{1}{{\stackrel{\sim }{r}}^2{sin}^2\theta }){\stackrel{\sim }{B}}_{\phi },\hfill \end{array}$$(B.13)