---
author:
- 'Kyle Augustson$^1$, Stéphane Mathis$^{1,2}$, Antoine Strugarek$^{3,1}$'
bibliography:
- 'tayler.bib'
title: On the Stability of a General Magnetic Field Topology in Stellar Radiative Zones
---

Motivation
==========

The radiative core of main-sequence low-mass stars and the radiative envelope of main-sequence massive stars likely host a fossil magnetic field [@neiner15; @braithwaite15].  This field is a remnant of the field built during the star’s birth and subsequently reinforced during convective phases of its evolution toward the main-sequence [@alecian13].  In particular, massive stars with an observed magnetic field typically possess a non-axisymmetric oblique magnetic dipole or a similarly simple magnetic field geometry [@moss90; @walder12; @wade16].  If a comparison is drawn between the stably-stratified regions of massive and low-mass stars [@strugarek11], given their hydrodynamic similarity, such non-axisymmetric magnetic fields may also exist within these regions for low-mass stars.  Fossil magnetic fields have also been proposed as an important source of angular momentum transport and mixing across the Hertzsprung-Russell diagram [e.g., @gough98; @heger05; @mathis05].  So, constraining the stability of a large class of magnetic fields within the convectively-stable, radiative regions is important for characterizing their influence on the transport of angular momentum over evolutionary timescales, understanding their topology that is observed at the surfaces of intermediate and high mass stars, and their consequences for the local stellar environment [e.g., @petit12]. 

As an example of how such fossil fields can form, the process of freezing out the magnetic field as the star evolves along the pre-main-sequence is depicted in Figure \[fig:relax\], where gravitational contraction decreases the radius of the star.  As the star slowly collapses, the gradual increase of the density and temperature deep within the star tends to lower the opacity, which eventually leads transition from convective heat transport to diffusive heat conduction at the edge of the core.  During these phases, rotationally-constrained convective motions will generate the magnetic field.  In contrast, once convection has halted in the stably-stratified layers, the field will undergo a slow Ohmic decay if the field has a stable configuration or a fast Alfvénic decay if it is unstable.  One way to distinguish which of these decay paths the magnetic field will take is to assess its stability to small (linear) displacements of fluid elements.  If the growth rate of those small perturbations is real and positive, the magnetic field undergoes the Tayler instability. 

![ image](Figure1.pdf){width="0.85\linewidth"}

The Tayler Instability
======================

The stability of axisymmetric magnetic field configurations within a quiescent, stably-stratified medium have been understood for quite some time, with @tayler73 addressing toroidal field configurations and @markey73 poloidal configurations.  The typical local instabilities arising in those systems are shown in Figure \[fig:tayler\], where there are three situations shown.  The equilibrium situation occurs when the magnetic field has no associated current (e.g., if it is potential), or if the Tayler stability criterion are met.  Yet for sufficiently large currents or sufficiently strong Lorentz forces, two other instabilities can be excited: the axisymmetric $m=0$ varicose instability, or the $m=1$ kink instability.  The latter of which grows most rapidly when excited.  Furthermore, such analyses indicate that only certain mixed (poloidal and toroidal) configurations of axisymmetric magnetic fields are stable within the radiative regions of stars [@tayler80; @braithwaite09].  The presence of rotation modifies the stability characteristics of these axisymmetric systems in that it tends to further stabilize them through the Coriolis force [@pitts85].  The precise form of the equilibrium states of the mixed-morphology magnetic fields has been considered extensively in both non-rotating and rotating systems [@prendergast56; @braithwaite06; @braithwaite08; @duez10a; @duez10b; @duez11; @braithwaite13; @emeriau15]. 

Generalizing the Tayler Instability
===================================

The class of local stability analysis established by @tayler73 can be generalized to global-scale geometries as can be found in @goossens80.  However, axisymmetric magnetic fields are not the final story in the study of the Tayler instability.  Rather, the analysis can be extended to configurations with both non-axisymmetric magnetic fields and differential rotation as shall be shown in an upcoming paper [@augustson16].  The resulting stability criteria are assessed here.  Such criteria help to restrict the number of magnetic field configurations that are possible within the stable regions of low-mass stars, thereby limiting the routes of angular momentum transport in the radiative interior and means of interaction with the dynamo-generated magnetic fields established in their overlying convective layers.\
The linearized equation of motion under the Cowling approximation [@cowling41] for a fluid element in a general, but non-rotating, coordinate system is

