Evolution of a magnetic field in a differentially rotating radiative zone
^{1} Université de Toulouse, UPSOMP, Institut de Recherche en Astrophysique et Planétologie, 31028 Toulouse Cedex 4, France
email: mathieu.gaurat@irap.omp.eu; laurene.jouve@irap.omp.eu; francois.lignieres@irap.omp.eu
^{2} CNRS, Institut de Recherche en Astrophysique et Planétologie, 14 avenue Édouard Belin, 31400 Toulouse, France
^{3} MaxPlanckInstitut für Sonnensystemforschung, JustusvonLiebigWeg 3, 37077 Göttingen, Germany
email: gastine@mps.mpg.de
Received: 18 March 2015
Accepted: 24 June 2015
Context. Recent spectropolarimetric surveys of mainsequence intermediatemass stars have exhibited a dichotomy in the distribution of the observed magnetic field between the kG dipoles of Ap/Bp stars and the subGauss magnetism of Vega and Sirius.
Aims. We would like to test whether this dichotomy is linked to the stability versus instability of largescale magnetic configurations in differentially rotating radiative zones.
Methods. We computed the axisymmetric magnetic field obtained from the evolution of a dipolar field threading a differentially rotating shell. A full parameter study including various density profiles and initial and boundary conditions was performed with a 2D numerical code. We then focused on the ratio between the toroidal and poloidal components of the magnetic field and discuss the stability of the configurations dominated by the toroidal component using local stability criteria and insights from recent 3D numerical simulations.
Results. The numerical results and a simple model show that the ratio between the toroidal and the poloidal magnetic fields is highest after an Alfvén crossing time of the initial poloidal field. For high density contrasts, this ratio converges towards an asymptotic value that can thus be extrapolated to realistic stellar cases. We then consider the stability of the magnetic configurations to nonaxisymmetric perturbations and find that configurations dominated by the toroidal component are likely to be unstable if the shear strength is significantly higher than the poloidal Alfvén frequency. An expression for the critical poloidal field below which magnetic fields are likely to be unstable is found and is compared to the lower bound of Ap/Bp magnetic fields.
Key words: stars: magnetic field / stars: rotation / stars: interiors / magnetohydrodynamics (MHD) / methods: numerical
© ESO, 2015
1. Introduction
In the past decades, the study of the interplay between magnetic fields and differential rotation in radiative zones has mainly been driven by constraints on stellar rotation. This interaction has in particular been invoked to explain the flat rotation profile of the solar radiative zone (see for example Mestel & Weiss 1987; Charbonneau & MacGregor 1993) or more recently, the slow rotation of the core of subgiants and giants revealed by asteroseismology (Deheuvels et al. 2014; Cantiello et al. 2014).
The evolution of the angular momentum in the presence of a magnetic field has mainly been investigated assuming axisymmetry and neglecting meridional circulation and turbulence. This was the case of the numerical simulations by Charbonneau & MacGregor (1992) that considered the spinup of a radiative stellar zone threaded by an initial poloidal field. The authors pointed out the role of the field geometry, more specifically, the existence of field lines anchored in the core, in reaching or failing to reach a solidbody rotation as a final state. They determined the timescale necessary to reach this state. This timescale is controlled by the damping of the Alfvén waves that are excited by the backreaction of the Lorentz force that follows the windingup of the poloidal field by the differential rotation. Charbonneau & MacGregor (1993), Rüdiger & Kitchatinov (1996), and Spada et al. (2010) studied the spindown of the radiative zone of the Sun under the same physical assumptions. The effect of meridional flows has been addressed in Mestel et al. (1988), who restricted their analysis to the axisymmetric case. The impact of nonaxisymmetric field components was then studied by Moss et al. (1990), Moss (1992) and more recently by Wei & Goodman (2015). It was shown that a state of solidbody rotation is reached on a timescale shorter than the global diffusive time, regardless of the existence of field lines anchored in the core. However, if the field strength is small compared to the differential rotation strength, the misaligned magnetic field can be axisymmetrised before solidbody rotation is reached.
These studies did not include the effect of possible nonaxisymmetric instabilities of the magnetic configurations. However, on the one hand, such instabilities are expected if the windingup of the poloidal field by differential rotation produces magnetic configurations that are dominated by the toroidal component. Purely toroidal field configurations can indeed be unstable to various kinds of instabilities, such as the Tayler instability, the magnetorotational instability (MRI), or the buoyancy instability (see for a review Spruit 1999). On the other hand, magnetic fields with toroidal and poloidal components of the same order of magnitude can remain stable, as was found in recent numerical simulations (e.g. Braithwaite & Nordlund 2006; Braithwaite 2007). For the angular momentum transport, the occurrence of such instabilities is a crucial issue since the development of a nonaxisymmetric instability could profoundly modify the redistribution of angular momentum. Estimates of the transport resulting from the Tayler instability have been proposed in Spruit (2002) on phenomenological grounds, while recent numerical simulations by Rüdiger et al. (2015) quantified the transport induced by the socalled azimuthalMRI in the simplified setup of a fluid of constant density confined between two rotating cylinders.
Beyond the problem of angular momentum transport, nonaxisymmetric instabilities could also affect the value of the magnetic field measured at the surface of stars. Notably, spectropolarimetry is an observational technique that allows measuring the lineofsight component of the magnetic field vector integrated over the whole visible surface. This averaged field could be strongly decreased as one goes from a wellorganized large scale magnetic field to a destabilized configuration where opposite polarities cancel each other, especially if the typical length scale left by the instability is much smaller than the stellar radius. This effect has been invoked to explain the magnetic dichotomy observed among intermediatemass stars between the strong mostly dipolar Ap/Bp magnetic fields and the ultraweak magnetic fields of other intermediatemass stars (Aurière et al. 2007; Lignières et al. 2014). In this scenario, the observed lower bound of Ap/Bp magnetic fields (~300 G for the dipolar strength) corresponds to the limit between stable and unstable magnetic fields. On phenomenological grounds, Aurière et al. (2007) approximated the maximum toroidal field produced by the windingup of the poloidal field to (where R is the radius of the star, Ω is the rotation rate and ρ is the density) and assumed that the magnetic configuration is unstable as soon as the toroidal field dominates the poloidal field at the stellar surface. They thus derived an orderofmagnitude estimate of the critical field that separates stable and unstable magnetic fields. This estimate appears to match the observed value for a typical Ap/Bp star.
Our ultimate goal is to investigate the occurrence and the effect of nonaxisymmetric instabilities of the magnetic field induced by the windingup of a poloidal field in a differentially rotating stellar radiative zone. In the present paper, we use axisymmetric simulations to explore the different magnetic configurations obtained for various differential rotation profiles, density stratifications, and boundary conditions. In this systematic parameter study, we focus on the ratio between the toroidal and poloidal magnetic field as a key parameter governing the stability of the magnetic configuration. We discuss the occurrence of instability in fields that are dominated by their toroidal component, for which we use local stability criteria for purely toroidal fields as well as results from recent 3D numerical simulations by (Jouve et al. 2015) performed for particular configurations. The advantage of these 2D axisymmetric numerical simulations over 3D simulations is to give access to a much broader parameter range and in particular to explore asymptotic behaviours in the regimes of low magnetic diffusivities and viscosities and large density stratifications that are unaccessible to 3D models.
In Sect. 2, the mathematical formulation of our problem is described. We set down the basic assumptions that we adopted and then introduce the equations governing the joint evolution of the differential rotation and the magnetic field. In Sect. 3, the results of the numerical simulations are presented. In Sect. 4, simple models that describe the evolution of the toroidal to poloidal field ratio are presented and compared to the numerical results. In Sect. 5, we discuss the stability of the magnetic configurations obtained by the simulations.
2. Mathematical formulation
2.1. Simplifying assumptions and governing equations
We consider the axisymmetric evolution of an initially poloidal field submitted to differential rotation in a spherical shell. Beyond the assumption that the flow remains axisymmetric over time, we also neglect the meridional circulation. The dynamics is therefore not affected by buoyancy effects and is only governed by the induction and the azimuthal momentum equations. The stellar equilibrium structure comes into play through the density stratification. We also consider timescales shorter than the magnetic diffusion time on a radius length scale. Under the assumptions of axisymmetry and without meridional circulation, the initial poloidal field only evolves through Ohmic dissipation and is thus considered as constant in time.
The governing equations are thus where the rotation rate Ω(r,θ,t) is related to the azimuthal velocity by u_{ϕ}(r,θ,t) = rsinθ Ω(r,θ,t) and the magnetic field reads (3)B_{ϕ}(r,θ,t) being its azimuthal component and B_{p}(r,θ), the timeindependent poloidal component. In these equations, the magnetic diffusivity η and the dynamic viscosity μ are uniform, while the density ρ is given by a polytropic model of a star with a polytropic index n = 3 (Rieutord et al. 2005).
2.2. From three to two dimensionless numbers
A possible and straightforward nondimensionalization of these equations is obtained by taking the stellar radius R as the reference length scale, 1/Ω_{0} as the reference timescale (where Ω_{0} is the surface rotation rate at the equator), B_{0} the surface poloidal magnetic field taken at the equator as a magnetic field unit, and ρ_{0} the surface density as a density unit. The dynamical Eqs. (1) and (2) are then controlled by three different dimensionless numbers, namely the Reynolds number Re = Ω_{0}R^{2}/ν_{0}, where ν_{0} = μ/ρ_{0} is the kinematic viscosity at the surface, the Elsasser number , and the magnetic Prandtl number P_{m} = ν_{0}/η. However, it is also possible to reduce the number of dimensionless parameters to two by instead employing the Alfvén timescale as the reference timescale and as a magnetic field unit for the toroidal component. In addition to the reduction of the number of dimensionless parameters allowed by this choice, it is motivated by the fact that B_{ϕ} is generated by the Ωeffect and hence related to Ω_{0} in a differentially rotating radiative zone. B_{0} remains the poloidal magnetic field unit and Ω_{0} the angular rate unit for this nondimensionalization. The ratio of the toroidal to poloidal magnetic fields is then expressed as (4)where t_{Ω} = 1/Ω_{0}, and where tildes indicate dimensionless quantities. With this nondimensionalization, Eqs. (1) and (2) read L_{u}, the Lundquist number, is defined as follows: (7)This reduction from three to two dimensionless numbers greatly facilitates the exploration of the parameter space.
2.3. Initial and boundary conditions
The initial magnetic field is a dipole (see white lines in Fig. 1): As stated above, the poloidal component is timeindependent, while the toroidal component will evolve in time. The boundary conditions are a perfect conductor at the bottom boundary and an insulating medium at the outer radius. This translates into (10)where R_{c} is the radius of the inner boundary. We here explore different initial and boundary conditions for the rotation. Three different radial differential rotations and two cylindrical profiles were considered:
where ϖ = rsinθ is the distance from the rotation axis. Two types of boundary conditions are used at r = R_{c}, namely a stressfree or a fixed rotation, while a stressfree boundary condition is always used at r = R. These different cases are summarized in Tables 3 and 4 and illustrated in Fig. 1 which displays a meridional cut of a radial (left frame) and a cylindrical (right frame) differential rotation profile.
From our reference density profile corresponding to an n = 3 polytrope, we create different profiles with ratios ρ_{c}/ρ_{0} ranging from 1 to 10^{6} (where ρ_{c} is the density at the radius r = R_{c}) by cutting the polytropic profile at different radii. The created profiles are then expanded on the whole computational domain. Figure 2 shows the radial dependence of the dimensionless Alfvén velocity / for the different density profiles used in our simulations. The corresponding time for an Alfvén wave to propagate from r = R_{c} to r = R along the equator would vary from 0.25 t_{Ap} for ρ_{c}/ρ_{0} = 1 to 22 t_{Ap} for ρ_{c}/ρ_{0} = 10^{6}.
Fig. 1 Meridional cut of for the radial profile defined by Eq. (12) (left) and for the cylindrical profile defined by Eq. (15) (right). The dotted black lines represent the isocontours of . The white lines represent the poloidal magnetic field lines and the red line the one on which Alfvén waves propagate the slowest. The radius at the inner boundary is r = R_{c} = 0.3 R. 

