© ESO, 2015

1. Introduction

Rotation affects the structure of stars, inducing fluid motions such as circulation, and can enhance the mixing of chemical elements. In stably stratified radiative interiors, however, fluid motions are largely restricted to isobaric surfaces, and strong turbulent fluid motions parallel to the isobars smooth velocity gradients in shells resulting in a shellular-like rotation (Zahn 1992). In order to smooth the radial velocity gradients, hydrodynamic processes appear not to be able to transport enough angular momentum between neighboring shells. Consequently, it is thought that some process involving magnetic fields must account for any angular momentum transport in the radial direction.

Even a relatively weak magnetic field is sufficient to couple different parts of the star and maintain a state of nearly uniform rotation. For the interior of the Sun, for example, a field of less than 1 gauss would be able to transmit the torque exerted by the solar wind through the interior (Mestel 1953). Measurements of the solar core’s rotation rate (Chaplin et al. 2001) show a very uniform rotation. This uniform rotation may be due to a magnetic field, but this field’s origin, configuration, and strength are not known. By analogy with the magnetic A stars, one might speculate that a “fossil” magnetic field could exist in the core of the Sun. Since no significant net field is seen at the surface (averaged over the solar cycle), the radial component of such a fossil would, however, have to be weak – on the order of a gauss or less. A field weaker than this will quickly wrap up into a predominantly toroidal field, under the action of the remaining differential rotation in the core. Eventually an instability will set in, limiting the growth of this toroidal field.

Near the end of the stellar lifetime, progenitors of white dwarfs and supernovae go through stages where the core contracts and spins up while the envelope expands and spins down. The degree of coupling between core and envelope by a magnetic field in this stage will affect the rotation rates of the stellar leftovers, pulsars, and white dwarfs. It is then fair to ask whether the rotation rate of the remnant is determined by the initial rotation of the progenitor, or if a secondary process must be responsible (Spruit & Phinney 1998; Spruit 1998).

Tayler (1973) showed that any purely toroidal field should be unstable in at least some part(s) of the star. He derived stability conditions that are local in the meridional plane but global in the azimuthal direction; indeed, it is the low azimuthal m modes that become unstable first. Physically, one can think of the axisymmetric mode m = 0 in the following way. A purely toroidal field B_φ = B_φ(r,θ), which must be axisymmetric in order to satisfy ∇·B = 0, can be thought of as a collection of discrete neighboring circular flux tubes. Since all stars have a very high plasma-β, magnetically-induced motions are restricted to be almost incompressible, ∇·v ≈ 0. Now, any flux tube whose volume is held constant can reduce its energy by becoming shorter and fatter, which in the case of a flux tube encircling the axis of symmetry means contracting towards the axis. Obviously there is more energy to be released if the field is stronger, so that if we have two flux tubes (of equal volume) lying at different distances from the symmetry axis, if the outer tube contains a greater magnetic flux than the inner tube, interchanging the two will release energy. In a stellar radiative zone, the stable stratification largely suppresses motion in the radial direction, v·g ≈ 0. However, to interchange flux tubes some radial motion is necessary; it is impossible for flux tubes confined to the same spherical surface to slip past one another. Given that the unstable m = 0 mode is axisymmetric, rotation, and thus Coriolis force, has no net effect on the stability conditions or growth rate, as it simply causes the annulus to rotate around the axis of symmetry. The non-axisymmetric modes are also driven by the reduction of energy by letting regions with strong magnetic field, and therefore a strong inwards-pointing curvature force, move inwards, and can be thought of in similar terms to the fluting instability in sunspots.

The growth timescale of the instability is on the order of the time taken for an Alfvén wave to travel around the star on a toroidal field line. This timescale is very short compared to the evolutionary timescale of the star, e.g. about 10 yr in a main-sequence star with a field of 1 kG. This form of instability is likely to be the first to set in as the field strength of the toroidal field is increased (Spruit 1999) but it is still uncertain what happens next. According to a scenario developed by Spruit (2002) the instability could lead to a regeneration of the poloidal field which will eventually be transformed into a toroidal field as a result of the remaining differential rotation, setting up a hydromagnetic dynamo loop. The balance between wrapping-up by differential rotation on the one hand and the destruction of the toroidal field by Tayler instability on the other determines the strength and configuration of the field. The rate at which angular momentum and chemical elements are transported through the stellar interior depends on the resulting field strength, configuration and the details of the instability. This process has been implemented as a sub-grid model in stellar evolution calculations by Woosley et al. (2002), Heger et al. (2003) and Maeder & Meynet (2003) and has been numerically reproduced in, with a very ideal configuration, simulations by Braithwaite (2006a). More recently the process has apparently been seen in simulations by Rüdiger et al. (2014)¹.

In the long-term, we would like to determine the type of magnetic field that is maintained by differential rotation in a stably-stratified star, under the action of magnetic instabilities, and to develop a quantitative theory for the transport of angular momentum and chemical elements by magnetic fields in stars. This will require extensive simulations which include differential rotation and have a spherical geometry (Boldt & Braithwaite, in prep.) In the meantime we aim to improve our understanding of the Tayler instability. Progress on this front has been made with local simulations by Braithwaite (2006b), and has also been investigated in global simulations by Lander & Jones (2011). More recently Bonanno & Urpin (2013a,c,b) have looked at various properties of the Tayler instability, in particular the effect of rotation, and find that rotation is essentially unable to suppress the instability; one of our aims in this paper is to confirm this result with a different method and to ascertain the growth rate of the instability in the presence of fast rotation.

Another good reason to study the Tayler instability is to gain a better understanding of the range of magnetostatic equilibria available in radiative and other non-convective stars and zones. Since the discovery by Babcock (1947) of a strong magnetic field in the A star 78 Vir, strong, large-scale fields are now known to exist in a subset of early-type stars as well as white dwarfs and neutron stars (see e.g. Braithwaite & Spruit 2015, for a review). These fields, known in the literature as fossil fields, are the relics of some previous stage of evolution, and have found their way into a stable equilibrium, evolving only on the very long diffusive timescale. An arbitrary initial magnetic field is known to relax in a non-convective star into such an equilibrium (Braithwaite & Spruit 2004; Braithwaite & Nordlund 2006), and there are apparently a large range of equilibria available, including both approximately axisymmetric “twisted torus” configurations as well as more complex non-axisymmetric configurations (Braithwaite 2008). Amongst the axisymmetric equilibria, it is an interesting question to ask in what ratio the poloidal and toroidal field components must be present to ensure stability, since both are unstable on their own. It seems that a relatively weak poloidal field can stabilize a toroidal field (Braithwaite 2009; Akgün et al. 2013) but the detail of how this works is still uncertain. One reason why the ratio of the two is particularly interesting is that a star with a predominantly poloidal field is oblate in shape, but if the field is predominantly toroidal, the star is prolate. Fast-spinning highly-magnetized neutron stars (the so-called millisecond magnetars) with predominantly toroidal magnetic fields might flip over until the magnetic and rotation axes are orthogonal, the damping of free precession reducing the kinetic energy by maximizing the moment of inertia about the rotation axis. Such stars would emit gravitational waves potentially detectable with the next generation of detectors.

In this paper we use numerical methods to investigate the linear properties of the Tayler instability, much of which is still uncertain. We provide a useful check on the analytical results, which may have missed some aspects of the mechanism. We use similar methods to those used by Braithwaite (2006b) to investigate the Tayler instability, but we investigate different aspects of the instability. In the next section, we look at what is known about the instability, before describing our numerical model in Sect. 3, presenting our results in Sect. 4 and summarizing in Sect. 5.

2. Stability of toroidal fields: analytic results

Tayler (1973) derived the stability conditions of purely toroidal fields in radiative stellar interiors, ignoring the effects of viscosity and of thermal and, magnetic diffusion. The occurrence of instability for the adiabatic case is independent of the field strength but only depends on the field configuration. However, when the instability is already present the growth rate depends on some measure of the field strength and the stratification. The field strength is characterized by the Alfvén frequency, which can be thought as the inverse of the timescale taken by an Alfvén wave to go around the star along a toroidal field line, defined as $\begin{array}{ccc}{\mathit{\omega }}_{\mathrm{A}}\mathrm{=}\frac{{\mathit{v}}_{\mathrm{A}}}{{\mathit{R}}_{\mathrm{\ast }}}\mathrm{=}\frac{{\mathit{B}}_{\mathit{\phi }}}{{\mathit{R}}_{\mathrm{\ast }}\sqrt{\mathrm{4}\mathit{\pi \rho }}}\mathit{,}& & \end{array}$(1)where ${\mathit{v}}_{\mathrm{A}}\mathrm{=}{\mathit{B}}_{\mathit{\phi }}\mathit{/}\sqrt{\mathrm{4}\mathit{\pi \rho }}$ is the Alfvén velocity in a toroidal field line and R_∗ is the radius of the star. The stratification is characterized by the Brunt-Väisälä or the buoyancy frequency, N, and radiative zones are always in the regime $\begin{array}{ccc}\mathit{N}\mathrm{\gg }{\mathit{\omega }}_{\mathrm{A}\mathit{,\phi }}\mathit{.}& & \end{array}$(2)This is the same as saying that the thermal energy density is much bigger than the magnetic energy density, i.e. a large plasma- β = P/ 8πB² parameter, which is true in all stars.