$$\begin{aligned}
    \rho\frac{\partial^2 {\boldsymbol{\xi}}}{\partial t^2} &= \frac{1}{4\pi}\left[\left({\boldsymbol{\nabla}\!\!\boldsymbol{\times}\!\!}{\delta{\boldsymbol{B}}}\right){\!\!\boldsymbol{\times}\!\!}{\boldsymbol{B}}+\left({\boldsymbol{\nabla}\!\!\boldsymbol{\times}\!\!}{{\boldsymbol{B}}}\right){\!\!\boldsymbol{\times}\!\!}\delta{\boldsymbol{B}}\right]-\delta\rho{\boldsymbol{\nabla}}{\Phi}-{\boldsymbol{\nabla}}{\delta P}, \label{eqn:leom}
  \end{aligned}$$

where ${\boldsymbol{\xi}}$ is the displacement, ${\boldsymbol{B}}$ the magnetic field, $\rho$ the density, $P$ the pressure, $\Phi$ the gravitational potential.  The Eulerian perturbations $\delta$ of those quantities follow directly from the continuity, pressure, and induction equations as

$$\begin{aligned}
    \delta\rho &= -{\boldsymbol{\nabla}\!{\!\boldsymbol{\cdot}\!}}{\left(\rho{\boldsymbol{\xi}}\right)},\\
    \delta P &= -{\boldsymbol{\xi}}{\!\boldsymbol{\cdot}\!\!\boldsymbol{\nabla}}P - \gamma P {\boldsymbol{\nabla}\!{\!\boldsymbol{\cdot}\!}}{{\boldsymbol{\xi}}},\\
    \delta {\boldsymbol{B}}&= {\boldsymbol{\nabla}\!\!\boldsymbol{\times}\!\!}{\left({\boldsymbol{\xi}}{\!\!\boldsymbol{\times}\!\!}{\boldsymbol{B}}\right)}. 
  \end{aligned}$$

Therefore, one has that

$$\begin{aligned}
    \rho\frac{\partial^2{\boldsymbol{\xi}}}{\partial t^2} &= \mathcal{F}[{\boldsymbol{\xi}}; \rho, P, \Phi, {\boldsymbol{B}}] \nonumber \\
    &= {\boldsymbol{\nabla}\!{\!\boldsymbol{\cdot}\!}}{\left(\rho{\boldsymbol{\xi}}\right)}{\boldsymbol{\nabla}}{\Phi}+{\boldsymbol{\nabla}}{\left[{\boldsymbol{\xi}}{\!\boldsymbol{\cdot}\!\!\boldsymbol{\nabla}}P + \gamma P {\boldsymbol{\nabla}\!{\!\boldsymbol{\cdot}\!}}{{\boldsymbol{\xi}}}\right]} \label{eqn:full} \\
    &+\left[\left({\boldsymbol{\nabla}\!\!\boldsymbol{\times}\!\!}{{\boldsymbol{\nabla}\!\!\boldsymbol{\times}\!\!}{\left({\boldsymbol{\xi}}{\!\!\boldsymbol{\times}\!\!}{\boldsymbol{B}}\right)}}\right){\!\!\boldsymbol{\times}\!\!}{\boldsymbol{B}}+\left({\boldsymbol{\nabla}\!\!\boldsymbol{\times}\!\!}{{\boldsymbol{B}}}\right){\!\!\boldsymbol{\times}\!\!}\left({\boldsymbol{\nabla}\!\!\boldsymbol{\times}\!\!}{\left({\boldsymbol{\xi}}{\!\!\boldsymbol{\times}\!\!}{\boldsymbol{B}}\right)}\right)\right], \nonumber
  \end{aligned}$$

where $\gamma$ is the ratio of specific heats. 

As was shown in @bernstein58 and to decide upon the stability of this system, one can consider simple solutions of the form ${\boldsymbol{\xi}}= Re\left[{\ensuremath{\pmb{\psi}}}(\mathbf{x}) \exp{\left(i\omega t\right)}\right]$, for which the equation of motion yields $-\omega^2 \rho{\ensuremath{\pmb{\psi}}}= \mathcal{F}\left[{\ensuremath{\pmb{\psi}}}\right]$. This is not general since one has not yet proven that these basis functions form a complete set on the Hilbert space for the Eulerian system.  Yet it can be shown that the vector function $\mathcal{F}$ is self-adjoint, a proof of which will be reserved for the upcoming paper [@augustson16].  With a properly defined inner product for the solutions ${\boldsymbol{\xi}}$, one can see that the dispersion relationship for a general displacement in an arbitrary coordinate system is given by

