© ESO, 2015

1. Introduction

Dynamical processes associated with the presence of a magnetic field in a stellar radiative zone have been the subject of various studies in the past decades. The main questions addressed by these studies are:

• What is the origin of the magnetic field observed at the surface ofroughly 10% of intermediate-mass and massive stars that are incapable of producing a convective dynamo as we know it from cooler stars?

• What is the effect of magnetic fields on the transport of angular momentum under the general conditions existing in a stellar radiative zone?

Interest in the first question has recently revived thanks to large spectropolarimetric surveys performed with the Narval and ESPaDOnS instruments. The chemically peculiar Ap-Bp stars, representing around 10% of the intermediate-mass star population, were confirmed to host invariable dipolar-like magnetic fields with intensities ranging from approximately 300 G to 30 kG (Donati & Landstreet 2009). Another class of magnetism has now been revealed in A-type stars thanks to recent improvements in the detection limit of those spectropolarimeters; it is called Vega-like magnetism, after the bright star in which it was first discovered (Lignières et al. 2009; Petit et al. 2011). The structure of the magnetic field in those stars is different from the strong dipole observed at the surface of Ap stars. Vega exhibits a surface magnetic field organised at much smaller scales and typical amplitudes for the line of sight field of less than 1 G. In these intermediate-mass stars and for both types of magnetism, the origin of the observed magnetic field is still unknown (e.g. Lignières et al. 2014). Various assumptions have been formulated in the literature fairly recently, some of which have already been ruled out. For example, it was argued that the surface magnetic field could be the consequence of the emergence of buoyant concentrations of field produced by a powerful convective dynamo acting in the core of these stars. Unfortunately, these fields are not likely to be able to pop up at the stellar surface within the lifetime of the star (MacGregor & Cassinelli 2003). Another hypothesis relies on the existence of a dynamo in the radiative envelope (Spruit 2002; Braithwaite 2006), but we still lack a definite answer to the question of whether such a process could really take place in stellar interiors (Zahn et al. 2007) and be efficient at producing the observed magnetic field. The most popular hypothesis is thus the one of a fossil field which would be the remnant of an early phase of the stellar evolution. In this case, the dichotomy between the strong dipolar fields of Ap stars and the sub-Gauss magnetism of the other A stars does not come up naturally and thus opens new questions on the origin of the magnetism of these stars. Tentative explanations for this dichotomy have involved the existence of an instability of a magnetic field wound up by a sufficiently long-lasting differential rotation (Aurière et al. 2007). Assuming a continuous initial distribution of fossil fields, Aurière et al. (2007) proposed that below a given threshold magnetic fields would become unstable and that small scales would consequently be created in the horizontal direction. The observed field (i.e. the disc-averaged field) would in this case become very small because of cancellation effects between the opposite polarities of the destabilised field configuration. The unstable cases would then be responsible for the Vega-like magnetism, whereas the stable cases (when the initial field is strong and differential rotation is rapidly quenched) would produce the strong dipolar structures of Ap stars.

The question of the angular momentum transport by magnetic fields in stellar radiative zone is a more general question that has also regained interest thanks to recent asteroseismic observations. The generality of this question arises because most stars possess a radiative zone somewhere in their interior. Initially, the problem of angular momentum transport was addressed in the context of the Sun’s radiative core which was shown by helioseismic inversions to rotate almost rigidly (Schou et al. 1998). Hydrodynamical transport processes such as meridional circulation or hydrodynamical instabilities proved to be inefficient at imposing such a low degree of differential rotation in this region (see Maeder & Meynet 2000, for a review on transport processes induced by rotation). Evidence of processes efficient at transporting angular momentum has been obtained in very different types of stars such as red giant or pre-main-sequence solar-type stars. Indeed, asteroseismology applied to sub-giants and red-giants stars observed by the CoRoT and Kepler satellites has allowed the rotation rate of their core to be constrained using the splitting of mixed modes (Deheuvels et al. 2012, 2014; Beck et al. 2012; Mosser et al. 2012a,b). At the end of the main sequence, the core of these stars contracts and thus tends to spin-up. However, when angular momentum transport caused only by rotational mixing and circulation is considered, the core is found to rotate two to three orders of magnitude faster than what is observed (Eggenberger et al. 2012; Marques et al. 2013). An efficient extraction of angular momentum has thus occurred at some point during the evolution of these stars, which is not reproduced by purely hydrodynamical processes. However, this transport process does not lead to uniform rotation since the asteroseismic analysis shows that a significant level of differential rotation still exists between the core and the envelope (Deheuvels et al. 2012, 2014). A similar situation is encountered in pre-main sequence (PMS) solar-type stars which are expected to be fast rotators partly because they are contracting towards the main sequence. This process, together with the coupling with their circumstellar disc, should spin them up close to break-up velocity (Hartmann & Stauffer 1989). However, the measurement of rotation rates of PMS stars show that they rotate typically at only a tenth of their break-up velocity (Bouvier et al. 1986). Clearly, a powerful mechanism is acting to prevent young stars from spinning up during their early PMS evolution. Semi-empirical models have been developed to explain the rotational history of solar-type stars from birth to the early main-sequence. They include the extraction of angular momentum through the interaction with a disc and magnetised winds while the internal transport of angular momentum is modelled by a coupling timescale between the radiative core and the convective envelope (Gallet & Bouvier 2013). Both an efficient transport of angular momentum in the radiative zone and a certain level of differential rotation are required in these models to fit the rotational velocities observed at the ZAMS. In all these cases, magneto-hydrodynamical (MHD) processes could possibly provide an additional efficient way of extracting angular momentum from stellar interiors and could thus explain the observed rotation profile in those stars.

The case of a poloidal field embedded in a differentially rotating shell is a canonical configuration to study the interplay between magnetic field and differential rotation. It has been mostly studied under the assumption of axisymmetry. In this case, the angular momentum evolution is governed by Ohmic diffusion of Alfvén waves propagating along the field lines. The phase mixing between Alfvén waves of adjacent field lines creates strong transverse gradients and controls the damping time of these waves (Charbonneau & MacGregor 1992; Spada et al. 2010). This picture is however too simplistic as various non-axisymmetric instabilities are likely to develop (see review by Spruit 1999), leading to a profound modification of the efficiency of the angular momentum transport as well as of the magnetic configuration. Under the conditions existing in a stellar radiative region, three main types of instabilities are likely to occur, which differ by the source of free energy they derive from. As reviewed by Spruit (1999), gravitational energy can be released through Parker-like instabilities (Parker 1966), magnetic free energy through purely magnetic instabilities (like the Tayler instability, Tayler 1973) while shear instabilities (like the magneto-rotational instability) depend on the free kinetic energy in the differential rotation (Balbus & Hawley 1991). In the stellar context, most studies have been dedicated to the first two types of instabilities. The buoyancy instability has been particularly used to explain the formation of active regions at the surface of cool stars. Unstable regions of toroidal field are thought to rise through the convective envelope from their region of creation (the tachocline) to the stellar surface to create bipolar magnetic structures (Choudhuri & Gilman 1987; Caligari et al. 1995; Emonet & Moreno-Insertis 1998; Fan 2008; Jouve & Brun 2009; Favier et al. 2012; Jouve et al. 2013). The Parker instability has also been invoked as a possible regeneration process for the poloidal magnetic field in the Sun and thus thought to take a large part in the dynamo loop in solar-like stars (Cline et al. 2003). The Tayler instability (TI) has been mostly considered to assess the stability of magnetic equilibrium solutions in stellar radiative environments. Indeed, Tayler (1973), Wright (1973) and Markey & Tayler (1973) were the first to show that a purely poloidal or purely toroidal magnetic field would be unstable under all general circumstances and that a stable equilibrium would thus have to be of mixed poloidal/toroidal type. This idea was then investigated more thoroughly including additional physical processes like rotation, and confirmed both by analytical studies and numerical models (Braithwaite & Nordlund 2006; Braithwaite 2007; Bonanno & Urpin 2008, 2013; Duez & Mathis 2010). The magneto-rotational instability (MRI) has been extensively analysed to understand the dynamics of accretion discs, in particular the transport of angular momentum induced by such an instability (see Fromang 2013, for a complete lecture on the subject). It has also been used in the context of planetary interiors (Petitdemange et al. 2008, 2013) where a modified version of the MRI was shown to be able to operate in the magnetospheric force balance relevant to the Earth’s core. In stellar physics, the MRI of a purely poloidal or purely toroidal field was studied using local stability criteria (Menou et al. 2004; Masada et al. 2006) or using global analysis in spherical (Arlt et al. 2003) or cylindrical (Ogilvie & Pringle 1996; Rüdiger et al. 2014, 2015) geometries, with a particularfocus on the induced transport of angular momentum in the last two articles. However, considering the level of differential rotation and the presence of strong magnetic fields in various situations encountered in stellar radiative zones, this mechanism could have attracted more attention. The main reason for favouring the TI over the MRI has been the strong stratification of stellar radiative zones, which has been thought to affect the MRI more severely (Spruit 1999).

We consider here the initial value problem of a purely poloidal field in a differentially rotating unstratified spherical shell. The origin of the magnetic field and the differential rotation is not addressed, they are assumed to be the result of physical processes at play during the star’s history. The differential rotation acting on the poloidal field will naturally produce a mixed poloidal/toroidal configuration which will interact with the rotation. The initial differential rotation is not enforced by the boundary conditions, so the final state is known to be solid body rotation. In a companion paper (Gaurat et al., in prep.), the evolution of the axisymmetric system was studied in 2D and tentative predictions were formulated on the stability of the induced magnetic configurations. In this work where 3D numerical simulations are performed, we wish to address the two questions quoted at the beginning of this introduction, with a particular focus on the following more specific interrogations:

• What are the stability conditions for an initial poloidal fieldwound up by differential rotation?

• What is the nature of the instability if there is any?

• What is the impact of the instability on the efficiency of the angular momentum transport?

• What is the impact of the instability on the large scale magnetic field?

• What would be the consequences for the observed magnetic fields and how does this relate to the properties of the magnetism of A-type stars?

After presenting the numerical setup in Sect. 2, we treat separately the case of the instability of the initial purely poloidal field, without rotation or shear (Sect. 3). This enables us to show the limited impact the poloidal field instability will have on the timescales considered in the next sections. Section 4 deals with the instability of the wound-up field, with various initial poloidal field configurations. In this section, a particular focus was made on the transport of angular momentum induced by the instability. We then interpret the results of this work in the context of the magnetic dichotomy observed in A-type stars in Sect. 5. A brief summary of the consequences of this work is then given as a conclusion in Sect. 6.

2. Numerical model

2.1. Governing equations

Assuming constant kinematic viscosity ν, magnetic diffusivity λ and density ρ₀, the system of equations for an incompressible magnetised fluid reads