Following Tayler, we make the analysis easier by using cylindrical coordinates (ϖ,φ,z). We now consider an axisymmetric, azimuthal toroidal field independent of height z, B = B(ϖ)e_φ, and look at the region around the axis. The important parameter is the radial field gradient p: $\begin{array}{ccc}\mathit{p}\mathrm{\equiv }\frac{\mathrm{d}\ln \mathit{B}}{\mathrm{d}\ln \mathit{\varpi }} \mathrm{as} \mathit{\varpi }\mathrm{\to }\mathrm{0.}& & \end{array}$(3)Now, if we perturb a toroidal field in equilibrium, the resulting displacement can be described by $\begin{array}{ccc}\mathit{\xi }\mathrm{~}{\mathrm{e}}^{\mathrm{i}\mathrm{(}\mathit{m\phi }\mathrm{+}\mathit{kz}\mathrm{)}\mathrm{+}\mathit{\sigma t}}\mathit{,}& & \end{array}$(4)where the variation in the ϖ direction is not included as its length scale is much larger. Although displacements are predominantly horizontal because of the vertical stratification, the incompressible continuity equation (whose use is justified by the very high plasma-β) shows that a horizontal displacement ξ_h is necessarily accompanied by a radial displacement on the order of ξ_r ~ ξ_hl_r/l_h, where l_h and l_r are horizontal and radial length-scales, respectively, so that near the axis l_r = 1 /k. The shape of some of the unstable modes is shown in Fig. 1. Note that the m = 0 mode is also known as the sausage mode and the m = 1 mode as the kink mode.

Fig. 1 Physical form of the different unstable azimuthal m modes present in a plasma column with a purely toroidal field. Below each mode the cross section of the co-moving surface is shown along with the cross section of the equilibrium state. The only axisymmetric mode is m = 0 and the only mode that preserves the column’s cross section is the m = 1 mode.

In this figure we can observe why perturbations like (4) in a purely toroidal field are unstable. Regions where magnetic field is pinched have higher magnetic pressure forcing the fluid to displace opposite to the pressure gradient, leading to a runaway process. Figure 1 shows the equilibrium configuration and the shape of the azimuthal m = 0 and m = 1 unstable modes. In a stellar interior we often expect a vertical stratification, and vertical displacements experience a damping effect because they have to do work against the stratification. However, given that the m = 1 mode mostly drives azimuthal displacements, the damping effect of the stable stratification is very weak.

Tayler (1957) showed that for the instability to develop the field must satisfy $\begin{array}{ccc}\mathit{p}\mathit{>}\frac{{\mathit{m}}^{\mathrm{2}}}{\mathrm{2}}\mathrm{-}\mathrm{1}\hspace{0.17em}\hspace{0.17em}\mathrm{\left(}\mathit{m}\ne \mathrm{0}\mathrm{\right)}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{and}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathit{p}\mathit{>}\mathrm{1}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{(}\mathit{m}\mathrm{=}\mathrm{0}\mathrm{)}\mathit{.}& & \end{array}$(5)These instability conditions are valid only for non-viscous, non-rotating, non-diffusive, radiative stellar interiors. Examining this stability condition (5) we observe that for a p = 1 field, the m = 1 mode is the only unstable mode, the m = 0 and 2 modes are marginally stable and the other modes are stable. For a somewhat steeper radial field gradient of p = 2, modes m = 0,1,2 are expected to be unstable.

To study the growth of the instability we define the local toroidal² Alfvén frequency as $\begin{array}{ccc}{\mathit{\omega }}_{\mathrm{A}\mathit{,\phi }}\mathrm{\equiv }{\mathit{v}}_{\mathrm{A}\mathit{,\phi }}\mathit{/}\mathit{\varpi }\mathrm{=}\frac{{\mathit{B}}_{\mathit{\phi }}}{\mathit{\varpi }\sqrt{\mathrm{4}\mathit{\pi \rho }}}\mathit{,}& & \end{array}$(6)where v_A,φ is the toroidal Alfvén speed and ϖ the distance to the axis, we expect the growth rate of the instability to be $\begin{array}{ccc}\mathit{\sigma }\mathrm{~}{\mathit{\omega }}_{\mathrm{A}\mathit{,\phi }}\mathit{.}& & \end{array}$(7)Given that toroidal magnetic fields in stars are buried under the surface, there is no direct observational evidence of their configuration, therefore we have no knowledge of the radial magnetic field gradient. Even less understood is the configuration of a magnetic field inside a differentially rotating star. However we know that towards the magnetic axis the field must tend to zero, and to avoid singularities the radial gradient p ≥ 1 for ϖ → 0, and for this reason in this paper we start by studying the p = 1 case. Since it is expected that there is a qualitative threshold at p = 3/2 (see Sect. 2.1), we then also study the case of a p = 2 field gradient.

The unstable growth rate in Eq. (7) is valid only for purely toroidal, non-rotating, non-stratified and non-diffusive cases. We now discuss the behavior of the Tayler instability for situations where these very ideal conditions are relaxed and discuss what happens when magnetic diffusion is included, rotation is present and the magnetic field configuration is not purely toroidal anymore. By doing so, the stability conditions and growth rates are modified, and furthermore growth rates develop length scale dependencies. A stable stratification has major effects on the behavior of the instability, this is however outside the scope of this paper and is already discussed in Braithwaite (2006b).

2.1. Effects of rotation

All stars are rotating and a big fraction rotates at very high speeds (Collins 1963; Ramírez-Agudelo et al. 2013).

2.1.1. Rotation and magnetic axes aligned

