Interplay between magnetic fields and differential rotation in a stably stratified stellar radiative zone

The present study aims at studying the flow and field produced by a stellar radiative zone which is initially made to rotate differentially in the presence of a large-scale poloidal magnetic field threading the whole domain. We focus both on the axisymmetric configurations produced by the initial winding-up of the magnetic field lines and on the possible instabilities of those configurations. We aim in particular at assessing the role of the stratification at stabilising the system. We perform 2D and 3D global Boussinesq numerical simulations started from an initial radial or cylindrical differential rotation and a large-scale poloidal magnetic field. Under the conditions of a large rotation frequency compared to the Alfv\'en frequency, a magnetic configuration strongly dominated by its toroidal component is built. The parameters of the simulations are chosen to respect the ordering of time scales of a typical stellar radiative zone. In this framework, the axisymmetric evolution is studied by varying the relative effects of the thermal diffusion, the Br\"unt-V\"ais\"al\"a frequency, the rotation and the initial poloidal field strength. We find that the axisymmetric state only depends on $t_{es}/t_{Ap}$, the ratio between the Eddington-Sweet circulation time scale and the Alfv\'en time scale. A scale analysis of the Boussinesq equations allows us to recover this result. In the cylindrical case, a magneto-rotational instability develops when the thermal diffusivity is sufficiently high to enable the favored wavenumbers to be insensitive to the effects of the stable stratification. In the radial case, the magneto-rotational instability is driven by the latitudinal shear created by the back-reaction of the Lorentz force on the flow. Increasing the level of stratification then leaves the growth rate of the instability mainly unaffected while its horizontal length scale grows.


Introduction
Considerable progress has been made recently about the knowledge 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). There is also now 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 to 10% of these stars do exhibit a strong surface magnetic field above 300G 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 2 classes of magnetism could exist in intermediatemass 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 2 types of magnetism and the absence of stars possessing fields with amplitudes Jouve Lignières Gaurat: Magnetic fields in stably stratified radiative zones between approximately 1G and 300G 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, very likely to be 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 only non-axisymmetric hydro or MHD instabilities which, in conjonction 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). It is for now still debated if 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/expansion phenomena (Deheuvels et al., 2012(Deheuvels et al., , 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). A more efficient transport of angular momentum seems then to be at play in those stellar radiative zones and magnetic fields are seriously considered as interesting candidates to play this role. In particular, the transport by travelling 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 unclear conclusions 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 has been published recently by Rüdiger et al. (2018). The main instabilities described in that review are the currentdriven 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/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 analysis were local in radius but global in the horizontal directions and took into account the effects of stratification, focusing in particular 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 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/toroidal magnetic field configuration, heat conductivity, viscosity and magnetic diffusivity) is then let free to evolve into potentially unstable equilibria. Our recent numerical studies 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. (2019) 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 similar to what is presented in this paper, the woundup 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 so 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 yet 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 Brünt-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 backreacts 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 Brünt-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. Sect.3 then discusses the axisymmetric joint evolution of the magnetic field and the flow. We then inves-tigate the stability of this axisymmetric configuration in Sect.4 and finally conclude in Sect.5.

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 let 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 Sections 2.2 and 2.3 and the numerical method is finally briefly described in Sect. 2.4.

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 where v is the velocity field, B is the magnetic field, Ω 0 is the rotation rate at the rotation axis, T (r, θ, t) = T (r)+T 1 (r, θ, t) is the temperature field with T (r) the temperature of the reference state and T 1 its fluctuation, ρ is the uniform density of the reference state, p 1 is the pressure fluctuation, gravity is proportional to 1/r 2 , α 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 Ap = d √ 4πρ/B 0 as the time unit where the surface radial magnetic field at the poles B 0 is the poloidal magnetic field unit, dΩ 0 √ 4πρ as the toroidal magnetic field unit, V Ap = d/t Ap as the meridional circulation unit, dΩ 0 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 Ω 2 0 ρ as the pressure unit. The full set of governing equations of the problem is given in appendix A, namely the equations for the 3 components of the velocity field, for the 3 components of the magnetic field and for the temperature field.
Five dimensionless numbers appear in the set of equations: 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/Ω 0 is the ratio between the Brünt-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 2 /η and finally the Prandtl numbers quantify the ratio of diffusivities or the ratio of diffusive time-scales where t ν = d 2 /ν is the viscous time-scale and t κ = d 2 /κ 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 2 = Lo/(P mLu). We choose in this study to fix the values of 2 dimensionless numbers, namely the Lundquist number Lu and the magnetic Prandtl number P m. We then focus on the effects of the 3 other parameters: the Lorentz number Lo, the ratio N/Ω 0 and the Prandtl number P r. 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.