$\begin{array}{ccc}& & {\nabla }\mathrm{·}{u}\mathrm{=}\mathrm{0}\\ & & {\mathit{\rho }}_{\mathrm{0}}\left(\frac{\mathit{\partial }{u}}{\mathit{\partial t}}\mathrm{+}{u}\mathrm{·}{\nabla }{u}\right)\mathrm{=}\mathrm{-}{\nabla }\mathit{p}\mathrm{+}\frac{\mathrm{1}}{{\mathit{\mu }}_{\mathrm{0}}}\mathrm{(}{\nabla }\mathrm{\times }{B}\mathrm{)}\mathrm{\times }{B}\mathrm{+}{\mathit{\rho }}_{\mathrm{0}}\mathit{\nu }{\Delta }{u}\\ & & \end{array}$(1)The equations are then non-dimensionalised using the viscous diffusion time d²/ν as a reference timescale (where d = r_o − r_i is the thickness of the spherical shell), ν/d as a velocity scale, ρ₀ as a density scale and $\sqrt{{\mathit{\rho }}_{\mathrm{0}}{\mathit{\mu }}_{\mathrm{0}}}\mathit{\lambda }\mathit{/}\mathit{d}$ as a typical scale for the magnetic field. The system of Eq. (1) can then be rewritten as $\begin{array}{ccc}& & \frac{\mathit{\partial }{u}}{\mathit{\partial t}}\mathrm{+}{u}\mathrm{·}{\nabla }{u}\mathrm{=}\mathrm{-}{\nabla }\mathit{p}\mathrm{+}\frac{\mathrm{1}}{\mathrm{P}{\mathrm{m}}^{\mathrm{2}}}\mathrm{(}{\nabla }\mathrm{\times }{B}\mathrm{)}\mathrm{\times }{B}\mathrm{+}{\Delta }{u}\mathit{,}\\ & & \end{array}$(2)The initial seed field imposed at the beginning of the simulation is then given in units of $\frac{{\mathit{B}}_{\mathrm{0}}\mathit{d}}{\sqrt{{\mathit{\rho }}_{\mathrm{0}}{\mathit{\mu }}_{\mathrm{0}}}\mathit{\lambda }}\mathrm{=}\frac{{\mathit{t}}_{\mathit{\lambda }}}{{\mathit{t}}_{\mathrm{ap}}}\mathrm{=}\mathrm{L}{\mathrm{u}}_{\mathrm{p}}\mathit{,}$(3)where B₀ is a measure of the surface magnetic field at the poles. ${\mathit{t}}_{\mathrm{ap}}\mathrm{=}\sqrt{{\mathit{\rho }}_{\mathrm{0}}{\mathit{\mu }}_{\mathrm{0}}}\mathit{d}\mathit{/}{\mathit{B}}_{\mathrm{0}}$ is the Alfvén timescale based on the poloidal field and t_λ = d²/λ is the magnetic diffusion timescale. $\mathrm{L}{\mathrm{u}}_{\mathrm{p}}\mathrm{=}{\mathit{B}}_{\mathrm{0}}\mathit{d}\mathit{/}\sqrt{{\mathit{\rho }}_{\mathrm{0}}{\mathit{\mu }}_{\mathrm{0}}}\mathit{\lambda }$ is the Lundquist number based on the poloidal magnetic field and Pm = ν/λ is the magnetic Prandtl number.

If, in addition to an initial seed, we also impose a differential rotation profile, then the maximal value of the rotation rate Ω₀, located on the axis of rotation, is proportional to the Reynolds number Re. The initial velocity is thus given in units of $\frac{{\mathrm{\Omega }}_{\mathrm{0}}{\mathit{d}}^{\mathrm{2}}}{\mathit{\nu }}\mathrm{=}\mathrm{Re}\mathit{.}$(4)This hydrodynamical problem is therefore controlled by three dimensionless parameters, namely $\mathrm{Re}\mathrm{=}\frac{{\mathrm{\Omega }}_{\mathrm{0}}{\mathit{d}}^{\mathrm{2}}}{\mathit{\nu }}\mathrm{;} \mathrm{L}{\mathrm{u}}_{\mathrm{p}}\mathrm{=}\frac{{\mathit{B}}_{\mathrm{0}}\mathit{d}}{\sqrt{{\mathit{\rho }}_{\mathrm{0}}{\mathit{\mu }}_{\mathrm{0}}}\mathit{\lambda }}\mathrm{;} \mathrm{Pm}\mathrm{=}\frac{\mathit{\nu }}{\mathit{\lambda }}\mathrm{·}$(5)In all our calculations, the value of the magnetic Prandtl number Pm has been set to 1. The study of the influence of this parameter is postponed to future works.

2.2. Numerical method and boundary conditions

The numerical simulations were computed with the code MagIC (Wicht 2002) which solves the MHD equations in a spherical shell using a poloidal toroidal decomposition for the velocity and the magnetic fields, $\begin{array}{ccc}{u}& \mathrm{=}& \mathrm{\nabla }\mathrm{\times }\mathrm{\nabla }\mathrm{\times }\hspace{0.17em}\mathrm{\left(}\mathit{W}\hspace{0.17em}{e}\mathrm{r}\mathrm{\right)}\mathrm{+}\mathrm{\nabla }\mathrm{\times }\hspace{0.17em}\mathrm{\left(}\mathit{Z}\hspace{0.17em}{e}\mathrm{r}\mathrm{\right)}\mathit{,}\\ {B}& \mathrm{=}& \mathrm{\nabla }\mathrm{\times }\mathrm{\nabla }\mathrm{\times }\hspace{0.17em}\mathrm{\left(}\mathit{C}\hspace{0.17em}{e}\mathrm{r}\mathrm{\right)}\mathrm{+}\mathrm{\nabla }\mathrm{\times }\hspace{0.17em}\mathrm{\left(}\mathit{D}\hspace{0.17em}{e}\mathrm{r}\mathrm{\right)}\mathit{,}\end{array}$where \begin{lxirformule}$W$\end{lxirformule} (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). Typical numerical resolutions employed in this study range from (N_r = 65, l_max = 170) for the more diffusive cases (i.e. low Lu_p, low Re) to (N_r = 97, l_max = 426) for the less diffusive ones. The spherical shell considered extends from longitude φ = 0 to φ = 2π, from colatitude θ = 0 to θ = π and from radius r = r_i = η/ (1 − η) to r = r_o = 1/(1 − η) where η = r_i/r_o is the aspect ratio of the shell, equal to 0.3 in all our calculations.

For the boundary conditions, we assume stress-free boundary conditions and non-penetrating radial velocity at top and bottom.

The magnetic field is matched to an external potential field at the top and the bottom boundaries: ${B}\mathrm{=}{\nabla }\mathrm{\Phi } \mathrm{\to } \mathrm{\Delta \Phi }\mathrm{=}\mathrm{0} \mathrm{for} \mathit{r}\mathrm{=}\mathrm{[}{\mathit{r}}_{\mathit{i}}\mathit{,}{\mathit{r}}_{\mathrm{o}}\mathrm{]}\mathit{.}$(8)We note here that applying such potential top and bottom boundary conditions to an internal dipolar field is not physical and thus leads to a mismatch between the interior and exterior of the domain. This mismatch results in the creation of a strong current sheet at the base of the domain which needs to be carefully treated. This is why the purely dipolar case for the initial poloidal field will not be considered here.

2.3. Initial conditions

In all our calculations, we will study the stability of a magnetic field configuration which results from the winding up of an initial poloidal field by differential rotation.

For the initial poloidal field, we choose a simple configuration which has a similar structure as a pure dipole (a l = 1, m = 0 configuration) but with an azimuthal current (the only non-zero component) that is constant in radius. The field lines at low latitudes are closed inside the domain and open up above around 30° in latitude. The advantage of this configuration is that contrary to a dipolar structure, the mismatch between the interior of the domain and the potential boundaries does not create strong current sheets likely to be unstable and to be difficult to handle numerically. The expression for the initial poloidal field is as follows: $\begin{array}{ccc}{\mathit{B}}_{\mathit{r}}& \mathrm{=}& \mathrm{3}\cos \mathit{\theta }\hspace{0.17em}\frac{\mathit{r}\mathrm{-}\mathrm{4}{\mathit{r}}_{\mathrm{o}}\mathit{/}\mathrm{3}\mathrm{+}{\mathit{r}}_{\mathit{i}}^{\mathrm{4}}\mathit{/}\mathrm{3}{\mathit{r}}^{\mathrm{3}}}{{\mathit{r}}_{\mathrm{o}}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{(}{\mathit{r}}_{\mathit{i}}\mathit{/}{\mathit{r}}_{\mathrm{o}}{\mathrm{)}}^{\mathrm{4}}\mathrm{)}}\mathit{,}\\ {\mathit{B}}_{\mathit{\theta }}& \mathrm{=}& \mathrm{-}\frac{\mathrm{3}\sin \mathit{\theta }}{\mathrm{2}}\hspace{0.17em}\frac{\mathrm{3}\mathit{r}\mathrm{-}\mathrm{8}{\mathit{r}}_{\mathrm{o}}\mathit{/}\mathrm{3}\mathrm{-}{\mathit{r}}_{\mathit{i}}^{\mathrm{4}}\mathit{/}\mathrm{3}{\mathit{r}}^{\mathrm{3}}}{{\mathit{r}}_{\mathrm{o}}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{(}{\mathit{r}}_{\mathit{i}}\mathit{/}{\mathit{r}}_{\mathrm{o}}{\mathrm{)}}^{\mathrm{4}}\mathrm{)}}\mathit{,}\\ {\mathit{B}}_{\mathit{\phi }}& \mathrm{=}& \end{array}$(9)This poloidal magnetic field is normalized such that the surface polar value is equal to 1.

The profile adopted for the initial differential rotation only depends on the dimensionless cylindrical radius s = rsinθ. It is normalised at 1 on the rotation axis. It is the same profile as the one used by Arlt & Rüdiger (2011b) which is consistent with a system rotating relatively fast (tending to satisfy the Taylor-Proudman theorem) and which possesses a certain degree of differential rotation, without being linearly unstable to hydrodynamical instabilities: $\mathrm{\Omega }\mathrm{=}\frac{\mathrm{1}}{\sqrt{\mathrm{1}\mathrm{+}{\mathit{s}}^{\mathrm{4}}}}\mathrm{·}$(10)Assuming such a profile, we adopt a setup where the shear parameter q = | ∂lnΩ /∂lns | ≈ 1, which constitutes an additional assumption in our model where the shear strength is approximately equal to the rotation rate itself. We note that when a radial (or shellular) profile is used for the differential rotation, a strong meridional circulation is created which quickly advects the poloidal field and strongly modifies its configuration before the differential rotation has time to build a significant toroidal field. This situation is not convenient since we are interested in the possibility of a wound-up field to become unstable. With a cylindrical profile, the amplitude of the meridional flow is extremely weak and the main process acting on the poloidal field will thus be the winding-up by the differential rotation.

Fig. 1 Azimuthal averages: initial rotation profile represented as isocontours of Ω (colours) and initial magnetic field represented by the poloidal field lines (dashed contours). Ω is given in units of Re.

Figure 1 shows the initial velocity and magnetic field profiles used in the simulations. The magnetic field lines (dashed lines) are superimposed to the isocontours of Ω = v_φ/s represented in colour. This type of representation helps to anticipate that the strongest winding-up of the magnetic field lines (in other words the strongest Ω-effect, measured by the intensity of the term (sB_p·∇Ω) in the induction equation) will occur where the angle between the isocontours of Ω and the magnetic field lines is maximum. From this figure, we can argue that a lobe of toroidal field will be predominantly created at mid-latitudes in each hemisphere, confined quite close to the top of the domain and that the magnetic field will change sign at the equator where the Ω-effect vanishes.

Fig. 2 Lup = 200: temporal evolution of the magnetic energy contained in the first azimuthal wavenumbers (left), volume rendering showing isosurfaces of the fluctuating part of the radial magnetic field ${\mathit{B}}_{\mathit{r}}^{\mathrm{\prime }}$ together with magnetic field lines during the linear phase of the instability (middle) and cut of the total radial magnetic field at the top of our domain when the instability has reached a saturated state.

Fig. 3 Volume rendering of the instability for Lup = 400 (left) and Lup = 1000 (right). Surfaces of constant ${\mathit{B}}_{\mathit{r}}^{\mathrm{\prime }}$ are represented, as well as some magnetic field lines inside the domain.

3. Instability of the initial poloidal field, without rotation and shear

3.1. Instability criterion

In 1973, Markey & Tayler showed that, if a star contains a poloidal magnetic field which possesses field lines closed within the star, hydromagnetic instabilities with a growth rate comparable to the Alfvén timescale would be triggered. Flowers & Ruderman (1977) then showed that a poloidal field with all field lines crossing the stellar surface was equally subject to magnetic instabilities. In the case where some field lines are closed in the domain, the magnetic field loops exert pressure on each other and tend to slip out of the equilibrium position (or in other words, be pushed out by pressure from the neighbours, Braithwaite 2007). Since the movement of these loops is restricted in the radial direction by the stable stratification of the star, they move in a direction parallel to the magnetic axis.