Acheson (1978), Pitts & Tayler (1985) and Spruit (1999) showed that both the stability conditions and the unstable growth rate are affected by fast rotation. When the rotation frequency is higher than the Alfvén frequency (Ω_∥>ω_A,φ) and the rotation and magnetic axes are aligned, the stability condition is modified from (5) to $\begin{array}{ccc}\mathit{p}\mathit{>}\frac{{\mathit{m}}^{\mathrm{2}}}{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{\left(}\mathit{m}\ne \mathrm{0}\mathrm{\right)} \mathrm{and} \mathit{p}\mathit{>}\mathrm{1}\mathrm{(}\mathit{m}\mathrm{=}\mathrm{0}\mathrm{)}& & \end{array}$(8)so that the condition of stability against all m modes is p< 3/2 rather than p< −1/2 as before, and furthermore, the growth rate when any unstable mode is present is damped by a factor of ω_A,φ/ Ω_∥ (see Acheson 1978; Spruit 1999).

2.1.2. Rotation and magnetic axes not aligned

It is not always the case that the magnetic axis, the rotation axis and the gravity are aligned, as in the cases looked at above. Figure 2 (upper panel) shows the poloidal field lines in an axisymmetric configuration, and a toroidal-field-like configuration is observed around the neutral field line around the star, where the poloidal field vanishes. This configuration suffers from the same kind of instability (Tayler instability when considered locally), the difference being that the magnetic axis (in this case the neutral line) is perpendicular to gravity, but because of the switch in geometry, it is now the poloidal field that drives the instability, instead of the toroidal field as before; during local considerations it therefore makes sense to switch the definitions of the two words. Also note that in the presence of rotation, the magnetic and rotation axis wont be alligned somewhere in the star. Kink-like instabilities appear with displacements perpendicular to the direction of stratification, v·g ≈ 0, as appreciated in Fig. 2 (lower panel). Displacements are mainly in the direction of the arrows, but, as discussed in Sect. 2, there are also necessarily smaller displacements in the direction parallel to the neutral line. The ratio between the displacements in these two directions should be equal to the ratio of length scales in the same two directions, as we can see from the constraint ∇·v ≈ 0, a consequence of the high plasma-β.

The Coriolis force affects motions perpendicular to the rotation axis; in the case where the magnetic and rotation axes are aligned, unstable displacements are mainly perpendicular to the rotation and magnetic axes, and stabilization in the case −1/2 <p< 3/2 is achieved when Ω_∥ ≥ ω_A,φ. When the magnetic and rotation axes are not aligned, the amplitude of the displacements perpendicular to the rotation axis decrease with increasing inclination angle, and the extreme case occurs when the rotation and magnetic axes are perpendicular to each other. Then naively we think that in order to stabilize the perpendicular rotator, a faster rotation rate is required, and as we assumed that the fluid is almost incompressible, ∇·v ≈ 0, we derive the stability condition $\begin{array}{ccc}{\mathrm{\Omega }}_{\mathrm{\perp }}\mathrm{\ge }{\mathit{\omega }}_{\mathrm{A}\mathit{,\phi }}{\mathit{\varpi }}_{\mathrm{0}}\mathit{k,}& & \end{array}$(9)where Ω_⊥ is the perpendicular rotation rate with respect to the magnetic field axis, ϖ₀ is the typical length scale of the motion perpendicular to the magnetic field and k is the wavenumber in the direction of the magnetic axis (here the neutral line) as defined in (4). Note that the new perpendicular rotating stability conditions now dependent on the wavenumber, unlike in the aligned rotator case.

Fig. 2 Top: mixed poloidal-toroidal field configuration. If the toroidal field is too weak or indeed absent, the poloidal field is subject to an instability. Bottom: left, an equilibrium configuration with poloidal field loops threaded by the neutral line, and on the right-hand-side the kink-like instability present around the neutral line. Unstable displacements are restricted to be approximately perpendicular to the stratification. Figure taken from Braithwaite (2007).

Fig. 3 Representations of the first two vertical n modes of the unstable m = 1 mode. Top: purely toroidal magnetic field configuration, where the orientation of the magnetic symmetry axis is perpendicular to the rotation axis . It also shows what an n = 1 and n = 2 disk looks like as the red and blue box, respectively. Bottom: slipping disks for the n = 1 (left) and n = 2 (right) cases. The magnitude of the deformation Δl⊥, and the deformation velocity Δv⊥, suffered by the disks is different depending on the unstable mode. Small n modes present large deformations, whereas high n modes show small deformations, Δl⊥(n = 1) > Δl⊥(n = 2) and Δv⊥(n = 1) > Δv⊥(n = 2).

This effect can be explained if we imagine the plasma column to be a series of disks subject to pressure on top and below. When the disks displace perpendicular to the magnetic axis they suffer a deformation fattening at one side and shrinking at the other. The magnitude of the deformation depends on the size of the disks, as shown in Fig. 3, where disk of different sizes correspond to different vertical n modes. This is simply a consequence of the constraint ∇·v ≈ 0, which can be written as kv_z ~ ϖ₀v_ϖ. The Coriolis force affects motion in the vertical direction; Greater deformations, in other words small n modes, feel a stronger stabilizing effect. Note that the Coriolis force also affects motions perpendicular to both the magnetic and rotation axis, but this motion, at least for the m = 1 mode, can be minimized by combining left- and right-handed spiral modes with the right phase.

2.2. Effects of magnetic diffusion

Until now we have only considered adiabatic perturbations, that is ignoring the effects of thermal and magnetic diffusion, as well as viscosity. We now relax the flux freezing condition and allow the field lines to slip through the plasma.

Pitts & Tayler (1985) and Spruit (1999) suggested that in the presence of any diffusive process, thermal or magnetic, the instability condition in fast rotators (8) would recover its original form given in condition (5) (see Spruit 1999, Appendix). The reason for this can be thought of in terms of energy dissipation: rotation tends to give stability because the Coriolis force, which acts perpendicularly to the velocity and has consequently no effect on the energies, drives the fluid back into its equilibrium position; diffusion damps these epicycle motions, making a return to the equilibrium energetically impossible; as magnetic energy is turned into heat the system moves irreversibly away from its original equilibrium.

Magnetic diffusion is scale dependent, and therefore when included in the analysis, the unstable growth rate of the Tayler instability becomes scale dependent as well (see e.g. Spruit 1999). For a given magnetic diffusivity, η, the diffusion rate is given by, ηk², where k is the wavenumber. At first glance one expects this to impose an upper limit for the unstable wavenumbers that are able to grow; perturbations at a larger wavenumbers should be smoothed away faster than they can grow. We have stability if $\begin{array}{ccc}\mathit{\eta }{\mathit{k}}^{\mathrm{2}}\mathit{>}\mathit{\sigma }\mathit{.}& & \end{array}$(10)Therefore for a given wavenumber we can define a critical diffusivity $\begin{array}{ccc}{\mathit{\eta }}_{\mathrm{crit}}\mathrm{=}{\mathit{\omega }}_{\mathrm{A}\mathit{,\phi }}\mathit{/}{\mathit{k}}^{\mathrm{2}}\mathit{.}& & \end{array}$(11)From now on, we will express the different diffusive environments as fractions of the critical magnetic diffusivity $\begin{array}{ccc}{\mathit{n}}_{\mathit{\eta }}\mathrm{=}\frac{\mathit{\eta }}{{\mathit{\eta }}_{\mathrm{crit}}}\mathit{,}& & \end{array}$(12)where one can always move back and forth between fraction of critical magnetic diffusivity and the wavenumber of interest with the following relation $\begin{array}{ccc}{\mathit{n}}_{\mathit{\eta }}\mathrm{=}{\mathit{k}}^{\mathrm{2}}\left(\frac{\mathit{\eta }}{{\mathit{\omega }}_{\mathrm{A}\mathit{,\phi }}}\right)\mathit{,}& & \end{array}$(13)where η is the magnetic diffusivity of the system and k = 1 /l_z is the vertical wavenumber.

2.3. Mixed poloidal-toroidal configurations

The addition of a poloidal field is expected to modify or suppress the instability. If one imagines the toroidal and poloidal fields to operate independently of each other (in some vague sense) then some of the energy liberated from the toroidal field must be used to bend the poloidal field lines. If the poloidal field is strong, the instability might be suppressed altogether. However, the effect will depend on the wavenumber of the mode; whilst all wavenumbers release the same energy from the toroidal field, higher wavenumbers have to do more work against the poloidal field. This gives rise to a wavenumber stability threshold.

In general, if an instability displaces a fluid element horizontally a distance ξ_h, it does so by means of a force per unit mass of F_h = σ²ξ_h where σ is the growth rate. In the case of the Tayler instability, σ ≈ ω_A,φ as discussed at the beginning of Sect. 2 (Eq. (1)). Imagine, in the cylindrical geometry used above, that there is a poloidal field in the z direction. The displacement caused by the instability drags the poloidal field with it, which reacts against this bending with a force per unit mass of ${\mathit{B}}_{\mathit{z}}{\mathit{B}}_{\mathit{\varpi }}\mathit{k}\mathit{/}\sqrt{\mathrm{4}\mathit{\pi \rho }}\mathrm{\approx }{\mathit{B}}_{\mathit{z}}^{\mathrm{2}}{\mathit{\xi }}_{\mathrm{h}}{\mathit{k}}^{\mathrm{2}}\mathit{/}\sqrt{\mathrm{4}\mathit{\pi \rho }}$, where k is the radial (vertical) wavenumber. For a perturbed displacement to grow it is required that the instability force is bigger than the restoring force by bending the poloidal field line, which gives a condition on the unstable wavenumbers $\begin{array}{ccc}\mathit{k}\mathit{<}\frac{\mathit{\sigma }}{{\mathit{v}}_{\mathrm{A}\mathit{,z}}}\mathit{,}& & \end{array}$(14)where ${\mathit{v}}_{\mathrm{A}\mathit{,z}}\mathrm{=}{\mathit{B}}_{\mathit{z}}\mathit{/}\sqrt{\mathrm{4}\mathit{\pi \rho }}$ is the vertical Alfvén speed. See Spruit (2002) and Braithwaite (2009) for a fuller discussion of this topic.

3. Model

To determine the effects of rotation and magnetic diffusion on the Tayler instability, local simulations of a radiative stellar interior are performed. For our local simulations we select a box containing the magnetic axis of the star, located at the center of the box, as shown in Fig. 4. Global simulations of the stellar radiative region are not suitable at this stage because of resolution limitations to accurately resolve small scales. However, global simulations might eventually be necessary to investigate the non-linear behavior in a differentially rotating star, or rather, global in two dimensions and local in the radial direction.

Fig. 4 Predominantly toroidal magnetic field configuration in the interior of a star. The box shows a zoomed-in region containing the toroidal field lines around the magnetic axis and bold arrows showing the direction of the parallel and perpendicular rotation axes.

The equation of state used is the ideal equation of state, in units such that the molar gas constant divided by the molecular mass is unity $\begin{array}{ccc}\mathit{P}\mathrm{=}\mathit{\rho }\mathit{T}\mathrm{and}\mathit{e}\mathrm{=}\frac{\mathit{P}}{\mathit{\gamma }\mathrm{-}\mathrm{1}}\mathrm{·}& & \end{array}$(15)To trace the evolution of the system we solve the standard magnetohydrodynamic (MHD) equations, momentum equation: $\begin{array}{ccc}\frac{\mathit{\partial }\mathit{\rho }{v}}{\mathit{\partial t}}\mathrm{=}\mathrm{-}{\nabla }\mathrm{·}\mathrm{\left(}\mathit{\rho }{vv}\mathrm{\right)}\mathrm{-}{\nabla }\mathit{P}\mathrm{+}\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}{\nabla }\mathrm{\times }{B}\mathrm{\times }{B}\mathrm{-}{g}\mathit{\rho }\mathrm{+}\mathrm{2}\mathit{\rho }{v}\mathrm{\times }\mathrm{\Omega }\mathit{,}& & \end{array}$(16)continuity equation, $\begin{array}{ccc}\frac{\mathit{\partial \rho }}{\mathit{\partial t}}\mathrm{=}\mathrm{-}\mathrm{\nabla }\mathrm{·}\mathit{\rho }{v}\mathit{,}& & \end{array}$(17)induction equation, $\begin{array}{ccc}\frac{\mathit{\partial }{B}}{\mathit{\partial t}}\mathrm{=}\mathrm{\nabla }\mathrm{\times }\left[\mathrm{(}{v}\mathrm{\times }{B}\mathrm{)}\mathrm{-}\mathit{\eta }\mathrm{\nabla }\mathrm{\times }{B}\right]\mathit{,}& & \end{array}$(18)and the specific internal energy: $\begin{array}{ccc}\frac{\mathit{\partial e}}{\mathit{\partial t}}\mathrm{=}\mathrm{-}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{e}{v}\mathrm{\right)}\mathrm{-}\mathit{P}\mathrm{\nabla }\mathrm{·}{v}& & \end{array}$(19)where ρ is the density, v the velocity, B the magnetic field, η the magnetic diffusivity, e the internal energy, g the gravity, P gas pressure, and T gas temperature. The ratio of specific heats γ is set to 5/3. We note that in the momentum Eq. (16) we neglected viscosity and the Coriolis term is explicitly included. In the induction Eq. (18) magnetic diffusion η is spatially constant. The energy equation neglects viscous and Ohmic heating and thermal conduction, since we are not interested at this stage in viscous and thermal effects. We do not include the effects of gravity in this study.