Open with DEXTER 
Fig. 2 at the equator as a function of for different density contrasts ρ_{c}/ρ_{0}. A n = 3 polytropic profile was cut at different radii to produce these various ρ_{c}/ρ_{0}. 

Open with DEXTER 
2.4. Numerical model
To solve Eqs. (5) and (6) with the boundary and initial conditions, we used the twodimensional axisymmetric code STELEM (Charbonneau & MacGregor 1992; Jouve et al. 2008). This code uses a finiteelement method in space and a thirdorder RungeKutta scheme in time. The computational domain is limited to the annulus r ∈ [R_{c},R], θ ∈ [0,π], where R_{c} = 0.3 R was chosen for all our simulations. P_{m} was fixed to 1 in all our simulations while the Lundquist number L_{u} is varied from 10 to 10^{5}. Depending on the value of L_{u} and on the boundary and initial conditions, different spatial and temporal resolutions are required. We used a spatial mesh ranging from N_{r} × N_{θ} = 128 × 128 to 512 × 512 nodes depending on the cases. The integration time was adjusted to satisfy the CourantFriedrichsLewy criterion. It varies from dt = 10^{3}t_{Ap} to dt = 10^{8}t_{Ap}.
With the initial conditions used in this study, our problem has an equatorial symmetry. This is clearly visible in the meridional cut of Fig. 1 for example. However, for comparison to future studies with more complex symmetry, we chose to include both hemispheres in our computations and always represent the full meridional plane.
We note that we do not present any results for L_{u}> 245 with a stressfree boundary condition at r = R_{c}. In fact, with this boundary condition there is a Hartmann boundary layer at the inner radius. The thickness δ of this boundary layer is known to be (Dormy et al. 1998) (16)where n is the unit vector normal to the boundary. The number of mesh points is increased according to Eq. (16) to ensure that the Hartmann layer is much larger than the size of the numerical grid. However, for large L_{u}, this imposes prohibitive restrictions on the time step.
Fig. 3 The reference case: meridional cut of , and for different times. t_{max}, the time at which B_{ϕ}/B_{p} is maximal, is equal to 1.85 t_{Ap} in this simulation. The colourbar at the bottom of each panel indicates the minimal and maximal values of each variable. The parameters of this simulation are L_{u} = 10^{2}, a density contrast of ρ_{c}/ρ_{0} = 7 × 10^{3}, an initial radial differential rotation given by Eq. (12) and a fixed rotation at the inner radius. 