To look at the onset of the poloidal field instability, we start from the purely poloidal field whose profile is given by Expression 9 and shown in Fig. 1 and the simulation is run without rotation or shear. Since some field lines are closed inside the domain and since a non-zero toroidal current initially exists, this configuration is a good candidate for non-axisymmetric instabilities to develop. It is expected in this case that high azimuthal wavenumbers will be excited and that the most unstable m will depend on the magnetic diffusivity or, in non-dimensionalised parameters, on the Lundquist number Lu_p.

3.2. Characteristics of the instability

The first calculation is performed for a Lundquist number Lu_p = 200. The initial perturbation is applied at all scales equivalently (white noise) on the poloidal magnetic field with a typical amplitude of about 10^-5. The energy in the non-axisymmetric modes is thus initially around 10^-10. At time t ≈ 0.8t_ap, an instability starts to grow.

Figure 2 shows some features of the poloidal field instability found in this configuration for a particular Lundquist number Lu_p = 200. The left panel shows the temporal evolution of the poloidal magnetic energy contained in the various azimuthal wavenumbers, from m = 0 to m = 11. In this situation, the most unstable wavenumbers are m = 5 and m = 6, as also illustrated in the mid-panel where a volume rendering of the fluctuating component of the radial magnetic field is shown. We find that the instability starts to develop close to the equator and at mid-radius, where the poloidal configuration possesses the strongest toroidal current. The last panel of Fig. 2 shows the total surface radial field when the instability has reached its non-linear phase. We clearly find that the initial simple magnetic field configuration is completely destroyed by the instability, leaving instead a complex non-axisymmetric field with rather small-scale structures very extended in latitude. This is very close to what Braithwaite & Spruit (2006) reported in their models of the instability of a poloidal field in a neutron star. They initiated their models with a different initial field configuration, inspired from the work of Flowers & Ruderman (1977) where the field is uniform in the domain and matches a potential field outside. They studied the instability in a compressible, stably-stratified atmosphere (when we are incompressible) but the typical structure of the radial magnetic field at the top of their domain is similar to what is found here in our simulations.

The effects of magnetic diffusion were then investigated by varying the Lundquist number Lu_p. Three additional cases were calculated, with Lu_p = 100, Lu_p = 400 and Lu_p = 1000. All cases were unstable with respect to the poloidal field instability. The surfaces of constant ${\mathit{B}}_{\mathit{r}}^{\mathrm{\prime }}$ during the linear phase of the instability are shown in Fig. 3 for Lu_p = 400 and Lu_p = 1000. It is clear from this figure that the wavenumber with the highest growth rate (or the favoured wavenumber) critically depends on the Lundquist number. Indeed, m = 4 is favoured for Lu_p = 100 (not shown), m = 10 for Lu_p = 400 and m = 14 − 15 for Lu_p = 1000, the length scale of the instability thus being strongly dependent on the magnetic diffusivity. The growth rate of the instability was measured in those various cases and it is found that, in agreement with theoretical studies, it is equal to the Alfvén frequency at the location of the instability, namely where the field possesses the strongest currents, around mid-radius and at the equator. In all cases, the non-linear outcome of the instability is the complete destruction of the large-scale axisymmetric structure into a complex non-axisymmetric field configuration.

We note here that for the poloidal field instability, the unstable modes tend to concentrate where the current is the strongest, which is consistent with the fact that this instability consists in releasing the magnetic energy contained in this non-potential configuration. As a consequence, a purely dipolar field, which does not possess currents, should not be unstable. However, as stated before, imposing a dipolar field inside the domain and a potential field outside will necessarily create a mismatch at one of the boundaries. This mismatch results in a strong current sheet near the base of the domain which can be the seed for the development of the poloidal field instability. This is probably what Zahn et al. (2007) found when they argue that the high-m instability of their initial dipole (which should be current-free and thus stable) is a classical poloidal field instability. It is indeed clear in their results that the instability develops at the base of their domain where a current sheet is artificially created because of the mismatch between their dipolar field and the potential field boundary conditions. We chose in this work to concentrate on an initial field configuration which would not lead to the creation of artificially strong current sheets close to the boundaries, in order to avoid possible numerical instabilities.

We thus now have to take into account that such a poloidal field configuration can and will become unstable on an Alfvén crossing time in our system. When solid-body rotation is added, Braithwaite (2007) showed that no significant effects were seen on the linear phase of the instability but that the non-linear outcome would be modified. In the cases we wish to study, the initial poloidal field configuration will be quickly modified by the shearing by differential rotation but our configurations may still be affected by the poloidal field instability on timescales of a few poloidal Alfvén times. However, we will see that the instabilities that are found in a case of a wound-up field grow on a much shorter timescale and the effects of the poloidal field instability on our results will thus be very limited.

4. Shearing the initial poloidal field

Now that we have checked the properties of the purely poloidal field instability in our simulations, we wish to study the stability of a magnetic field resulting from the action of differential rotation on this initial poloidal field. We will first analyse the large-scale axisymmetric evolution, with a particular focus on the configuration of the generated toroidal field.

4.1. Evolution of the unperturbed system

The initial evolution is a linear growth of the toroidal magnetic field intensity created by the winding-up of poloidal field lines by differential rotation (the so-called Ω-effect). We quickly get two lobes of toroidal field of opposite sign, one in each hemisphere. This antisymmetry is simply due to the initial configuration of our poloidal field and our differential rotation which produces an Ω-effect of opposite sign in both hemispheres. This phase is well understood by looking at the induction equation in an ideal MHD case in a kinematic regime (i.e. when the velocity is considered constant in time): $\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{B}}_{\mathit{\phi }}}{\mathit{\partial t}}\mathrm{=}\mathit{s}\hspace{0.17em}{B}\mathrm{p}\mathrm{·}\mathrm{\nabla }\mathrm{\Omega }& & \end{array}$The resulting toroidal field configuration is a relatively thin region close to the surface whose intensity grows in time as B_φ = sB_p·∇Ω t as long as ∇Ω and B_p are constant, which is approximately the case during this winding-up phase.

Fig. 4 Temporal evolution of the toroidal magnetic (red) and kinetic (blue) energies for an initial poloidal field strength of Lup = 60 (solid lines), Lup = 100 (dashed lines) and Lup = 200 (dash-dotted lines). In all cases, Re = 2 × 104.

An important parameter of our simulations is the ratio between the toroidal Alfvén frequency and the rotation rate ${\mathit{\omega }}_{{\mathrm{A}}_{\mathit{\phi }}}\mathit{/}\mathrm{\Omega }\mathrm{=}{\mathit{B}}_{\mathit{\phi }}\mathit{/}\mathrm{(}\mathit{r}\hspace{0.17em}\mathrm{\Omega }\sqrt{{\mathit{\mu }}_{\mathrm{0}}{\mathit{\rho }}_{\mathrm{0}}}\mathrm{)}$. During the kinematic phase where the profiles of both B_p and Ω are assumed to be constant, we thus expect that this ratio can be written as: $\begin{array}{ccc}\frac{{\mathit{\omega }}_{{\mathrm{A}}_{\mathit{\phi }}}}{\mathrm{\Omega }}\mathrm{(}\mathit{r,\theta ,t}\mathit{/}{\mathit{t}}_{\mathrm{ap}}\mathrm{)}\mathrm{=}\frac{\mathit{d}\hspace{0.17em}\hspace{0.17em}{\mathit{B}}_{\mathrm{p}}\mathrm{(}\mathit{r,\theta ,t}\mathrm{=}\mathrm{0}\mathrm{)}\mathrm{·}\mathrm{\nabla \Omega }\mathrm{(}\mathit{r,\theta ,t}\mathrm{=}\mathrm{0}\mathrm{)}}{{\mathit{B}}_{\mathrm{0}}\hspace{0.17em}\hspace{0.17em}\mathrm{\Omega }\mathrm{(}\mathit{r,\theta ,t}\mathrm{=}\mathrm{0}\mathrm{)}}\hspace{0.17em}\hspace{0.17em}\frac{\mathit{t}}{{\mathit{t}}_{\mathrm{ap}}}& & \end{array}$and is independent of the poloidal Lundquist number Lu_p and on the Reynolds number Re. This ratio will thus only depend on the initial and boundary conditions on the poloidal field and the rotation (see Gaurat et al., in prep.). Since the geometries of the initial poloidal field and differential rotation and the boundary conditions are kept the same in all our cases (except when the influence of the initial poloidal field configuration is investigated), this ratio will only depend on t/t_ap in all our calculations. This property will be of prime interest for the understanding of the development of possible instabilities. Figure 4 shows the evolution of the toroidal magnetic and kinetic energies as a function of t/t_ap, for three different values of the initial poloidal field strength Lu_p and for Re = 2 × 10⁴. It is indeed found that in the kinematic phase, until the toroidal field reaches its maximum value, the evolution of the toroidal magnetic and kinetic energies is similar for all Lu_p. The maximum value reached by the toroidal field is a bit reduced in the case where Lu_p = 60 since diffusive effects are not negligible here.

Fig. 5 Characteristics of the instability for Re = 2 × 104 and Lup = 60. Panels a)–d) are snapshots of, respectively, the azimuthally averaged toroidal field on which the poloidal field lines are superimposed, a meridional cut of ${\mathit{B}}_{\mathit{r}}^{\mathrm{\prime }}$ and a volume rendering showing surfaces of constant ${\mathit{B}}_{\mathit{r}}^{\mathrm{\prime }}$ and magnetic field lines at t = 0.33tap. Panel b) shows the temporal evolution of the toroidal magnetic energy contained in the first 12 azimuthal wavenumbers.

After this kinematic phase, the magnetic tension becomes sufficiently strong to be able to back-react on the differential rotation. A maximum of toroidal magnetic field is then reached at a time which is measured here to be approximately equal to 0.67t_ap for all initial poloidal field strength considered. It is natural that the backreaction occurs at a time which is approximately equal to an Alfvén wave travel time based on the background poloidal magnetic field since the dynamical evolution is entirely governed by the dynamics of Alfvén waves (see Gaurat et al., in prep.).

When the toroidal magnetic field reaches its maximum, the Ω-effect changes sign in the domain, thus locally creating a toroidal magnetic field of opposite sign. As seen in Fig. 4, the subsequent evolution thus consists in oscillations of B_φ and Ω with a period again approximately equal to the period of an Alfvén wave travelling along the poloidal field lines. Periodic oscillations are then damped and a final state of uniform rotation and zero toroidal magnetic field is reached on a time typical of magnetic diffusion acting on the length scale of Alfvén waves. This typical time can be rather short as small scales are created by the phase-mixing mechanism (Gaurat et al., in prep.; Charbonneau & MacGregor 1992). Indeed, Alfvén waves propagating at different speeds on neighbouring surfaces progressively build-up strong gradients of field between those surfaces, thus creating a typical small scale l. We estimate the phase-mixing timescale to be the magnetic diffusion time based on this length scale: t_pm = l²/λ. The ratio between the phase-mixing timescale and the Alfvén timescale t_ap is thus given by: t_pm/t_ap = t_pm/t_λ × t_λ/t_ap = l²/d² × Lu_p. In our cases where Lu_p ≈ 100, we only need to produce scales 3 to 4 times smaller than the size of the domain to get a damping time of the periodic oscillations of a few Alfvén times, as found here and as visible in Fig. 4.

4.2. Evidence for a magnetic instability

