A&A 478, 1-8 (2008)
DOI: 10.1051/0004-6361:20077172

Stability of toroidal magnetic fields in rotating stellar radiation zones

L. L. Kitchatinov1,2 - G. Rüdiger1


1 - Astrophysikalisches Institut Potsdam, An der Sternwarte 16, 14482, Potsdam, Germany
2 - Institute for Solar-Terrestrial Physics, PO Box 291, Irkutsk 664033, Russia

Received 26 January 2007 / Accepted 14 October 2007

Abstract
Aims. Two questions are addressed as to how strong magnetic fields can be stored in rotating stellar radiation zones without being subjected to pinch-type instabilities and how much radial mixing is produced if the fields are unstable.
Methods. Linear equations are derived for weak disturbances of magnetic and velocity fields, which are global in horizontal dimensions but short-scaled in radius. The linear formulation includes the 2D theory of stability to strictly horizontal disturbances as a special limit. The eigenvalue problem for the derived equations is solved numerically to evaluate the stability of toroidal field patterns with one or two latitudinal belts under the influence of rigid rotation.
Results. Radial displacements are essential for magnetic instability. It does not exist in the 2D case of strictly horizontal disturbances. Only stable (magnetically modified) r-modes are found in this case. The instability recovers in 3D. The minimum field strength $B_{\rm min}$ for onset of the instability and radial scales of the most rapidly growing modes are strongly influenced by finite diffusion, the scales are indefinitely short if diffusion is neglected. The most rapidly growing modes for the Sun have radial scales of about 1 Mm. In the upper part of the solar radiation zone, $B_{\rm min} \simeq 600$ G. The toroidal field can exceed this value only marginally, for otherwise the radial mixing produced by the instability would be too strong to be compatible with the observed lithium abundance.

Key words: instabilities - magnetohydrodynamics (MHD) - stars: interiors - stars: magnetic fields - Sun: magnetic fields

1 Introduction

The old problem of the hydromagnetic stability of stellar radiation zones attains renewed interest in relation to Ap star magnetism (Braithwaite & Spruit 2004), the solar tachocline (Gilman 2005), and the transport processes in stably stratified stellar interiors (Barnes et al. 1999).

Which fields the stellar radiative cores can possess is rather uncertain. The resistive decay in radiation zones is so slow that primordial fields can be stored there (Cowling 1945). Whether the fields of 105 G that can influence g-modes of solar oscillations (Rashba et al. 2006) or whether still stronger fields that can cause neutrino oscillations (Burgess et al. 2003) can indeed survive inside the Sun mainly depends on their stability.

Among several instabilities to which the fields can be subjected (Acheson 1978), the current-driven (pinch-type) instability of toroidal fields (Tayler 1973) is probably the most relevant one because it proceeds via almost horizontal displacements. The radial motions in radiative cores are strongly suppressed by buoyancy. Watson (1981) estimated the ratio of radial (ur) to horizontal ($u_\theta$) velocities of slow (subsonic) motions in rotating stars as $
u_r/u_\theta \sim \Omega^2/N^2 ,
$where $\Omega $ is the basic angular velocity and N the buoyancy (Brunt-Väisälä) frequency. This ratio is small in radiative cores (Fig. 1). If the radial velocities are completely neglected, the stability analysis can be done in a 2D approximation with purely horizontal displacements (Watson 1981; Gilman & Fox 1997). Some radial motion still is, however, excited. How much radial mixing the Tayler (1973) instability produces is not well-known so far. As the radial mixing is relevant to the transport of light elements, a theory of the mixing compared with the observed abundances can help to restrict the amplitudes of the internal magnetic fields (Barnes et al. 1999). In the present paper the vertical mixing produced by the Tayler instability is estimated and then used to evaluate the upper limit on the magnetic field amplitude.

In this paper, the equations governing the linear stability of toroidal magnetic fields in a differentially rotating radiation zone are derived. They are solved for two latitudinal profiles of the toroidal field with one and with two belts in latitude but only for rigid rotation. The unstable modes are expected to have the longest possible horizontal scales (Spruit 1999). Accordingly, our equations are global in horizontal dimensions. They are, however, local in radius, i.e. the radial scales are assumed to be short, $kr \gg
1$ (k is radial wave number). The computations confirm that the most rapidly growing modes have $kr \sim 10^3$, but they are global in latitude. The derived equations reproduce the 2D approximation as a special limit. An exactly solvable case shows that the Tayler instability is missing in the 2D approximation. It recovers, however, in the 3D case. Finite diffusion is found important for the instability. It prevents the preference for the indefinitely short radial scales found for ideal MHD. The disturbances that first become unstable as the field strength grows have small but finite radial scales. For those scales, the minimum field producing the instability is strongly reduced by allowing for finite thermal conductivity. This field amplitude is about 600 G for the upper part of the solar radiation zone. Considering the transport of chemical species by the Tayler instability, we find that the field strength can only be slightly above this marginal value. Otherwise, the intensity of radial mixing would not be compatible with the observed lithium abundance.

  \begin{figure}
\par\includegraphics[width=8.8cm, height=6.0cm]{7172f1.eps}
\end{figure} Figure 1: The buoyancy frequency (2) in the upper part of the radiative core of the Sun after the solar structure model by Stix & Skaley (1990). The convection zone of the model includes the overshoot layer so that $N/\Omega _\odot $ is large immediately beneath the convection zone.
Open with DEXTER