3.1. Numerical code

We use the stagger code, a parallel 3D MHD code originally developed by (Nordlund & Gaalsgard 1995), and continuously improved over the years by its user community. The Stagger-code uses a sixth-order explicit finite-difference scheme for numerical derivatives and the corresponding fifth- order interpolation scheme. The solution of the hydrodynamic equations is advanced in time using an explicit third-order Runge-Kutta integration method (Williamson 1980). The code operates in a staggered, Eulerian, Cartesian mesh, which seem to be a better implementation than cylindrical and spherical coordinates, as a Cartesian mesh avoids the coordinate singularity in the axis present in cylindrical and spherical grid-based codes. Another advantage of a Cartesian based code is reflected in the simplicity of the equations with respect to another coordinate system. Cartesian coordinates also have their own disadvantages such as the noise introduced by fitting a cylindrical shape, e.g. torus, in squared boxes, or the extra computational time cost due to the un-used corners of the box.

For the simulations of purely toroidal fields a resolution of 48² in the azimuthal direction and 32 in the vertical direction was used, for the mixed poloidal-toroidal fields the resolution was 48³. Test cases with a resolution of 96³ were performed in order to compare the results with the low resolution simulations and no significant differences were observed. Periodic boundaries are used in all simulations; see below.

The code also contains highly localized diffusivities that prevent the growth of wavelengths close to the Nyquist frequency. We however decided to switch off this artificial diffusive terms in order to avoid non-physical effects. This was possible given the clean initial conditions, where the development and growth of the Tayler instability occurs at very early times when no artificial stabilization is required.

3.2. Initial conditions

The initial conditions correspond to the equilibrium configuration of a plasma containing a magnetic field of the form: $\begin{array}{ccc}{B}\mathrm{=}{\mathit{B}}_{\mathit{\phi }}\mathrm{\left(}\mathit{\varpi }\mathrm{\right)}{{e}}_{\mathit{\phi }}\mathrm{+}{\mathit{B}}_{\mathit{z}}{{e}}_{\mathit{z}}\mathit{.}& & \end{array}$(20)The magnetic configuration in this paper has no radial component, B_ϖ = 0, and is a combination of an azimuthal, B_φ, and a vertical, B_z, components with no dependence on height z or angle φ.

The azimuthal component is an axisymmetric torus defined by $\begin{array}{ccc}{\mathit{B}}_{\mathit{\phi }}\mathrm{\left(}\mathit{\varpi }\mathrm{\right)}\mathrm{=}{\mathit{B}}_{\mathrm{0}\mathit{,\phi }}{\left(\frac{\mathit{\varpi }}{{\mathit{\varpi }}_{\mathrm{0}}}\right)}^{\mathit{p}}{\mathrm{e}}^{\mathrm{-}\mathrm{(}\mathit{\varpi }\mathit{/}{\mathit{\varpi }}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}\mathit{,}& & \end{array}$(21)where ϖ₀ is a constant length, set to one quarter of the horizontal size of the box, p is the radial field gradient, set to either 1 or 2 in the simulations, and B_0,φ is a constant. The exponential term exp(ϖ²/${\mathit{\varpi }}_{\mathrm{0}}^{\mathrm{2}}\mathrm{)}$ shapes the field such that it behaves as B ∝ ϖ^p near the axis and vanishes farther away so as to avoid problems at the boundaries. The local toroidal Alfvén frequency, as defined in (6), is ${\mathit{\omega }}_{\mathit{A}\mathit{,\phi }}\mathrm{=}{\mathit{B}}_{\mathit{\phi }}\mathit{/}\mathit{\varpi }\sqrt{\mathrm{4}\mathit{\pi \rho }}$. The maximum initial field strength in the box corresponds to a minimum plasma beta parameter of β = 10, i.e. thermal pressure dominates over magnetic pressure.

The initial poloidal field is implemented as a uniform field in the vertical direction $\begin{array}{ccc}{\mathit{B}}_{\mathit{z}}\mathrm{=}{\mathit{B}}_{\mathrm{0}\mathit{,}\mathit{z}}& & \end{array}$(22)for this field configuration the vertical Alfvén frequency is also constant everywhere, ω_A,z = v_A,z(2π/L_z), where ${\mathit{v}}_{\mathrm{A}\mathit{,z}}\mathrm{=}{\mathit{B}}_{\mathit{z}}\mathit{/}\sqrt{\mathrm{4}\mathit{\pi \rho }}$ is the vertical Alfvén velocity and L_z is the vertical extent of the box.

Fig. 5 a) Initial purely toroidal magnetic configuration. The region of interest is between 0 ≤ ϖ ≲ ϖ0 given that the field strength falls off at greater ϖ. Colored areas represent surfaces of equal distance to the axis. Field lines are also plotted; their density is related to the increase in the field strength proportional to the distance from the axis. b) Mixed poloidal-toroidal field configuration. The field lines show a spiral morphology, for this particular case (with p = 1) two spirals turns are completed in the box, i.e. toroidal Alfvén frequency is double the vertical Alfvén frequency ωA,φ = 2ωA,z.

A visual example of the field configuration, for purely toroidal and mixed poloidal-toroidal fields, is shown in Fig. 5.

The initial conditions for the thermodynamic variables are determined by the equilibrium solutions found by setting ∂/∂t = v = η = 0 in Eqs. (16) to (19), with the initial magnetic field given by Eqs. ((20), (22) and (21)). No height dependence is considered, therefore pressure, density and temperature depend only on the distance from the axis. For a full derivation of these stability configurations see Braithwaite (2006b; Sect. 3.3).

Four different cases for two different field configurations are simulated:

• 1

  rotating diffusive star,

• 2

  perpendicular rotator,

• 3

  rotating star with a quadratic radial field gradient: p = 2, and

• 4

  mixed poloidal-toroidal field.

• Case 1:

  rotating diffusive star with linear radial field gradient,p = 1. From (21) we see that the toroidal Alfvén frequency varies as ω_A ∝ ϖ^p−1 at low ϖ, so that it is approximately independent of location in the inner part of the computational box. The magnetic and rotation axes are aligned. We explore the parameter space of rotation-only, magnetic diffusion-only and the simultaneous action of rotation and diffusion.