Open with DEXTER 
3. Description of the results
The evolution of the magnetic field and the rotation rate is first qualitatively described pointing out the successive phases and the associated physical mechanisms. We then focus on the ratio between the toroidal and the poloidal magnetic field as a key parameter in the study of the stability of a magnetic field. The influence of the Lundquist number, the density contrast, the initial and the boundary conditions on the highest value reached by this ratio is investigated in detail.
3.1. Evolution of the magnetic field and the differential rotation in a typical case
Figure 3 shows the evolution of the toroidal magnetic field B_{ϕ}, the ratio between B_{ϕ} and the poloidal magnetic field B_{p}, and the rotation rate Ω for a typical case (hereafter called the reference case). The parameters of this simulation are L_{u} = 10^{2}, a density contrast ρ_{c}/ρ_{0} = 7 × 10^{3}, an initial radial differential rotation given by Eq. (12) and a fixed rotation at the inner radius. Different phases in the evolution of B_{ϕ} and Ω can be distinguished. The first is the windingup phase, where the shearing of the poloidal field by the differential rotation, the socalled Ωeffect, generates B_{ϕ} without any backreaction on the differential rotation (panel a). Since B_{ϕ} is generated by the term B_{p}·∇Ω, the fact that B_{ϕ} is antisymmetric with respect to the equator directly relates to the properties of symmetry of Ω and B_{p}. After this phase, the Lorentz force backreacts on the differential rotation, leading to the propagation of Alfvén waves along the poloidal magnetic field lines in both directions. In this simulation, we mainly observe an outward propagation from the internal radius because the initial perturbation of the system is concentrated close to the bottom boundary.
In panel d, a region of opposite sign is visible on B_{ϕ} close to the core. It is due to the reflection of Alfvén waves on the equator, where conditions on B_{ϕ} and Ω allow for such reflections. Indeed, the ability for a boundary to reflect Alfvén waves depends on the characteristics of B_{ϕ} and Ω at this boundary. In our simulations, the equator behaves as a boundary where B_{ϕ} = 0 (insulating condition) and ∂Ω /∂θ = 0 (stressfree condition). This pair of boundary conditions induces a perfect reflection of Alfvén waves at the equator (Schaeffer et al. 2012). The subsequent evolution is characterized by the propagation and the reflection of the waves either at the equator, at the surface, or at the bottom boundary.
The variable Alfvén velocity and the different distances that the Alfvén waves have to travel before reflexion lead to a phase shift between waves on neighbouring magnetic field lines. As we can see in panel e, strong gradients of B_{ϕ} and Ω are then created between these waves. The effect of magnetic and viscous diffusion is thus increased. This is the phasemixing mechanism as described in Ionson (1978), Heyvaerts & Priest (1983), or Parker (1991). The oscillations of B_{ϕ} then decrease in amplitude, as shown in panel f. After this dissipative phase controlled by the phasemixing mechanism, the system evolves towards a steady state known as Ferraro’s law of isorotation (Ferraro 1937), in which Ω is constant along poloidal field lines. In our case, the only possible steady state compatible with the boundary conditions and Ferraro’s law is uniform rotation. This state of uniform rotation is reached on timescales of the order of a few diffusive timescales which are much longer than the timescales we are interested in here. It is thus not shown in Fig. 3.
3.2. Evolution of the ratio between the toroidal and the poloidal magnetic field
Since nonaxisymmetric instabilities are expected to occur when the magnetic configuration is dominated by the toroidal component, we now focus on the highest value of B_{ϕ}/B_{p} achieved during the windingup process. In the cases we considered, this maximum is always reached before the phasemixing episode takes place (panel c for the simulation shown in Fig. 3). Figure 4 shows the evolution of the spatial maximal value of B_{ϕ}/B_{p} for the reference case. This quantity is defined as the highest value reached in the whole domain. At first, we observe a nearly linear growth of the spatial maximal value of B_{ϕ}/B_{p}. As can be seen in panels a−c of Fig. 3, during this phase the location of the maximum moves as an Alfvén wave along a poloidal field line. The B_{ϕ}/B_{p} ratio reaches its highest value in time and space, max(B_{ϕ}/B_{p}), at a time denoted t_{max}. At later times in Fig. 4, the spatial maximal value of B_{ϕ}/B_{p} oscillates due to reflections of Alfvén waves at the boundaries and rapidly decreases towards 0 as the system evolves towards a steady state through the dissipation induced by the phasemixing mechanism. The temporal evolution presented in Fig. 4 is typical for most of our simulations. We note that (r_{max},θ_{max}), the location at which max(B_{ϕ}/B_{p}) is reached, is visible in panel c of Fig. 3 at midradius and low latitude. In this zone, the poloidal magnetic field lines do not reach the external radius, they are closed within the shell (see Fig. 1).
To illustrate the physical mechanisms that give rise to max(B_{ϕ}/B_{p}), one can look at the evolution of B_{ϕ}/B_{p} and Ω along the poloidal field line on which max(B_{ϕ}/B_{p}) is reached. Figure 5 displays this evolution for the parameters of the reference case, except that L_{u} = 10^{4} instead of 10^{2}. Top panels show the growth and decrease of B_{ϕ}/B_{p} between 0 and 2 t_{max}, max(B_{ϕ}/B_{p}) being reached at time t_{max} and location s_{max}. The bottom panels show that Ω behaves like a standing wave rather than a traveling wave. This is due to the boundary conditions used in the simulation and the spatial scale of the initial differential rotation, which is comparable to the length of the field line. B_{ϕ} then also behaves as a standing wave, but this is not clearly visible in the top panels where B_{ϕ}/B_{p} is represented. Indeed, since B_{p} is not uniform along the field line, the standing wave character of B_{ϕ} is not obvious in B_{ϕ}/B_{p}. We also note that at approximately t = t_{max}, the gradient of Ω projected onto the poloidal field line, ∂Ω /∂s = e_{s}·∇Ω where s is the curvilinear coordinate on the field line, changes sign. According to the induction Eq. (5) in which the effects of magnetic diffusion are neglected (L_{u} ≫ 1), ∂(B_{ϕ}/B_{p}) /∂t changes sign together with ∂Ω /∂s, and we indeed observe in Fig. 5 the decrease of B_{ϕ}/B_{p} after t_{max}.
Fig. 4 Evolution of the spatial maximal value of for the reference case. The horizontal dashed line indicates and the vertical dashed line indicates . 

Open with DEXTER 
3.3. Summary of the parameter study
The quantity max(B_{ϕ}/B_{p}) defined in the previous section determines whether the magnetic configuration is locally dominated by the toroidal or the poloidal field, and t_{max} gives the time at which max(B_{ϕ}/B_{p}) is reached. Tables 1−4 summarize the values of and obtained by varying the Lundquist number L_{u}, the density contrast ρ_{c}/ρ_{0}, the initial differential rotation profile, and the boundary condition for Ω at r = R_{c}, respectively. In the following, these results are discussed, starting with the effect of L_{u}. The two last columns show the predictions of models that are presented in the next section.
3.4. Influence of the diffusivities
Figure 6 displays max(B_{ϕ}/B_{p}) and t_{max} as a function of L_{u} for the data listed in Table 1. Both max(B_{ϕ}/B_{p}) and t_{max} increase with L_{u}. We expect max(B_{ϕ}/B_{p}) to increase when the effects of diffusion on the toroidal magnetic field decrease and thus when L_{u} increases. To explain why t_{max} increases with L_{u}, we note that the diffusion term in the induction equation is expected to be negative (positive) near local maxima (minima) of B_{ϕ}. Thus, diffusion contributes to a faster change of sign of ∂B_{ϕ}/∂t and hence to a smaller t_{max}.
More interestingly, we see in this figure that an asymptotic value is reached for L_{u} ≳ 10^{4}. The same asymptotic behaviour with L_{u} is observed for all cases considered. The fact that the magnetic and viscous diffusions do not affect the maximum of B_{ϕ}/B_{p} above a certain L_{u} is not surprising since this maximum is reached on a timescale of the order of an Alfvén time, this timescale becoming much shorter than the diffusive timescale as L_{u} increases. This asymptotic behaviour constitutes an interesting feature when applications to realistic stellar radiative zones are considered.
Fig. 5 (top panels) and (bottom panels) as a function of s/L for different times. s is the curvilinear coordinate of the poloidal field line in which max(B_{ϕ}/B_{p}) is reached, and L the length of this field line. s = 0 locates the inner radius and s = L the equator. The simulation is the same as in Fig. 3 with the exception of L_{u} = 10^{4}. The vertical dashed line indicates s_{max} the location of max(B_{ϕ}/B_{p}) on the poloidal field line. 