2 The model

2.1 Background state and basic assumptions

The stability of the rotating radiation zone of a star containing magnetic field is considered. The field is assumed axisymmetric and purely toroidal, so it can be written as


 \begin{displaymath}{\vec B} = {\vec e}_\phi r\sin\theta\sqrt{\mu_0\rho}\
\Omega_{\rm A}\left( r,\theta\right)
\end{displaymath} (1)

in terms of the Alfvén angular frequency $\Omega _{\rm A}$. In this equation, $\rho$ is the density, $r,\theta ,\phi$ are the usual spherical coordinates, and  ${\vec e}_\phi$ is the azimuthal unit vector. Equation (1) automatically ensures that the toroidal field component vanishes - as it must - at the rotation axis. Centrifugal and magnetic forces are assumed weak compared to gravity, $g/r \gg \Omega^2$ and $g/r \gg \Omega^2_{\rm A}$. Deviation of the fluid stratification from spherical symmetry can thus be neglected.

The stabilizing effect of a subadiabatic stratification of the radiative core is characterized by the buoyancy frequency,


 \begin{displaymath}N^2 = \frac{g}{C_{\rm p}}\frac{\partial s}{\partial r} ,
\end{displaymath} (2)

where $s = C_{\rm v}\ln \left( P/\rho^\gamma\right)$ is the specific entropy of ideal gas. The ratio $N/\Omega$ is high in the radiative cores of not too rapidly rotating stars like the Sun (see Fig. 1).

The larger N, the more the radial displacements are opposed by the buoyancy force. Radial velocities should therefore be low. They are often neglected in stability analysis, which might be dangerous as certain instabilities may even be suppressed by the neglect (we shall see later how this indeed happens for the Tayler instability).

The consequences of neglecting the radial velocity perturbations, u'r, can be seen from the expression for (divergence-free) velocity, ${\vec u}'$, in spherical geometry in terms of the two scalar potentials for the poloidal, $P_{\rm v}$, and toroidal, $T_{\rm v}$, flows

 
                              $\displaystyle {\vec u}'$ = $\displaystyle \frac{{\vec e}_r}{r^2}\hat{\cal L}P_{\rm v}
- \frac{{\vec e}_\the...
...frac{1}{\sin\theta}\frac{\partial^2 P_{\rm v}}
{\partial r\partial\phi}\right),$  
$\displaystyle \hat{\cal L}$ = $\displaystyle \frac{1}{\sin\theta}\frac{\partial}{\partial\theta}
\sin\theta\fr...
...tial}{\partial\theta} +
\frac{1}{\sin^2\theta}\frac{\partial^2}{\partial\phi^2}$ (3)