$$\begin{aligned}
    \omega^2 &= -\frac{\langle{\ensuremath{\pmb{\psi}}},\mathcal{F}\left[{\ensuremath{\pmb{\psi}}}\right]\rangle}{\langle\rho{\ensuremath{\pmb{\psi}}},{\ensuremath{\pmb{\psi}}}\rangle}
             = -\int {\ensuremath{\pmb{\psi}}}^{*}{\!\boldsymbol{\cdot}\!}\mathcal{F} dV \left[\int \rho {\ensuremath{\pmb{\psi}}}^{*}{\!\boldsymbol{\cdot}\!}{\ensuremath{\pmb{\psi}}}dV\right]^{-1},  \nonumber \\
             &= 2\Delta W \left[\int \rho {\ensuremath{\pmb{\psi}}}^{*}{\!\boldsymbol{\cdot}\!}{\ensuremath{\pmb{\psi}}}dV\right]^{-1},
  \end{aligned}$$

where the integral is taken over the region of interest, which for stars are their convectively-stable zones. 

The displacement is unstable if the change in the potential energy of the system ($\Delta W$) is negative.  In general, one finds that this energy can be split into three parts as $\Delta W = \Delta W_L + \Delta W_B + \Delta W_P$, with $\Delta W_L$ being the work due to Lorentz forces, $\Delta W_B$ being the work due to buoyancy, and with $\Delta W_P$ being the pressure work.  To find general classes of magnetic fields that are stable in a given radiative region, one needs to be able to find an expression for the work that can be minimized.  This is possible within the context of separable coordinate systems.  For this work, the spherical coordinate system is used.  Therefore, for compactness and expedience, ${\ensuremath{\pmb{\psi}}}$, $\rho$, $P$, $\Phi$, and ${\boldsymbol{B}}$ are projected onto the spherical spin vector harmonics (SVH).  The SVH are a complete orthonormal set of vector functions that are formed from specific combinations of the spherical harmonics and their derivatives [@varshalovich88].  In particular, they correspond to the joint eigenstates of the angular momentum and spin-1 operators.  Relative to the RST basis [@rieutord87], vector operations such as the dot and cross products are simpler to perform on SVH-projected vector-valued functions.  An explicit representation of the SVH in terms of scalar spherical harmonics is as follows:

$$\begin{aligned}
    \mathbf{Y}_{\ell,1}^m &= \frac{1}{\sqrt{\left(\ell+1\right)\left(2\ell+1\right)}}\left[-\left(\ell+1\right){\hat{\mathbf{r}}}+ r{\boldsymbol{\nabla}}{}\right] Y_\ell^m, \label{eqn:svhp1}\\
    \mathbf{Y}_{\ell,0}^m  &= \frac{-i r}{\sqrt{\ell\left(\ell+1\right)}}{\hat{\mathbf{r}}}{\!\!\boldsymbol{\times}\!\!}{\boldsymbol{\nabla}}{Y_\ell^m}, \label{eqn:svhp0}\\
    \mathbf{Y}_{\ell,-1}^m &= \frac{1}{\sqrt{\ell\left(2\ell+1\right)}}\left[\ell{\hat{\mathbf{r}}}+ r{\boldsymbol{\nabla}}{}\right] Y_\ell^m \label{eqn:svhm1},
  \end{aligned}$$

![ Tayler instabilities in a cylindrical current channel.  From left to right: shows an equilibrium configuration for an azimuthal magnetic field, a varicose (m=0) instability, and a kink type (m=1) instability.  Field lines are marked with arrows.[]{data-label="fig:tayler"}](Figure2.pdf){width="\linewidth"}

where the second lower index $\nu$ on the $\mathbf{Y}_{\ell,\nu}^m$ indicates the corresponding spin-1 basis vector.  Indeed, when the work integrands are expanded on the SVH basis, one can show that

$$\begin{aligned}
    &\Delta W_L = \!\!\int_{r_b}^{r_t} \!\!  dr r^2\sum_{i=1}^4\sum_{\substack{\ell,m,\nu,\mu,\lambda \\ \ell_i, m_i,
    \nu_i}} w^{\mu,\lambda}_{\substack{\ell,m,\nu \\ \ell_i, m_i, \nu_i}},
  \end{aligned}$$