Open with DEXTER 
Influence of the density contrast ρ_{C}/ρ_{0} on and .
3.5. Influence of the density profile
The dependence on the density contrast was studied numerically by varying ρ_{c}/ρ_{0} between 1 and 10^{6} (see Table 2). For a given configuration (L_{u} = 10^{5}, an initial differential rotation given by Eq. (12) and Ω fixed at the inner radius), the two top frames of Fig. 7 show that cosθ_{max} and respectively increases and decreases towards an asymptotic value as ρ_{c}/ρ_{0} increases. Thus, the location of the maximum B_{ϕ}/B_{p} remains unchanged above a certain ρ_{c}/ρ_{0}. As B_{ϕ} is produced through the propagation of Aflvén waves along field lines, the field line where max(B_{ϕ}/B_{p}) is reached is the one that goes through θ_{max} and r_{max}. This field line will remain the same above a certain ρ_{c}/ρ_{0} ratio. The same asymptotic behaviour with ρ_{c}/ρ_{0} is observed for the other simulations not listed in Table 2, performed for different L_{u} and initial differential rotations. For all the cases considered, the asymptotic value of is comprised between 0.6 and 0.8, while cosθ_{max} is comprised between 0.3 and 0.4. Hence, max(B_{ϕ}/B_{p}) is always reached in a zone located between the middle and the top of the radiative zone, close to the equator.
Influence of the initial differential rotation on and .
Influence of the boundary condition of Ω at the inner radius on and .
Table 2 shows that and monotonically increase with ρ_{c}/ρ_{0}. If we consider and , we find that an asymptotic limit is reached for large ρ_{c}/ρ_{0}. This is shown in the two bottom frames of Fig. 7. In other words, and as ρ_{c}/ρ_{0} tends towards realistic stellar values, where C1 and C2 are independent of the density contrast. For the different simulations performed, the value of C1 varies between 3 × 10^{3} and 2 × 10^{1} and C2 between 10^{2} and 3 × 10^{2}. The variations of these constants are mainly due to the initial differential rotation considered and, to a smaller extent, to the boundary condition adopted at the inner radius. This asymptotic behaviour is expected if the physical mechanism that gives rise to max(B_{ϕ}/B_{p}) is confined in the stellar interior since above a certain ρ_{c}/ρ_{0} ratio the density profile is no longer modified in the internal layers. Indeed, max(B_{ϕ}/B_{p}) is reached along a field line that remains confined in these internal layers.
3.6. Influence of the initial and boundary conditions
The initial differential rotation can have a strong effect on max(B_{ϕ}/B_{p}), but it only slightly modifies t_{max}. Indeed, as observed in Table 3, the value of max(B_{ϕ}/B_{p}) varies from 18.5 t_{Ap}/t_{Ω} for the 1 /r^{2} profile (Eq. (11)), which is the strongest shear we can use that is hydrodynamically stable, to 0.31t_{Ap}/t_{Ω} for the cylindrical profile given by Eq. (13). Generally, the stronger the shear, the larger max(B_{ϕ}/B_{p}). But we also find that the cylindrical differential rotation induces less B_{ϕ} than the radial one because the isocontours of Ω are more aligned with the poloidal magnetic field lines in the cylindrical case, causing the Ωeffect to be less efficient. The value of t_{max} is approximately independent of the radial differential rotation profile. This indicates that the time to reverse the shear does not seem to depend on its intensity. However, t_{max} varies from 0.96 t_{Ap} to 1.35 t_{Ap} for the two cylindrical profiles considered.
Our variables of interest max(B_{ϕ}/B_{p}) and t_{max} are less influenced by the boundary conditions on Ω than by all the other parameters considered in this study. As illustrated in Table 4, for radial differential rotation profiles, max(B_{ϕ}/B_{p}) and t_{max} are always slightly larger for a fixed value than for a stressfree boundary condition, 1.01 to 1.66 and 1.03 to 1.37 times higher, respectively, depending on the cases. For the cylindrical profiles, max(B_{ϕ}/B_{p}) and t_{max} reach the same values, independently of the boundary conditions.
The boundary conditions other than the one on Ω at r = R_{c} only have a limited impact on max(B_{ϕ}/B_{p}) and t_{max}. We thus do not show their influence here. If longer timescales were considered (and final steady states), they would play a significant role, however.
Fig. 6 max(B_{ϕ}/B_{p}) (top) and t_{max} (bottom) as a function of L_{u} with the physical ingredients used for the reference case. 

Open with DEXTER 
Fig. 7 cosθ_{max}, (top frames), and (bottom frames) as a function of ρ_{c}/ρ_{0}. The parameters of the simulations are L_{u} = 10^{5}, a profile of Ω defined by Eq. (12) and Ω fixed at the inner radius. 

Open with DEXTER 
4. Simple models
In this section two models are presented that provide an estimate of max(B_{ϕ}/B_{p}) and t_{max}. Their predictions are compared with the results of the simulations. The first model is a local model that uses the same phenomenological assumptions as Aurière et al. (2007). The second model consists of estimating the poloidal field line on which max(B_{ϕ}/B_{p}) is reached and then of solving a 1D Alfvén wave equation along this field line.
4.1. Local model of Aurière et al. (2007)
Starting from the ideal induction equation for the toroidal component (17)Aurière et al. (2007) assumed that the toroidal field backreacts on the differential rotation after a time τ_{max} = ℓ/v_{Ap}, corresponding to the period of an Alfvén wave of wavelength ℓ =  ∇Ω  /Ω, the length scale of the initial differential rotation. Then, by integrating Eq. (17) from t = 0 to t = τ_{max} and assuming that the differential rotation is constant up to that time, we obtain (18)Taking the spatial maximum of B_{ϕ}/B_{p} in Eq. (18), max(B_{ϕ}/B_{p}),r_{max},θ_{max} and then τ_{max} can be determined.
Of the assumptions of this model, neglecting magnetic diffusion appears justified since the numerical results have shown that max(B_{ϕ}/B_{p}) becomes independent of L_{u} when L_{u} ≳ 10^{4}−10^{5}. Relating τ_{max} to the length scale of the initial differential rotation ℓ leads to a simple expression, but then B_{ϕ}/B_{p} taken at τ_{max} does not depend on the initial gradient of Ω and in particular on its projection onto poloidal field lines. For example, if Ferraro’s law is verified initially, B_{ϕ} should remain equal to zero, but the model incorrectly predicts that induction takes place. Another problem is that with a local model, the effects of the boundary conditions are not taken into account.
4.2. Model of standing Alfvén waves on a poloidal magnetic field line
Since the initial perturbations of Ω and B_{ϕ} propagate as Alfvén waves along B_{p}, a more accurate model would be to solve a onedimensional wave equation along the poloidal field line where max(B_{ϕ}/B_{p}) is reached. As we expect this field line to be the one for which the Ωeffect lasts the longest, we assume that max(B_{ϕ}/B_{p}) is reached on the field line that maximizes the Alfvén propagating time over the line length L (taken from the inner radius to the equator or the surface) . The field line obtained using this condition for a case with ρ_{c}/ρ_{0} = 7 × 10^{3} is indicated in red in Fig. (1).
Neglecting the magnetic and viscous diffusions and assuming that the scale of variation of rsinθ and B_{p} along B_{p} is larger than the scale of ∂Ω /∂s (Mestel & Weiss 1987), the 1D wave equation to be solved on this field line is (19)As already noticed, Ω behaves approximately as a standing wave for the initial differential rotation considered here (see Fig. 5). We thus search for normalmode solutions f(s)exp(iωt) of the wave equation. In addition, we consider the fundamental mode as an approximate solution of the problem because the initial perturbation is dominated by its largescale component. The equation for f(s) is (20)A simple and explicit solution is obtained assuming that v_{Ap} is a constant equal to the mean Alfvén velocity along the poloidal field line. We then obtain B_{ϕ}/B_{p}(s,t = t_{max}) by integrating the induction equation up to t_{max}(21)for a fixed boundary condition at s = 0 and a stressfree boundary condition at s = L, with and where ΔΩ = Ω(s = 0,t = 0)−Ω(s = L,t = 0) and g(s) = r(s)sinθ(s). (22)for a stressfree boundary condition at s = 0 and s = L and with .
A more accurate solution can be obtained using the WentzelKramersBrillouin (WKB) approximation (23)for a fixed boundary condition at s = 0 and a stressfree boundary condition at s = L, with t_{max} = τ(L) and where . (24)for a stressfree boundary condition at s = 0 and s = L and with .
We now compare this WKB solution and the model of Aurière et al. (2007), hereafter called SW model for “standing wave” and Aurière model, with the results of the numerical simulations.
4.3. Comparison of the models with the numerical simulations
Fig. 8 as a function of s/L for different times. s is the curvilinear coordinate of the poloidal field line on which max(B_{ϕ}/B_{p}) is reached, and L the length of this field line. s = 0 locates the inner radius and s = L the equator. The parameters of the simulation are L_{u} = 245, ρ_{c}/ρ_{0} = 7 × 10^{3}, the differential rotation profile given by Eq. (12) and a stressfree boundary condition. 