The preceding system with Re = 2 × 10⁴ and Lu_p = 60 is now perturbed, at an instant very close to the initial time (around t = 6 × 10^-3t_ap). A white noise perturbation is applied to the poloidal magnetic field to all combinations of l,m except l = 0,m = 0 with the same typical amplitude ${\mathit{b}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{=}{\mathrm{10}}^{-5}$. The system first continues a very similar evolution compared to what we get in the unperturbed case. The toroidal magnetic field grows linearly with time and the poloidal field and differential rotation stay relatively constant. However, as soon as the toroidal magnetic field undergoes a situation where the maximum value of ω_Aφ/ Ω reaches the value of about 0.2, a magnetic instability starts to grow.

The instability is illustrated in Fig. 5. The first panel a shows the averaged toroidal magnetic field generated in this case in coloured levels, while the dashed black lines correspond to the poloidal field lines. Panel b shows the temporal evolution of the toroidal magnetic energy distribution between the first 12 azimuthal wavenumbers, including the axisymmetric mode m = 0. Panel c shows a meridional cut of the fluctuating radial magnetic field (i.e. the total radial field to which the axisymmetric component is subtracted) and panel d is a 3D rendering showing surfaces of constant ${\mathit{B}}_{\mathit{r}}^{\mathrm{\prime }}$ and a few magnetic field lines. Panel b clearly shows the exponential growth of different wavenumbers m, with m = 4,m = 5 and m = 6 being initially favoured. Around t = 0.1t_ap when the instability starts to grow and at the location of the development of the instability, ω_Aφ/ Ω ≈ 0.2 and B_φ/B_p is around 30. We find here that the typical growth rate of the most unstable wavenumbers m = 4,m = 5 and m = 6 is approximately equal to the toroidal Alfvén frequency, i.e. σ/ω_Aφ ≈ 1. From panels c and d which are snapshots taken during the linear phase of the instability, we clearly see that this instability is confined in a region of maximum toroidal field, i.e. close to the surface at a latitude around 40°. The lengthscale of the instability in the latitudinal direction is rather short compared to the typical scale of the background magnetic field, we indeed find typically 8 nodes in the θ-direction inside the envelope of toroidal magnetic field located between latitudes of 20° and 50°.

In previous works with very similar setups (Arlt & Rüdiger 2011a,b; Bonanno & Urpin 2013), such an instability has been attributed to the Tayler instability of the toroidal field. However, several arguments suggest that this is not a Tayler instability that is found here, but a magnetic instability of a different nature. First, the instability does not seem to be mostly located on the gradients of the magnetic field, which would be expected in the case of a Tayler instability. We now come back to the criterion first derived by Goossens & Tayler (1980) in spherical geometry for the Tayler instability of a purely toroidal field of the form ${\mathit{B}}_{\mathit{\phi }}\mathrm{=}\mathit{f}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathit{P}}_{\mathit{l}}^{\mathrm{1}}\mathrm{\left(}\cos \mathit{\theta }\mathrm{\right)}$, where ${\mathit{P}}_{\mathit{l}}^{\mathrm{1}}$ is the associated Legendre polynomial of degree l. Since we work here with spherical coordinates, it is more directly applicable to our cases than the criterion derived by Tayler (1973) in cylindrical geometry. For non-axisymmetric modes (i.e. m ≠ 0), the instability occurs when $\frac{{\mathit{B}}_{\mathit{\phi }}^{\mathrm{2}}}{\mathrm{4}\mathit{\pi }{\mathit{r}}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathit{\theta }}\left({\mathit{m}}^{\mathrm{2}}\mathrm{-}\mathrm{2}{\cos }^{\mathrm{2}}\mathit{\theta }\mathrm{-}\mathrm{2}\sin \mathit{\theta }\cos \mathit{\theta }\frac{\mathrm{1}}{{\mathit{B}}_{\mathit{\phi }}}\frac{\mathit{\partial }{\mathit{B}}_{\mathit{\phi }}}{\mathit{\partial \theta }}\right)\mathit{<}\mathrm{0.}$(11)This criterion clearly indicates that the Tayler instability will first develop where the latitudinal gradients of magnetic field are maximum. This is not what is found here. In addition, this criterion can be used to evaluate the most unstable azimuthal wavenumber. For any configuration, the m = 1 mode is favoured but higher wavenumbers may also be excited if the gradient of B_φ is sufficiently large. We computed this criterion for our particular configuration (even if it cannot be represented on a particular Legendre polynomial) to give us some hints on the possible values of m which could be unstable in our case. Applying Goossens’ criterion, we find that the configuration is stable for m ≥ 4. As higher values of m are excited in our 3D simulation, the Tayler instability cannot solely account for the instabilities observed here. Finally, the growth rate of the instability does not match the theoretical prediction for the Tayler instability when rotation dominates (i.e. Ω >ω_Aφ). As stated by Szklarski & Arlt (2013), and as also found in our simulations, the maximum ratio ω_Aφ/ Ω never reaches values above 0.4, and our instability starts to grow at ω_Aφ/ Ω = 0.2. We note that, in the limit of small magnetic diffusivity, this maximum ratio is independent of the initial poloidal field strength and is only related to the efficiency of the Ω-effect which is itself only related to the initial geometry of the poloidal field and of the differential rotation. In cases where Ω ≫ ω_Aφ, the TI growth rate is reduced by a factor of ω_Aφ/ Ω compared to the non-rotating case, where the growth rate is ω_Aφ. This means that in our calculations, the growth rate should be σ = ω_Aφ/ Ω × ω_Aφ = 0.2 ω_Aφ, when we find a growth rate of ω_Aφ. Szklarski & Arlt (2013) claim that no instability is found in any of their self-consistent calculations. We suspect that this is due to their stable stratification which probably inhibited the instability we found in our unstratified cases.

4.3. Analysis of the instability

We argue here that this is not the Tayler instability which is found in those configurations but the magneto-rotational instability of a magnetic field in the presence of a differential rotation. The configurations of magnetic field obtained in our simulations are rather complex and do not possess typical profiles for which the conditions of stability are known. In this section, we thus proceed to a local linear stability analysis, using the dispersion relation derived by Ogilvie (2007) for the case of a purely toroidal field with differential rotation, both possessing arbitrary profiles. The results of this analysis are then compared with the properties of the instability observed in the numerical simulation.