• Case 2:

  perpendicular rotator with linear radial field gradient, p = 1. Rotation and magnetic axes are perpendicular. We explore the rotation rate parameter space.

• Case 3:

  rotating star with a quadratic radial field gradient, p = 2. The Alfvén frequency increases linearly with the distance to the axis. For this case rotation and magnetic axes are aligned. We explore the rotation rate parameter space.

• Case 4:

  mixed poloidal-toroidal configuration. The poloidal field is constant in the box, therefore the vertical Alfvén frequency is also constant. The toroidal field has a linear radial field gradient, p = 1, thus the azimuthal Alfvén frequency is also constant in the inner part of the box. We explore the poloidal-to-toroidal field strength ratio parameter space.

Given that the initial configurations correspond to equilibrium solutions, an initial perturbation is included to trigger instability.

3.3. Initial perturbation

The initial conditions in Sect. 3.2 correspond to a perfect hydrostatic equilibrium, in order to see any instability grow we seed tiny perturbations in the velocity field. Given the cylindrical symmetry of this problem, we can perturb the different azimuthal, m, and vertical, n, modes independently. The shape and description of the perturbation used here are the same as used by Braithwaite (2006b, Sect. 3.5). The perturbation has its maximum at ϖ = 0 and decreases with increasing ϖ, almost vanishing at ϖ = ϖ₀. There is no perturbation to the vertical velocity component.

The initial perturbation has the form $\begin{array}{ccc}{v}& \mathrm{=}& \sum _{\mathit{m}\mathrm{=}\mathrm{0}\mathit{,n}}{\mathit{V}}_{\mathrm{0}}^{\mathit{n}}\frac{\mathit{\varpi }}{{\mathit{\varpi }}_{\mathrm{0}}}\exp \left(\mathrm{-}\frac{\mathrm{5}}{\mathrm{2}}\frac{{\mathit{\varpi }}^{\mathrm{2}}}{{\mathit{\varpi }}_{\mathrm{0}}^{\mathrm{2}}}\right)\cos \mathrm{\left(}\mathit{kz}\mathrm{\right)}{{e}}_{\mathit{\varpi }}\\ & & \end{array}$The amplitude coefficients of the perturbed velocity field are normalized to the maximum Alfvén velocity in the box: ${\mathit{V}}_{\mathit{m}}^{\mathit{n}}\mathit{/}{\mathit{v}}_{\mathrm{A}\mathit{,}\max }\mathrm{=}\mathrm{1}\mathrm{\times }{\mathrm{10}}^{-8}$. Such a value is large enough to trigger the Tayler instability and small enough to let us follow the linear growth phase for several growth timescales.

3.4. Boundary conditions

The box represents a small volume of a radiative zone in a stellar interior. Figure 4 is a cartoon of the box location in star with respect to the magnetic field axis.

The difference in length scales between the instability and the radius of the star are huge. The instability is many orders of magnitude smaller³ than the radius of the star (l_r ≪ R_∗), these instabilities are local in the radial direction and global in the azimuthal direction, therefore different shells of the star are in effect disconnected from each other.

Periodic boundaries are the best option for our simulation, given that we want to avoid the selection of special locations in the stellar interior and to prevent additional effects such as current sheets or shear.

4. Results

First, test cases without rotation or magnetic diffusion were simulated with field gradient p = 1. The Tayler instability is properly traced for a number of growth timescales. More physical ingredients are built up each step in order to approach our simple model with the conditions found in a radiative stellar interior.

For the analysis of the simulations we need a formalism for measuring the amplitude of the various azimuthal m and vertical n modes, keeping in mind that n is the dimensionless vertical wavenumber; the n = 1 mode corresponds to a wavelength equal to the height of the computational box. We extract information from two fields: the change in the magnetic field δB = B−B(t = 0), and the velocity field v. If we imagine the field to be frozen into the fluid, then δB reflects the displacement ξ, the integral of the velocity v.

We can now decompose the φ component of these fields using Fourier series in order to extract the amplitudes of the different azimuthal m modes. For instance with the velocity field $\begin{array}{ccc}{v}\mathrm{\left(}\mathit{\varpi ,\phi ,}\mathit{z}\mathit{,}\hspace{0.17em}\mathit{t}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{{A}}_{\mathrm{0}}\mathrm{\left(}\mathit{\varpi ,}\mathit{z}\mathit{,t}\mathrm{\right)}\mathrm{+}\sum _{\mathit{m}\mathrm{=}\mathrm{1}}^{\mathit{m}\mathrm{=}\mathrm{\infty }}\hspace{0.17em}\mathrm{\Re }\mathrm{(}{{A}}_{\mathit{m}}\mathrm{(}\mathit{\varpi ,}\mathit{z}\mathit{,}{\mathrm{t}}^{\mathrm{)}}\hspace{0.17em}{\mathrm{e}}^{\mathit{im}\mathit{\phi }}\mathrm{)}\mathit{,}& & \end{array}$(23)where A₀ and A_m are the Fourier coefficients given by $\begin{array}{ccc}{{A}}_{\mathrm{0}}\mathrm{\left(}\mathit{\varpi ,}\mathit{z}\mathit{,t}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathit{\pi }}{\mathrm{\int }}_{\mathrm{-}\mathit{\pi }}^{\mathit{\pi }}{v}\mathrm{\left(}\mathit{\varpi ,\phi ,z,t}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\phi }& & \end{array}$and $\begin{array}{ccc}{{A}}_{\mathit{m}}\mathrm{\left(}\mathit{\varpi ,}\mathit{z}\mathit{,t}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathit{\pi }}{\mathrm{\int }}_{\mathrm{-}\mathit{\pi }}^{\mathit{\pi }}{v}\mathrm{\left(}\mathit{\varpi ,\phi ,z,t}\mathrm{\right)}\hspace{0.17em}\cos \mathrm{\left(}\mathit{m\phi }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\phi }\mathit{.}& & \end{array}$We then compute the amplitude of each m mode as $\begin{array}{ccc}\mathit{v}{\mathit{\_}\mathrm{Amp}}_{\mathit{m}}\mathrm{\left(}\mathrm{t}\mathrm{\right)}\mathrm{=}{\left(\frac{{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{\varpi }}_{\max }}\mathrm{\int }\begin{array}{c}\mathit{L}z\\ \mathrm{0}\end{array}{\mathrm{A}}_{\mathrm{m}}^{\mathrm{\ast }}\hspace{0.17em}{\mathrm{A}}_{\mathrm{m}}\hspace{0.17em}\mathrm{2}\mathit{\pi \varpi }\hspace{0.17em}\mathrm{d}\mathit{\varpi }\mathrm{d}\mathit{z}}{{\mathit{\varpi }}_{\max }^{\mathrm{2}}{\mathit{L}}_{\mathit{z}}\mathit{\pi }}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}& & \end{array}$(24)where ϖ_max is set to ϖ₀/ 2, as the magnetic field strength in the region between ϖ = 0 to ϖ₀/ 2 is approximately proportional to the distance from the axis to the power p, i.e. B_φ ∝ ϖ^p. A similar procedure can be applied to the magnetic fields to derive the magnetic displacement amplitude, δB_Amp_m(t).

Fig. 6 Amplitude of the first four m modes as a function of time, extracted from the magnetic field δB_Ampm(t), and normalized to the wavelength corresponding to the n = 1 vertical mode. The azimuthal m = 1 mode, solid red line, is the only unstable mode.

Figure 6 shows the evolution of the amplitudes of the first three m modes as a function of time. The m = 1 mode is the only unstable mode, as predicted above in (5). When the amplitude of the unstable m = 1 mode is comparable to the wavelength, δB_Amp_m ≈ λ_{n = 1}, the system goes into the non-linear regime. At this point energy can be transferred between modes and m ≠ 1 modes begin to grow. This is observed to happen around t/τ_A ≈ 10 (where τ_A = 1 /ω_A) in Fig. 6.

The displacement amplitude computed using Eq. (24) contains the contribution of all vertical wavelengths. For further analysis, it will be necessary to measure the amplitudes of each individual vertical n modes separately. A second Fourier decomposition, similar to the one applied in Eq. (23), is performed in the vertical direction, $\begin{array}{ccc}\mathit{\delta }{v}\mathrm{\left(}\mathit{\varpi ,\phi ,z,t}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{{C}}_{\mathit{m,}\mathrm{0}}\hspace{0.17em}\mathrm{\left(}\mathit{\varpi ,t}\mathrm{\right)}\mathrm{+}\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathit{n}\mathrm{=}\mathrm{\infty }}\mathrm{\Re }\mathrm{(}{{C}}_{\mathit{m,n}}\mathrm{(}\mathit{\varpi ,}{\mathit{t}}^{\mathrm{)}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{i}\mathit{nz}}\mathrm{)}\mathit{,}& & \end{array}$(25)where C_m,0 and C_m,n are the new Fourier coefficients given by: $\begin{array}{ccc}{{C}}_{\mathit{m,}\mathrm{0}}\mathrm{=}\frac{\mathrm{1}}{\mathit{\pi }{\mathit{L}}_{\mathit{z}}}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{L}}_{\mathit{z}}}{\mathrm{\int }}_{\mathrm{-}\mathit{\pi }}^{\mathit{\pi }}\mathit{\delta }{v}\hspace{0.17em}\cos \mathrm{\left(}\mathit{m\phi }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\phi }\hspace{0.17em}\mathrm{d}\mathit{z}& & \end{array}$and $\begin{array}{ccc}{{C}}_{\mathit{m,n}}\mathrm{=}\frac{\mathrm{1}}{\mathit{\pi }{\mathit{L}}_{\mathit{z}}}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{L}}_{\mathit{z}}}{\mathrm{\int }}_{\mathrm{-}\mathit{\pi }}^{\mathit{\pi }}\mathit{\delta }{v}\hspace{0.17em}\cos \mathrm{\left(}\mathit{m\phi }\mathrm{\right)}\hspace{0.17em}\cos \mathrm{\left(}\mathit{nz}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\phi }\hspace{0.17em}\mathrm{d}\mathit{z}& & \end{array}$and are used to compute the displacement amplitude for each individual pair of azimuthal m and vertical n modes, v_Amp_m,n(t), similar to the procedure applied in Eq. (24).