Open with DEXTER 
The values of max(B_{ϕ}/B_{p}) and t_{max} given by the models are compared with the results of the simulations in Tables 1−4.
First, we observe that the Aurière model and the SW model reproduce the asymptotic behaviour of max(B_{ϕ}/B_{p}) when ρ_{c}/ρ_{0} increases. Indeed, Table 2 indicates that max(B_{ϕ}/B_{p}) is given with a mean relative error of 13% by the Aurière model against 24% for the SW model. In this particular case, the Aurière model agrees better with the simulations than the SW model. However, in the other simulations that are not listed in this paper it is generally the opposite trend. The SW model provides very accurate predictions of t_{max} with a mean relative error of only 4%, while the Aurière model predicts values of t_{max} several orders of magnitude higher than the simulations (we have thus chosen not to list them in the tables).
As expected, the effect of the initial differential rotation is much better reproduced by the SW model than by the Aurière model, with a mean relative error on max(B_{ϕ}/B_{p}) of 26% and 280%, respectively (see Table 3). While the Aurière model does not take into account the effects of the boundary conditions, the SW model estimates max(B_{ϕ}/B_{p}) with a mean relative error of only 16% for the different boundary conditions considered (see Table 4). Finally, Fig. 8, for the stressfree BC case, and Fig. 5, for the fixed BC, indeed show that the evolution of the perturbation is not far from the standing waves prescribed by the SW model.
We conclude that the SW model provides a useful approximation for the maximum ratio B_{ϕ}/B_{p} in a differentially rotating star with a dipolar field, provided that the length scale of the initial differential rotation is not too small compared to the stellar radius.
5. Towards the instability of the magnetic field
We now turn to the discussion of possible instabilities of the magnetic configurations found in our simulations. Figure 9 shows the configuration obtained at time t_{max} in a typical simulation. A 3D rendering of the field lines is shown, coloured by the total magnetic intensity. We clearly see in this figure that we have a complex structure with mixed poloidal and toroidal components, even if the toroidal field dominates at midlatitudes. Moreover, differential rotation is still present in the spherical domain at this time of the simulation and may strongly influence the stability conditions of this complex magnetic field.
One way to study the stability of such a complex magnetic configuration without approximation is to perform 3D numerical simulations where the full set of MHD equations is solved. This approach has been followed in a companion paper (Jouve et al. 2015), although the magnetic configurations studied in 3D have been limited to cases of uniform density, cylindrical differential rotation, and low L_{u} values. Jouve et al. (2015) found that a magnetic instability is triggered and destroys the largescale magnetic field if the ratio of the poloidal Alfvén time t_{Ap} to the rotation timescale t_{Ω} is sufficiently high. In unstable cases, an enhanced transport of angular momentum due to the turbulence induced by the instability is found, confirming the results of Rüdiger et al. (2015) in cylindrical geometry. This fully 3D study also served to test the capacity of approximate local criteria to determine the stability of these complex configurations. In particular, it was found that the dispersion relation of Ogilvie (2007) could predict the nature of instability and provide reasonable estimates of its growth rates. Here, we used these local criteria to predict the stability of the configuration obtained through our axisymmetric simulations. As compared to Jouve et al. (2015), we therefore study the stability of a wider range of configurations, including the effects of the density stratification, of a radial differential rotation profile, and of stable stratification.
5.1. Local stability analysis under simplifying assumptions
Fig. 9 3D view of the magnetic configuration obtained with our 2D simulations at t_{max}. It is computed for L_{u} = 10^{2}, a density contrast ρ_{c}/ρ_{0} = 7 × 10^{3}, the radial differential rotation profile defined by Eq. (12) and Ω fixed at the inner radius. 

Open with DEXTER 
From numerous previous studies, we know that a magnetic field is more likely to be unstable if it is dominated by one of its components, namely when it is either purely poloidal (Markey & Tayler 1973; Flowers & Ruderman 1977) or purely toroidal (Tayler 1973; Pitts & Tayler 1985). In the previous sections, we have established in which conditions the initial differential rotation produces a configuration dominated by the toroidal field. For the stability analysis, we only consider the situations where the toroidal magnetic field is so dominant that the poloidal component can be neglected. We return to this assumption in the next subsection.
A major difficulty in performing a linear stability analysis of our magnetic configurations is that the axisymmetric background field evolves in time. Therefore, an important quantity to consider here is the ratio between the growth time of possible instabilities and the evolution timescale of the background magnetic field. From now on, we assume that this ratio is low and thus that the background state can be assumed to be steady. This assumption is discussed in the next subsection.
The local magnitude of the differential rotation is measured by qΩ, where q = ∂lnΩ /∂lnϖ is the shear parameter and ϖ is the cylindrical radius. For the different profiles considered (see Eqs. (11) to (15)), the shear parameter is always of order unity, meaning that the differential rotation is of the order of the local rotation rate.
Under these assumptions, we may proceed to a local linear stability analysis, using the dispersion relation derived by Ogilvie (2007) for the case of a purely toroidal field under the influence of a differential rotation, both possessing arbitrary profiles. Two instabilities may be expected in this situation: the Tayler instability (TI), which derives from free magnetic energy, and the MRI, whose source of energy is the free kinetic energy of the differential rotation. The steady and axisymmetric basic state consists of a differential rotation profile Ω(ϖ,z) and a purely toroidal magnetic field B_{ϕ}(ϖ,z), where z is the coordinate along the rotation axis. If we consider a perturbation such that the displacement is in the direction e, we obtain the following form of the dispersion relation: (25)where ω is the frequency of the perturbation, m is the azimuthal wavenumber, v_{Aϕ} = B_{ϕ}/ is the toroidal Alfvén velocity, and N is the BruntVäisälä frequency taken as uniform in the whole domain. As noted previously, the buoyancy force does not affect the evolution of the axisymmetric magnetic field since the meridional flows are neglected. It influences the stability of the magnetic configurations through the value of the N/ Ω ratio, however. As noted by Jouve et al. (2015), another important quantity of the magnetic stability is the ratio of the rotation rate Ω to the toroidal Alfvén frequency ω_{Aϕ} = v_{Aϕ}/r.
Fig. 10 Contours of the toroidal magnetic field (white lines) and of the growth rate of the instability (coloured contours) for the most unstable perturbation (m = 7 for the top frame and m = 1 for the bottom frame) which is in the direction e = e_{ϖ}. The parameters of the simulations are N = 0 and, respectively for the top and the bottom frames Ω /ω_{Aϕ} ≈ 10 and Ω /ω_{Aϕ} ≈ 1 at the maximum of the toroidal field. 