with,

$$\begin{aligned}
    &\sum_{i=1}^4 w^{\mu,\lambda}_{\substack{\ell,m,\nu \\ \ell_i, m_i, \nu_i}}\!\!\!=\nonumber \\ 
    &{_{\lambda,\mu}L_{\substack{\ell_1,\ell_2\\ \nu_1, \nu_2}}^{m_1, m_2}}\!\!\left(\psi_{\ell_3,\nu_3}^{m_3}B_{\ell_4,\nu_4}^{m_4}\right)\psi_{\ell,\nu}^{*m} 
    \mathcal{J}_{\substack{\ell_3, m_3, \nu_3 \\ \ell_4, m_4, \nu_4}}^{\ell_2, m_2, \lambda}
    \mathcal{J}_{\substack{\ell_1, m_1, \nu_1 \\ \ell_2, m_2, \nu_2}}^{\ell, m, \nu},
  \end{aligned}$$

where each $\ell$ ranges from zero to infinity, each $m$ ranges between $-\ell$ and $\ell$, and where $\nu$, $\mu$, and $\lambda$ range between $-1$ and $1$. The integral is taken between radii $r_b$ and $r_t$, which demark the bottom and top boundaries of the radiatively stable region.  The $\mathcal{J}$ coefficients arise from the projection of the cross products of the SVH basis vectors in Equation \[eqn:full\] back onto the basis.  The $L$ symbol is a function of radius that arises from the Lorentz force, and thus it is a second-order differential operator involving the radial functions of the displacement and the magnetic field as

$$\begin{aligned}
    &_{\lambda,\mu}L_{\substack{\ell_1,\ell_2 \\ \nu_1, \nu_2}}^{m_1, m_2}\!\left(\psi_{\ell_3,\nu_3}^{m_3} B_{\ell_4,\nu_4}^{m_4}\right)\!  = \nonumber \\
    &\!-\frac{B_{\ell_1, \nu_1}^{m_1}}{4\pi}\left[E_{\nu_2,\lambda}^{\ell_2, m_2}\frac{\partial^2}{\partial r^2} 
      + \frac{F_{\nu_2,\lambda}^{\ell_2, m_2}}{r}\frac{\partial}{\partial r}
      + \frac{G_{\nu_2,\lambda}^{\ell_2, m_2}}{r^2}\right]\psi_{\ell_3,\nu_3}^{m_3} B_{\ell_4,\nu_4}^{m_4}\mathcal{I}_{\lambda,\mu} \nonumber\\
    &+\frac{1}{4\pi}\left[D_{\nu_1,\lambda}^{\ell_1,m_1} \frac{\partial B_{\ell_1,\lambda}^{m_1}}{\partial r} 
      + C_{\nu_1,\lambda}^{\ell_1,m_1}\frac{B_{\ell_1,\lambda}^{m_1}}{r}\right]\times \nonumber \\
    &\qquad\qquad\qquad\left[D_{\nu_2,\mu}^{\ell_2,m_2}\frac{\partial}{\partial r}+\frac{C_{\nu_2,\mu}^{\ell_2,m_2}}{r}\right]\psi_{\ell_3,\nu_3}^{m_3}B_{\ell_4,\nu_4}^{m_4},\end{aligned}$$

where the coefficient matrices $C$, $D$, $E$, $F$, and $G$ describe the projection of the curl and double curl operators onto the spin vector harmonic basis, and $\mathcal{I}$ is the unit tensor.  Assuming that the star is spherically symmetric, namely that the gradient of the gravitational potential is only in the radial direction, then one has that $g=-\partial_r{\Phi}$. So, tackling the buoyancy work integral, it can be seen that