4.1. Case 1 – rotating diffusive star

We first simulated the effects of just rotation and just diffusion on their own, for the constant-Alfvén-frequency case, p = 1. Rotation is parallel to the magnetic axis. We compare our results to those presented by Braithwaite (2006b) and show that we obtained similar results.

The rotation-only case has a threshold for the rotation rate equal to the Alfvén frequency, Ω_∥ = ω_A,φ. Rotation rates faster than that completely suppress the Tayler instability, presumably imposing some epicyclic motion of the perturbed fluid elements around an equilibrium point.

Fig. 7 Amplitude δB_Ampm,n(t) of the m = 1, n = 1 mode as a function of time for various rotation rates. Blue lines correspond to rotation rates below the Alfvén frequency, and red lines correspond to rotation rates above the Alfvén frequency.

Figure 7 shows the amplitude of the m = 1, n = 1 mode in simulations with different Ω_∥, ranging from Ω_∥/ω_A,φ = 0−10. It is clearly observed how the system is unstable for Ω_∥/ω_A,φ< 1 (blue lines), and is stabilized for rotation rates Ω_∥/ω_A,φ> 1 (red lines).

In contrast, the diffusion-only case shows no threshold beyond which the Tayler instability is completely stabilized.

Fig. 8 Growth rate of the Tayler instability for the m = 1 and n = 1 mode as a function of the relative magnetic diffusivity η/ηcrit measured from the velocity field, v_Ampm,n(t).

Figure 8 shows the measured unstable growth rates in terms of n_η. It is observed that even at diffusivities as high as n_η = 50; unstable growth is still present.

For the simultaneous action of both rotation and magnetic diffusion a total of 56 simulations were conducted. Figure 9 summarizes the measured growth rates normalized to the Alfvén frequency, σ/ω_A,φ.

We consider a system to be stable if no unstable growth is observed for a hundred Alfvén crossing timescales, 100τ_A, corresponding to the white region in Fig. 9. We see that at very low magnetic diffusivities, left side of the figure, the system behaves similarly to the non-diffusive case, the system is stable for rotation rates Ω_∥>ω_A,φ. Now if we move from the stable region in a straight horizontal line towards the right side of the figure we go towards higher magnetic diffusive environments and penetrate a region of unstable growth. This is a very interesting result as one would prima facie expect that two stabilizing processes, such as rotation and magnetic diffusion, to collaborate in suppressing the instability, but instead the presence of magnetic diffusion disturbs the stabilization introduced by rotation.

Presumably this happens because the Coriolis force pushes a fluid element onto an epicyclic equilibrium motion; magnetic diffusion constantly removes energy from the magnetic field such that there is no longer sufficient energy to return the fluid element to an equilibrium epicycle, and it will slowly spiral outwards.

The dashed yellow line in Fig. 9 is a relation between n_Ω∥ and n_η. If we assume that the boundary dividing the stable and unstable regimes can be separated by a linear relation we obtain $\begin{array}{ccc}{\mathit{n}}_{{\mathrm{\Omega }}_{\mathrm{\parallel }}}\mathrm{\approx }\mathrm{2}{\mathit{n}}_{\mathit{\eta }}\mathrm{+}\mathrm{1.}& & \end{array}$(26)It is not clear why the threshold should be here. The factor of 2 may come from the 2 in the Coriolis force itself.

Fig. 9 Contours of mean growth rate plotted as a function of normalized rotation rate nΩ∥ = Ω∥/ωA,φ and the relative magnetic diffusion nη = η/ηcrit. The white region represents stability, for at least 100τA. The color contours correspond to increasing growth rates of the instability, where black corresponds to small and blue to fast growth rates. The yellow dashed line corresponds to the relation nΩ∥ = 2nη + 1 separating the stable and unstable regions.

Alternatively the stability condition can be expressed as $\begin{array}{ccc}{\mathit{\omega }}_{\mathrm{A}\mathit{,\phi }}\mathit{<}{\mathrm{\Omega }}_{\mathrm{\parallel }}\mathrm{-}\mathit{\eta }{\mathit{k}}^{\mathrm{2}}\mathit{.}& & \end{array}$(27)This is the stability condition for fast rotating, diffusive systems as a function of the perturbed wavenumber k. The growth rate of the instability is in general a function of the Alfvén frequency, rotation rate, and magnetic diffusion rate, σ = σ(ω_A,Ω_∥,ηk²) for the non-stratified case.

4.2. Case 2 – perpendicular rotator

We simulated a rotator perpendicular to the magnetic field axis of symmetry. The magnetic field is oriented in the z-axis, , and the rotation axis is oriented along the ϖ(φ = 0)-direction Ω = Ω_⊥(φ = 0)e_ϖ. As in the previous section, the field gradient is p = 1. We expect rotation to have a different effect on different vertical wavelengths, as given by Eq. (9). Figure 10 shows the growth rates of the first three vertical n modes (all m = 1) as a function of the rotation rate normalized to the critical rotation for each n mode. We can observe how the growth rate of the instability is damped with increasing rotation rate and furthermore how this damping affects different vertical modes. It is however important to remember that the horizontal axis of the plot is normalized to the critical rotation rate for each vertical mode. This means that given a rotation rate of e.g. Ω_⊥ = 4Ω_{⊥ ,crit}(n = 1) = 20ω_A,φ, the n = 1 mode is damped by ~90%, n = 2 mode by ~45%, and n = 3 mode by ~25%. In agreement with our prediction that perpendicular rotation damps different vertical modes with a different strength.

Fig. 10 Growth rate of the Tayler instability, normalized to the Alfvén frequency, as a function of the perpendicular rotation rate, normalized to the predicted critical rotation rate for each n mode. The lines approximately on top of each other show the dependence of the damping on wavenumber is as predicted.

Even if we did not see complete stability in our simulation, Fig. 10 suggests that there is stabilization for each vertical mode at a rotation rates of Ω_⊥ ≈ 5Ω_{⊥ ,crit}(n). Apparently our stability condition (9) is missing a factor of order unity. This is not surprising, given that the damping effect should depend on the precise arrangement of the displacements within the mode.

4.3. Case 3 – quadratic radial field gradient