Open with DEXTER 
For a given configuration B_{ϕ}(ϖ,z,t_{max}) and Ω(ϖ,z,t_{max}), corresponding to a Ω(ϖ,z,t_{max}) /ω_{Aϕ}(ϖ,z,t_{max}) ratio, obtained with the 2D numerical simulations and a value of N/ Ω, the dispersion relation provides the growth rate σ = ℑ(ω(ϖ,z)) of a local perturbation characterized by m and e. Two cases are illustrated in Fig. 10. In this figure, only the northern hemisphere is shown since the results are symmetric with respect to the equator. The contours of the toroidal magnetic field are drawn (white lines) together with the contours of the growth rate of the instability (coloured contours) for the most unstable perturbation, which is in the direction e = e_{ϖ}. The only difference in the initial setup of these two simulations is the differential rotation profile (Eq. (13) was used for the case in the top panel and Eq. (12) for the bottom panel). As a consequence of a different windingup between these two cases, the ratio Ω /ω_{Aϕ} reaches very different values. As shown in Jouve et al. (2015), the nature of the instability that is likely to be triggered in these magnetic configurations depends on this Ω /ω_{Aϕ} frequency ratio. For the case in the top panel Ω /ω_{Aϕ} ≈ 10 at the maximum of toroidal field, the differential rotation thus plays a dominant role, and the most vigorous instability is the MRI, the most unstable mode having an azimuthal wave number m = 7. In contrast, for the other case where Ω /ω_{Aϕ} ≈ 1 at the maximum of toroidal field, the currentdriven Tayler instability is triggered and favours the m = 1 mode. The location of the instability is also quite different in both cases. In the MRI case, the most unstable region is concentrated at rather low latitudes and very extended in radius. In the TI case, however, the unstable zones are mainly located on strong gradients of the field close to the bottom of the domain at latitudes 30° and 60°. In all cases studied, the value of Ω /ω_{Aϕ} is found to vary between about 1 and 10 depending on the initial and boundary conditions. We thus expect our axisymmetric configurations to be subject to an instability, the nature of which only depends on this Ω /ω_{Aϕ} ratio. Both the MRI and the TI are found to be likely to exist in our situations. The location of the instability, the growth rates as well as the typical lengthscale of the most unstable mode will then differ and the consequences on the observed surface field may be very different between the TI and the MRI cases.
We showed in Sect. 3.5 that when the density contrast is increased between the top and bottom of our domain, both the highest value of B_{ϕ}/B_{p} and its location tend to an asymptotic value. For the stability conditions, we find that for the case where the m = 1 Tayler instability is favoured (bottom panel of Fig. 10), increasing the density contrast does not significantly vary the nature, location, and growth rate of the instability. Indeed, as ρ_{c}/ρ_{0} is increased from 7 × 10^{3} (case shown in the figure) to 2 × 10^{4} and to 10^{6}, the whole magnetic configuration is left mostly unchanged and the value of Ω /ω_{Aϕ} increases from 1 to 1.25 and to 1.54. In these three cases, the m = 1 mode is found to be the most unstable to the TI and is always located at the base of the domain in two distinct regions in latitude (i.e. at 30° and 60° where the gradients are the strongest). The highest growth rate is always around σ ≈ 7−9 Ω_{Aϕ}, where and B_{ϕm} is the maximal toroidal field. This is an interesting feature of this work since it implies that our results can be extrapolated to realistic values of the density contrast not only for the location and value of max(B_{ϕ}/B_{p}), but also for the instabilities that are likely to be triggered in the stellar interior.
The effect of the stratification has also been explored by varying the parameter N/ Ω_{0}, where Ω_{0} is the rotation rate at the surface. For the case where the Tayler instability is favoured, increasing N/ Ω_{0} from 0 (case shown in the bottom panel of Fig. 10) to 8 and to 16 decreases the maximum growth rate of the instability from 9.2 Ω_{Aϕ} to 8.2 Ω_{Aϕ} and to 5.2 Ω_{Aϕ}. The location of the instability remains unchanged. For the case where the MRI is favoured, the highest growth rate strongly decreases, the most favoured wave number m decreases when N/ Ω_{0} increases, and the instability totally disappears for N/ Ω_{0} ≥ 2. This agrees with Spruit (1999), who stated that the MRI is more affected by the stable stratification than the TI. It also agrees with Masada et al. (2006), who showed that the most unstable azimuthal wave number m decreases when N/ Ω_{0} increases. Nevertheless, with a realistic stable stratification (N/ Ω_{0} ≃ 10), the most unstable perturbations are no longer in the direction e = e_{ϖ}, but rather in the direction e = e_{θ}. Indeed, such horizontal perturbations are not affected by the stratification, and magnetic instabilities are still present in these cases. For the same case as the one shown in the bottom panel of Fig. 10, the m = 1 TI still possesses a significant growth rate σ ≃ 9 Ω_{Aϕ} in the direction e = e_{θ}. The MRI can still be triggered for e = e_{θ}, but with a growth rate σ ≃ 0.4 Ω_{Aϕ} much lower than in the unstratified cases, where σ ≃ 2.6 Ω_{Aϕ} (case shown in the top panel of Fig. 10). We thus still expect nonaxisymmetric instabilities to develop in the stably stratified cases, the main difference being that when Ω /ω_{Aϕ} > 1 (MRI regime), the growth rate is significantly reduced.
5.2. Discussion of the stability condition
One of the simplifying assumptions in the previous section was to consider that the instability would grow on a steady axisymmetric background state. However, both the background magnetic field and the differential rotation always oscillate on the short timescales considered in our study.
From Fig. 4, this oscillation frequency can be estimated to be 1/2 t_{max}. To freely develop, an instability with a growth rate σ should therefore fulfil the condition (26)As discussed in the previous subsection, if Ω /ω_{Aϕ} > 1, the MRI is favoured. Its maximum growth rate is related to the shear parameter and the rotation rate in the following way (see for example Balbus & Hawley 1991): (27)If Ω /ω_{Aϕ} ≲ 1, the TI is favoured and the growth rate reads (28)From the standing Alfvén wave model presented in Sect. 4.2, the toroidal Alfvén frequency ω_{Aϕ} and t_{max} are given by Eqs. (23) and (24), or by Eqs. (21) and (22) in a simplified case. The quantity ΔΩ, which appears in these equations, can be expressed in terms of the rotation rate Ω and the shear parameter q: ΔΩ ≈ ϖ∂Ω /∂ϖ = q Ω. The condition for the development of a nonaxisymmetric instability, Eq. (26), can then be rewritten as (29)which is valid for both the MRI and the TI cases.
In practice, effects ignored in this estimation will make this condition more restrictive. First, as shown in the previous subsection, the stable stratification is expected to significantly decrease the instability growth rate in the MRI case. Second, to derive (29), we assumed that the toroidal component dominates the poloidal component. This assumption does not impose an additional constraint since, according to the SW model, condition (29) also ensures that the toroidal component dominates. Nevertheless, even in the presence of a small poloidal component, a more conservative instability condition might be necessary (see for the TI case Linton et al. 1996; Braithwaite 2009). Finally, the 3D simulations by Jouve et al. (2015) indicate that a condition is more appropriate than to account for the effect of the oscillating background state. Indeed, in the regime perturbations do grow exponentially, but are killed by the background field reversal before they reach the level of energy of the axisymmetric field. Taking these three effects into account, the condition for the nonaxisymmetric instabilities therefore is rather . In the particular case of the 3D simulations by Jouve et al. (2015), the condition is .
5.3. Discussion of the critical surface field
When written in terms of quantities observables at the stellar surface, the condition for instability reads (30)where, as we just discussed, the factor C takes into account effects on the instability ignored in the simple condition (29). In the case studied by 3D numerical simulations (where ρ is uniform, ΔΩ/Ω_{0} ~ 1 and the initial field is dipolar), it is found to be ~10^{2}. This expression defines a surface critical field, B_{crit}, below which nonaxisymmetric instabilities can be triggered by differential rotation. This field is proportional to and depends on internal properties, such as differential rotation and field geometry, through the ratio ΔΩ/Ω_{0} and . The actual value of ΔΩ/Ω_{0} will depend on the mechanism that enforces differential rotation in the star and might be much lower than the ΔΩ/Ω_{0} ~ 10^{1}−1 values we considered here as initial conditions.
Following Aurière et al. (2007), the magnetic dichotomy of intermediatemass stars would be due to the development of nonaxisymmetric instabilities separating stable strong field configurations, the Ap/Bp stars, from unstable weaker field configurations whose surface average field becomes very weak after the destabilization. Accordingly, the lower limit of Ap/Bp magnetic fields would be given by a stability condition such as (30) (the present stability condition is global and thus differs from the local condition found by Aurière et al. (2007)). As we know that the observed lower bound of Ap/Bp dipolar fields is ≃300 G and that for a typical Ap star (R = 3 R_{⊙}, P = 5 days, log g = 4 dex and T_{eff} = 10^{4} K), the condition for instability (30) would match observations if . A density profile calculated for such a star thanks to the stellar evolutionary code MESA yields , which implies that ΔΩ/Ω_{0} ≳ 5 × 10^{3}. However, in making this estimate, we assumed that any nonaxisymmetric instability taking place deep inside the star strongly affects the observed surface field. This is not necessarily true, in which case higher ΔΩ/Ω_{0} would be compatible with the Aurière scenario.
Given the present uncertainties on the stability condition (see previous subsection), the level of the differential rotation and the effect of an internal instability on the surface field, the quantitative comparison between the theoretical critical field and the observed lower bound of Ap/Bp fields is still rather approximative. More realistic 3D simulations will help, although the density contrasts encountered in a stellar envelope remain beyond reach for such simulations. It will be more relevant to compare the dependence of the critical field on the surface rotation rate with observations. Existing data indeed point towards a linear dependency of the lower bound of Ap magnetic fields (Lignières et al. 2014) compatible with the stability condition (30).
6. Conclusion
We performed an extensive study of the evolution of a magnetic field in a differentially rotating radiative zone. Among other results, we provided a systematic estimate of the ratio between the toroidal and the poloidal magnetic field. We tested the influence of the diffusivities and of the density profile and characterized the impact of different boundary conditions and differential rotation profiles. The large parametric study that 2D numerical simulations made possible has allowed finding asymptotic regimes for the diffusion and the density contrast. This in turn enabled us to extrapolate our numerical results to realistic stellar diffusivities and density contrasts. A simple onedimensional standing wave model was also found to provide reliable estimates of max(B_{ϕ}/B_{p}) and the time at which this maximum is reached.
We then discussed the stability of the magnetic configurations dominated by the toroidal component using a local dispersion relation derived by Ogilvie (2007) and results of full 3D numerical simulations performed in a simplified case (Jouve et al. 2015). We argued that these magnetic configurations are likely to be subject to a magnetorotational instability or to the Tayler instability, the nature of the instability depending solely on the ratio of the rotation rate to the toroidal Alfvén frequency Ω /ω_{Aϕ}. Such nonaxisymmetric instabilities are expected to modify the largescale axisymmetric configuration when the surface poloidal field is weaker than a threshold value B_{crit}. In the context of the Aurière scenario to explain the dichotomy of intermediatemass star magnetism, our estimate of B_{crit} is compatible with the lower bound of Ap/Bp magnetic fields, although a precise quantitative comparison remains difficult at this stage.
We neglected meridional flows in our simulations. The effect of a prescribed circulation has been considered by Mestel et al. (1988) and Moss et al. (1990). They found that it is negligible as long as the circulation velocity is very low compared to the Alfvén speed, which is the case for typical EddingtonSweettype circulations. Another hypothesis in this work is that of axisymmetry, when it is known from observations that Ap/Bp stars are generally oblique rotators. Nonaxisymmetric components of the initial poloidal field could have an influence on the stability of the magnetic field. Moss et al. (1990), Moss (1992), and Wei & Goodman (2015) have shown that if the magnetic field is symmetrized before the differential rotation decays. In such extreme cases, the magnetic configurations are predicted to be unstable by (30). In the opposite regime, , the strong magnetic field quickly suppresses the differential rotation and the magnetic field remains inclined with respect to the rotation axis. These situations are similar to axisymmetric cases in which magnetic configurations are also predicted to be stable. In intermediate cases, we do not know the effect of the nonaxisymmetric components of the field on the stability conditions. This question needs to be addressed.
To proceed, we need 3D numerical simulations that include the effects of the stable stratification as well as the mechanism that forces the differential rotation. The premainsequence contraction phase can in principle generate the required differential rotation (Lignières et al. 2014). While realistic diffusivities and density contrasts are beyond what can be reached by 3D simulations, we expect that the dependency of B_{crit} with the rotation rate will not strongly depend on these parameters. The relation between B_{crit} and the surface rotation will then be compared to observational constraints on the lower bound of Ap/Bp magnetic fields. Another application of such simulations is to consider the role of nonaxisymmetric instabilities on the efficiency of the angular momentum transport in stellar radiative zones.
Acknowledgments
The authors acknowledge financial support from the Programme National de Physique Stellaire (PNPS) of CNRS/ INSU and from the Agence Nationale pour la Recherche (ANR) through the IMAGINE project. T.G. is supported by the Special Priority Program 1488 PlanetMag of the German Science Foundation. The authors wish to thank Sébastien Deheuvels for providing us with a density profile from a stellar evolutionary model.
References
 Aurière, M., Wade, G.A., Silvester, J., Lignières, F., et al. 2007, A&A, 475, 1053 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Balbus, S.A., & Hawley, J.F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Braithwaite, J. 2007, A&A, 469, 275 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Braithwaite, J. 2009, MNRAS, 397, 763 [NASA ADS] [CrossRef] [Google Scholar]
 Braithwaite, J., & Nordlund, Å. 2006, A&A, 450, 1077 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cantiello, M., Mankovich, C., Bildsten, L., et al. 2014, ApJ, 788, 93 [NASA ADS] [CrossRef] [Google Scholar]
 Charbonneau, P., & MacGregor, K.B. 1992, ApJ, 387, 639 [NASA ADS] [CrossRef] [Google Scholar]
 Charbonneau, P., & MacGregor, K.B. 1993, ApJ, 417, 762 [NASA ADS] [CrossRef] [Google Scholar]
 Deheuvels, S., Dogan, G., Goupil, M., et al. 2014, A&A, 564, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dormy, E., Cardin, P., & Jault, D. 1998, Earth Planet. Sci. Lett., 160, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Ferraro, V.C.A. 1937, MNRAS, 97, 458 [NASA ADS] [CrossRef] [Google Scholar]
 Flowers, E., & Ruderman, M.A. 1977, ApJ, 215, 302 [NASA ADS] [CrossRef] [Google Scholar]
 Heyvaerts, J., & Priest, E.R. 1983, A&A, 117, 220 [NASA ADS] [Google Scholar]
 Ionson, J.A. 1978, ApJ, 226, 650 [NASA ADS] [CrossRef] [Google Scholar]
 Jouve, L., Brun, A.S., Arlt, R., Brandenburg, A., et al. 2008, A&A, 483, 949 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jouve, L., Gastine, T., & Lignières, F. 2015, A&A, 575, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lignières, F., Petit, P., Aurière, M., Wade, G.A., & Böhm, T. 2014, IAU Symp., 302, 338 [NASA ADS] [Google Scholar]
 Linton, M.G., Longcope, D.W., & Fisher, G.H. 1996, ApJ, 469, 954 [NASA ADS] [CrossRef] [Google Scholar]
 Markey, P., & Tayler, R.J. 1973, MNRAS, 163, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Masada, Y., Sano, T., & Takabe, H. 2006, ApJ, 641, 447 [NASA ADS] [CrossRef] [Google Scholar]
 Mestel, L., & Weiss, N.O. 1987, MNRAS, 226, 123 [NASA ADS] [CrossRef] [Google Scholar]
 Mestel, L., Tayler, R. J., & Moss, D. L. 1988, MNRAS, 231, 873 [NASA ADS] [CrossRef] [Google Scholar]
 Moss, D. 1992, MNRAS, 257, 593 [NASA ADS] [Google Scholar]
 Moss, D. L., Mestel, L., & Tayler, R. J. 1990, MNRAS, 245, 550 [NASA ADS] [Google Scholar]
 Ogilvie, G. 2007, The Solar Tachocline, eds. D. W. Hughes, R. Rosner, & N. O. Weiss (Cambridge, UK: Cambridge University Press), 299 [Google Scholar]
 Parker, E.N. 1991, ApJ, 376, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Pitts, E., & Tayler, R.J. 1985, MNRAS, 216, 139 [NASA ADS] [CrossRef] [Google Scholar]
 Rieutord, M., Dintrans, B., Lignières, F., Corbard, T., & Pichon, B. 2005, in SF2A2005: Semaine de l’Astrophysique Française, eds. F. Casoli, T. Contini, J. M. Hameury, & L. Pagani, 759 [Google Scholar]
 Rüdiger, G., & Kitchatinov, L.L. 1996, ApJ, 466, 1078 [NASA ADS] [CrossRef] [Google Scholar]
 Rüdiger, G., Gellert, M., Spada, F., & Tereshin, I. 2015, A&A, 573, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schaeffer, N., Jault, D., Cardin, P., & Drouard, M. 2012, Geophys. J. Int., 191, 508 [CrossRef] [Google Scholar]
 Spada, F., Lanzafame, A.C., & Lanza, A.F. 2010, MNRAS, 404, 641 [NASA ADS] [CrossRef] [Google Scholar]
 Spruit, H.C. 1999, A&A, 349, 189 [NASA ADS] [Google Scholar]
 Spruit, H.C. 2002, A&A, 381 [Google Scholar]
 Tayler, R.J. 1973, MNRAS, 161, 365 [NASA ADS] [CrossRef] [Google Scholar]
 Wei, X., & Goodman, J. 2015, ApJ, 806, 50 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1 Meridional cut of for the radial profile defined by Eq. (12) (left) and for the cylindrical profile defined by Eq. (15) (right). The dotted black lines represent the isocontours of . The white lines represent the poloidal magnetic field lines and the red line the one on which Alfvén waves propagate the slowest. The radius at the inner boundary is r = R_{c} = 0.3 R. 