(cf. Chandrasekhar 1961). Here and in the following, disturbances are marked by dashes. For zero radial velocity, the whole class of disturbances corresponding to the poloidal flow vanishes. Then the remaining toroidal flow can only produce toroidal magnetic disturbances so that, in the expression for the magnetic field perturbations (${\vec B}'$),
 
$\displaystyle {\vec B}' = \frac{{\vec e}_r}{r^2}\hat{\cal L}P_{\rm m}
- \frac{{...
...rac{1}{\sin\theta}\frac{\partial^2 P_{\rm m}}
{\partial r\partial\phi}\right) ,$     (4)

the poloidal field potential becomes zero, $P_{\rm m} = 0$.

It can be seen from Eq. (3) that the horizontal part of the poloidal flow can remain unchanged when the radial velocity is reduced if the radial scale of the flow is reduced proportionally. The disturbances can thus avoid the stabilizing effect of the stratification by decreasing their radial scale instead of losing their poloidal parts.

Our stability analysis is local in the radial dimension; i.e. we use Fourier modes $\exp ({\rm i}kr)$ with $kr \gg
1$. It will be confirmed a posteriori that the most unstable modes do indeed prefer short radial scales. The analysis remains, however, global in horizontal dimensions. Our formulation is very close to the thin layer approximation by Cally (2003) and Miesch & Gilman (2004), though we do not impose sharp boundaries and use other dependent variables. We also include finite diffusion and focus on marginal stability instead of growth rates for highly supercritical fields.

The instabilities of toroidal fields or differential rotation proceed via non compressive disturbances. The characteristic growth rates of the instabilities are low compared to p-modes frequencies. In the short-wave approximation, the velocity field can be assumed divergence-free, ${\rm div}~ {\vec u}' = 0$. Note that even a slow motion in a stratified fluid can be divergent if its spatial scale in the radial direction is not small compared to scale height. In the short-wave approximation (in radius), the divergency can, however, be neglected.

  
2.2 Equations

We start from the linearized equations for small velocity perturbations, i.e.
 
$\displaystyle \frac{\partial{\vec u}'}{\partial t}
+ \left({\vec u}\cdot\nabla\...
...ght){\vec B}\right)
=-\left(\frac{1}{\rho}\nabla P\right)' + \nu\Delta{\vec u}'$     (5)

magnetic field,

 \begin{displaymath}\frac{\partial{\vec B}'}{\partial t} =
\nabla\times\left( {\...
...+ {\vec u}'\times{\vec B} - \eta\nabla\times{\vec B}'\right) ,
\end{displaymath} (6)

and entropy,

 \begin{displaymath}\frac{\partial s'}{\partial t} + {\vec u}\cdot \nabla s' +
{\vec u}'\cdot\nabla s = \frac{C_{\rm p}\chi}{T}\Delta T' .
\end{displaymath} (7)

The background flow is rotation with angular velocity $\Omega
(r,\theta)$ and the mean magnetic field is the toroidal one of the form (1). Perturbations of velocity and magnetic field are expressed in terms of scalar potentials according to Eqs. (3) and (4). The identities
                        $\displaystyle r \left({\vec r}\cdot\nabla\times{\vec B}'\right)$ = $\displaystyle \hat{\cal L} T_{\rm m},\ \ \
r \left({\vec r}\cdot{\vec B}'\right) = \hat{\cal L}P_{\rm m},$  
$\displaystyle r \left({\vec r}\cdot\nabla\times{\vec u}'\right)$ = $\displaystyle \hat{\cal L} T_{\rm v},\ \ \
r^3\left({\vec r}\cdot\nabla\times\n...
...\hat{\cal{L}} + r^2\frac{\partial^2}{\partial r^2}\right)
\hat{\cal L}P_{\rm v}$  

are used to reformulate the equations in terms of the potentials (using the potentials helps to reduce the equation system). The radial component of Eq. (6) then gives the equation for the poloidal magnetic field and the radial components of curled Eqs. (5) and (6) give the equations for toroidal flow and toroidal magnetic field, respectively[*].

The perturbations are considered as Fourier modes in time, azimuth, and radius in the form of ${\rm exp}\left({\rm i}(-\omega t
+ m\phi + kr)\right)$. For an instability, the eigenvalue $\omega$should possess a positive imaginary part. Only the highest order terms in kr for the same variable were kept in the same equation.

When deriving the poloidal flow equation, the pressure term was transformed as

\begin{displaymath}{\vec r}\cdot\nabla\times\nabla\times\left(\frac{1}{\rho}
\n...
...left(\frac{1}{\rho^2}
(\nabla\rho)\times(\nabla P)\right)' =
\end{displaymath}


\begin{displaymath}\frac{\vec r}{C_{\rm p}}\cdot\nabla\times\left(\frac{1}{\rho}...
...imes\nabla s'\right) = \frac{g}{rC_{\rm p}}
\hat{\cal L}s' .
\end{displaymath}

This transformation neglects the pressure disturbances in the baroclinic forcing. More precisely, local thermal disturbances are assumed to occur at constant pressure so that $\rho'/\rho = -T'/T$or $s' = -C_{\rm p}\rho'/\rho$. This assumption is again justified by the incompressible nature of the perturbations. Another interpretation of this assumption is given by Acheson (1978), who assumed zero disturbances of total (including magnetic) pressure to involve magnetic buoyancy instability. In our derivations, assumptions of constant total or only gas pressure are identical because the effect of magnetic buoyancy appears in higher order in (kr)-1 compared to the terms kept in the equations to follow. Note that the pressure disturbances are neglected in the barocline term alone and the geostrophic balance is still possible.

In order to use normalized variables, the time is measured in units of $\Omega_0^{-1}$ ($\Omega_0$ is a characteristic angular velocity), the velocities are scaled in units of $r\Omega_0$, the normalized eigenvalue ( $\hat\omega$) is measured in units of $\Omega_0$, and the other normalized variables are

 
                                   A = $\displaystyle \frac{k}{\Omega_0 r^2\sqrt{\mu_0\rho}}\ P_{\rm m},\ \
B = \frac{1...
...ga_0 r^2\sqrt{\mu_0\rho}}\ T_{\rm m},\ \
V = \frac{k}{\Omega_0 r^2}\ P_{\rm v},$  
W = $\displaystyle \frac{1}{\Omega_0 r^2}\ T_{\rm v},\ \
S = \frac{{\rm i} k r g}{C_...
...Omega}{\Omega_0},\ \ \ \
\hat\Omega_{\rm A} = \frac{\Omega_{\rm A}}{\Omega_0} .$ (8)

Introducing the factors kr in the normalizations of poloidal potentials makes them equal by order of magnitude to the toroidal potentials, while in Eqs. (3) and (4) it was $P_{\rm v} \ll rT_{\rm v}$, $P_{\rm m} \ll
rT_{\rm m}$.

The equation for the poloidal flow then reads

 
                              $\displaystyle \hat\omega\left(\hat{\cal L}V\right)$ = $\displaystyle -\hat{\lambda}^2\left(\hat{\cal L}S\right)
- {\rm i}\frac{\epsilon_\nu}{\hat{\lambda}^2}\left(\hat{\cal L}V\right)$  
    $\displaystyle - 2\mu\hat\Omega\left(\hat{\cal L}W\right)
- 2\left(1-\mu^2\right...
...u}\frac{\partial W}{\partial\mu}
-2 m^2\frac{\partial\hat\Omega}{\partial\mu} W$  
    $\displaystyle + 2\mu\hat\Omega_{\rm A}\left(\hat{\cal L}B\right)
+ 2\left(1-\mu...
...\partial B}{\partial\mu}
+2 m^2\frac{\partial\hat\Omega_{\rm A}}{\partial\mu} B$  
    $\displaystyle - m\hat\Omega_{\rm A}\left(\hat{\cal L}A\right)
- 2m\frac{\partia...
...)\frac{\partial\hat\Omega_{\rm A}}
{\partial\mu}
\frac{\partial A}{\partial\mu}$  
    $\displaystyle + m\hat\Omega\left(\hat{\cal L}V\right)
+ 2m\frac{\partial\left(\...
...^2\right)\frac{\partial\hat\Omega}{\partial\mu}
\frac{\partial V}{\partial\mu},$ (9)

where $\mu = \cos\theta$, and

 \begin{displaymath}\hat\lambda = \frac{N}{\Omega_0 k r}
\end{displaymath} (10)

can be understood as special normalization for radial wavelength. The first term on the righthand side describes the stabilizing effect of the stratification. It vanishes for small $\hat\lambda$. Apart from this stabilizing buoyancy term, the wavelength is only present in diffusive terms. The second term of the RHS includes the action of finite viscosity,

 \begin{displaymath}\epsilon_\nu = \frac{\nu N^2}{\Omega_0^3 r^2}.
\end{displaymath} (11)

Similarly, we use

 \begin{displaymath}\epsilon_\eta = \frac{\eta N^2}{\Omega_0^3 r^2},\ \ \ \ \ \ \
\epsilon_\chi = \frac{\chi N^2}{\Omega_0^3 r^2}
\end{displaymath} (12)

for the diffusive parameters $\eta$ and $\chi$ below (note that the buoyancy frequency N drops off the diffusive terms of the resulting equations). The second and the following lines of Eq. (9) describe the influences of the basic rotation and the toroidal field. Only latitudinal derivatives of $\Omega $ and $\Omega _{\rm A}$ appear. All radial derivatives are absorbed by disturbances that vary on much shorter radial scales than $\Omega $or $\Omega _{\rm A}$. The complete system of five equations also includes the equation for toroidal flow
 
                           $\displaystyle \hat\omega\left(\hat{\cal L}W\right)$ = $\displaystyle - {\rm i}\frac{\epsilon_\nu}{\hat\lambda^2}\left(\hat{\cal L}W\ri...
...Omega\left(\hat{\cal L}W\right)
- m\hat\Omega_{\rm A}\left(\hat{\cal L}B\right)$  
    $\displaystyle - mW\frac{\partial^2}{\partial\mu^2}
\left(\left(1-\mu^2\right)\h...
...{\partial^2}{\partial\mu^2}
\left(\left(1-\mu^2\right)\hat\Omega_{\rm A}\right)$  
    $\displaystyle + \left(\hat{\cal L}V\right)\frac{\partial}{\partial\mu}
\left(\l...
...frac{\partial}{\partial\mu}
\left(\left(1-\mu^2\right)\hat\Omega_{\rm A}\right)$  
    $\displaystyle + \left(\frac{\partial}{\partial\mu}\left(\left(1-\mu^2\right)^2
...
...}\right) -
2\left(1-\mu^2\right)\hat\Omega\right)\frac{\partial V}{\partial\mu}$  
    $\displaystyle - \left(\frac{\partial}{\partial\mu}\left(\left(1-\mu^2\right)^2
...
...
2\left(1-\mu^2\right)\hat\Omega_{\rm A}\right)
\frac{\partial A}{\partial\mu},$ (13)

the equation for toroidal magnetic field
 
                           $\displaystyle \hat\omega\left(\hat{\cal L}B\right)$ = $\displaystyle - {\rm i}\frac{\epsilon_\eta}{\hat\lambda^2}\left(\hat{\cal L}B\r...
...al L}\left(\hat\Omega B\right)
- m\hat{\cal L}\left(\hat\Omega_{\rm A} W\right)$  
    $\displaystyle - m^2\frac{\partial\hat\Omega}{\partial\mu} A
- \frac{\partial}{\...
...)^2\frac{\partial\hat\Omega}{\partial\mu}
\frac{\partial A}{\partial\mu}\right)$  
    $\displaystyle + m^2\frac{\partial\hat\Omega_{\rm A}}{\partial\mu} V
+ \frac{\pa...
...partial\hat\Omega_{\rm A}}{\partial\mu}
\frac{\partial V}{\partial\mu}\right) ,$ (14)

the poloidal field equation

 \begin{displaymath}\hat\omega\left(\hat{\cal L}A\right) =
- {\rm i}\frac{\epsil...
...al L}A\right)
- m\hat\Omega_{\rm A}\left(\hat{\cal L}V\right)
\end{displaymath} (15)

and the entropy equation,

 \begin{displaymath}\hat\omega S = - {\rm i}\frac{\epsilon_\chi}{\hat\lambda^2}
S + m\hat\Omega S + \hat{\cal L}V .
\end{displaymath} (16)

In the simplest case of $\Omega=\Omega_{\rm A}=0$ and vanishing diffusivities, the only nontrivial solution of the above equations are short (in radius) gravity waves with

 \begin{displaymath}\omega = \frac{N}{kr}\sqrt{l(l+1)},\ \ l=1,2,...
\end{displaymath} (17)

Instabilities are thus only possible because of magnetic field or nonuniform rotation. It should be kept in mind that the above equations are only valid for $kr \gg
1$. We shall see that the maximum growth rates do indeed belong to the short radial scales.

  
2.3 2D approximation ( $\hat\lambda\gg 1$)

The ratio of $N/\Omega_0$ in radiation zones is so high (Fig. 1) that $\hat\lambda$ (10) can also be high in spite of $kr \gg
1$. Equation (9) in the leading order in $\hat\lambda$ then gives S = 0. Next, Eqs. (16) and (15) successively yield V=0 and A=0. Diffusive terms can also be neglected. The equation system reduces to two coupled equations for toroidal magnetic field and toroidal flow (Gilman & Fox 1997),
 
                          $\displaystyle \hat\omega\left(\hat{\cal L}W\right)$ = $\displaystyle m\hat\Omega\left(\hat{\cal L}W\right)
- m\hat\Omega_{\rm A}\left(\hat{\cal L}B\right)$  
    $\displaystyle - mW\frac{\partial^2}{\partial\mu^2}
\left(\left(1-\mu^2\right)\h...
...\partial^2}{\partial\mu^2}
\left(\left(1-\mu^2\right)\hat\Omega_{\rm A}\right),$  
$\displaystyle \hat\omega\left(\hat{\cal L}B\right)$ = $\displaystyle m\hat{\cal L}\left(\hat\Omega B\right)
- m\hat{\cal L}\left(\hat\Omega_{\rm A} W\right).$ (18)

Here the wave number k drops out and the equations describe 2D disturbances within decoupled spherical shells. The condition for validity of the 2D approximation is, therefore, $N/(\Omega k r) \gg
1$.

2.4 Equatorial symmetry of fluctuating and background fields

In the present paper the interaction of toroidal magnetic fields and the basic rotation is formulated only for rigid rotation. This is to single out the instabilities that are magnetic by origin (uniform rotation has a minimum energy for given angular momentum so that rotational energy cannot feed any instability). The results for the joint stability of differential rotation and magnetic field will be the subject of a separate publication.

For the toroidal field, two simple geometries are considered. First, the quantity $\Omega _{\rm A}$ is taken as a constant so that the toroidal field has only one belt symmetric with respect to the equator. Second, two magnetic belts are considered with equatorial antisymmetry, i.e. with a node of $B_\phi$ at the equator,

 \begin{displaymath}\Omega_{\rm A} = b\ \Omega \cos\theta ,
\end{displaymath} (19)

where b is the normalized field strength, b=1 means equipartition of rotational and magnetic energies.

Two types of equatorial symmetry are allowed by the equations for the linear disturbances in Sect. 2.2. Different dependent variables of the same symmetry type can, however, have different equatorial symmetries. As the toroidal flow disturbances are present in all instabilities allowed by our equations, the symmetry notations are related to the potential W of that flow. Standard notations, Sm and Am, are used for the modes with symmetric and antisymmetric W about the equator, respectively, and m is the azimuthal wave number. The symmetry of magnetic disturbances with the same symmetry notation differ, however, between the cases of one and two belts of background toroidal field. For the case of one belt, $\Omega_{\rm A} = {\rm const}$, it is

\begin{displaymath}{\rm Sm\ modes:}\ \ W,\ B\ \ {\rm symmetric}, \ \ V,\ A,\ S\ \ \
{\rm antisymmetric},
\end{displaymath}


\begin{displaymath}{\rm Am\ modes:}\ \ W,\ B\ \ {\rm antisymmetric}, \ \ V,\ A,\ S\ \ \
{\rm symmetric}.
\end{displaymath}

For two belts of toroidal fields of Eq. (19) the symmetry notations imply

\begin{displaymath}{\rm Sm\ modes:}\ \ W,\ A\ \ {\rm symmetric}, \ \ V,\ B,\ S\ \ \
{\rm antisymmetric},
\end{displaymath}


\begin{displaymath}{\rm Am\ modes:}\ \ W,\ A\ \ {\rm antisymmetric}, \ \ V,\ B,\ S\ \ \
{\rm symmetric}.
\end{displaymath}

In this latter case, Sm modes have vector field ${\vec B}'$ which is mirror-symmetric about equatorial plane (symmetric ${B_\phi}'$, Br' and antisymmetric ${B_\theta}'$) and mirror-antisymmetric flow ${\vec u}'$ (symmetric ${u_\theta}'$ and antisymmetric ${u_\phi}'$, ur') and the other way round for Am modes. The symmetry definition for two belts of toroidal field are the same as in Gilman et al. (2007).

3 Results and discussion

3.1 Nonexistence of 2D magnetic instabilities

The particular case where $\Omega $ and $\Omega _{\rm A}$ are both constant strongly simplifies the 2D Eqs. (18). The equations, then, attain constant coefficients and can be solved analytically. The solution in terms of Legendre polynomials $W,B\sim
P^m_l(\mu )$ gives the eigenfrequencies

 \begin{displaymath}\frac{\hat\omega}{m} = 1 - \frac{1}{l(l+1)} \pm
\sqrt{\hat\O...
...^2\left( 1 - \frac{2}{l(l+1)}\right)
+ \frac{1}{l^2(l+1)^2} }
\end{displaymath} (20)

of magnetically modified r-modes (Longuet-Higgins 1964; Papaloizou & Pringle 1978) describing stable horizontal patterns drifting in longitude. We shall see in Sect. 3.3.1 that the case of constant $\Omega $ and $\Omega _{\rm A}$shows Tayler instability in the 3D case.

The 2D approximation, however, misses the instability. Though the result is obtained for the particular case of constant $\Omega _{\rm A}$, it is most probably valid in general. Indeed, the toroidal field can be understood as consisting of closed flux tubes. Noncompressive disturbances conserve the volume and magnetic flux of the tubes so that the magnetic energy of a tube is proportional to the square of its length (Zeldovich et al. 1983). The 2D disturbances also conserve the area of a segment on a spherical surface encircled by a magnetic field line. The circular lines of background field have minimum length for any given encircled area. Any 2D disturbance increases the length, thereby increasing magnetic energy. Therefore, there is no possibility of feeding a 2D instability by magnetic energy release. A more formal proof of this statement can use the expression for magnetic disturbances

$\displaystyle {\vec B}' = {\vec\nabla}\times\left( {\vec\xi}\times{\vec
B}\right)$      

in terms of the displacement vector $\vec\xi$. It is easy to see that the contribution of ${\vec B}'$ to the magnetic energy density integrated over a spherical surface of constant r is positively definite for 2D displacements ${\vec\xi} =
{\vec\nabla}\times\left({\vec r} \psi\right)$ ($\psi$ is an arbitrary scalar function).

The non-existence of Tayler instability in 2D clearly shows that radial displacements are necessary for the instability. Therefore, radial mixing is unavoidable. Transport of chemical species by the Tayler instability will be estimated in Sect. 3.3.2.

The conclusion about the absence of 2D instabilities is of course restricted to the case of rigid rotation. Rotational shear can excite its own 2D instability (Watson 1981; Dziembowski & Kosovichev 1987; Garaud 2001) fed by rotational energy release, and magnetic field can catalyze the process (Gilman & Fox 1997).

Gilman et al. (2007) find that 3D instabilities in the solar tachocline with strong toroidal fields differ only slightly from their 2D counterparts. This is probably because of large differential rotation in the tachocline. We shell see that the results for rigid rotation differ strongly between the 2D and 3D cases. The amount of differential rotation required for the difference to fade is still unknown.

3.2 Ideal fluids, $\chi = \eta = \nu = 0$

We proceed by discussing numerical solutions for Eqs. (9), (13)-(16) for 3D disturbances. Consider first the case of ideal fluids with vanishing diffusion. Constant $\Omega _{\rm A}$ give the simplest realization of the Tayler instability. The instability criterion for nonaxisymmetric disturbances (Goossens et al. 1981) reads

 \begin{displaymath}\frac{\partial}{\partial\mu}\ln\left(B_\phi^2\right) <
\frac...
... - m^2}{\mu\left(1-\mu^2\right)} ,\ \ \
{\rm for}\ \mu\geq 0.
\end{displaymath} (21)

It is satisfied with constant $\Omega _{\rm A}$ for m=1 and close to the pole. The axisymmetric instability with m=0 does not exist in this case.

From Sect. 3.1 we know that the instability can appear only in 3D formulation when radial displacements are allowed. Figure 2 gives the resulting stability map. Two features are remarkable: (i) rotation suppresses the instability so that subeqiupartition fields, $\Omega _{\rm A} < \Omega $, are stable, and (ii) the unstable disturbances for the smallest field producing the instability have indefinitely short radial scales, $\hat\lambda \rightarrow 0$. Buoyancy frequency, N, enters the equations for linear disturbances of Sect. 2.2 in combination (10) with radial wave number, k. Therefore, the disturbances can avoid the stabilizing effect of stratification by reducing their radial scale. This is the strategy the Tayler instability follows. Another possibility would be zero entropy disturbances (the $\hat\lambda$-parameter (10) enters the equations as a factor of entropy perturbation). Entropy perturbations can be avoided with 2D disturbances only. Tayler instability cannot use this opportunity because it has no 2D counterpart.

  \begin{figure}
\par\includegraphics[width=7.0cm]{7172f2.eps}
\end{figure} Figure 2: The stability map for constant $\Omega _{\rm A}$ and zero diffusivities. There is no instability for weak fields with $\Omega _{\rm A} < \Omega $. The threshold field strength for the instability increases with increasing vertical wavelength $\hat\lambda$.
Open with DEXTER

Rotational quenching of Tayler instability is a controversial issue. Spruit (1999) finds that rotation modifies the instability but does not switch it off. Cally (2003), however, concluded that the polar kink instability (m=1 mode of Tayler instability) is totally suppressed if $\Omega _{\rm A}$ is smaller than $\Omega $near the pole. Simulations of Braithwait (2006) also show that the instability does not develop when the toroidal field is below a critical value roughly equal to the equipartition level of $\Omega_{\rm A}\simeq\Omega$.


  \begin{figure}
\par\includegraphics[width=7.0cm]{7172f3.eps}
\end{figure} Figure 3: The small $\hat\lambda$ part of the stability map for two belts of toroidal magnetic field (19). Rotation does not suppress instability in this case. Parameter b is the polar value of the ratio $\Omega _{\rm A}/\Omega $.
Open with DEXTER

The controversy can be resolved by observing that the rotational quenching is sensitive to details of the toroidal field profile. Already, Pitts & Tayler (1985) suggested that suppression of the instability by rotation may be exceptionally strong when the ratio of Alfvén velocity to rotation velocity is uniform, $\Omega_{\rm A}/\Omega = {\rm const.}$, but rotation is unlikely to lead to complete stability for the general configuration of the toroidal field. Computations for the case of two belts of the toroidal field of Eq. (19) confirm this expectation. The stability map of Fig. 3 shows that the disturbances with sufficiently small radial scales remain unstable when the (normalized) toroidal field b is small[*].

The high sensitivity of the Tayler instability in ideal fluids to details of the $\Omega_{\rm A}(\mu )$ profile is mainly of academic interest. Figures 2 and 3 show that the smallest b producing the instability corresponds to indefinitely small radial scales. Finite diffusion, therefore, must be included. We shall see that the strong difference between the cases of one and two belts of toroidal field fades when diffusion is allowed for.

3.3 Finite diffusion

From now on the values

 \begin{displaymath}\epsilon_\nu = 2\times 10^{-10}, \ \ \ \ \ \ \ \
\epsilon_\eta = 4\times 10^{-8}, \ \ \ \ \ \ \ \
\epsilon_\chi = 10^{-4},
\end{displaymath} (22)

which are characteristic of the upper part of the solar radiative core, are used for the diffusion parameters (11) and (12). The relation $\chi\gg\eta\gg\nu$ of Eq. (22) is typical of stellar radiation zones.

For ideal fluids the Tayler instability operates with extremely small radial scales. When there is no stabilizing stratification, however, finite vertical scales are preferred (Arlt et al. 2007). The thermal conductivity decreases the stabilizing effect of the stratification and reduces the critical field strengths for the instability if the (normalized) radial scale $\hat\lambda$ is not too small.

  \begin{figure}
\par\includegraphics[width=7.0cm]{7172f4.eps}
\end{figure} Figure 4: The same as in Fig. 2 but with the finite diffusivities (22). The critical field strengths for the onset of the instability are strongly reduced compared to ideal MHD.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=7.0cm]{7172f5.eps}
\end{figure} Figure 5: The normalized growth rate $\hat\sigma$ as function of the magnetic field amplitude for the S1 mode and for $\hat\lambda = 0.6$ (where the line of Fig. 4 has a minimum). The dotted lines for strong, $\Omega _{\rm A} > \Omega $, and weak, $\Omega _{\rm A} < \Omega $, fields give approximations by the linear, $\hat\sigma \sim \hat\Omega_{\rm A}$, and parabolic, $\hat\sigma \sim \hat\Omega^2_{\rm A}$, laws, respectively.
Open with DEXTER

  
3.3.1 Fields with equatorial symmetry
The reduction of the critical field for onset of the instability is exceptionally strong for the case of uniform $\Omega _{\rm A}$. The stability map for diffusion parameters of Eq. (22) is shown in Fig. 4. The minimum field producing the instability is reduced by more than two orders of magnitude compared to ideal fluids (Fig. 2). The characteristic minima in Fig. 4 correspond to small, $\hat\lambda \la 1$, but finite vertical scales.

For strong fields ( $\Omega _{\rm A} > \Omega $) the basic rotation and diffusion are not important and there only characteristic frequency to scale the growth rates is $\Omega _{\rm A}$. The dependence of Fig. 5 does indeed approach the $\sigma \sim
\Omega_{\rm A}$ relation ($\sigma$ is the growth rate) in the strong-field limit. The growth rate drops by almost four orders of magnitude when $\hat\Omega_{\rm A}$ is reduced below 1. In the weak-field regime it is $\sigma \sim \Omega^2_{\rm A}/\Omega$(Spruit 1999). The growth rates for weak fields ( $\Omega _{\rm A} < \Omega $), where the instability exists due to finite diffusion, are rather small. Nevertheless, the growth rate of, e.g., $\hat\sigma = 10^{-5}$ means the e-folding time of about 1000 yr for the Sun, which is very short compared to evolutionary time scales. We shall see that the instability of weak fields is actually too vigorous to be compatible with observed lithium abundance.

The structures of the unstable modes differ between strong and weak field regimes. For the weak fields, the growing disturbance pattern is distributed over the entire sphere. For strong fields, it is much more concentrated at the poles (Fig. 6) but still remains global in latitude.

  \begin{figure}
\par\includegraphics[width=8.0cm]{7172f6_1.eps} \includegraphics[width=8.0cm]{7172f6_2.eps}
\end{figure} Figure 6: Toroidal field lines of the most rapidly growing eigenmodes of Fig. 5 for weak background toroidal field, $\Omega _{\rm A} = 0.1\Omega $ ( upper panel), and strong field, $\Omega _{\rm A} = 10\Omega $ ( lower panel).
Open with DEXTER

  
3.3.2 Fields with equatorial antisymmetry
Now a toroidal field with two belts of opposite polarity (19) is considered with the diffusivity set (22). Such a field geometry can result from the action of differential rotation on dipolar poloidal fields. Figure 7 shows the corresponding stability map. It is now similar to the map of Fig. 4 for one field belt. Strong dependence on the $\Omega_{\rm A}(\mu )$-profile is no longer observed. Axisymmetric instability does not exist for profile (19). Instability is again only found for the kink mode of m=1.
  \begin{figure}
\par\includegraphics[width=7.0cm]{7172f7.eps}
\end{figure} Figure 7: Stability map for the field model (19) with two belts and equatorial antisymmetry.
Open with DEXTER

In physical values, one finds for the upper part of solar radiation zone (with $\rho\simeq$ 0.2 g/cm3)

 \begin{displaymath}B_\phi \simeq 10^5 b \ \ \ \ {\rm G}
\end{displaymath} (23)

and

 \begin{displaymath}\lambda \simeq 10\hat\lambda\ {\rm Mm},
\end{displaymath} (24)

so that from Fig. 7 follows a critical magnetic field for the instability slightly below 600 G. The mode that first becomes unstable if the field exceeds this critical value has a vertical wavelength between 1 and 2 Mm. For larger field strengths, of course, there is a range of unstable wavelengths. The maximum growth rates appears, however, at wavelengths $\la$1 Mm (Fig. 8).
  \begin{figure}
\mbox{
\includegraphics[width=4.1cm]{7172f8_1.eps} \hspace{0.2cm}
\includegraphics[width=4.3cm]{7172f8_2.eps} }
\end{figure} Figure 8: Growth rates in units of $\Omega $ for b=0.01 ( left) and b=0.1 ( right). The highest rates exist for $\hat\lambda \;\buildrel
<\over{\scriptstyle\sim}\; 0.1$ independent of the magnetic field amplitude, all for most rapidly growing S1 modes.
Open with DEXTER

The flow field of the instability also mixes chemical species in radial direction. Such an instability can thus be relevant to the radial transport of the light elements (Barnes et al. 1999). The effective diffusivity, $D_{\rm T} \simeq u'\ell$ (u' and $\ell$ are rms velocity and correlation length in radial direction) can be roughly estimated from our linear computations assuming that the instability saturates when mixing frequency $u'/\ell$approaches the growth rate $\sigma \simeq u'/\ell$ and $\ell \simeq
\lambda/2$. With Eq. (24), this yields

 \begin{displaymath}D_{\rm T} \simeq 7 \times 10^9\ \hat\sigma\ \
{\rm cm}^2~{\rm s}^{-1} ,
\end{displaymath} (25)

where $\hat\sigma$ is the normalized growth rate given in Fig. 9. For the range of 0.01 < b <0.2 where the plot is closely approximated by the parabolic law $\hat\sigma = 0.1b^2$, Eq. (25) can be rewritten in terms of $B_\phi$(Eq. (23)) as

 \begin{displaymath}D_{\rm T} \simeq 7\times 10^4 \left(\frac{B_\phi}{1\ {\rm kG}}
\right)^2\ \ \ \ {\rm cm}^2~{\rm s}^{-1}.
\end{displaymath} (26)

As is known, diffusivities in excess of 104 cm2 s-1 in the upper radiative core are not compatible with the observed solar lithium abundance. Hence, the toroidal field amplitude can only slightly exceed the marginal value of about 600 G. In our (simplified) formulation, the observed solar lithium abundance seems to exclude[*] any concept of hydromagnetic dynamos driven by Tayler instability in the upper radiation zone of the Sun. We should not forget, however, that the superrotation ( $\partial \Omega/\partial r >0$) at the bottom of the convection zone in the equatorial region acts to stabilize. A mild stabilizing effect can also be expected from the compositional gradient due to gravitational settling of chemical species (Richard et al. 1996; Vauclair 1998). The critical field amplitudes for Tayler instability may, thus, be somewhat higher than the computed 600 G. Even this value is, however, above the toroidal field strengths obtained in magnetic models for the laminar solar tachocline (cf. Rüdiger & Kitchatinov 2007). The axisymmetric tachocline models are, therefore, expected to be stable to nonaxisymmetric disturbances.
  \begin{figure}
\par\includegraphics[width=7.0cm]{7172f9.eps}
\end{figure} Figure 9: The growth rate as a function of the toroidal field amplitude for $\hat{\lambda} = 0.1$. The dotted line shows the parabolic approximation $\hat\sigma = 0.1b^2$. The scale on the right gives the radial diffusivity of chemical species estimated after (25).
Open with DEXTER

4 Summary

The linear stability of toroidal magnetic field in rotating stellar radiation zones is analyzed assuming that the vertical scale of the fluctuations is short compared to the local radius ("short-wave approximation''). The analysis is global in horizontal dimensions. Stability computations confirm that the most rapidly growing perturbations have short radial scales: $kr \sim 10^3$.

We have shown that the pinch-type instability of toroidal fields require nonvanishing radial displacements. The instability does not exist in a 2D approximation with zero radial velocities. The maximum amplitude of stable toroidal magnetic fields for the Sun we found is about 600 G. This value results only for rigid rotation. It will most probably increase if the stabilizing influence of the positive radial gradient of $\Omega $ in the equatorial region of the tachocline is included in the model.

The field strength in the upper part of the solar radiative interior can only marginally exceed the resulting critical values. Otherwise the instability would produce too strong a radial mixing of light elements. After our results, all the axisymmetric hydromagnetic models of the solar tachocline (Rüdiger & Kitchatinov 1997; MacGregor & Charbonneau 1999; Garaud 2002; Sule et al. 2005; Kitchatinov & Rüdiger 2006) have stable toroidal fields. On the contrary, the strong fields $\sim$105 G, which are able to modify g-modes, or even stronger fields that may be relevant to neutrino oscillations (Burgess et al. 2003) are strongly unstable with e-folding times shorter than one rotation period.

The 3D computations of joint instabilities of toroidal fields and differential rotation (Gilman & Fox 1997; Cally 2003; Rüdiger et al. 2007) will be discussed in a separate paper. Another tempting extension is to include of the poloidal field. The field can be important in view of the very short vertical scales of the unstable modes.

Acknowledgements
This work was supported by the Deutsche Forschungsgemeinschaft and by the Russian Foundation for Basic Research (project 05-02-04015).

References

 

Copyright ESO 2008