Initial and boundary conditions
In this work, we focus on initial conditions which will produce a large-scale magnetic field, likely to be unstable to MHD instabilities under certain circumstances. To do so, we start from a poloidal field which will be acted upon by an initial differential rotation. The winding-up of the initial poloidal field by the differential rotation will naturally produce 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: With this choice of normalization, B 0 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 v(r, θ, t = 0) = v ϕ (r, θ, t = 0) e ϕ = r sin θ Ω(r, θ, t = 0) e ϕ (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 ∇ 2 T = 0. Fixed values are imposed for the temperature at both boundaries : Finally, to ensure that the flow is stable with respect to convection, the Brünt-Väisälä frequency N must be real which means that ∆T > 0.

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 choose 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: where Ω 0 is the rotation rate at the equator at r = r o and where c 1 = 0.980 and c 2 = 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 − Ω 0 )/Ω 0 ≈ 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 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.

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 (https://github.com/magic-sph/magic) which solves the MHD equations in a spherical shell using a poloidal toroidal decomposition for the mass flux and the magnetic fields: 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 where η = r i /r o is the aspect ratio of the shell, equal to 0.3 in all our calculations.

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 perform 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/Ω 0 and P r = ν/κ, while the two others Lu and P m = ν/η are fixed to Lu = 50 and P m = ν/η = 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.

Physical parameters in stellar radiative zones and the parameter range of our numerical simulations
Parameters such as the Brünt 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/Ω 0 is typically much larger than one, except for stars rotating near break-up velocity for which N/Ω 0 ∼ 1, meanwhile the Prandtl number Pr is always much smaller than 1. For example, stellar structure models of a 3M star indicate that, during the main sequence, N ∼ 1−2×10 −3 s −1 away from the convective core (Talon & Charbonnel, 2008) while Pr ≤ 4 × 10 −6 (Garaud et al., 2015). For rotation periods between 1 and 2.7 days, N/Ω 0 is then comprised between ∼ 10 and 75. For comparison, this ratio is much higher, N/Ω 0 ∼ 300, in the radiative zone of the Sun. As we shall see below, the product Pr(N/Ω 0 ) 2 is also a relevant pa-rameter 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 ) 2 ≤ 0.02 for a 3M main sequence star. The situation is different in the solar radiative zone where Pr(N/Ω 0 ) 2 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 2 Ω) ∼ 10 −14 , 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/Ω 0 ) 2 1. It corresponds to the following time scale order : /Ω 0 and t N = 1/N . As shown in table 1, the simulations performed respect that order.
The timescale associated with the initial poloidal field is t Ap = d √ 4πρ/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 Ap = B p / √ 4πρ is dominated by the density increase. For example, at a radius r = R /3, the Lorentz number would be 2.7 × 10 −3 assuming a density ratio of 10 8 and a dipolar-like radial increase B p ∝ 1/r 3 . 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 2 rms /(4πρr 2 c Ω 2 ) ∼ 2.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 P r (cases 2, 3 and 4), then P r and N/Ω 0 while keeping P rN 2 /Ω 2 0 fixed (cases 2, 5, 6 and R9 for the radial case). We also decrease the value of P rN 2 /Ω 2 0 (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 Fig. 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) (left panel) and the radial differential rotation defined by Eq. (14) (right 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. This evolution of the magnetic field, namely a nearlinear 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 (for example 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 for example Ionson, 1978;Spruit, 1999) finally leads to uniform rotation as nothing enforces differential rotation in our set-up. 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.
Rad 2 6.25 × 10 −2 2 × 10 −5 10 −3 2.5 × 10 −1 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, i.e. 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Ω 0 √ 4πρ. 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.
Beside these spatial distributions, the ratio of the azimuthal Alfvén frequency, ω Aϕ = B ϕ /r √ 4πρ 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) . 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 Figs. 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 backreacts 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.

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 backreaction 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  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/Ω 0 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 2 /κ)(N/Ω 0 ) 2 (see for example Spiegel & Zahn (1992)) 1 . 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 Ap /t es = P m/(LuP rN 2 /Ω 2 0 ) varies between 8×10 −4 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 P r and N/Ω 0 . A striking feature is that, after a transient phase, most of our axisymmetric solutions are controlled by  Fig. 2. Evolution of the ratio between the toroidal and poloidal magnetic energies for cases C1 (left) and R1 (right). E ϕ seems much stronger in the cylindrical case but it's because B ϕ is more extended and located closer to the surface where B p is much smaller. If we look locally, we also find that B ϕ /B p is stronger in the cylindrical case but mostly because of the location of B ϕ (where B p is small). the product P r N 2 /Ω 2 0 , rather than by the two parameters P r and N/Ω 0 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 P r N 2 /Ω 2 0 = 0.25 but different P r and N/Ω 0 , namely P r = 10 −2 , N/Ω 0 = 5 (blue), P r = 10 −3 , N/Ω 0 = 15.8 (green) and P r = 10 −4 , N/Ω 0 = 50 (red). The two cases P r = 10 −3 and P r = 10 −4 are even indistinguishable on this plot after about t = 10 −3 t Ap . Figure 4 also displays the energy evolution of two simulations at P rN 2 /Ω 2 0 = 0.04, the P r = 10 −3 (cyan) and P r = 10 −4 (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 P r N 2 /Ω 2 0 . On a short timescale however, flows having the same P rN 2 /Ω 2 0 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 P r = 10 −2 case (run R2). The thermal damping is so efficient in the case P r = 10 −4 (for which thermal diffusivity is 100 times larger than for the case P r = 10 −2 ) that the oscillations of the kinetic energy are not even visible. The invariance of the solution with P rN 2 /Ω 2 0 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 P r N 2 /Ω 2 0 = 0.25. Then, not only the kinetic and magnetic energies evolution are close for identical value of PrN 2 /Ω 2 0 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 PrN 2 /Ω 2 0 = 0.04 (two left panels) and PrN 2 /Ω 2 0 = 0.25 (two right panels), in both cases for 2 different values of Pr, namely 1.6 × 10 −3 and 1.6 × 10 −5 for the left panels (cases R7 and R8) and 10 −2 and 10 −4 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.
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 PrN 2 /Ω 2 0 , Lu and P m if time is scaled by t Ap and B ϕ by d∆Ω √ 4πρ. The Lorentz number Lo only appears as a scaling factor of the ratio The simulations that verify the required time ordering do show this PrN 2 /Ω 2 0 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 tκ t Ap = LuPr P m = 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 PrN 2 /Ω 2 0 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 PrN 2 /Ω 2 0 = 0.25 but two different Lorentz numbers (runs R1 and R2). The ratio indeed increased by a factor ≈ 2.5 2 , as Lo was reduced from 2.5 × 10 −3 (right panel of Fig.2) to Lo = 10 −3 ) (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 PrN 2 /Ω 2 0 instead of the three Lo, N/Ω 0 and Pr.
We thus consider how the rotation and magnetic configurations depend on the values of P rN 2 /Ω 2 0 . As shown on the right panel of Fig. 4, the maximum ratio of toroidal to poloidal magnetic energy is 1500 for P rN 2 /Ω 2 0 = 0.04 while it is close to 4000 for P rN 2 /Ω 2 0 = 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: ∂Ω ∂z = 0, where 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 P r N 2 /Ω 2 0 and that stable stratification favors the creation of magnetic field configurations more strongly dominated by their toroidal component.

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 P rN 2 /Ω 2 0 = 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 −3 for the cylindrical case and Lo = 10 −3 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).

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  Emagφ/Emagp Fig. 4. Temporal evolution of the kinetic energy (on a long timescale, left and zoomed in, mid-panel) and magnetic energy ratio (right) for 5 different cases: 3 cases with P rN 2 /Ω 2 0 = 0.25 (blue: R2, green: R5, red: R6 and black: R9) and 2 cases with P rN 2 /Ω 2 0 = 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 P rN 2 /Ω 2 0 , as long as P r is small enough, but different for different values of P rN 2 /Ω 2 0 . In particular, the amount of toroidal field produced is much less in the case where P rN 2 /Ω 2 0 = 0.04.
Ω/Ω 0  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 −2 Ω for the m = 4 mode in the cylindrical case and 2.5×10 −2 Ω for the m = 1 mode in the radial case.
We now emphasize the differences between the 2 cases. First, the time at which the instability starts to grow is quite different. Indeed, in the cylindrical case, the axisym-  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, i.e. at t = 0.15t ap for C2 and t = 0.35t ap for R2. The colorbar applies to the rotation rate. metric 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 the case R2 and C2 lies in the structure of the unstable modes. From the mid 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: 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 , 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 dom-inated 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) are 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 (ω Aϕ /Ω) 2 E −1 k P m 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 (ω Aϕ /Ω) 2 E −1 k P m when the instability develops in the radial case is of the order of 2 × 10 3 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 −2 ). We now wish to investigate the effect of varying the thermal diffusion on the instability.

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 P r, keeping the same value for N/Ω 0 . 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 P r = 0.1 and P r = 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 P r = 10 −2 , it is clear that for the cylindrical case, the instability is largely suppressed by the increase of the stable stratification effect. In particular, for P r = 1, the axisymmetric solution becomes completely stable to any non-axisymmetric perturbation. In other words, the preferentially radial displacements that were unstable at P r = 10 −2 are inhibited at P r = 1. The transition between P r = 10 −2 and P r = 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 : With the dimensionless parameters used in our calculations, the critical lengthscale reads : Ω N The computation of this quantity gives a value ranging from 4% to 0.4% when P r goes from 10 −2 to P r = 1, for these cases where P rN 2 /Ω 2 0 = 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 P r = 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 (left 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 P r = 10 −2 case. As we show below, this comes with the fact that the unstable displacements become more and more horizontal. On figure 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 P rN 2 /Ω 2 0 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 lenghtscale) 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 P rN 2 /Ω 2 0 = 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.
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 (i.e. increasing) the vertical wavenumber accordingly. The possibility to adapt the vertical lengthscale to get the same growth rate also exists when the shear instability 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.

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 Section 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 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 P m = 1 and P r 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 ω = ω/Ω 0 .
We solve numerically that polynomial equation 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 P rN 2 /Ω 2 0 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 Figure 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 σ/Ω 0 , 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 ≈ 0.1, compared to 2.5×10 −2 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.
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 ≈ 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, i.e. 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, i.e. with P r = 10 −1 instead of P r = 10 −2 . 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 equation C.4, we see that all the terms involving P rN 2 /Ω 2 0 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 P rN 2 /Ω 2 0 vanishes. Of course, we are not really in this limit here but there can be a factor of at least 10 between the poloidal wavenumbers such that the effect of the stable stratification becomes very weak on the value of the growth rates.
On Figure 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/Ω 0 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, theoreti-  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.5t ap (during the linear phase of the instability) for cases R3 (P r = 0.1, top) and R4 (P r = 1, bottom). We clearly see that the perturbation direction is orthogonal to gravity, especially for the most stratified case R4. Note that the background flow and field are very similar in both cases.
cally, 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.
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 −1 compared to 5 × 10 −2 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 ∂Ω ∂θ is set to 0. When the Prandtl number is now in-   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.
creased to P r = 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.

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 , 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 Figure 13, both for the cylindrical (left panel) and the radial case (right panel). 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 nonaxisymmetric components.
The conclusion here is similar to the unstratified case  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, i.e. 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 −3 while in the radial case, it must be even weaker, of the order of 10 −3 . Both values are compatible with the values expected to be found in stellar interiors (see Sec.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.

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 P m, the axisymmetric evolution depends only on one dimensionless parameter P rN 2 /Ω 2 0 measuring the level of stratification, instead of the 3 independent parameters L o , P r and N/Ω 0 . 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 P rN 2 /Ω 2 0 ) 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 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.
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, i.e. 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 −3 , 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;Rüdiger et al., 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.