$$\begin{aligned}
    &\Delta W_B = \!\!\!\!\sum_{\substack{\ell, m, \ell_1 \\ \ell_2, m_2,\nu_2}} \!\!\!\int_{r_b}^{r_t} \!\!  dr \frac{\left(-1\right)^{m_2} g r^2}{2\ell+1}
    \left(\sqrt{\ell+1}\psi_{\ell,1}^{*m}- \sqrt{\ell}\psi_{\ell,-1}^{*m}\right) \nonumber \\
    &\left[\sqrt{\ell+1}\left(\frac{\partial}{\partial r}+\frac{\ell+2}{r}\right)\mathcal{K}_{\ell_1;\substack{\ell,1\\ \ell_2,\nu_2}}^{-m, m_2}\right.\nonumber\\
    &\qquad\qquad\left.-\sqrt{\ell}\left(\frac{\partial}{\partial r}-\frac{\ell-1}{r}\right)\mathcal{K}_{\ell_1;\substack{\ell,-1\\ \ell_2,\nu_2}}^{-m, m_2}\right]
    \rho_{\ell_1}^{m_2-m} \psi_{\ell_2,\nu_2}^{m_2}. 
  \end{aligned}$$

Similarly, the pressure work integral can be identified as

$$\begin{aligned}
    &\Delta W_P = 
      \!\!\int_{r_b}^{r_t} \!\!  dr r^2\!\!\!\sum_{\substack{\ell,m,\ell_1 \\ \ell_2,m_2}}\frac{1}{\sqrt{\left(2\ell+1\right)\left(2\ell_2+1\right)}}\nonumber\\
    &\left[\sqrt{\ell+1}\left(\frac{\partial}{\partial r} +\frac{\ell+2}{r}\right) \psi_{\ell,1}^{*m+m_2}\right.  \nonumber \\
    &\qquad\qquad+\!\left.\sqrt{\ell}\left(-\frac{\partial}{\partial r} +\frac{\ell-1}{r}\right)\psi_{\ell,-1}^{*m+m_2} \right]\nonumber\\
    &\left\{\sum_{\nu_1}\Bigg(\psi_{\ell_1,\nu_1}^{m}\left[\!\sqrt{\ell_2+1}\left(\frac{\partial }{\partial r} 
      - \left.\frac{\ell_2}{r}\right) \mathcal{K}_{\ell;\substack{\ell_1,\nu_1\\ \ell_2,1}}^{m, m_2} \right.  \right.  \right.\nonumber\\
    &\left. -  \sqrt{\ell_2}\left(\frac{\partial }{\partial r} + \frac{\ell_2+1}{r}\right) 
      \mathcal{K}_{\ell;\substack{\ell_1,\nu_1\\\ell_2,-1}}^{m, m_2}\right] P_{\ell_2}^{m_2}\Bigg)\nonumber\\
    &+\left(-1\right)^{m+m_2}\gamma \mathcal{H}_{\substack{\ell,\ell_1,\ell_2\\ m, m_2}}P_{\ell_1}^m
    \left[\sqrt{\ell_2+1}\left(\frac{\partial}{\partial r} +\frac{\ell_2+2}{r}\right) \psi_{\ell_2,1}^{m_2}\right.  \nonumber \\
    &\left.+\sqrt{\ell_2}\left(-\frac{\partial}{\partial r}+\frac{\ell_2-1}{r}\right)\psi_{\ell_2,-1}^{m_2} \right]\Bigg\}. 
  \end{aligned}$$

Here, $\mathcal{H}$ and $\mathcal{K}$ are coefficients related to the 3-j and 6-j symbols that arise from integrals over products of SVH that are then either projected onto the scalar spherical harmonics, which, along with $\mathcal{J}$, are closely related to those defined in @varshalovich88 and @strugarek13. 

Conclusions
===========

With the expanded form $\Delta W$ in hand, one can then find conditions under which the system with a chosen general magnetic field is linearly stable or unstable to an arbitrary displacement by integrating the terms with radial derivatives of ${\ensuremath{\pmb{\psi}}}$ by parts and then explicitly minimizing the radial integrals with respect to ${\ensuremath{\pmb{\psi}}}$. This will be demonstrated more completely in an upcoming paper [@augustson16].  As applied to stellar radiative zones, this will permit the determination of the stability of certain classes of magnetic fields that have a broad range of non-axisymmetric components.  In particular, it may be possible to assess why the magnetic field configuration where the magnetic axis of symmetry is oblique to the rotation axis of the star is the most commonly observed. 

Acknowledgments {#acknowledgments .unnumbered}
===============

[K. C. Augustson and Stéphane Mathis acknowledge support from the ERC SPIRE 647383 grant.  A. Strugarek acknowledges support from the Canadian Institute of Theoretical Astrophysics (National Fellow) and from the Canadian Natural Sciences and Engineering Research Council.]{}