Open with DEXTER  
In the text 
Fig. 2 at the equator as a function of for different density contrasts ρ_{c}/ρ_{0}. A n = 3 polytropic profile was cut at different radii to produce these various ρ_{c}/ρ_{0}. 

Open with DEXTER  
In the text 
Fig. 3 The reference case: meridional cut of , and for different times. t_{max}, the time at which B_{ϕ}/B_{p} is maximal, is equal to 1.85 t_{Ap} in this simulation. The colourbar at the bottom of each panel indicates the minimal and maximal values of each variable. The parameters of this simulation are L_{u} = 10^{2}, a density contrast of ρ_{c}/ρ_{0} = 7 × 10^{3}, an initial radial differential rotation given by Eq. (12) and a fixed rotation at the inner radius. 

Open with DEXTER  
In the text 
Fig. 4 Evolution of the spatial maximal value of for the reference case. The horizontal dashed line indicates and the vertical dashed line indicates . 

Open with DEXTER  
In the text 
Fig. 5 (top panels) and (bottom panels) as a function of s/L for different times. s is the curvilinear coordinate of the poloidal field line in which max(B_{ϕ}/B_{p}) is reached, and L the length of this field line. s = 0 locates the inner radius and s = L the equator. The simulation is the same as in Fig. 3 with the exception of L_{u} = 10^{4}. The vertical dashed line indicates s_{max} the location of max(B_{ϕ}/B_{p}) on the poloidal field line. 