In his work, Ogilvie (2007) develops the linear stability analysis of a purely toroidal magnetic field with differential rotation in an ideal MHD case. Instead of using a classical normal mode expansion for the disturbances, he considers a solution localised in the meridional plane (on the single point s₀,z₀) in a Gaussian envelope within which the displacement has a plane-wave form. The solutions are thus expressed as: $\begin{array}{ccc}\mathit{\xi }\mathrm{=}\mathrm{\Re }\mathrm{[}\begin{array}{c}\mathrm{˜}\\ \mathit{\xi }\end{array}\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}\exp {\mathrm{(}\mathrm{i}\mathit{m\phi }\mathrm{-}\mathrm{i}\mathit{\omega t}\mathrm{)}}^{\mathrm{]}}& & \end{array}$where $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{\xi }\end{array}$ is further expanded as: $\begin{array}{ccc}\begin{array}{c}\mathrm{˜}\\ \mathit{\xi }\end{array}\mathrm{=}\mathit{\xi ̂}\exp \mathrm{(}\mathrm{i}{\mathit{k}}_{\mathit{s}}\mathit{s}\mathrm{+}\mathrm{i}{\mathit{k}}_{\mathit{z}}\mathit{z}\mathrm{)}\exp \left(\mathrm{-}\frac{\mathrm{(}\mathit{s}\mathrm{-}{\mathit{s}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{l}}^{\mathrm{2}}}\mathrm{-}\frac{\mathrm{(}\mathit{z}\mathrm{-}{\mathit{z}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{l}}^{\mathrm{2}}}\right)& & \end{array}$where is a constant complex vector, l is the localisation scale and k = k_se_s + k_ze_z is a real wavevector. Instability sets in when ℑ(ω) > 0.

This technique was also used to show the existence of linear instabilities in the context of accretion discs (Terquem & Papaloizou 1996). This analysis is relevant in the limits where the spatial scale of the instability is small compared to the scale of the envelope which is itself small compared to the typical scale of variation of the basic state: $\begin{array}{ccc}\mathrm{1}\mathit{/}\mathit{k}\mathrm{\ll }\mathit{l}\mathrm{\ll }\mathit{L}& & \end{array}$where L is the scale of the gradients of the background state.

In our case here, the scale of the instability is indeed small compared to the envelope in which it is contained, but it is less clear if the localisation scale is small compared to the gradients of the basic state. However, we consider the disturbances which are found in our simulations to be sufficiently localised to be able to apply this linear stability analysis and gain some insight from it. We nevertheless keep in mind that the results may somewhat differ from what is found in the simulations since we are not exactly in the asymptotic limit of interest. Three other assumptions of this stability analysis are not strictly fulfilled in our simulations, namely the ideal MHD, the purely toroidal field, and the time-independence of the background state.

Fig. 6 Maximal growth rate of the instability for the magnetic field and the differential rotation of the case with Re = 2 × 104 and Lup = 60 as a function of Ω0/ωAφ0 for different values of m (left) and as a function of the wavenumber m for Ω0/ωAφ0 = 5 (right).

4.3.1. Local linear stability results

We first adapt the local dispersion relation derived by Ogilvie (2007) in a general case to our situation where the flow is incompressible and without gravity. We are then left with the two possible instabilities of interest: the TI and the MRI. As mentioned above, the steady and axisymmetric basic state consists in a differential rotation profile Ω₀Ω(s,z) where Ω₀ is a measure of the rotation rate at the surface and a purely toroidal magnetic field B₀B(s,z) where B₀ is the maximum field strength. The dispersion relation of Ogilvie (2007) is non-dimensionalised using d for the length scale and B₀ and Ω₀ for the typical scales of magnetic field and rotation rate, respectively. This produces a unique dimensionless parameter Ω₀/ω_Aφ0 where the Alfvén frequency ${\mathit{\omega }}_{{\mathrm{A}}_{\mathit{\phi }\mathrm{0}}}\mathrm{=}{\mathit{B}}_{\mathrm{0}}\mathit{/}\mathit{d}\sqrt{{\mathit{\mu }}_{\mathrm{0}}{\mathit{\rho }}_{\mathrm{0}}}$. It is quite noticeable here that the influence of the wavevector k will only enter the dispersion relation through its direction, and not its norm. We thus only need to choose a particular direction for the perturbation, the dispersion relation being independent of the poloidal wavenumber itself.

If we consider the perturbation such that k is aligned with the rotation axis (i.e. a displacement in the direction e_s perpendicular to the axis since the perturbation is transverse in this analysis), we get the following form of the dimensionless dispersion relation, $\begin{array}{ccc}\left[{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}\frac{{\mathit{m}}^{\mathrm{2}}{\mathit{B}}^{\mathrm{2}}}{{\mathit{s}}^{\mathrm{2}}}\mathrm{-}\mathrm{2}{\left(\frac{{\mathrm{\Omega }}_{\mathrm{0}}}{{\mathit{\omega }}_{{\mathrm{A}}_{\mathit{\phi }\mathrm{0}}}}\right)}^{\mathrm{2}}\hspace{0.17em}\mathit{s}\hspace{0.17em}\mathrm{\Omega }\hspace{0.17em}{e}\mathit{s}\mathrm{·}\mathrm{\nabla }\mathrm{\Omega }\mathrm{+}\mathrm{2}\hspace{0.17em}\mathit{B}\hspace{0.17em}{e}\mathit{s}\mathrm{·}\mathrm{\nabla }\left(\frac{\mathit{B}}{\mathit{s}}\right)\right]& & \\ \mathrm{\times }\left[{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}\frac{{\mathit{m}}^{\mathrm{2}}{\mathit{B}}^{\mathrm{2}}}{{\mathit{s}}^{\mathrm{2}}}\right]\mathrm{=}{\left[\mathrm{2}\hspace{0.17em}\left(\frac{{\mathrm{\Omega }}_{\mathrm{0}}}{{\mathit{\omega }}_{{\mathrm{A}}_{\mathit{\phi }\mathrm{0}}}}\right)\hspace{0.17em}\mathit{\omega }\hspace{0.17em}\mathrm{\Omega }\mathrm{+}\frac{\mathrm{2}\mathit{m}{\mathit{B}}^{\mathrm{2}}}{{\mathit{s}}^{\mathrm{2}}}\right]}^{\mathrm{2}}& & \mathit{,}\end{array}$(12)where s is the cylindrical radius and m the azimuthal wavenumber considered. The ratio Ω₀/ω_Aφ0 will be varied in our analysis to investigate the effect of rotation and differential rotation on the possible instabilities and then fixed to the value of Ω₀/ω_Aφ0 ≈ 5 found in our simulations when the instability starts to appear to compare with the results of the 3D calculations. The profiles of magnetic field and of differential rotation are arbitrary in this formulation, we can thus choose to focus on the profiles we get from the simulations when the instability starts to kick in. This dispersion relation is an algebraic quartic equation in ω and is thus easily solved numerically. In all cases, the perturbation with k parallel to the rotation axis was found to be the most unstable. We thus show results only for this particular direction of the poloidal wavevector.

The first results of this local analysis are obtained by choosing the toroidal magnetic field configuration and the differential rotation we get in the 3D simulations in the case where Lu_p = 60 and Re = 2 × 10⁴, at approximately the time when the instability starts to grow, i.e. at t = 0.1t_ap (see Fig. 5). At this time, the maximum ratio between toroidal and poloidal magnetic field is already equal to approximately 30. We may thus assume that the effect of the poloidal magnetic field on the stability of the main concentration of toroidal magnetic field can be ignored in the analysis. As stated before, at the time when the instability starts to develop in the 3D simulation, the ratio Ω /ω_Aφ at the maximum of toroidal field is close to the value of 5.

Figure 6 shows the maximal growth rates of the magnetic instabilities in our configuration, as a function of the wavenumber m and of the ratio Ω₀/ω_Aφ0. The left panel represents the maximal growth rate (in units of the toroidal Alfvén frequency ω_Aφ0) as a function of Ω₀/ω_Aφ0 for different wavenumbers (different colours) and the right panel illustrates the dependence of the maximal growth rate on the azimuthal wavenumber for the particular value of Ω₀/ω_Aφ0 = 5 obtained in the 3D simulation. As we can see on the left panel, our configuration is always subject to a magnetic instability for a certain m independently of the value of Ω₀/ω_Aφ0. However, the nature of the instability changes as this ratio is varied. Indeed, when Ω₀/ω_Aφ0< 1, the dominant instability favours the m = 1 mode with a growth rate σ/ω_Aφ0 ≈ 2.26. In this regime, the m = 2 and m = 3 modes are also unstable and the growth rate quickly drops to 0 and below so that there is no instability for m ≥ 5. In this regime where the magnetic field dominates, the main instability is the Tayler instability which favours the m = 1 mode and requires strong gradients of field to develop.

Fig. 7 Azimuthally-averaged toroidal field and meridional cut of the fluctuating radial field during the growing phase of the instability for cases where Re = 2 × 104 and Lup = 60 (left panels) and with Re = 4 × 104 and Lup = 120 (right panels). In the northern hemisphere and for each snapshot, the contours of the growth rate for m = 5 and Ω0/ωAφ0 = 5 predicted by the linear analysis are superimposed to the colour contours.

At the other limit where Ω₀/ω_Aφ0> 4, the favoured mode always possesses a wavenumber equal to Ω₀/ω_Aφ0. This is true for m = 4,5 and 8 but also true for the other modes not represented here. From the figure, we also find that the growth rate of this most unstable mode is approximately equal to Ω₀/ 2ω_Aφ0. The oblique dashed line represented on the figure is the line of equation y = x/ 2, showing the maximum growth rate which can be reached for the favoured mode. This regime is consistent with the magneto-rotational instability regime.

Indeed, using a local linear analysis, Masada et al. (2006) estimate the maximal growth rate of the non-axisymmetric MRI to be σ_max = q Ω₀/ 2 for ${\mathit{m}}_{\max }\mathrm{=}\sqrt{\mathrm{-}{\mathit{q}}^{\mathrm{2}}\mathrm{+}\mathrm{4}\mathit{q}}\hspace{0.17em}{\mathrm{\Omega }}_{\mathrm{0}}\mathit{/}\mathrm{2}\hspace{0.17em}{\mathit{\omega }}_{{\mathrm{A}}_{\mathit{\phi }\mathrm{0}}}$, where we recall that the shear parameter q is defined as q = | ∂lnΩ /∂lns |. In our case where q ≈ 1, we thus indeed expect a maximum growth rate of about σ_max = Ω₀/ 2 and a m_max close to Ω₀/ω_Aφ0 (as also found in the global analysis of Ogilvie & Pringle 1996). To confirm that the instability found in this regime is the MRI, we artificially changed the direction of the gradient of Ω to make the rotation rate increase with the cylindrical radius s instead of decrease with s as we have in the simulations. In this case, the MRI disappears: all the modes are indeed stabilised for Ω₀/ω_Aφ0> 4.

For the intermediate regime where 1 < Ω₀/ω_Aφ0< 4, both instabilities coexist, with similar growth rates of the order of the Alfvén frequency. We note that in this regime, if the gradient of Ω is artificially modified to make the MRI disappear, the m = 3 mode is stabilised, the m = 2 mode is stabilised for Ω₀/ω_Aφ0> 2 and the m = 1 mode is stabilised for Ω₀/ω_Aφ0> 2.5. The TI is thus visible on the m = 1 mode for values of the rotation rate below Ω₀/ω_Aφ0 = 2.5. Previous studies (e.g. Spruit 1999) state that in the rotating case (but without differential rotation), the growth rate of the TI is reduced by approximately a factor of ω_Aφ0/ Ω compared to the non-rotating case. This is indeed what is found here in the regime Ω₀/ω_Aφ0< 2.5. However, for sufficiently high values of the rotation compared to the Alfvén speed (Ω₀/ω_Aφ0> 2.5), the TI does not seem to just have a reduced growth rate in our differentially rotating cases, it completely disappears. We find that the whole system is then stabilised when rotation dominates in our particular setup.

The right panel of Fig. 6 shows the growth rate of the instability on our configuration, as a function of the wavenumber m for Ω₀/ω_Aφ0 = 5, i.e. in the MRI regime. As expected from theoretical works on the MRI and as noted above, the growth rate increases with m until it reaches a maximum value of about q/ 2 × Ω₀. The growth rate then seems to saturate to values close to the maximum for higher values of m. This indicates that several modes are supposed to be unstable in our configuration with similar growth rates, with a maximum around m ≈ 5. However, we note that this linear analysis does not take into account the effect of diffusion which will probably act on high wavenumbers to significantly decrease their growth rates. We thus do not expect to get very high-m modes growing in our simulations.

4.3.2. Comparison to 3D results

We now return to the results of the 3D calculations and make a tentative comparison with what is found in the local linear analysis. In Fig. 5, panel b, the temporal evolution of each mode from m = 0 to m = 11 of the toroidal magnetic energy is shown. From this panel, it is rather clear that at t ≈ 0.1 − 0.15t_ap, the most unstable modes are the m = 3, m = 4 and m = 5 modes, which is consistent with what we expect at Ω /ω_Aφ ≈ 5. We can calculate the typical growth rate of those modes and normalise it with the maximum Alfvén frequency, in order to compare with the previous results. We find a maximum growth rate of about σ/ω_Aφ ≈ 0.8 in this case, as stated in Sect. 4.2. When the initial poloidal field strength is increased, the same unstable modes are found since the typical value of Ω /ω_Aφ is always around 5 at instability onset. A growth rate of σ/ω_Aφ0 ≈ 1 is then found for Lu_p = 100 and Lu_p = 200. The rotation rate was also varied in our simulations and the same instability is again found in more rapidly rotating cases. Again, for Re = 4 × 10⁴ and Lu_p = 60,120,200 and 400, the maximum growth rate is found to be around σ/ω_Aφ ≈ 1 and the instability always starts to grow when Ω /ω_Aφ ≈ 5. For all our cases where the shear parameter is fixed, we thus find a typical growth rate which only depends on the rotation rate Ω, such that: σ/ Ω ≈ 0.2.

This is in agreement with what was found in the local linear analysis where the maximal growth rate was dependent on the rotation rate and the shear parameter only: we had a maximal growth rate of σ/ Ω₀ = q/ 2 where Ω₀ is a measure of the rotation rate at the surface. Even if the dependency of the growth rate on the rotation rate is recovered in our 3D calculations, we find growth rates which are typically 2 to 3 times smaller than what we expect from the local linear analysis, since q ≈ 1 in our simulations. There could be several reasons for this discrepancy. First, Ohmic diffusion is not taken into account in the linear analysis and may affect the growth rate of our instability which possesses rather small scales in the horizontal direction. Another reason could be linked to the fact that the background magnetic field on which the instability develops is not steady. The strength of the field, as well as its spatial configuration are modified during the growth of the instability. This evolution of the equilibrium field is of course not taken into account in the linear analysis and could significantly alter the instability. In particular, since the instability develops during the kinematic phase of the evolution of the axisymmetric system, the ratio Ω /ω_Aφ decreases since the amplitude of the toroidal field grows linearly. As a consequence, the most unstable mode of azimuthal wavenumber m ≈ Ω /ω_Aφ at a given time is not the same at a later time: the m of the most unstable mode tends to decrease as time evolves. This is indeed what is seen in Fig. 5 panel b where the modes of lower azimuthal wavenumbers seem to take over (to be of highest amplitude) as the background field gets amplified until the m = 2 and then the m = 1 modes become dominant at the saturation of the instability after t ≈ 0.5t_ap.

Fig. 8 Instability of an initial dipole with Re = 2 × 104 and Lup = 1000 (to have the same initial poloidal energy as the reference case with Lup = 100). Panel a): initial velocity (colours) and magnetic field lines (contours). Panel c): volume rendering of the instability. Panels b) and d): same as Fig. 7 for Ω0/ωAφ0 = 5 and m = 4 (panel b)) and m = 7 (panel d)). Velocity is given in units of Re.

Figure 7 shows the location of the unstable regions of magnetic field in a simulation with Re = 2 × 10⁴ and Lu_p = 60 (left panel) and with Re = 4 × 10⁴ and Lu_p = 120 (right panel). The azimuthally-averaged toroidal field and a meridional cut of the fluctuating radial field during the growing phase of the instability are shown, together with isocontours of the growth rate of the instability predicted by the linear theory for m = 5 and Ω₀/ω_Aφ0 = 5. In both cases, the location of the instability in the 3D simulation agrees quite well with the location of the maximum growth rate predicted by the linear stability analysis. The instability is confined in a thin layer close to the surface in these cases, where the toroidal field also concentrates. We note that the toroidal field gets more squeezed against the top boundary in the right panel case where the diffusive effects are lower (the viscosity and magnetic diffusivity were reduced by a factor of 2 compared to the case on the left panel), showing that some effects of diffusion are still visible on the axisymmetric evolution. Since the diffusion effects are decreased, smaller scales are able to exist, as can also be seen on the meridional structure of the fluctuating part of the radial field. In both cases anyway, it is the MRI which is found to be the dominant instability for our magnetic structure.

4.4. Other poloidal field configurations

Other initial poloidal fields were also tested. A classical dipolar field (with the bottom boundary initial mismatch) and a confined dipole were introduced with a similar level of magnetic energy compared to the previous case with Lu_p = 100 where the instability was clearly visible and where the effects of dissipation were limited. The expressions of B_r and B_θ are given below for the dipolar field and the confined dipole (respectively Eqs. (13) and (14)):

Fig. 9 Instability of an initial confined dipole with Re = 2 × 104 and Lup = 300 (to have the same initial poloidal energy as the reference case with Lup = 100). Same as Fig. 8 except that panel d) shows contours of the growth rate of the m = 8 mode.