We saw above that in the case where p = 1, rotation can completely stabilize the magnetic field, at least where there is no diffusion. We now investigate the case with a radial field gradient of p = 2, which is predicted to behave differently. The Alfvén frequency grows proportional to distance from the axis ω_A,φ(ϖ) ∝ ϖ. For this reason we expect that the instability grows faster at larger radii than close to the axis. We introduce a new quantity, δB_Amp_m,n(ϖ,t), corresponding to the displacement amplitude, for a given pair of azimuthal and vertical modes, as a function of time and distance to the axis ϖ: $\begin{array}{ccc}\mathit{\delta B}{\mathit{\_}\mathrm{Amp}}_{\mathit{m,n}}\mathrm{\left(}\mathit{\varpi ,t}\mathrm{\right)}\mathrm{=}{\left(\frac{{{C}}_{\mathit{m,n}}^{\mathrm{\ast }}\mathrm{\left(}\mathit{\varpi ,t}\mathrm{\right)}{{C}}_{\mathit{m,n}}\mathrm{\left(}\mathit{\varpi ,t}\mathrm{\right)}}{{\mathit{\varpi }}^{\mathrm{2}}{\mathit{L}}_{\mathit{z}}\mathit{\pi }}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{·}& & \end{array}$(28)This new quantity allows us to study the growth rate of the instability at different distances from the axis.

Fig. 11 Growth rate σ of the Tayler instability as a function of the solid-body rotation rate Ω∥, in the p = 2 case. Normalized local growth rates are given for three different distances from the axis, σ(ϖ). The dashed line corresponds to the predicted growth rate of the Tayler instability for slow rotators, σ = ωA,φ. The dotted-dashed line is the predicted behavior of the instability for fast rotators, $\mathit{\sigma }\mathrm{=}{\mathit{\omega }}_{\mathrm{A}\mathit{,\phi }}^{\mathrm{2}}\mathit{/}\mathrm{2}{\mathrm{\Omega }}_{\mathrm{\parallel }}$.

Local growth rates, σ(ϖ), are computed from the local amplitude displacements, δB_Amp_m,n(ϖ,t). This quantity is then divided by the local Alfvén frequency, ω_A,φ(ϖ), and the normalized local growth rate is obtained, σ(ϖ) /ω_A,φ(ϖ).

Figure 11 shows the behavior of these normalized local growth rates as a function of rotation rate. The first main feature of the figure is that even for the very fast rotators, instability is still present. This is in agreement with the prediction that instability is always present in fast rotating stars containing a purely toroidal field with a radial field gradient p ≳ 3/2, obtained from the stability conditions (8).

There are two behaviors observed in this plot separated by a critical rotation rate at Ω_∥ ≈ ω_A,φ. In the slow rotation regime where Ω_∥/ω_A,φ ≲ 1, the growth rate σ ≈ ω_A,φ, whereas in the fast rotation regime where Ω_∥/ω_A,φ ≳ 1, we have instead $\mathit{\sigma }\mathrm{\approx }{\mathit{\omega }}_{\mathrm{A}\mathit{,\phi }}^{\mathrm{2}}\mathit{/}\mathrm{2}{\mathrm{\Omega }}_{\mathrm{\parallel }}$. Which is exactly what Pitts & Tayler (1985) and Spruit (1999) predicted a couple of decades ago⁴.

4.4. Case 4 – mixed poloidal-toroidal field

This section is intended to shed some light in the stability of mixed poloidal-toroidal magnetic field configurations in radiative stellar interiors. We focus our analysis on the behavior of the instability known to be present in purely toroidal fields, but now in mixed poloidal-toroidal field configurations. As discussed in Sect. 2.3 a poloidal field could damp, and if strong enough probably be able to suppress, the instability of a toroidal field. Tayler’s (1973) stability conditions however are not valid for mixed field configurations.

We note that the results presented in this section should be reviewed carefully because of the appearance of an unknown, spurious mode in the simulations growing always at the smallest scale in the vertical direction⁵. Such a mode was absent in purely poloidal and purely toroidal field simulations. We work under the hypothesis that any well resolved unstable mode growing faster than this unknown perturbation is real, and relevant for our analysis. The nature of the unknown mode will be studied in detail in a forthcoming paper and it must be established whether it is astrophysically relevant.

We first introduce some important definitions to avoid misunderstandings; the azimuthal Alfvén frequency, ω_A,φ = v_A,φ/ϖ, is calculated with the toroidal Alfvén velocity. In the initial conditions, the radial field gradient of the toroidal field is chosen for simplicity to grow linearly with the distance from the axis, p = 1, therefore the Alfvén frequency is constant throughout the volume of interest. The vertical Alfvén frequency ω_A,z is calculated with the vertical Alfvén velocity, v_A,z and is given by $\begin{array}{ccc}{\mathit{\omega }}_{\mathrm{A}\mathit{,z}}\mathrm{=}\frac{{\mathit{v}}_{\mathrm{A}\mathit{,z}}}{{\mathit{L}}_{\mathit{z}}}\mathrm{=}\frac{{\mathit{B}}_{\mathit{z}}}{{\mathit{L}}_{\mathit{z}}\sqrt{\mathrm{4}\mathit{\pi \rho }}}\mathrm{·}& & \end{array}$(29)The poloidal field strength is constant everywhere, B_z = B₀, thus the vertical Alfvén velocity and the vertical Alfvén frequency are constant.

For this field configuration field lines follow spirals as shown in Fig. 5b. The number of spirals a field line completes in the box depends on the relative strengths of the toroidal and poloidal fields. The more complete spirals a field line draws, the stronger the toroidal field is. We call this number of spirals in the box the pitch-number, P_n, which is given by $\begin{array}{ccc}{\mathit{P}}_{\mathit{n}}\mathrm{=}\frac{{\mathit{\omega }}_{\mathrm{A}\mathit{,\phi }}}{\mathrm{2}\mathit{\pi }{\mathit{\omega }}_{\mathrm{A}\mathit{,z}}}\mathrm{;}& & \end{array}$(30)re-writing the instability condition (14) in terms of the pitch-number and solving for the n = 1 mode, we obtain that instability is expected to be present for P_n> 1. We note that P_n is only independent of ϖ in the special case p = 1, which is what we use here.

Fig. 12 Growth rate of the first two azimuthal m modes as a function of the pitch-number. The solid black line corresponds to the m = 0 mode, stable at very high and low pitch-numbers but with a peak in unstable growth near the critical pitch–number, Pn ≈ 1. The red dashed line corresponds to the m = 1 mode, unstable for very high pitch–numbers and stabilized near the critical pitch-number, Pn ≈ 1.

Figure 12 shows the behavior of the first two azimuthal modes, m = 0 and m = 1, as a function of pitch–number. The m = 1 mode has a growth rate on the order of the Alfvén frequency for a predominantly toroidal field configuration, i.e. large pitch-number, and drops suddenly when the pitch–number reaches a critical value. This drop of the growth rate agrees with the stability conditions (14). We note that no complete suppression of the instability is observed, even for very small pitch-numbers. However, these small growth rates are the effect of non-linear interaction with the spurious modes.

The second interesting mode shown in Fig. 12 is the m = 0 mode. This mode was marginally stable for the purely toroidal field (5) and it is being perturbed in the mixed configuration. Growth of the m = 0 mode begins to appear with increasing poloidal field strength, reaching a maximum near the critical poloidal-to-toroidal field ratio, P_n ≈ 1.

The presence of this unstable mode, not predicted analytically, might be important for the hydromagnetic dynamo theory, given that this instability triggers axisymmetric vertical motions which it is often claimed are necessary for closure of the dynamo loop. We will investigate the importance of this mode in a future paper.

5. Discussion and summary

Solid-body rotation, magnetic diffusion and field configuration are expected to have a strong influence on the presence and evolution of hydro-magnetic instabilities in toroidal fields in radiative stellar interiors. We present a set of 3D MHD simulations of purely toroidal and of mixed poloidal-toroidal magnetic fields in non-convective stellar interiors. A combination of various radial field strength gradients (p = dlnB/ dlnϖ), solid body rotation rates (Ω), magnetic diffusivities (η), and poloidal-to-toroidal field strengths (P_n) is scoped.

We observe the presence of the kink type m = 1 Tayler instability in purely toroidal fields. We measure its growth rate, and compare our results with the various analytical results for non-rotating and non-diffusive (Tayler 1973), fast rotating (Acheson 1978; Pitts & Tayler 1985), highly diffusive (Spruit 1999), and rotating-diffusive (Pitts & Tayler 1985; Spruit 1999) radiative stellar interiors. We also study the development of instabilities in a mixed poloidal-toroidal field configuration.

Magnetic diffusion is not able to suppress the occurrence of the instability even at high diffusivity, in accordance with the results of Braithwaite (2006a).

Diffusion and rotation parallel to magnetic axis: if the magnetic and rotation axes are parallel, fast rotation (Ω_∥>ω_A,φ) is able to stabilize the configuration (see Fig. 7 ), where ω_A,φ is the azimuthal Alfvén frequency. This is in the case where the field gradient p = 1. However as soon as magnetic diffusion is included, the stability conditions now depend on a combination of the rotation rate, the magnetic field strength and the magnetic diffusion rate given by Eq. (27), Which suggests that very fast rotation can stabilize the system only if the magnetic diffusion rate does not subtract energy fast enough such that Coriolis force is unable maintain the fluid elements in stable epicyclic motions.

Perpendicular rotation: if the magnetic field axis and the rotation axis are perpendicular, the damping strength introduced by rotation is different for different vertical wavelengths. The damping effect is stronger for small vertical wavenumbers k compared to large wavenumbers because of the larger deformation occurring perpendicular to the rotation axis at smaller wavenumbers. It is found that the instability is suppressed if Ω_⊥ ≳ 5ω_A,φϖ₀k, where ϖ₀ is the length scale of the motions in the horizontal direction.

Radial field gradient: to analyze the effects of the radial field gradient and rotation in the instability we investigated two cases, the linear radial field gradient, p = 1 and the quadratic radial field gradient, p = 2. As described above, in the p = 1 case (without diffusion) stability is reached for fast rotation (meaning faster than the azimuthal Alfvén frequency). In this p = 2 case, very fast rotation is never able to stabilize the configuration. Linear growth rates behave broadly as predicted by Pitts & Tayler (1985), where for slow rotators the growth rate is roughly equal to the Alfvén frequency, whilst for fast rotators, the growth rate is damped by a factor of ω_A,φ/ 2Ω_∥. It is probably safe to assume that this is also the case for all values of p above 3/2, as originally predicted. We have therefore provided a check on the theory, useful in light of the fact that the originators of the analytic results were, we think it is fair to say, less than 100% certain they were correct.

Mixed poloidal-toroidal field: finally we simulated a mixed poloidal-toroidal field, and in particular the stabilizing influence of a uniform poloidal field. In the p = 1 case, where in the absence of a poloidal field, only the m = 1 mode is unstable, we find that the m = 1 mode is indeed suppressed by the poloidal field, and that the stability criterion is where it is predicted to be. Interestingly, we also find that the m = 0 mode is unstable in the vicinity of this threshold.

The presence of the unstable axisymmetric m = 0 mode in the mixed poloidal-toroidal field configuration may be quite important. In the case where the field gradient is such that it is marginally stable in the purely toroidal case, it becomes dominant in the vicinity of the threshold of stabilization by the poloidal field, which of course is precisely where the dynamo is likely to be. It opens the possibility for axisymmetric turbulent radial motions which affect the mixing and the regeneration of the poloidal field. It is often claimed by various members of the MHD community that a poloidal field cannot be renewed without the m = 0 mode; whether this claim is correct is, if the m = 0 mode is present and perhaps even dominant, no longer relevant.

Indeed, as far as the field gradient p is concerned, in the fast-rotating (which is expected to occur much more frequently in nature than the slowly rotating) adiabatic case we can summaries behavior in terms of the threshold: where p< 1 all modes are stable; where 1 <p< 3/2 just the m = 0 mode is unstable; where p> 3/2 the m = 1 mode is also unstable, with thresholds at higher values of p introducing instability at progressively higher values of m. However, these thresholds are probably of little practical importance. Whilst it is often said that the Tayler instability is dominated by the m = 1 mode because it is that mode which is in some sense the first to become unstable, this is misleading: only in the rather special circumstance of −1/2 <p< 1 and Ω_∥<ω_A is the m = 1 mode the only one present. Firstly, the second of those two conditions is rarely fulfilled in nature. Secondly, in the context of the fluctuating field envisaged in the Tayler-Spruit dynamo model, the field gradient p will have values all the way from p = −∞ to p = + ∞. It seems perfectly plausible that the dynamo could be dominated by the m = 0 mode, or indeed by m = 20. Roughly speaking, the instability is present where p is positive and absent where it is negative: the dynamo can therefore be expected to be intermittent in some sense.

Acknowledgments

The authors would like to thank Alfio Bonanno and Henk Spruit for fruitful discussions. The analysis and plots presented here were performed using VisIt (Visualization Tool) and Ipython notebook, and the bibliography was put together thanks to the Astrophysical Database ADS and Mendeley.

References

All Figures

	Fig. 1 Physical form of the different unstable azimuthal m modes present in a plasma column with a purely toroidal field. Below each mode the cross section of the co-moving surface is shown along with the cross section of the equilibrium state. The only axisymmetric mode is m = 0 and the only mode that preserves the column’s cross section is the m = 1 mode.
In the text

Fig. 2 Top: mixed poloidal-toroidal field configuration. If the toroidal field is too weak or indeed absent, the poloidal field is subject to an instability. Bottom: left, an equilibrium configuration with poloidal field loops threaded by the neutral line, and on the right-hand-side the kink-like instability present around the neutral line. Unstable displacements are restricted to be approximately perpendicular to the stratification. Figure taken from Braithwaite (2007).

In the text

Fig. 3 Representations of the first two vertical n modes of the unstable m = 1 mode. Top: purely toroidal magnetic field configuration, where the orientation of the magnetic symmetry axis is perpendicular to the rotation axis . It also shows what an n = 1 and n = 2 disk looks like as the red and blue box, respectively. Bottom: slipping disks for the n = 1 (left) and n = 2 (right) cases. The magnitude of the deformation Δl⊥, and the deformation velocity Δv⊥, suffered by the disks is different depending on the unstable mode. Small n modes present large deformations, whereas high n modes show small deformations, Δl⊥(n = 1) > Δl⊥(n = 2) and Δv⊥(n = 1) > Δv⊥(n = 2).

In the text

	Fig. 4 Predominantly toroidal magnetic field configuration in the interior of a star. The box shows a zoomed-in region containing the toroidal field lines around the magnetic axis and bold arrows showing the direction of the parallel and perpendicular rotation axes.
In the text

Fig. 5 a) Initial purely toroidal magnetic configuration. The region of interest is between 0 ≤ ϖ ≲ ϖ0 given that the field strength falls off at greater ϖ. Colored areas represent surfaces of equal distance to the axis. Field lines are also plotted; their density is related to the increase in the field strength proportional to the distance from the axis. b) Mixed poloidal-toroidal field configuration. The field lines show a spiral morphology, for this particular case (with p = 1) two spirals turns are completed in the box, i.e. toroidal Alfvén frequency is double the vertical Alfvén frequency ωA,φ = 2ωA,z.