Open with DEXTER  
In the text 
Fig. 6 max(B_{ϕ}/B_{p}) (top) and t_{max} (bottom) as a function of L_{u} with the physical ingredients used for the reference case. 

Open with DEXTER  
In the text 
Fig. 7 cosθ_{max}, (top frames), and (bottom frames) as a function of ρ_{c}/ρ_{0}. The parameters of the simulations are L_{u} = 10^{5}, a profile of Ω defined by Eq. (12) and Ω fixed at the inner radius. 

Open with DEXTER  
In the text 
Fig. 8 as a function of s/L for different times. s is the curvilinear coordinate of the poloidal field line on which max(B_{ϕ}/B_{p}) is reached, and L the length of this field line. s = 0 locates the inner radius and s = L the equator. The parameters of the simulation are L_{u} = 245, ρ_{c}/ρ_{0} = 7 × 10^{3}, the differential rotation profile given by Eq. (12) and a stressfree boundary condition. 

Open with DEXTER  
In the text 
Fig. 9 3D view of the magnetic configuration obtained with our 2D simulations at t_{max}. It is computed for L_{u} = 10^{2}, a density contrast ρ_{c}/ρ_{0} = 7 × 10^{3}, the radial differential rotation profile defined by Eq. (12) and Ω fixed at the inner radius. 

Open with DEXTER  
In the text 
Fig. 10 Contours of the toroidal magnetic field (white lines) and of the growth rate of the instability (coloured contours) for the most unstable perturbation (m = 7 for the top frame and m = 1 for the bottom frame) which is in the direction e = e_{ϖ}. The parameters of the simulations are N = 0 and, respectively for the top and the bottom frames Ω /ω_{Aϕ} ≈ 10 and Ω /ω_{Aϕ} ≈ 1 at the maximum of the toroidal field. 

Open with DEXTER  
In the text 