The classical formula for the dipole is used, normalised at the poles at the bottom boundary in this case: $\begin{array}{ccc}{\mathit{B}}_{\mathit{r}}& \mathrm{=}& \cos \mathit{\theta }\hspace{0.17em}\frac{{\mathit{r}}_{\mathit{i}}^{\mathrm{3}}}{{\mathit{r}}^{\mathrm{3}}}\mathit{,}\\ {\mathit{B}}_{\mathit{\theta }}& \mathrm{=}& \end{array}$(13)For the confined dipole, the field vanishes at the top and bottom boundaries, the field is normalised at the poles at r = 2r_ir_o/ (r_i + r_o): $\begin{array}{ccc}{\mathit{B}}_{\mathit{r}}& \mathrm{=}& \mathrm{16}\cos \mathit{\theta }\hspace{0.17em}\frac{{\mathit{r}}_{\mathit{i}}^{\mathrm{2}}\hspace{0.17em}{\mathit{r}}_{\mathrm{o}}^{\mathrm{2}}\hspace{0.17em}\mathrm{(}\mathit{r}\mathrm{-}{\mathit{r}}_{\mathit{i}}{\mathrm{)}}^{\mathrm{2}}\hspace{0.17em}\mathrm{(}{\mathit{r}}_{\mathrm{o}}\mathrm{-}\mathit{r}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{4}}}\mathit{,}\\ {\mathit{B}}_{\mathit{\theta }}& \mathrm{=}& \end{array}$(14)Similar instabilities are found in those cases but show different characteristics. Similarly to the previous case, we compare the results provided by the local linear analysis with the characteristics of the instability found in the global 3D calculations. Figures 8 and 9 illustrate the results for both initial poloidal configurations. The shearing of the initial dipole by the differential rotation produces two lobes of toroidal field of opposite polarity, as in the previous case but much less confined in a thin layer close to the surface. Again, because of the initial poloidal field and differential rotation geometries, the maximum Ω-effect is now located at mid-radius, close to the equator where we indeed find a maximum of toroidal field. This magnetic configuration is unstable to a magnetic instability which starts to develop at a time where Ω /ω_Aφ ≈ 5. This is quite remarkable that the typical value of Ω /ω_Aφ for which the instability is triggered is the same for this different initial poloidal field configuration. We note that this value is not reached at the same time t/t_ap as with the configuration of the previous case though, since this time will depend on the initial field geometry, as stated in Sect. 4.1. The most unstable mode is a m = 4 mode which typically grows on a toroidal Alfvén crossing time. Interestingly, the field seems to be unstable in two distinct regions, one located around the maximum of toroidal field inside the domain and the other one located close to the surface where strong radial gradients of B_φ exist, as seen on the right panels of Fig. 8. We note that the poloidal wavevector k tends to be parallel to e_θ for the unstable modes located close to the surface and aligned with the rotation axis for the unstable modes of the interior. In this situation, the linear analysis predicts a maximum growth rate for the azimuthal wavenumber m = 7 but the growth rate has almost already reached its maximal value of σ/ω_Aφ ≈ 2.5 for m = 4−5. As seen in Fig. 8, the location of the maximum growth rate differs significantly when the m = 4 or m = 7 modes are considered, contrary to the case with the reference poloidal field configuration where all unstable modes were approximately located in the same region of strong toroidal field (see Fig. 7). We indeed find here that higher azimuthal wavenumbers m tend to be unstable close to the surface where the radial gradient of B_φ is strong. On the contrary, lower values of m concentrate inside the domain where the latitudinal gradient dominates. Two regions of instability thus also appear in the local analysis, the higher m modes being more sensitive to the strong near-surface radial gradient of B_φ. However, it is less clear in the full 3D calculations if the unstable modes located close to the surface possess higher wavenumbers. Indeed, the m = 4 mode seems to dominate almost everywhere in the unstable regions.

Fig. 10 Top panel: same as Fig. 4 but for an unperturbed case which remained stable (dashed lines) and a perturbed case where the instability was found (solid lines). Bottom panels: azimuthally-averaged rotation rate for the unperturbed and perturbed cases at time t = 0.75tap (left) and t = 1.5tap (right).

Figure 9 shows the results of both the 3D calculations and the local linear analysis for the confined dipole. Again, a magnetic instability typically grows in the 3D simulations on a toroidal Alfvén time, with preferred m = 1,2 and 3 modes, i.e. lower than in the two previous cases. This can be seen on the volume rendering where the isocontours of the fluctuating part of the radial field are shown and look quite axisymmetric. The unstable modes are in this case located at low latitudes where the gradient of toroidal field with respect to the cylindrical radius s is the strongest. In this case, the linear analysis predicts a growth rate which quickly reaches a value around σ/ω_Aφ ≈ 2.5 for m ≈ 5 but which keeps on slightly increasing for higher values of m. However, higher m modes tend to be unstable on thinner and thinner regions around the negative gradient of | B_φ | with respect to s, as can be seen on panel d of Fig. 9 where the contours of the growth rate for m = 8 are represented. These very thin regions where high wavenumbers are supposed to be the most unstable may have reduced growth rates in the 3D calculations since they may be subject to Ohmic diffusion which is not taken into account in the local linear analysis. Moreover, in this case, the instability starts to develop when the maximum ratio between toroidal and poloidal fields is only around B_φ/B_p ≈ 8. As a consequence, we may not be able to neglect the effect of the poloidal component on the development of the instability of the magnetic field. In particular, additional magnetic tension could act to stabilise the high-wavenumber modes and thus cause the low m modes to be favoured compared to the purely toroidal case.

In both initial configurations, we confirm that the instability is of magneto-rotational type, since it vanishes for a positive gradient of Ω with respect to s. As in the reference case, those configurations are also found to be unstable to the TI when ω_Aφ dominates over Ω, which is not the case when the instability sets in in the simulation. In addition, the TI unstable modes are located on the regions where | B_φ | increases with s, i.e. much deeper inside the spherical shell than the unstable regions shown here. In the magnetic configurations studied in this work, we have shown that the MRI dominates over the TI. One of the key parameters which puts the configuration in a MRI-dominated regime is the Alfvén frequency to rotation rate ratio. In all cases shown here where an initial poloidal field is sheared by a cylindrical differential rotation, the induced toroidal field is always dominated by rotation. As a consequence, the growth rate of the MRI is always higher than the growth rate of the TI, independently of the initial field configuration. We suspect however that a different profile for the differential rotation, a density gradient through the domain or a stable stratification may have an impact on this balance.

4.5. Effect on the transport of angular momentum

We focus in this section on the transport of angular momentum induced by the development of the instability. It is particularly interesting to focus on the effect of MHD processes on the degree of differential rotation they allow. Indeed, the stellar magnetic field has often been invoked as a good candidate for forcing a flat rotation profile in a stellar radiative envelope where other transport processes such as meridional circulation or turbulence associated with hydrodynamical instabilities were not efficient enough to reduce the differential rotation. We expect here a form of turbulence associated with the magnetic instability and we may get an efficient transport of angular momentum which we intend to measure in our 3D calculations.

In Figure 10 the temporal evolution of the toroidal magnetic and kinetic energies are represented on the top panel, in a case with Re = 2 × 10⁴ and Lu_p = 60 when a non-axisymmetric perturbation was applied to the system (solid lines) and without perturbation (dashed lines). These two cases enable to quantify the effect of the instability on the redistribution of angular momentum and on the evolution of the wound-up magnetic field, as compared to the purely axisymmetric evolution. Between t = 0 and t ≈ 0.5t_ap, the evolution is rather similar between the two cases, the toroidal field takes its energy from the differential rotation and the magnetic energy thus increases while the kinetic energy decreases at the same rate. At t ≈ 0.555t_ap, a maximum of toroidal field is reached in the unstable case, while it is reached at t = 0.674t_ap in the axisymmetric unperturbed case. As a consequence, the maximum magnetic energy is larger by about 8% compared to the case where the instability developed.

If we look at the temporal evolution of the magnetic energy in the various azimuthal wavenumbers m of Fig. 5, we see that t = 0.555t_ap approximately corresponds to the time at which the instability non-linearly saturates and where the energy contained in the non-axisymmetric modes becomes comparable to the energy of the m = 0 mode. The maximum value of the toroidal field here is thus limited by the non-linear interactions between modes occurring because of the development of the instability. We saw that in the axisymmetric case, the subsequent evolution of the magnetic field consists in damped oscillations with periods of about an Alfvén time. In the unstable case, the oscillations are also observed but their amplitude is significantly reduced and the final state of zero toroidal field is reached much faster.

Concerning the kinetic energy, a similar behaviour is found. In the unstable case, by the time the toroidal magnetic energy reaches its maximum, the kinetic energy has already almost reached its minimal value. The oscillations due to travelling Alfvén waves are also of much smaller amplitude than in the unperturbed case. As we said before, the final state consists of a solid body rotation. In the case where the instability developed, we can almost argue that this state of solid body rotation is reached on a timescale much closer to the Alfvén timescale than to the diffusive timescale, i.e. 60 times faster for this case with Lu_p = 60.

On the bottom panels of Fig. 10, the averaged rotation rate in a meridional plane is shown for 2 different times t = 0.75t_ap (left) and t = 1.5t_ap (right) in the perturbed and unperturbed cases. The colour table is different for the two times but kept the same at each time for the stable and unstable cases. At time t = 0.75t_ap when the kinetic energy has already reached its minimal value in the perturbed case and is still decreasing in the unperturbed case, the shell in the unperturbed case is still differentially rotating when a state of solid body rotation is already almost reached in the other case. For example, at the equator, the bottom boundary rotates 33% faster than the top in the unperturbed case at t = 0.75t_ap. For the case where the instability acted to redistribute the angular momentum, this differential rotation is reduced to 0.01% at the equator. Higher or lower values of the rotation rate still exist in this case in a thin layer close to the rotation axis where the toroidal field vanishes and where the instability was probably inefficient at transporting angular momentum since the unstable modes are primarily located around the strongest concentrations of B_φ, as seen in the previous section (see Fig. 7). A very efficient transport of angular momentum has thus occurred in the unstable case in the bulk of the domain due to the strong Maxwell and Reynolds stresses induced by the instability. At time t = 1.5t_ap represented on the bottom right panels of Fig. 10, a change of direction of the gradient of Ω is clearly visible in both cases but again a strong reduction of the differential rotation has taken place in the unstable case. The top boundary now rotates 36% faster than the bottom in the unperturbed case and only 12% faster in the other case. We thus find here a typical decay time for our differential rotation of the order of the poloidal Alfvén time t_decay ≈ t_ap. We have yet not intended to perform a detailed parameter study to determine a scaling law of this decay rate with the various simulation parameters.

In a recent work, Rüdiger et al. (2015) discusses the transport properties of an unstable toroidal magnetic field subject to the so-called azimuthal magneto-rotational instability (AMRI). For a rotation profile proportional to 1 /s and 1 /s² and in cylindrical geometry, they find that the turbulent eddy viscosity ν_t associated with the AMRI scales as the rotation rate, but with a strong dependence on Pm. Their value of the turbulent viscosity is in fact a maximum value since for each Reynolds number, they keep only the Lundquist number which maximises ν_t, thus leading to a scaling independent of the magnetic field strength. In their particular setup and for Pm = 1, they find a minimal decay time of their differential rotation of t_decay ≈ 200 t_Ω, the coefficient of proportionality 200 being established numerically. The same kind of parametric study was not performed here but in typical cases, we find a decay time of the order of the poloidal Alfvén time t_decay ≈ t_ap. In our notations, we have t_ap = (Re/Lu_p)t_Ω where the ratio Re/Lu_p was varied between 50 to 1000 in our simulations. We thus find a decay time which is indeed of the same order as that found by Rüdiger et al. (2015). Similar values were also established by Arlt et al. (2003) with again another setup of a purely axial field becoming unstable with respect to the MRI. In any case, this typical time to flatten the rotation profile is much faster than the magnetic diffusion time at the scale of the radius of the star. More importantly, the decay time in the unstable cases is also found to be faster than the phase-mixing time which is a diffusion time associated with the small scales created by oscillations on neighbouring magnetic surfaces getting out of phase. We thus get to the conclusion that this particular instability provides an interesting additional mechanism to transport angular momentum in stellar radiative zones.