In the text

	Fig. 6 Amplitude of the first four m modes as a function of time, extracted from the magnetic field δB_Ampm(t), and normalized to the wavelength corresponding to the n = 1 vertical mode. The azimuthal m = 1 mode, solid red line, is the only unstable mode.
In the text

	Fig. 7 Amplitude δB_Ampm,n(t) of the m = 1, n = 1 mode as a function of time for various rotation rates. Blue lines correspond to rotation rates below the Alfvén frequency, and red lines correspond to rotation rates above the Alfvén frequency.
In the text

	Fig. 8 Growth rate of the Tayler instability for the m = 1 and n = 1 mode as a function of the relative magnetic diffusivity η/ηcrit measured from the velocity field, v_Ampm,n(t).
In the text

Fig. 9 Contours of mean growth rate plotted as a function of normalized rotation rate nΩ∥ = Ω∥/ωA,φ and the relative magnetic diffusion nη = η/ηcrit. The white region represents stability, for at least 100τA. The color contours correspond to increasing growth rates of the instability, where black corresponds to small and blue to fast growth rates. The yellow dashed line corresponds to the relation nΩ∥ = 2nη + 1 separating the stable and unstable regions.

In the text

	Fig. 10 Growth rate of the Tayler instability, normalized to the Alfvén frequency, as a function of the perpendicular rotation rate, normalized to the predicted critical rotation rate for each n mode. The lines approximately on top of each other show the dependence of the damping on wavenumber is as predicted.
In the text

Fig. 11 Growth rate σ of the Tayler instability as a function of the solid-body rotation rate Ω∥, in the p = 2 case. Normalized local growth rates are given for three different distances from the axis, σ(ϖ). The dashed line corresponds to the predicted growth rate of the Tayler instability for slow rotators, σ = ωA,φ. The dotted-dashed line is the predicted behavior of the instability for fast rotators, $\mathit{\sigma }\mathrm{=}{\mathit{\omega }}_{\mathrm{A}\mathit{,\phi }}^{\mathrm{2}}\mathit{/}\mathrm{2}{\mathrm{\Omega }}_{\mathrm{\parallel }}$.

In the text

Fig. 12 Growth rate of the first two azimuthal m modes as a function of the pitch-number. The solid black line corresponds to the m = 0 mode, stable at very high and low pitch-numbers but with a peak in unstable growth near the critical pitch–number, Pn ≈ 1. The red dashed line corresponds to the m = 1 mode, unstable for very high pitch–numbers and stabilized near the critical pitch-number, Pn ≈ 1.

In the text