Fig. 11 Temporal evolution of the poloidal magnetic energy contained in the first 12 azimuthal modes for Re = 2 × 104 with Lup = 60 (left), Lup = 200 (middle) and Lup = 400 (right), corresponding to Lo = 3 × 10-3,10-2 and 2 × 10-2, respectively. In each panel, the vertical dash-dotted line represents the time at which the maximum value of the toroidal magnetic field is reached in the domain.

5. Possible application to A-type stars

One motivation for this type of studies is the application to the magnetic fields observed at the surface of intermediate-mass stars, possessing a radiative envelope. As said in the introduction, spectropolarimetric studies have lead to a distinction between two populations in those stars. The first population are the Ap stars which show a strong dipolar-like magnetic field (above 300G). The other population is formed by the non-Ap stars which do not exhibit strong magnetic fields. The following scenario was proposed by Aurière et al. (2007) to explain this dichotomy. If we consider an initial poloidal field and differential rotation in a radiative envelope, two possibilities would then occur:

• The initial poloidal field is strong compared to the initialdifferential rotation, i.e. the winding-up time is long compared tothe typical time of the back-reaction on the shear. In this case, thedifferential rotation is quickly quenched by the magnetic field andthe configuration relaxes to a stable equilibrium. This constitutesthe magnetic field observed in Ap stars.

• The initial poloidal field is weak compared to the initial differential rotation, i.e. the winding-up time is short compared to the back-reaction timescale. In this case, a strong toroidal field is built in a differentially rotating star. A magnetic instability of some type develops and creates small scales in the horizontal direction. Opposite polarities are created on the stellar surface and produce a weak line-of-sight-averaged field. The stars that underwent this instability would constitute the observed sub-Gauss magnetism.

Aurière et al. (2007) estimated the winding-up time to be approximately 1 /qΩ and the back-reaction timescale equal to l/v_ap where v_ap is the poloidal Alfvén speed and where the shear length scale $\mathit{l}\mathrm{=}\mathrm{(}\frac{\mathit{\partial }\ln \mathrm{\Omega }}{\mathit{\partial s}}{\mathrm{)}}^{-1}\mathrm{=}\mathit{s}\mathit{/}\mathit{q}$. The ratio of the two timescales, estimated at the stellar surface, will then be exactly the ratio t_Ω/t_ap. The first case above thus corresponds to t_ap ≪ t_Ω and the second case to t_ap ≫ t_Ω. At the limit t_ap = t_Ω, a critical field strength can be calculated: ${\mathit{B}}_{{\mathrm{0}}_{\mathrm{crit}}}\mathrm{=}{\mathrm{\Omega }}_{\mathrm{0}}\mathit{d}\sqrt{{\mathit{\rho }}_{\mathrm{0}}{\mathit{\mu }}_{\mathrm{0}}}$. For a typical main sequence Ap star (with a period of 5 days, a radius of 3 R_⊙ and a temperature T_eff = 10 000 K), a critical field of around 300 G is found, in agreement with the observational constraints.

5.1. A critical number to distinguish between stable and unstable cases

In this work, we are exactly in the right setup to test this scenario by varying the dimensionless number Lu_p/ Re = t_Ω/t_ap in our simulations (where Pm = 1), starting with the initial poloidal field and the cylindrical differential rotation described in Sect. 2.3. We note here that the parameter Lu_p/ Re measuring the ratio of the rotation timescale t_Ω to the Alfvén timescale t_ap was introduced in the geophysical community as the Lorentz number Lo (e.g. Christensen & Aubert 2006). We shall use this notation now to ease the description of our calculations. We shall keep in mind that those calculations are performed assuming an incompressible fluid in an unstratified atmosphere, which is obviously not realistic for a stellar radiative zone. The same kind of study in a stably stratified spherical shell will be the subject of future works. However, the interaction of the magnetic field with the differential rotation and the development of a magnetic instability are processes which are well taken into account in this work. We thus wish here to take a first step towards testing the above scenario in a more realistic situation.

From Sect. 4.3.2, the growth rate of the instability in our numerical simulations is equal to σ = 0.2Ω ≈ 0.1Ω₀ since in all cases, the rotation rate at the location of the instability Ω is typically half the value of the rotation rate on the axis Ω₀. This gives a typical timescale for the growth of the instability of t_inst = 1 /σ = 10/Ω₀ = 10 t_Ω. We also saw that the lifetime of the toroidal field configuration was approximately the Alfvén crossing time based on the poloidal field t_ap. Indeed, after t_ap, regions of opposite polarities for the toroidal field are created by the modified gradient of differential rotation. As a consequence, the instability is able to develop and reach the level of the axisymmetric mode only if t_inst/t_ap = 10 t_Ω/t_ap ≪ 1, i.e. t_Ω/t_ap = Lo ≪ 0.1.

In Fig. 11, the temporal evolution of the poloidal energy contained in the first azimuthal wavenumbers is shown for 3 different cases with Lo = 3 × 10^-3,10^-2 and 2 × 10^-2. It is rather clear from the figure that the ratio between the dynamical timescales in our system, has a very strong impact on the level at which the instability saturates. In particular, for a sufficiently small value of Lo (for example for Lo = 3 × 10^-3), the instability has enough time to grow and reach the level of energy of the axisymmetric mode before the toroidal field gets to its maximal value. In fact, it is even the instability which saturates the field and prevents the toroidal field from reaching its highest possible value (obtained without perturbing the system, see Fig.10). On the contrary, when Lo is increased, the typical growth time for the instability becomes comparable to the lifetime of the background field on which it grows. We find here that even when the timescale of the instability is 10 times smaller than the poloidal Alfvén crossing time (t_Ω/t_ap = 10^-2), the level of energy of the axisymmetric mode is not reached by the non-axisymmetric modes and the configuration will remain mostly axisymmetric. Similar cases were calculated with a rotation rate twice as strong and thus a differential rotation also twice as strong, keeping the same ratio between the dynamical timescales. It is also found in this case that t_Ω/t_ap = 10^-2 seems to be the critical value below which the axisymmetric structure will be severely altered by the instability and above which the system will relax into a quasi-axisymmetric equilibrium, which may not be that of the initial purely poloidal field as we will discuss in the next section.

5.2. How does the surface magnetic field look like?

Now that we have seen that the ratio t_Ω/t_ap was crucial to determining the level of saturation of the magnetic instability found in our simulations, we focus on the resulting field which would be observed at the stellar surface for an initial value of t_Ω/t_ap above and below the critical value quoted previously. We are thus here interested in the structure of the surface radial field obtained in a particular setup where a poloidal field of arbitrary strength initially coexist with a cylindrical differential rotation of arbitrary strength in the stellar interior. We expect a dramatic change in this radial field structure between cases where the instability fully develops and cases where the instability saturates at a low value.

Fig. 12 Radial magnetic field close to the surface at time t = tap for the first 2 cases shown in the previous figure for Re = 2 × 104 (top panels) and for cases at Re = 4 × 104 (bottom panels) and with Lo = 3 × 10-3 (left) and Lo = 10-2 (right). The left panels correspond to the cases where the instability reaches the level of energy of the axisymmetric mode.

Figure 12 represents the radial component of the magnetic field close to the top of our computational domain for 4 different cases. The 2 left panels correspond to cases with Lo = 3 × 10^-3 and the 2 right panels to cases with Lo = 10^-2. The top panels were computed with Re = 2 × 10⁴ and the bottom ones with Re = 4 × 10⁴. This figure shows that the magnetic configuration dramatically differs after one Alfvén time when Lo is increased from 3 × 10^-3 to 10^-2. As expected, for a low enough value, the non-axisymmetric magnetic instability rapidly develops during the winding-up phase of the simulation. For both initial Re shown here, the resulting surface field after an Alfvén time is organised at rather small scales, strongly non-axisymmetric and the initial poloidal configuration has been completely modified. Since regions of opposite polarities exist in each hemisphere, the averaged intensity measured along the line of sight would be significantly smaller than the initial poloidal field strength introduced. On the contrary, when Lo is increased, the MRI still grows but saturates orders of magnitude below the axisymmetric mode. The resulting surface field at t = t_ap is predominantly axisymmetric and possesses strong polar caps reminiscent of the initial poloidal field considered. However, we note that this surface field does not match the initial poloidal field configuration. Indeed, after an Alfvén time, the field has already relaxed to a mixed toroidal/poloidal structure characterised by one twisted torus per hemisphere close to the poles. This can be seen in the figure by the belts of B_r of opposite sign lying at high latitudes for both Reynolds number cases. This mixed poloidal/toroidal form is a classical stable equilibrium found both analytically (Wright 1973) or semi-analytically (Duez & Mathis 2010) and numerically (Braithwaite & Nordlund 2006).

As stated in Sect. 4.1, the ratio B_φ/B_p is proportional to d${\mathrm{\Omega }}_{\mathrm{0}}\sqrt{{\mathit{\mu }}_{\mathrm{0}}{\mathit{\rho }}_{\mathrm{0}}}\mathit{/}{\mathit{B}}_{\mathrm{0}}\mathrm{=}{\mathit{t}}_{\mathrm{ap}}\mathit{/}{\mathit{t}}_{\mathrm{\Omega }}$ in the kinematic phase of the evolution, during which the instability grows. As a consequence, when Lo = t_Ω/t_ap increases, the ratio between toroidal and poloidal fields decreases. For kink-type instabilities such as the TI, the poloidal field component stabilises the azimuthal field and the vertical wavelength tends to be proportional to the ratio B_p/B_φ (e.g. Linton et al. 1996). We could have situations (like in a stratified zone) where the vertical wavelength of possible instabilities is limited to small values (compared to the horizontal scale). In this case, for a fixed toroidal field, the poloidal field is thus stabilising both because the field line pitch decreases (i.e. B_p/B_φ too large) and/or because the wavelength of the possible instability is not allowed to exist in the domain. However, for the MRI found here, it is yet unclear how stabilising would a large poloidal field component be. If we assume that poloidal fields also have a stabilising effect on the MRI observed here, the higher Lo cases could also experiment a reduction in the growth rate of the instability because of a stronger impact of the poloidal field. While the ratio B_φ/B_p gets closer to 1, the field thus tends to quickly relax into a mixed poloidal/toroidal configuration. An example of such a configuration is shown in Fig. 13 where magnetic field lines with diameters proportional to the total magnetic energy are represented. This snapshot corresponds to the bottom right panel of Fig. 12 where hints for a twisted torus close to the poles were observed. On the 3D representation of the field lines, it is rather clear that twisted structures have formed in both hemispheres, although the toroidal component is still quite dominant. As stated before, if the field strength of our initial poloidal field was further increased, the ratio between toroidal and poloidal fields would decrease and the twist of the field lines (and consequently their pitch angle) would increase, producing twisted tori with a stronger poloidal component. It is worth noting here that the internal magnetic field configuration in the cases where the instability fully developed is in fact also rather similar to one twisted torus per hemisphere, even if the toroidal component is stronger compared to the poloidal one, for the reasons just stated above. As a consequence, even if the observed surface radial field appears drastically different between the stable and unstable cases as shown in Fig. 12, the “hidden” internal configuration may always relax after a few Alfvén times into a mixed poloidal/toroidal configuration similar to what is shown in Fig. 13.

Fig. 13 Magnetic field lines representation for the case at Re = 4 × 104 and Lo = 10-2 at t = tap. The diameter of the field lines is proportional to the magnetic energy to emphasise the structure of the magnetic field where it is maximum.

As already stated in Sect. 4.5, the turbulence induced by the magnetic instability could produce an additional transport of angular momentum which causes the final state of solid-body rotation to be reached quicker in the perturbed cases than in the axisymmetric ones. Here, the instability develops in all cases since they are all perturbed, but we have cases where the instability has fully developed to reach the level of the axisymmetric modes and other cases where the growth time of the instability was too long to reach a significant energy. By studying the differential rotation profile in each of these cases, we find that the angular momentum was efficiently transported again only in the case where the instability fully developed. A state of solid-body rotation is then quickly reached only when the quantity t_Ω/t_ap remains relatively small.

In agreement with the scenario at the beginning of this section, the ratio t_Ω/t_ap seems to be crucial to distinguish between cases exhibiting a large-scale axisymmetric surface magnetic field and cases where a much more complex non-axisymmetric structure dominates. A critical value of this ratio around 10^-2 is found to separate the two populations. We thus find that for our initial configuration and Pm = 1, the critical poloidal field below which the non-axisymmetric modes dominate is: $\begin{array}{ccc}{\mathit{B}}_{{\mathrm{0}}_{\mathrm{crit}}}\mathrm{=}{\mathrm{10}}^{-2}{\mathrm{\Omega }}_{\mathrm{0}}\mathit{d}\sqrt{{\mathit{\rho }}_{\mathrm{0}}{\mathit{\mu }}_{\mathrm{0}}}\mathit{.}& & \end{array}$This would indicate that the measured large-scale field in Ap stars should be at least equal to this value and that the minimum field observed in those stars should therefore be proportional to the rotation rate and to the square root of the surface density. The dependence on rotation of the minimum field in Ap stars has indeed already been observed (Lignières et al. 2014) and is still being investigated through detailed spectropolarimetric surveys of additional target stars. However, we have to note here that all the calculations were computed using a constant density, without entropy stratification and with a fixed profile for the differential rotation. The results may differ in particular on the nature of the instability and on the critical value of Lo. Additional progress is needed to improve the realism of our simulations in order to compare with spectropolarimetric observations, such efforts are underway.

6. Summary and conclusions

In this work, the interaction between differential rotation and magnetic fields in a simplified model of a stellar radiative zone has been investigated. Starting from an initial poloidal field and a cylindrical differential rotation, the magnetic field is first wound up to produce a strong toroidal component which backreacts on the differential rotation after a time close to an Alfvén crossing time. During this winding-up phase, a non-axisymmetric magnetic instability develops, which may in some cases completely modify the axisymmetric structure and which constitutes an efficient additional transport of angular momentum. This instability of a wound-up field has been thought to exist and discussed in the literature for a long time (e.g Spruit 1999) and has recently been studied numerically (Braithwaite 2006; Arlt & Rüdiger 2011b,a; Szklarski & Arlt 2013). However, it has very often been interpreted as a kink-type Tayler instability (TI). In this work and in our particular setup, we show that the instability is in fact the magneto-rotational instability (MRI) of a toroidal field in the presence of a differential rotation. In the context of stellar physics, there were only a few studies (e.g. Arlt et al. 2003; Masada et al. 2006) dedicated to the MRI. In Arlt et al. (2003) and Rüdiger et al. (2014, 2015), a particular focus on the transport of angular momentum induced by the MRI was made. They showed that Maxwell stresses due to magnetic field perturbations could rapidly lead to a state of solid-body rotation. Masada et al. (2006) studied the development of magnetic instabilities in proto-neutron stars. They also showed that for a sufficiently high shear parameter, the non-axisymmetric MRI would dominate over the TI. The crucial parameter which distinguishes the TI and the MRI regimes is the ratio between the rotation rate and the toroidal Alfvén frequency Ω /ω_Aφ. In a setup where a poloidal field is wound up by a cylindrical differential rotation with a shear parameter q = | ∂lnΩ /∂lns | ≈ 1 in an unstratified region, this ratio is always around 5 when the instability starts to grow. This points towards a situation falling in the MRI regime where azimuthal wavenumbers of approximately Ω /ω_Aφ are the most unstable (and not m = 1 like in the TI), with growth rates proportional to the rotation rate, as predicted by theory. The analysis of the instability is subtle since the background equilibrium field on which the instability develops is itself evolving with a characteristic timescale equal to the poloidal Alfvén time. As a consequence, the instability freely develops and strongly modifies the axisymmetric structure only when its growth rate (proportional to Ω) is large compared to the poloidal Alfvén frequency (proportional to B_p) or in other words only when a weak poloidal field is wound up in a rapidly rotating shell. In the opposite case, the field relaxes into an axisymmetric mixed poloidal/toroidal configuration.

The results suggest that, when differential rotation is present, slowly rotating radiative envelopes possessing a strong poloidal field would tend to present a stable axisymmetric surface field whereas weaker poloidal fields embedded in more rapidly rotating zones would become unstable and create non-axisymmetric surface fields with several polarity inversions in each hemisphere. This scenario could potentially explain the magnetic dichotomy between strongly magnetic slowly rotating Ap stars and the other weakly magnetised (or non-magnetic) normal A stars. Our calculations point towards a linear dependence between the minimum field strength expected in Ap stars and their rotation rate. This relation has already been partially confirmed through detailed spectropolarimetric observations (Lignières et al. 2014). In this scenario, a source of differential rotation needs to be identified. The contraction occurring during the PMS phase might provide the required forcing to maintain a significant differential rotation.

We do not exclude the development of the Tayler instability (or other magnetic instabilities) for other setups which could be relevant for stellar radiative zones. In our configurations, the MRI always dominates over the TI but this may be because of the lack of stable stratification and of an overestimate of the shear parameter to be used, which is not constrained by observations. As discussed in Masada et al. (2006), a weak N/ Ω ratio (where N is the Brunt-Väisälä frequency) together with a large shear parameter q ≈ 1 could be relevant for proto-neutron stars. In this case, the MRI could indeed be dominant. However, for parameters possibly more relevant to main-sequence stars (N/ Ω = 10³ and q< 0.1), the TI would dominate for Ω = 5ω_Aφ. This TI regime would then probably produce a different structure of the surface magnetic field. However, in the TI regime, we still expect to recover the magnetic dichotomy: weak non-axisymmetric surface fields would result from the instability of a toroidal field-dominated configuration and strong axisymmetric stable fields would exist for poloidal and toroidal components of similar strengths.

Entropy and density stratifications have not been taken into account in this work and are of course needed for a correct modelling of the objects we are interested in. The effects of the density stratification have already been discussed above. As far as entropy stratification is concerned, we expect both instabilities to be affected by a positive value of the Brunt-Väisälä frequency. It may then cause the instability found in our unstratified cases to be significantly reduced or limited to much smaller length scales although strong latitudinal gradients of differential rotation and magnetic fields could still be efficient at developing the MRI or the TI. Moreover, a stable stratification may naturally favour a differential rotation profile constant on concentric spheres instead of cylinders as was studied here (Zahn 1992). In the unstratified case, such a profile leads to a strong meridional circulation which makes the analysis of the stability of the wound-up toroidal field difficult. In the stratified case, a solution with a reduced amplitude of the meridional flow should exist, which should ease the study of different rotation profiles. We also note that we only studied a spin-up problem here, where the differential rotation was not forced by any external torque. The situation might be different if a forcing on the differential rotation was applied in the calculations, to take into account the presence of a convective core for example or a contraction/expansion of the shell. Additional numerical simulations are thus needed to efficiently tackle this problem and to see if non-axisymmetric instabilities can still survive and produce similar outcomes in more realistic situations.

Acknowledgments

T.G. is supported by the Special Priority Program 1488 PlanetMag of the German Science Foundation. 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. CPU time was provided by the HPC resources of CALMIP under the allocation 2014A-P1118. The authors wish to thank François Rincon and Gordon Ogilvie for very fruitful discussions and the anonymous referee for helpful comments.

References

All Figures

	Fig. 1 Azimuthal averages: initial rotation profile represented as isocontours of Ω (colours) and initial magnetic field represented by the poloidal field lines (dashed contours). Ω is given in units of Re.
In the text

Fig. 2 Lup = 200: temporal evolution of the magnetic energy contained in the first azimuthal wavenumbers (left), volume rendering showing isosurfaces of the fluctuating part of the radial magnetic field ${\mathit{B}}_{\mathit{r}}^{\mathrm{\prime }}$ together with magnetic field lines during the linear phase of the instability (middle) and cut of the total radial magnetic field at the top of our domain when the instability has reached a saturated state.

In the text

	Fig. 3 Volume rendering of the instability for Lup = 400 (left) and Lup = 1000 (right). Surfaces of constant ${\mathit{B}}_{\mathit{r}}^{\mathrm{\prime }}$ are represented, as well as some magnetic field lines inside the domain.
In the text

	Fig. 4 Temporal evolution of the toroidal magnetic (red) and kinetic (blue) energies for an initial poloidal field strength of Lup = 60 (solid lines), Lup = 100 (dashed lines) and Lup = 200 (dash-dotted lines). In all cases, Re = 2 × 104.
In the text

Fig. 5 Characteristics of the instability for Re = 2 × 104 and Lup = 60. Panels a)–d) are snapshots of, respectively, the azimuthally averaged toroidal field on which the poloidal field lines are superimposed, a meridional cut of ${\mathit{B}}_{\mathit{r}}^{\mathrm{\prime }}$ and a volume rendering showing surfaces of constant ${\mathit{B}}_{\mathit{r}}^{\mathrm{\prime }}$ and magnetic field lines at t = 0.33tap. Panel b) shows the temporal evolution of the toroidal magnetic energy contained in the first 12 azimuthal wavenumbers.

In the text

	Fig. 6 Maximal growth rate of the instability for the magnetic field and the differential rotation of the case with Re = 2 × 104 and Lup = 60 as a function of Ω0/ωAφ0 for different values of m (left) and as a function of the wavenumber m for Ω0/ωAφ0 = 5 (right).
In the text

Fig. 7 Azimuthally-averaged toroidal field and meridional cut of the fluctuating radial field during the growing phase of the instability for cases where Re = 2 × 104 and Lup = 60 (left panels) and with Re = 4 × 104 and Lup = 120 (right panels). In the northern hemisphere and for each snapshot, the contours of the growth rate for m = 5 and Ω0/ωAφ0 = 5 predicted by the linear analysis are superimposed to the colour contours.

In the text

Fig. 8 Instability of an initial dipole with Re = 2 × 104 and Lup = 1000 (to have the same initial poloidal energy as the reference case with Lup = 100). Panel a): initial velocity (colours) and magnetic field lines (contours). Panel c): volume rendering of the instability. Panels b) and d): same as Fig. 7 for Ω0/ωAφ0 = 5 and m = 4 (panel b)) and m = 7 (panel d)). Velocity is given in units of Re.

In the text

	Fig. 9 Instability of an initial confined dipole with Re = 2 × 104 and Lup = 300 (to have the same initial poloidal energy as the reference case with Lup = 100). Same as Fig. 8 except that panel d) shows contours of the growth rate of the m = 8 mode.
In the text

	Fig. 10 Top panel: same as Fig. 4 but for an unperturbed case which remained stable (dashed lines) and a perturbed case where the instability was found (solid lines). Bottom panels: azimuthally-averaged rotation rate for the unperturbed and perturbed cases at time t = 0.75tap (left) and t = 1.5tap (right).
In the text

	Fig. 11 Temporal evolution of the poloidal magnetic energy contained in the first 12 azimuthal modes for Re = 2 × 104 with Lup = 60 (left), Lup = 200 (middle) and Lup = 400 (right), corresponding to Lo = 3 × 10-3,10-2 and 2 × 10-2, respectively. In each panel, the vertical dash-dotted line represents the time at which the maximum value of the toroidal magnetic field is reached in the domain.
In the text

	Fig. 12 Radial magnetic field close to the surface at time t = tap for the first 2 cases shown in the previous figure for Re = 2 × 104 (top panels) and for cases at Re = 4 × 104 (bottom panels) and with Lo = 3 × 10-3 (left) and Lo = 10-2 (right). The left panels correspond to the cases where the instability reaches the level of energy of the axisymmetric mode.
In the text

	Fig. 13 Magnetic field lines representation for the case at Re = 4 × 104 and Lo = 10-2 at t = tap. The diameter of the field lines is proportional to the magnetic energy to emphasise the structure of the magnetic field where it is maximum.
In the text