Angular momentum transport in massive stars and natal neutron star rotation rates

TL;DR

This paper investigates how magnetic and rotational instabilities in massive stars influence angular momentum transport, predicting neutron star spin rates that align with observations and explaining the rarity of rapidly rotating young neutron stars.

Contribution

It demonstrates that Tayler instability-driven angular momentum transport significantly reduces core spin, leading to more accurate predictions of neutron star rotation periods.

Findings

Tayler instability effectively transports angular momentum in massive stars.

Predicted neutron star periods are 50-200 ms, consistent with observations.

Additional processes may increase initial neutron star spin rates.

Abstract

The internal rotational dynamics of massive stars are poorly understood. If angular momentum (AM) transport between the core and the envelope is inefficient, the large core AM upon core-collapse will produce rapidly rotating neutron stars (NSs). However, observations of low-mass stars suggest an efficient AM transport mechanism is at work, which could drastically reduce NS spin rates. Here we study the effects of the baroclinic instability and the magnetic Tayler instability in differentially rotating radiative zones. Although the baroclinic instability may occur, the Tayler instability is likely to be more effective for AM transport. We implement Tayler torques as prescribed by Fuller et al. 2019 into models of massive stars, finding they remove the vast majority of the core's AM as it contracts between the main sequence and helium-burning phases of evolution. If core AM is conserved during core-collapse, we predict natal NS rotation periods of $P_{\rm NS} \approx 50-200 \, {\rm ms}$, suggesting these torques help explain the relatively slow rotation rates of most young NSs, and the rarity of rapidly rotating engine-driven supernovae. Stochastic spin-up via waves just before core-collapse, asymmetric explosions, and various binary evolution scenarios may increase the initial rotation rates of many NSs.

Figures

Tables

Table 1. Table 1: List of variables used in main text and appendices

$\overrightarrow{b}$	Perturbed magnetic field
$\overrightarrow{B}$	Magnetic field
$B_0$	Background magnetic field magnitude
$c_P$	Specific heat at constant pressure
$c_{\text{s}}$	Adiabatic sound speed
$c_T$	Isothermal sound speed
$c_V$	Specific heat at constant volume
$\hat{e}$	Unit vector in $\overrightarrow{e}$ direction
$f$	Coriolis frequency
$F$	Angular overlap integral
$g$	Gravitational acceleration
$G$	Gravitational constant
$H$	Pressure scale height
$H_{\lambda }$	Hough function
$I_{\text{NS}}$	Neutron star moment of inertia
$j$	Specific angular momentum
$J$	Angular momentum
$k$	Total wave number
$(k_r,l,m)$	Wave number in $(r,\theta ,\varphi )$ direction
$(k_x,k_y,k_z)$	Wave number in $(x,y,z)$ direction
$k_{⟂}$	Wave number in horizontal direction
$\overrightarrow{L}$	Perturbed Lorentz force
$\ell$	Wave length
${\ell }_r$	Wave length in radial direction
${\ell }_h$	Wave length in horizontal direction
$m$	Mass
$m_{\text{p}}$	Proton mass
$M$	Total mass
$M_{\odot }$	Solar mass
$N$	Brunt-Väisälä frequency
$N_{\text{eff}}$	Effective Brunt-Väisälä frequency
$N_T$	Thermal part of Brunt-Väisälä frequency
$N_{\theta }$	Baroclinic characteristic frequency
$N_{\mu }$	Composition part of Brunt-Väisälä frequency
$P$	Pressure
$P_{\text{NS}}$	Rotation period of neutron star
$q$	Dimensionless shear
$(r,\theta ,\varphi )$	Spherical coordinates
$r_c$	Baroclinic critical radius
$S$	Entropy
$t$	Time
$T$	Temperature
$\overrightarrow{v}$	Velocity
$(x,y,z)$	Local Cartesian coordinates
$\alpha$	Torque parameter
$\gamma$	Growth rate
${\Gamma }_1$	Adiabatic index
$\delta Q$	Eulerian perturbation of $Q$
$\eta$	Magnetic diffusivity
$\kappa$	Thermal diffusivity
$\lambda$	Eigenvalue of Hough function
$\mu$	Composition
$\overrightarrow{\xi }$	Lagrangian displacement vector
$\rho$	Density
$\nu$	Rotation parameter
${\nu }_{\text{AM}}$	Effective viscosity
$\omega$	Angular oscillation frequency
${\omega }_A$	Alfvén frequency
${\omega }_{\text{c}}$	Critical Alfvén frequency
$\Omega \left(\overrightarrow{\Omega }\right)$	Angular rotation frequency (vector)

Equations

$$r\sin\theta\bigg{(}\cos\theta\frac{\partial}{\partial r}-\frac{\sin\theta}{r}\frac{\partial}{\partial\theta}\bigg{)}\Omega^{2}=-\frac{1}{\rho^{2}}(\nabla\rho\times\nabla P)_{\phi}$$

$$N_{\theta }^2\equiv -\frac{g}{\rho r}\frac{\partial \rho }{\partial \theta }$$

$$N_{\theta }^2=2q{\Omega }^2\sin \theta \cos \theta ,$$

$$q=\frac{d\ln \Omega }{d\ln r}$$

$$\delta Q\propto e^{-i\omega t}.$$

$$-{\omega }^2\vec{\xi }=-\frac{1}{\rho }\nabla \delta P-\frac{g}{\rho }\delta \rho \hat{r}+2i\omega \vec{\Omega }\times \vec{\xi }$$

$$\delta \rho +\nabla \cdot (\rho \vec{\xi })=0.$$

$$\frac{\delta \rho }{\rho }=\frac{1}{{\widetilde{c_{\mathrm{s}}}}^2}\frac{\delta P}{\rho }+{\xi }_r\frac{{\tilde{N}}^2}{g}+{\xi }_{\theta }\frac{{\widetilde{N_{\theta }}}^2}{g}$$

$$\delta Q\propto\delta Q(\theta)\exp\bigg{[}i\bigg{(}\int k_{r}dr+m\phi-\omega t\bigg{)}\bigg{]}\,.$$

$$\gamma =\mathrm{Im}(\omega )=-\frac{1}{2}\frac{N_T^2}{N^2}\kappa k_r^2\pm \frac{1}{2}\mathrm{sign}(k_r)\frac{N_{\theta }^{2}F}{(N^2\lambda {)}^{\frac{1}{2}}}$$

$$k_r^2\simeq \frac{\lambda N^2}{r^2{\omega }^2}$$

$$q{\Omega }^2{\Omega }_{\max }^2\gtrsim N_T^{2}N\frac{\kappa }{r^2}.$$

$$\gamma =\mathrm{Im}(\omega )=-\frac{N_T^2\lambda }{2\kappa k_r^4{r}^2}.$$

$$x$$

$$y$$

$$z$$

$$\delta Q\propto\exp\big{[}i(k_{x}x+k_{y}y+k_{z}z-\omega t)\big{]}\,.$$

$$f\equiv 2\Omega \cos \theta .$$

$$\beta \equiv \frac{\partial f}{\partial y}=\frac{2\Omega \sin \theta }{r_0}.$$

$$f{v}_y=\frac{1}{\rho }\frac{\partial }{\partial x}\delta P,$$

$$f{v}_x=-\frac{1}{\rho }\frac{\partial }{\partial y}\delta P,$$

$$g\delta \rho =-\frac{\partial }{\partial z}\delta P,$$

$$\frac{\partial v_x}{\partial x}+\frac{\partial v_y}{\partial y}+\frac{\partial v_z}{\partial z}=0.$$

$$\gamma =\mathrm{Im}(\omega )=-\frac{N_T^2}{N^2}\kappa k^2+\frac{N_{\theta }^2}{2\Omega \cos \theta }\frac{k_x}{k_z}.$$

$$\delta Q\propto\exp\big{[}i(k_{r}r+m\phi+l\theta-\omega t)\big{]}\,.$$

$$g\frac{\delta\rho}{\rho}=\bigg{[}N_{T}^{2}\bigg{(}1+i\frac{\kappa k_{r}^{2}}{\omega}\bigg{)}^{-1}+N_{\mu}^{2}\bigg{]}\xi_{r}\,,$$

$$-i\omega \vec{b}=\nabla \times (-i\omega \vec{\xi }\times \vec{B})-\eta k_r^2\vec{b},$$

$$-{\omega }^2\vec{\xi }=-\frac{1}{\rho }\nabla \delta P-g\frac{\delta \rho }{\rho }+2i\omega \Omega \hat{z}\times \vec{\xi }+\vec{L}=0,$$

$$L_r$$

$$L_{\theta }$$

$$L_{\varphi }$$

$${\omega }_{\mathrm{A}}=\frac{B_{\varphi }}{\sqrt{4\pi \rho r^2}}$$

$$\begin{split}&\tilde{\omega}^{6}-\tilde{\omega}^{4}(4\tilde{\Omega}^{2}\cos^{2}\theta+A_{t}+A_{\mu}+2m^{2}+2hk+h^{2})\\ -&\tilde{\omega}^{3}\cdot 8\tilde{\Omega}m+\tilde{\omega}^{2}\big{[}m^{2}A_{t}+m^{2}A_{\mu}+m^{2}-4m^{2}\cos^{2}\theta\\ +&2hk(4\tilde{\Omega}^{2}\cos^{2}\theta+A_{\mu}+m^{2})+h^{2}(4\tilde{\Omega}^{2}\cos^{2}\theta+A_{t}+A_{\mu})\big{]}\\ +&\tilde{\omega}\cdot 8mkh\tilde{\Omega}\cos^{2}\theta-hkm^{2}A_{\mu}\\ +&i\tilde{\omega}^{5}(2h+k)-i\tilde{\omega}^{3}\big{[}k(4\tilde{\Omega}^{2}\cos^{2}\theta+A_{\mu}+2m^{2})\\ +&2h(4\tilde{\Omega}^{2}\cos^{2}\theta+A_{t}+A_{\mu}+m^{2})+h^{2}k\big{]}\\ -&i\tilde{\omega}^{2}\cdot 8(k+h)\tilde{\Omega}\cos^{2}\theta m\\ +&i\tilde{\omega}\big{[}k(m^{2}A_{\mu}+m^{2}-4m^{2}\cos^{2}\theta)\\ +&hm^{2}(A_{t}+A_{\mu})+h^{2}k(4\tilde{\Omega}^{2}\cos^{2}\theta+A_{\mu})\big{]}\\ =&\,0\,.\end{split}$$

Full text

1 Introduction

Internal rotation rates of massive stars are weakly constrained and poorly understood. While we have a basic understanding of massive stellar evolution and the compact objects that are produced upon core-collapse, we do not even have a zeroth order understanding of their internal rotation rates. This ignorance prevents the development of a detailed understanding of massive star evolution, as rotational effects (e.g., rotational mixing, Maeder & Meynet 2000) can have a substantial impact on a star’s evolution. More importantly, the core rotation rate of a massive star is crucial for determining the outcome of core-collapse. A variety of transients such as gamma-ray bursts (GRBS, Woosley 1993; MacFadyen & Woosley 1999; Metzger et al. 2011) and superluminous supernovae (Kasen & Woosley, 2009; Nicholl et al., 2017) are thought to be engine-powered events that tap into the huge rotational energy reservoir ($E_{\rm rot}\gtrsim 10^{52}$ erg) that can be supplied by a rapidly rotating neutron star (NS) or black hole (BH). Theoretical predictions for core rotation rates depend on the unknown efficiency of AM transport. Zero AM transport implies that nearly all compact objects will be born maximally rotating (Heger et al., 2000), while instantaneous AM tranpsort (i.e., rigid stellar rotation) implies compact objects will be essentially non-rotating. Both models conflict with observations, indicating that strong but imperfect AM coupling mechanisms are at work.

Observational constraints on core rotation arise primarily from compact object rotation rates. While some young neutron stars (e.g., the Crab pulsar) have natal spin periods of $P_{0}\sim 20\,{\rm ms}$ (Kaspi & Helfand, 2002), the majority of NSs have longer natal spin periods of $P_{0}\sim 50-100\,{\rm ms}$ (Faucher-Giguère & Kaspi, 2006; Popov et al., 2010; Popov & Turolla, 2012; Gullón et al., 2014), and some young NSs have natal spin periods as long as $P_{0}\sim 400{\rm ms}$ (Gotthelf et al., 2013). Spin rates of accreting BHs in X-ray binaries (XRBs) can be estimated from accretion disk modeling, and current estimates suggest a broad range of spin-rates ($0.1\lesssim a\lesssim 1$) (Miller & Miller, 2015). Many of these spin rates may have been increased by post-BH formation mass transfer, while some may reflect natal spin rates, but these can be affected by binary interactions (Valsecchi et al., 2010; Qin et al., 2018b, a). BH mergers detected by LIGO (Abbott et al., 2016, 2017a, 2017c, 2017b) are now providing the first spin rates of non-accreting BHs, and the low values of $\chi_{\rm eff}$ measured so far may indicate low BH natal spins (The LIGO Scientific Collaboration et al., 2018; Zaldarriaga et al., 2018).

A crucial piece of observational evidence arises from asteroseismic rotation rates measured in low-mass ($M\lesssim 3\,M_{\odot}$) stars. These measurements have been performed for main-sequence, red-giant branch, horizontal branch, and white dwarf stars (Beck et al., 2012; Mosser et al., 2012; Deheuvels et al., 2014, 2015; Benomar et al., 2015; Hermes et al., 2017; Gehan et al., 2018). Their message is clear: stellar cores and compact objects rotate orders of magnitude slower than they would in the absence of AM transport, and the rotation rates are slower than predicted by most AM transport models (Cantiello et al., 2014). Hence, massive stellar cores may also rotate slower than prior predictions. Perhaps the most realistic predictions come from Heger et al. (2005) and Wheeler et al. (2015) upon inclusion of the Tayler-Spruit dynamo (Spruit, 2002) and/or magneto-rotational instability (MRI) into stellar models, and they predict typical NS rotation rates $P_{\rm NS}\sim 10\,{\rm ms}$, at the fast end of the distribution inferred from radio pulsars. However, the same physics implemented in low-mass stars predict rotation rates roughly an order of magnitude faster than observed in red giant cores and white dwarfs (Cantiello et al., 2014).

Recently, Fuller et al. (2019) argued that the Tayler instability (e.g., Spruit 1999) can grow to larger saturated amplitudes and provide more efficient AM than predicted by Spruit (2002), and they provided an updated prescription for the “Tayler torques" arising from the instability. Models incorporating this prescription produced a good match with asteroseismic measurements of stars on the main sequence, red giant branch, red clump, and white dwarf cooling track. In this work, we extend the models of Fuller et al. (2019) to the massive star regime ($M\gtrsim 10\,M_{\odot})$, and we predict the AM contained in the core as it approaches core-collapse, and the corresponding rotation rate of a NS or BH if AM is conserved during collapse. Incorporating Tayler torques into massive stellar models, we find significantly slower core rotation rates (by roughly one order of magnitude) than prior predictions. We also examine the previously neglected baroclinic instability, and we provide a more rigorous derivation of the dispersion relation for Tayler instabilities.

In Sections 2 and 3, we investigate the baroclinic and Tayler instabilities. Section 4 describes our stellar models and AM transport results. In Section 5, we discuss the implications of our results for massive stellar evolution and massive stellar death, and we conclude in Section 6. Most of the detailed calculations regarding baroclinic and Tayler instabilities are presented in the appendices.

In what follows we refer to many variables whose definitions can be found in Table 1.

2 Baroclinic Instability

A well-known hydrodynamic instability that might be important for AM transport is the baroclinic instability, which has only been examined in a few previous astrophysical works (e.g., Tassoul & Tassoul, 1982; Fujimoto, 1987, 1988; Kitchatinov, 2014), but has been extensively studied in Earth atmospheric and oceanic contexts (see review in Pedlosky, 1992). Except in the case of a cylindrical rotation profile, differential rotation generally displaces the surfaces of constant density from isobars, creating a baroclinic stellar structure with a density/entropy gradient along isobars. The baroclinic instability is sourced by the potential energy released by latitudinal exchange between high and low density fluid elements, and is analogous to convection in the horizontal direction driven by the latitudinal entropy gradient. Despite its ubiquity, the baroclinic instability has frequently been neglected in astrophysics due to a series of papers claiming that it does not usually occur in stars (Knobloch & Spruit, 1982, 1983; Spruit et al., 1983; Spruit & Knobloch, 1984; Zahn & Zinn-Justin, 1993). Below we show why the conclusion of these works are mostly erroneous, but why baroclinic instabilities often remain unimportant.

Under the assumption of an axisymmetric background, baroclinicity is related to differential rotation by (Kitchatinov, 2014)

[TABLE]

where $(r,\theta,\phi)$ are spherical coordinates, $\Omega$ is the angular frequency, $\rho$ is the density and $P$ is the pressure. In our analysis we ignore the effect of centrifugal distortion, such that $\nabla P$ is in the radial direction (accounting for centrifugal distortion, one can simply redefine the radial coordinate to be perpendicular to isobars). Baroclinic instability occurs when $\partial\rho/\partial\theta\neq 0$ which can occur when there is differential rotation. We define

[TABLE]

to be the characteristic buoyancy frequency for baroclinic instability, where $g$ is the gravity. In our analysis we shall assume shellular rotation, ($\partial\Omega/\partial\theta=0$), such that equation (1) reduces to

[TABLE]

where

[TABLE]

is the dimensionless shear.

We start by analyzing the stability of oscillation modes, assuming the time dependence of each perturbation variable $\delta Q$ is

[TABLE]

where $\omega$ is the complex oscillation frequency. In the rotating frame, we have the perturbed momentum equation

[TABLE]

where $\vec{\xi}$ is the perturbed displacement and $\hat{r}$ is the unit vector in the radius direction, and $\delta P$ and $\delta\rho$ are the Eulerian pressure and density perturbations. We have made the Cowling approximation by neglecting the perturbation to the gravitational acceleration. We also have the perturbed continuity equation

[TABLE]

We show in appendix A that the energy equation with the baroclinic term and thermal diffusion is

[TABLE]

with all the notations made clear in appendix A.

The set of equations (6), (7) and (8) are difficult to solve, even for linear theory, because the Coriolis and baroclinic terms break the spherical symmetry of the problem. The growing modes of interest are often low frequency and nearly incompressible gravity modes and Rossby modes, so here we consider two approximations that are often used for such modes, namely, the traditional approximation and the geostrophic approximation.

2.1 Traditional Approximation

Here we perform a linear stability analysis for modes computed using the traditional approximation (e.g., Bildsten et al. 1996; Lee & Saio 1997), decomposing the spatial dependence of the mode as

[TABLE]

We perform a local analysis and make a WKB approximation $k_{r}\gg 1/H,1/r$, since the radial wavenumbers of low frequency modes are large due to the strong buoyancy force.

We show in appendix B.1 that when thermal diffusion is small (such that $k^{2}\kappa\ll\omega$), this approximation leads to the growth rate

[TABLE]

where $N$ is the Brunt-Väisälä frequency, $N_{T}$ is its thermal component, $\kappa$ is the thermal diffusivity, $\lambda$ is the horizontal eigenvalue of a Hough function, and $F$ is an angular overlap integral (see appendix B.1 for details). The first term in equation 10 is due to thermal diffusion and always causes damping in radiative regions with $N_{T}^{2}$ and $N^{2}$ positive. The second term is due to baroclinicity and causes damping or driving depending on the sign of $k_{r}$. Hence, only modes propagating in one direction will be driven. In the absence of thermal diffusion, typical growth rates are thus $\gamma\sim q\Omega^{2}/N$, remarkably independent of radial wavenumber $k_{r}$ (though with some dependence on horizontal wavenumber through the $F$ and $\lambda$ terms). This explains the numerical results of Kitchatinov (2014), who found a nearly constant growth rate as a function of $k_{r}$, except at large $k_{r}$ where thermal diffusion is important. The transition to damped modes occurs when $NN_{\theta}^{2}\lesssim N_{T}^{2}\kappa k_{r}^{2}$.

We pause to reconcile the result of equation 10 with previous works such as Spruit et al. (1983) and Zahn & Zinn-Justin (1993), which claim that baroclinic instability requires very large shear and is unlikely to occur in stars, even in the absence of diffusive effects. Prior claims relied on a theorem (Charney & Stern, 1962) from geophysics based on potential vorticity. This theorem involves global analysis, whereas equation 10 is derived based on a local analysis. If one integrates equation 10 for a standing mode to determine its global growth rate, the baroclinic growth term is equal to zero, because the inwardly and outwardly propagating component of the standing mode (i.e., the positive and negative signs of $k_{r}$) cancel each other out. Only non-WKB terms contribute to the global growth rate, and these terms require large shear for growth (they typically require $q\Omega^{2}\gtrsim N^{2}$, a criterion similar to the Richardson criterion for instability of shear flows).

Nonetheless, the theorem from Charney & Stern (1962) is not generally applicable, as is well-known in geophysics (see review in Pedlosky 1992), and it should not be considered a necessary condition for baroclinic instability. The subtle reason is that the potential vorticity argument is based on a global calculation in which boundary terms can be substituted after assuming no divergences within an integration domain. However, such divergences are ubiquitous at critical layers where the local frequency of a mode goes to zero, and its wavenumber diverges. The critical layers allow a growing mode propagating in one direction to be absorbed without propagating back to become damped. For this reason, growing baroclinic modes are ubiquitous in atmoshperic contexts and may generally exist in stars as well.

However, because critical layers are required for the excitation of baroclinic modes, only modes of low enough frequency can be excited. Consider a mode whose local frequency (measured in a corotating frame) at a critical layer is $\omega(r_{c})=0$. Then the mode’s local frequency at a different layer in the star, whose angular rotation rate is different by $\Delta\Omega$, is $\omega(r)=m\Delta\Omega$. Using the g-mode dispersion relation

[TABLE]

we find that mode growth requires $N_{\theta}^{2}(m\Delta\Omega)^{2}r^{2}\gtrsim\lambda^{3/2}N_{T}^{2}N\kappa\,$. The maximum possible value of $\Delta\Omega$ is simply the maximal rotation rate at any point in the stars, $\Omega_{\rm max}$. Since $m\leq\ell\lesssim\lambda^{1/2}$, and $N_{\theta}^{2}\sim q\Omega^{2}$ baroclinic instability requires

[TABLE]

In Sun-like stars and red giants, $\Omega/N$ is typically of order $10^{-3}$ or smaller. In many situations, we find that equation 12 is not satisfied anywhere in the radiative regions of stars, so critical layers will not exist for modes with frequencies high enough to avoid radiative damping. For such standing modes, baroclinic growth will typically not occur, so the baroclinic instability will not operate. An important exception is rapidly rotating stars, and stars with strong shear layers near the surface due to ongoing accretion (e.g., Piro & Bildsten 2007), where the baroclinic instability is more likely to occur.

When the thermal diffusion is large, we take the limit $\kappa\rightarrow\infty$, and show in Appendix B.1 that the growth rate reduces to

[TABLE]

This is always a damping mode so the instability does not exist at large thermal diffusivity. However, our analysis does not include compositional diffusivity or viscosity. Interestingly, when included, they give rise to new classes of instabilities such as the ABCD instability (Spruit et al., 1983). We find it unlikely that such instabilities contribute significantly to AM transport because they vastly over predict the differential rotation in the radiative core of the Sun.

2.2 Geostrophic Approximation

An alternative method to study the problem locally is to apply the geostrophic approximation (e.g., Spruit & Knobloch 1984). To do this, we restrict ourselves to low-frequency modes and neglect higher order terms of $\omega$. We perform the analysis in a thin box near a given point $(r_{0},\theta_{0},\phi_{0})$, establishing a local coordinate frame $(x,y,z)$ with the unit vectors in their directions $(\hat{e}_{x},\hat{e}_{y},\hat{e}_{z})$ corresponding to the original unit vectors $(\hat{e}_{\phi},-\hat{e}_{\theta},\hat{e}_{r})$ respectively (Figure 1). Then within the box we have

[TABLE]

For local modes, we can express the perturbation variables as

[TABLE]

Since we neglect the high order terms of $\omega$, the most significant restoring term in the horizontal direction is the Coriolis force, so we define a Coriolis frequency $f$ as

[TABLE]

This is known as the beta-plane approximation111 The name is from the definition of the Rossby parameter $\beta$

$\beta\equiv\frac{\partial f}{\partial y}=\frac{2\Omega\sin\theta}{r_{0}}\,.$

In our ‘thin box’ $f$ varies linearly with the polar angle, i.e. $f=f_{0}+\beta y$.

in geophysical fluid dynamics. (Note that in geophysics the latitude $\varphi$ is often used instead of the polar angle $\theta$ we use here, so the Coriolis parameter can be equivalently defined as $f=2\Omega\sin\varphi$.) Due to the strong stratification in the radial direction, the radial displacements are expected to be much smaller than the horizontal ones, and the radial component of the Coriolis force can be neglected. Then the equation of motion (6) reduces to

[TABLE]

We also use the incompressible approximation

[TABLE]

When the thermal diffusion is small but not negligible, we show in Appendix B.2 that, with the energy equation, the growth rate is

[TABLE]

where $k^{2}=k_{x}^{2}+k_{y}^{2}+k_{z}^{2}$ is the total wave number. This is similar to the results we found for the traditional approximation, but notably different because the growth rate now depends on both the sign and the value of $k_{z}$. The reason is the different dispersion relation between the gravito-inertial modes studied in the traditional approximation, and the Rossby modes studied in the geostrophic approximation. The growth rate of equation 23 is generally smaller than that of equation 10 by a factor $\sim k_{x}/k_{z}\ll 1$. Hence, thermal diffusion will also typically suppress the growth of Rossby modes growing via the baroclinic instability.

3 Tayler Instability

Due to their complexity, magnetohydrodynamic instabilities in stellar interiors have received less attention than hydrodynamic instabilities, but they can be more important for AM transport for two reasons. First, they are less quenched by diffusion because the magnetic diffusivity $\eta$ is much smaller than the thermal diffusivity $\kappa$ in stellar interiors. Second, hydrodynamic instabilities can only transport AM via a Reynolds stress $T\sim 4\pi\rho\langle v_{r}v_{\phi}\rangle$ that requires radial motion, which is suppressed by large buoyancy forces. In contrast, magnetic instabilities can produce a Maxwell stress $T\sim B_{r}B_{\phi}$ that does not necessarily require radial motion. As shown by Spruit (2002) , the magnetic Tayler instability (Tayler, 1973; Spruit, 1999) could play an especially important role in AM transport process.

Prior studies of the Tayler instability have been restricted to a polar region of the star where the geometry becomes Cartesian (e.g., Acheson & Gibbons 1978; Spruit 1999; Zahn et al. 2007). Denissenkov & Pinsonneault (2007) argued that the effective horizontal length scale of the instability must be much smaller than $r$ when one considers the spherical geometry of a star. 222The argument by Denissenkov & Pinsonneault 2007 confuses the perturbation wavelength $\ell_{r}\sim k_{r}^{-1}$ with the Lagrangian displacement $\xi_{r}$, and it is invalid for constraining the vertical wavelength $\ell_{r}\sim\omega_{\rm A}\ell_{h}/N$ or horizontal wavelength $\ell_{h}\sim r$. The arguments of Denissenkov & Pinsonneault (2007) are essentially arguments about non-linear geometric and Coriolis terms, which do not affect the linear instability, but may affect its non-linear evolution and saturation. Here we present the dispersion relation for Tayler modes at an arbitrary latitude in the star. We will find that the dispersion relation is unaltered apart from some factors of $\cos^{2}\theta$, and that order of magnitude estimates based on a polar analysis are valid.

We start from local oscillation analysis, assuming the perturbation variables $\delta Q$ to have dependence

[TABLE]

We neglect viscosity, but include thermal and Ohmic diffusion. We apply the same energy equation as above (though we now neglect baroclinicity), and make an incompressible approximation in the sense that sound speed $c_{\mathrm{s}}\rightarrow\infty$. The eigenvalue of Laplacian can be safely approximated as $k_{r}^{2}$ by WKB approximation, hence the energy equation reduces to

[TABLE]

where $\kappa$ is the thermal diffusivity, and $N_{\mu}$ is the composition part of Brunt-Väisälä frequency. Note that thermal diffusion suppresses the thermal component of the Brunt-Väisälä frequency, such that the compositional component is often more important in stars. From the linear MHD equations, we have the perturbed induction equation

[TABLE]

where $\vec{B}$ and $\vec{b}$ are the background and perturbed magnetic fields. We also have the perturbed equation of motion

[TABLE]

where $\vec{L}$ is the perturbed Lorentz force and $\hat{z}$ is the unit vector in $z$ direction. We show in Appendix C that

[TABLE]

where

[TABLE]

is the Alfvén frequency.

We show in appendix C.1 (with all the notations below made clear) that with some assumptions on the background field, these equations lead to the dispersion relation

[TABLE]

This is the local Tayler instability dispersion relation found by Zahn et al. (2007)333We believe there is a typo the dispersion relation of Zahn et al. (2007), originating from a typo in equation (A3) of Spruit (1999). The Alfvén velocity should be defined as $\omega_{\mathrm{A}}=V_{\mathrm{A}}/\varpi$, as illustrated from equation (2.12) in Acheson & Gibbons (1978). Spruit includes a redundant $m$ in his definition, leading to Zahn’s dispersion relation, different from ours by factors of $m$ in a few terms.

, but at an arbitrary latitude $\theta$. It differs only by factors of $\cos^{2}\theta$, which of course equal unity at the pole, and are of order unity at mid latitudes.

In Appendix C.2, we analyze the growing modes of this dispersion relation in the limits of rapid rotation ($\omega_{\rm A}\ll\Omega$), a magnetic field profile $B_{\phi}\propto\sin\theta$, and fast thermal diffusion, each of which are likely to be realized in many stars. The key conclusions are

Only $m=1$ modes can grow. Modes with $m=0$ or $m\geq 2$ are always damped. 2. 2.

Growing modes exist where $\theta<\pi/3$ or $\theta>2\pi/3$. They do not exist at the equator. 3. 3.

Modes only grow in the presence of finite (but not too large) magnetic diffusivity, and the fastest growing modes have growth rates

[TABLE]

These results confirm the findings of Spruit (1999, 2002); Zahn et al. (2007) at the pole and generalize them to arbitrary latitudes.

3.1 Non-linear Saturation

There is general agreement that equation 32 is a linear dispersion relation for Tayler modes, with the fastest growing modes having $m=1$ and growth rate comparable to equation 33. However, the non-linear evolution, and likely turbulent saturation of the instability, is not well understood or agreed upon. Similarly, the AM transport via Maxwell stresses in the saturated state remains unclear. Spruit (2002) proposed that the instability saturates when a turbulent damping rate is equal to the mode growth rate, and when the magnetic energy generation by winding radial field lines is equal to the magnetic energy dissipation by turbulence. Zahn et al. (2007) contested Spruit’s argument on the grounds that the dynamo loop proposed to regenerate the radial component of the magnetic field cannot operate.

Fuller et al. (2019) re-examined the non-linear evolution of the Tayler instability. They agreed with Zahn et al. (2007) about the dynamo loop closure problem, and used stability considerations to propose an alternative criterion for the strength of the radial component of the magnetic field. They also investigated the rate at which magnetic perturbations non-linearly cascade to smaller scales, equating this with the energy damping rate. Fuller et al. (2019) argued that Spruit (2002) greatly overestimated the damping rate by assuming that energy in the ordered background field $\vec{B}$ could be damped, when only energy in the disordered (i.e., varying on short length scales) perturbed field $\vec{b}$ cascades to smaller scales to be damped. Whereas Spruit (2002) assumed an energy damping rate $\dot{E}\sim\gamma B^{2}$, with $\gamma$ given by equation 33, Fuller et al. (2019) estimated an energy damping rate $\dot{E}\sim\gamma|b|^{2}$. The smaller energy damping rate of Fuller et al. (2019) allows the radial and horizontal magnetic field to grow to larger amplitudes, exherting stronger Maxwell stresses and causing more AM transport than prior predictions.

3.2 Tayler Torques

Fuller et al. (2019) estimated the torque density $T$ from Maxwell stresses in the saturated state of the Tayler instability, finding

[TABLE]

where $\alpha$ is a constant of order unity, defined via the saturated Alfvén frequency $\omega_{\rm A}=\alpha\Omega(q\Omega/N_{\rm eff})^{1/3}$. We use a fiducial value of $\alpha=1$, which approximately matches asteroseismic measurements of core rotation rates in low-mass stars. The effective stratification is typically $N_{\rm eff}\simeq N_{\mu}$ except where no signifcant composition gradient exists, as discussed in Fuller et al. (2019). We implement Tayler torques in our models via an effective viscosity operating on gradients of rotation frequency,

[TABLE]

As in Fuller et al. (2019), these toques are only implemented above a critical shear rate

[TABLE]

because magnetic diffusion eliminates the instability at smaller shears. We assume rotation constant on spheres in radiative regions due to the much more rapid AM transport in the horizontal direction relative to the radial direction. Convective regions of our models are effectively rigidly rotating due to the assumption of a very large effective turbulent viscosity in convective zones.

4 Results

4.1 Models

We predict internal rotation rates of massive stars using stellar models constructed with the MESA stellar evolution code (Paxton et al., 2011; Paxton et al., 2013, 2015, 2018). The models include the effective AM diffusivity of equation 35 when the shear exceeds equation 36, while the classic prescription for Tayler torques (Spruit, 2002) in MESA has been turned off, nearly identical to the low-mass models of Fuller et al. (2019). We study massive stars with initial masses ranging from $12-45\,M_{\odot}$, and roughly solar initial metallicity $Z=0.02$. Our models include mass loss via the “Dutch" prescription (with efficiency $\eta=0.5$), and moderate convective overshoot. To determine the dependence of remnant rotation rate on progenitor rotation rate, we examine models with initial equatorial rotational velocities of $50,150$ and $450\,\mathrm{km/s}$. Whereas typical massive stars have rotation velocities in the range $50-150$ km/s, 450 km/s represents an extremely rapidly rotating massive star (de Mink et al., 2013). We run our models from the zero age main sequence (ZAMS) to the moment when the silicon core of a star exceeds $1.5\,M_{\odot}$, i.e., just before the onset of core silicon burning. At this point, the star is only $\sim$days from core-collapse, and extrapolation of our models predicts no more significant AM transport (see the discussion below).

To gauge the efficiency of AM transport at different phases of evolution, we compute the specific AM contained within a mass coordinate $m$,

[TABLE]

with

[TABLE]

In the absence of AM transport, $j$ is constant at every mass coordinate within the star. However, as the core contracts during stellar evolution, its angular velocity increases relative to the envelope, creating shear within the star. Tayler torques act upon this shear, causing the core to spin slower by transporting AM from the core to the envelope. Phases of evolution where the core $j$ decreases the most are most important for determining its pre-collapse rotation rate. A useful diagnostic is the core’s specific AM at a mass coordinate $m_{1.5}=1.5\,M_{\odot}$, denoted $j(m_{1.5})=j_{1.5}$, which is a typical baryonic NS mass. In the absence of other effects, the core AM within this mass coordinate determines the natal spin rates of NSs via

[TABLE]

and we adopt a typical NS moment of inertia $I_{\rm NS}=10^{45}\,{\rm g}\,{\rm cm}^{2}$.

Figure 2 shows the specific AM $j_{1.5}$ of a $14\,M_{\odot}$ stellar model as it evolves, and the corresponding NS rotation rate $P_{\rm NS}$. The different colors correspond to different burning processes during stellar evolution, in the sense that when the central mass fraction of an element has decreased below $10^{-2}$, we define that burning phase to have ceased. Hence, the blue line denotes core hydrogen burning, the green line shows hydrogen shell burning and core helium burning, etc. Note the most significant AM loss occurs immediately after core hydrogen burning, where roughly 99% of the core’s AM is lost. The sudden AM loss occurs as the helium core contracts and the hydrogen envelope expands, generating shear that is damped by Tayler torques. Fuller et al. (2019) found very similar behavior in low-mass red giant stars, which lose most of their core AM just after hydrogen depletion, as the helium core contracts and the star moves up the red giant branch.

Figure 3 demonstrates why so much AM is lost during hydrogen shell-burning, which is primarily related to the ratio of $\Omega/N_{\rm eff}$ that determines the AM diffusivity of equation 35. On the main sequence, the value of $\Omega/N_{\rm eff}$ is large enough that AM transport is efficient and nearly rigid rotation is maintained. The core loses a small amount of AM due to the moderate contraction of the core and expansion of the envelope during core hydrogen burning. After the main sequence, the core contracts (and the envelope expands) by a large factor, generating internal shear and allowing Tayler torques to extract large amounts of AM from the core. When the helium core first starts to contract, the value $\Omega/N_{\rm eff}\sim 10^{-3}$ is not too small, allowing most of the core’s AM to be extracted. During core helium burning, the value $\Omega/N_{\rm eff}\sim 10^{-4}$ decreases, slowing AM loss from the core and forming the plateau at $\sim 10^{5}-10^{6}$ years in Figure 2. Large amounts of differential rotation are present within the star from this stage onward. After core helium burning, the carbon core contracts, again increasing the shear and allowing for more efficient core AM extraction. This produces the steeper fall off of $j_{1.5}$ at $\sim 10^{4}$ years in Figure 2. Once core carbon burning begins, the core’s AM content has nearly reached its final value, and only drops by another factor of $\sim$2 before core-collapse.

In Figure 4, we plot the value of $j_{1.5}$ for stars of several initial masses. In each model, the vast majority of the core’s AM is lost just after core hydrogen burning during helium core contraction. We see from Figure 4 that the AM extraction rate during helium core contraction is similar (or even larger) in massive stars. However, in more massive stars, the helium core does not have to contract as much before helium burning begins, such that less AM is extracted by the time the rotation profile reaches the critical shear level (equation 36), after which very little AM is transported. Correspondingly, the more extended helium-burning cores of massive stars contract by a larger factor between the end of helium burning and the beginning of oxygen burning. In these higher mass stars, more AM is extracted during carbon/oxygen core contraction, compared to that extracted during helium core contraction. The two effects somewhat balance, such that the values of $j_{1.5}$ in massive stars at the start of oxygen burning are only a few times larger than in lower mass stars. While the AM extraction rate during core carbon/oxygen burning is larger than it is during helium burning (i.e., the slope $d\log j_{1.5}/dt$ in Figure 4 is large during carbon/oxygen burning), the duration of these phases is very short. Hence, the fraction of core AM lost during carbon/oxygen burning is much smaller than that lost during hydrogen shell-burning.

The value of $j_{1.5}$ at the end of each track in Figures 2 and 4 corresponds to the expected natal NS rotation period $P_{\rm NS}$ shown by the right axis. Figure 5 shows these values of $P_{\rm NS}$ for each of our models, with different initial surface velocities marked by different colors. These models predict decreasing NS rotation periods with increasing stellar mass, but all of our models have $P_{\rm NS}$ in the range $50-200\,\mathrm{ms}$. As discussed above, the faster NS spins of high-mass models stems from the fact that their cores have a larger physical extent that allows them to retain more AM.

Remarkably, for most models, the final NS rotation periods have little dependence on the initial surface velocity or total AM of the star. This result is similar to that of Fuller et al. (2019) for low-mass stars, who found post-main sequence internal rotation rates scaled very weakly with initial rotation rate, due to the convergent evolution in core rotation rate caused by the strong scaling of equation 35 with rotation rate, and the fact that AM tends to be transported out of the core until the minimum shear of equation 36 is realized. The dispersion of final rotation periods for stars of the same mass in Figure 5 appears to stem from slightly different evolutionary histories (due to small numerical differences coupled with “chaotic” stellar evolution as described by Sukhbold et al. 2018), rather than being directly tied to the initial rotation rate. We find that models with higher resolution produce nearly identical results (suggesting our models are largely converged) apart from some small residual scatter similar to that shown in Figure 5.

Finally, we must check whether Tayler torques can operate as prescribed by equation 35 during all phases of evolution of our models. The instability can no longer reach its equilibrium and transport AM according to equation 35 if the Tayler mode growth time $t_{T}\sim\Omega/\omega_{\rm A}^{2}$ is longer than the AM transport timescale $t_{\rm AM}\sim r^{2}/\nu_{\rm AM}\sim N_{\rm eff}^{2}/\Omega^{3}$. Using the saturated Alfvén frequency from Fuller et al. (2019), $\omega_{\rm A}\sim\Omega(q\Omega/N_{\rm eff})^{1/3}$, estimating $N_{\rm eff}\sim\sqrt{G\rho}$, $j=\Omega r^{2}$, $\rho\sim mr^{-3}$ and using $q\sim 1$, we find $t_{\rm T}/t_{\rm AM}>1$ when $r_{\rm min}\lesssim j^{2}/(Gm)$. Evaluating this timescale at a mass coordinate $m=1.4\,M_{\odot}$ at the end of our runs yields $r_{\rm min}\sim 1\,{\rm cm}$, which is much smaller than the corresponding radial coordinate. Hence our models are always in an acceptable regime. We note that our models can reach a regime where $t_{\rm T}$ is longer than the evolution time, which can occur around the silicon burning phase. However, the AM transport time $t_{\rm AM}$ at this phase is much longer than the remaining lifetime of the star, so the effects of Tayler torques are negligible.

The pre-collapse AM distribution may affect the dynamics of the SN explosion process (Foglizzo, 2017). Figure 6 shows profiles of the angular frequency and specific AM within the 14 $M_{\odot}$ model shown in Figure 2. We have plotted these quantities at two times, when the silicon core and the iron core reach $\approx 1.4\,M_{\odot}$. We find a pre-collapse rotation frequency of $\Omega\sim 2\times 10^{-3}$ rad/s in the iron core, corresponding to a specific AM of $j\sim 10^{13}\,{\rm cm}^{2}\,{\rm s}^{-1}$. Sharp gradients in the rotation frequency are associated with composition gradients between different burning shells. Neither the rotation profile nor the AM profile are flat in the core. At the end of the star’s evolution, Tayler torques become inefficient at redistributing AM within radiative zones, but convective torques still enforce profiles with constant $\Omega$ and hence $j\propto r^{2}$. After convection subsides, shells change their radial coordinate and their rotation rate, but their specific AM is approximately conserved. The final rotation profile thus depends on the details of the star’s evolution, but the AM profile as a function of mass remains roughly constant.

5 Discussion

The AM transport predicted by our models has important implications for massive star evolution. During the main sequence, we find the core and envelope remain tightly coupled, with little internal shear. Hence, prior stellar evolution calculations exhibiting mixing due to shear instabilities (e.g., Meynet & Maeder 1997; Georgy et al. 2013) have likely overestimated the extent of mixing, though rotational mixing via Eddington-Sweet circulation can still occur. This mixing is typically most important on the main sequence, where both our model and the Tayler-Spruit dynamo (Spruit, 2002) predict nearly rigid rotation. Main sequence rotational mixing for such models are unlikely to be greatly altered, but post-main sequence mixing should be re-examined in light of the slower internal rotation rates we predict.

Our calculations have crucial implications for the supernova explosions that mark their deaths, and for the spins of compact object remnants left behind. For our simple models of single stars at solar metallicity, the prediction is clear: AM transport is efficient, massive stellar cores rotate slowly, and compact objects are likely to be born slowly rotating.

For most NSs, our models predict initial spin periods $P_{\rm NS}\sim 100-200\,{\rm ms}$ as shown in Figure 5, nearly independent of initial rotation velocity of the star. We caution that these predictions have some dependence on the Tayler instability saturation parameter $\alpha$, which Fuller et al. (2019) found to be $\alpha\approx 1$ based on calibration by asteroseismic measurements of low-mass stars. There is also dependence on the uncertain moment of inertia $I$ for typical NSs, but neither effect produces uncertaintly larger than a factor of $\sim\!2$. However, red giant cores and white dwarfs exhibit scatter by a factor of a few around the predicted rotation periods. It would not be surprising if a similar scatter exists for massive star cores, potentially producing NSs with natal rotation periods of tens to hundreds of milliseconds.

Interestingly, our predictions are near the initial spin rates estimated for typical pulsars (Faucher-Giguère & Kaspi, 2006; Popov et al., 2010; Popov & Turolla, 2012; Gullón et al., 2014), which likely have a broad distribution of tens to hundreds of milliseconds. Our predictions are also similar to the spins of central compact object (CCO) pulsars, which have $P_{\rm NS}\sim 100-400\,{\rm ms}$ (Halpern & Gotthelf, 2010; Gotthelf et al., 2013). Since the spin-down times of these “anti-magnetar” COOs are orders of magnitude larger than their ages, the observed spin period is nearly identical to their natal spin period, justifying a direct comparison with our predictions.

However, it is clear that some pulsars are born spinning faster than our estimates, e.g., the Crab pulsar which is estimated to have had a natal spin rate of $P_{\rm NS}\sim 20\,{\rm ms}$ (Kaspi & Helfand, 2002). We posit that most NSs have had their spin rates altered by stochastic AM transport processes during the final stages of massive star evolution or during the SN explosion (Spruit & Phinney, 1998). Fuller et al. (2015) considered the AM transported by internal gravity waves excited by the vigorous convection due to late nuclear burning stages. The stochastic influx of AM into the core of the star results in a Maxwellian distribution of NS spin frequencies with typical natal spin periods of $P_{\rm NS}\sim 50-100\,{\rm ms}$. Additionally, three-dimensional simulations of core-collapse SNe demonstrate asymmetric explosions that stochastically spin the NS to periods as short as $P_{\rm NS}\sim 20\,{\rm ms}$ (Müller et al., 2018a), though prompt explosions of low-mass or ultra-stripped progenitors produce much less spin-up (Müller et al., 2018b; Rantsiou et al., 2011). Given the slow core rotation predicted by our models, these stochastic mechanisms could be the dominant processes that determine natal NS rotation periods. In these scenarios, the rotation rate and spin-axis of a young NS is essentially uncorrelated with the spin of its progenitor star. Indeed, the spin-orbit misalignment of the more recently formed NS in the binary pulsar system PSR J0737-3039 is consistent with this hypothesis (Farr et al., 2011). The relatively slow spins of CCOs may indicate they originated from prompt explosions, with spins inherited from their progenitor star.

Very rapidly rotating NSs and BHs are almost certainly formed under rare circumstances, and such objects are thought to be the central engines of broad-lined Ic SNe, gamma-ray bursts, and superluminous SNe (e.g., Woosley 1993; Kasen & Bildsten 2010; Metzger et al. 2011). We believe that close binary evolution or homogeneous evolution (e.g., Yoon et al. 2006) may be required to form these rapidly rotating compact objects. Tidal spin-up of a helium star progenitor (e.g., Kushnir et al. 2016; Qin et al. 2018b) will allow for rapidly rotating BHs. Whether this scneario can allow for rapidly rotating NSs is unclear, as Tayler torques may be able to couple the inner core of the star with the overlying helium envelope, preventing rapid core rotation. Preliminary binary runs indicate this to be the case, but we defer a more detailed treatment of binary scenarios to future work.

The leading alternative to our models is that of Kissin & Thompson (2015, 2018). In stark contrast to our models, those models feature rigid rotation in radiative regions (enforced by fossil fields) and differential rotation in convective zones (resulting from inward convective AM pumping). Interestingly, these models also predict NS rotation periods of $\sim$hundreds of milliseconds for typical progenitors. But they allow for much more rapid NS rotation for certain stars (e.g., $25\,M_{\odot}$) with very thick surface/core convective zones, and very slow NS rotation for more massive stars that have lost most of their AM to winds. Distinguishing between our model and theirs will benefit from more detailed comparisons with data for low-mass stars.

6 Conclusion

We have examined angular momentum (AM) transport in massive stars via the baroclinic and Tayler instabilities. In contrast to many previous works (Knobloch & Spruit, 1982, 1983; Spruit & Knobloch, 1984; Zahn, 1992), and in agreement with other works (Tassoul & Tassoul, 1982; Fujimoto, 1987; Kitchatinov, 2014), we find that unstable baroclinic modes generally exist in differentially rotating radiative zones, when thermal diffusion is neglected. In most stars, however, these modes will be stabilized by thermal diffusion, such that baroclinic instabilities are unlikely to play an important role in AM redistribution.

We also analyze the Tayler instability at arbitrary latitudes inside of stars, finding that unstable modes likely do not exist near the stellar equator and growth rates peak near the pole of the star. Typical growth rates and wavenumbers are similar to those found in prior work (e.g., Spruit 1999, 2002; Zahn et al. 2007). We then apply the Tayler instability saturation criterion and AM transport prescription of Fuller et al. (2019) to massive stellar models using MESA. We find that after the main sequence, roughly 99% of core AM is lost during hydrogen shell-burning while the helium core is contracting and the star is expanding into a red supergiant. The combination of large shear and modest buoyancy frequency during this stage allows for efficient AM transport, preventing the spin-up of the contracting core as its AM is transferred to the expanding envelope. A smaller, but still significant, amount of AM is lost during helium burning and beyond. The small core spin rates and internal shears of our models imply that rotational mixing effects may be somewhat overestimated in stellar models not taking magnetic torques into account.

Assuming AM conservation during core-collapse and the subsequent supernova explosion, our models predict natal rotation periods of $P\sim 100-200$ ms for NSs born from progenitors of $12-20\,M_{\odot}$. We predict faster rotation for higher-mass stars, but even our $45\,M_{\odot}$ model would produce a relatively slowly rotating NS with $P\approx 50$ ms. Hence, our models predict rotation periods similar to those inferred for newborn pulsars (Faucher-Giguère & Kaspi, 2006; Popov et al., 2010; Popov & Turolla, 2012; Gullón et al., 2014), though perhaps near the slow end of the distribution. In light of the slow core rotation of our models, the spins of many NSs may be primarily determined by stochastic AM redistribution processes just before core-collapse (Fuller et al., 2015) or during the explosion process (Spruit & Phinney, 1998). Relatively rapidly rotating pulsars like the Crab pulsar ($P_{i}\sim 20\,{\rm ms}$, Kaspi & Helfand 2002) can be accounted for by core-collapse spin-up as found in some supernova simulations (Müller et al., 2017; Müller et al., 2018a).

Because AM is generally transported from the contracting helium core into the expanding hydrogen envelope, our models predict that many black holes will be born slowly rotating, which we investigate in a companion paper (Fuller & Ma, in prep). However, it is important to emphasize that the majority of massive stars evolve in interacting binaries (Sana et al., 2012), and hence the predictions of our single-star models are merely one piece of the puzzle. We suspect AM accretion via mass transfer and tidal spin-up in close binaries will be able to produce more rapidly rotating compact objects capable of powering various energetic transients, a subject we plan to investigate in future work.

7 Acknowledgements

We thank Yulin Gong for helpful comments on the manuscript, and Adam Jermyn and Anthony Piro for useful insights on the baroclinic instability. This research is funded in part by an Innovator Grant from The Rose Hills Foundation, the Sloan Foundation through grant FG-2018-10515, and by the Gordon and Betty Moore Foundation through Grant GBMF7392.

Appendix A Energy Equation with Thermal Diffusion and Baroclinic Term

Here we derive the energy equation for an ideal gas with thermal diffusion and baroclinicity. We work in a corotating frame such that the background velocity $v$ is zero, and we assume an axisymmetric background structure. The entropy for an ideal gas is

[TABLE]

and the Brunt-Väisälä frequency is

[TABLE]

where $\Gamma_{1}=c_{P}/c_{V}$. We assume that temperature diffusion is the only process that changes entropy, which is a good approximation for short wavelength disturbances deep inside a star, such that

[TABLE]

Taking the Eulerian perturbation and assuming time dependence of each perturbed variable $\delta Q$,

[TABLE]

we have

[TABLE]

To obtain equation 46 we have used a standard WKB approximation by setting $\nabla^{2}\delta T=-k^{2}\delta T$, where $k$ is the wavenumber of the perturbation, and dropping terms on the right hand side with weaker $k$ dependence.

Next, following Kitchatinov (2014), we neglect higher order centrifugal terms (valid in relatively slowly rotating stars) by defining our radial coordinate $r$ to be perpendicular to isobars $r$ so that

[TABLE]

Next, we assume rotation constant on spherical shells, due to efficient horizontal AM transport within radiative regions of a star. Defining the latitudinal Brunt-Väisälä frequency

[TABLE]

and using the baroclinicity relation for shellular rotation (e.g., Kitchatinov 2014), we have

[TABLE]

where

[TABLE]

is the dimensionless shear. The latitudinal Brunt-Väisälä frequency $N_{\theta}^{2}$ is a horizontal buoyancy force due to baroclinicity arrising from differential rotation.

Combining equations 49, 50, and 41, we have

[TABLE]

For an ideal gas

[TABLE]

The Lagrangian perturbation $\Delta\mu$ is related to the Eulerian perturbation $\delta\mu$ by

[TABLE]

where we have assumed $\mu$ has only $r$ dependence, again due to efficient latitudinal transport of composition and AM. Then

[TABLE]

together with (53), we have

[TABLE]

Assuming $c_{P}$ and $c_{V}$ to be constant (as they are for an ideal gas), we have

[TABLE]

Defining the compositional buoyancy frequency

[TABLE]

we have

[TABLE]

[TABLE]

Note that

[TABLE]

where $N_{T}$ is the thermal component of the buoyancy frequency. Then

[TABLE]

can be rewritten

[TABLE]

where

[TABLE]

Equation 66 is the energy equation with thermal diffusion. When thermal diffusion is not important, i.e. $\kappa\rightarrow 0$,

[TABLE]

where $c_{\text{s}}$ is the adiabatic sound speed. The energy equation reduces to

[TABLE]

which is the adiabatic energy equation with a baroclinic term. When thermal diffusion is strong, i.e. $\kappa\rightarrow\infty$

[TABLE]

where $c_{T}$ is the isothermal sound speed. In this limit, the energy equation reduces to

[TABLE]

The thermal and baroclinic terms are suppressed by thermal diffusion, and only the composition gradient is important. We note that the baroclinic term has vanished due to our assumption of composition being constant on isobars. If a latitudinal composition gradient can persist, equation 71 would contain a compositional baroclinic term.

Appendix B Growth Rate of the Baroclinic Instability

Because baroclinicity and Coriolis forces break the spherical symmetry of a star, calculating the frequencies and growth rates of global scale modes becomes more difficult. Here we derive the growth rate of the baroclinic instability using the traditional approximation and the geostrophic approximation, following the assumptions made in the main text.

B.1 Traditional Approximation

We start by defining perturbation variables to have time and spatial dependence

[TABLE]

The equation of motion (equation (6)) can be expressed in spherical coordinates

[TABLE]

[TABLE]

[TABLE]

Since the Coriolis force in the radial direction is expected to be much smaller than the buoyancy force, equation (73) reduces to

[TABLE]

We further expect the radial perturbation to be small, such that $\xi_{r}\ll\xi_{\theta}$, hence equation (75) reduces to

[TABLE]

These approximations are known as the traditional approximation (see e.g., Bildsten et al., 1996; Lee & Saio, 1997). It is straightforward to solve for $\xi_{\theta}$ and $\xi_{\phi}$ in terms of $\delta P$ from the above equations, yielding

[TABLE]

where we have defined $\nu=2\Omega/\omega$ and $\mu=\cos\theta$. The angular parts of the equation of motion then yield the relation

[TABLE]

where

[TABLE]

Now consider the continuity equation

[TABLE]

We use the incompressible approximation to neglect the first two terms in equation 82, but we keep the $\xi_{\theta}$ term in order to preserve its baroclinic effect. Then we have

[TABLE]

Inserting equation 80 and making a WKB approximation, we find

[TABLE]

Recalling the definition of $N_{\theta}^{2}$ in equation 50, and substituting the expression for $\xi_{\theta}$, we have

[TABLE]

where

[TABLE]

From the definition of $\hat{\mathcal{L}_{\mu}}$ and $\hat{\mathcal{L}_{y}}$, we have $\hat{\mathcal{L}_{\mu}}\sim\hat{\mathcal{L}_{y}}\sim 1$ for the modes with smallest angular wavenumber. For all but the fastest rotating stars, $N_{\theta}^{2}r/g\sim\Omega^{2}r/g\ll 1$, thus we can neglect the second term in the left-hand side of equation 85 to find

[TABLE]

Now, combining the energy equation (66) with the radial equation of motion (76), we have

[TABLE]

We then make the WKB approximation $\partial\delta P/\partial r\approx ik_{r}\delta P\approx ik\delta P\gg g\delta P/\tilde{c_{\mathrm{s}}}^{2}$. Expressing $\xi_{\theta}$ in terms of $\delta P$, and limiting ourselves to low-frequency oscillations with $\tilde{N}^{2}\gg\omega^{2}$, we find

[TABLE]

Now combining with equation 87, we cancel the variable $\xi_{r}$ to obtain

[TABLE]

For very low-frequency modes ($\omega\ll\Omega\ll N$) obeying the traditional approximation, $k_{r}\approx(\tilde{N}/\omega)k_{\perp}$, and $\hat{\mathcal{L}_{y}}\sim rk_{\perp}\omega/\Omega$, so the baroclinic term in equation 90 is of order $(q\Omega/\tilde{N})r^{2}k_{\perp}^{2}$, while the first term is of order $r^{2}k_{\perp}^{2}$. For higher frequency modes with $\omega>\Omega$, the baroclinic term is smaller by a factor $\Omega/\omega$. Hence we can treat the baroclinicity as a small perturbation. The unperturbed equation is

[TABLE]

The eigenfunctions of the operator $\hat{\mathcal{L}_{\mu}}$ are Hough functions $H_{\lambda}(\mu)$ satisfying

[TABLE]

To first order in the baroclinic perturbation, we have the dispersion relation

[TABLE]

where

[TABLE]

is a dimensionless number depending on $\nu$. We thus obtain the dispersion relation

[TABLE]

which is identical to the usual dispersion relation for Hough modes, but modified to include thermal diffusion and with a baroclinic driving term proportional to $N_{\theta}^{2}$.

To solve for the growth rate of baroclinic modes, we first consider the case with small thermal diffusivity, so to zeroth order (when thermal diffusion is neglected), $\tilde{N}^{2}\sim N^{2}$, $\tilde{N}_{\theta}^{2}\sim N_{\theta}^{2}$. Again using the fact that the baroclinic term is much smaller than the buoyancy term, we find

[TABLE]

Expanding equation 95 to first order in thermal diffusivity and baroclinicity,

[TABLE]

Substituting the real part of equation 96 for $\omega$

[TABLE]

hence

[TABLE]

So the growth rate is

[TABLE]

The first term is always a damping term caused by thermal diffusion, while the second term accounts for baroclinicity and can cause either growth or damping. For the typical case in stars where $q<0$ such that $N_{\theta}^{2}<0$, only baroclinic modes with negative $Fk_{r}$ can grow.

In the limit of rapid thermal diffusion, i.e. $\kappa\rightarrow\infty$, we have $\tilde{N}^{2}\sim N_{T}^{2}\omega/i\kappa k_{r}^{2}+N_{\mu}^{2}$ and $\tilde{N_{\theta}}^{2}\sim N_{\theta}^{2}\omega/i\kappa k_{r}^{2}$. Then equation 95 reduces to

[TABLE]

This is a quadratic equation of the form $\omega^{2}-B\omega-C^{2}=0$, and in the limit of rapid thermal diffusion $B^{2}\ll C$. Then we have

[TABLE]

So the growth rate is

[TABLE]

Baroclinic modes can still grow if $N_{\theta}^{2}rk_{r}\gtrsim N_{T}^{2}$, although their growth rate is strongly suppressed by the effect of thermal diffusion. Growth requires large wavenumbers $rk_{r}\gtrsim N_{T}^{2}/\Omega^{2}$, which is typically an enormous number in stellar interiors, so the growth rate will be very small. In fact, the required wavenumbers can be so large that molecular viscosity can become important. Hence, in many situations where thermal diffusion is important, the baroclinic instability will be quenched or grow too slowly to be relevant.

B.2 Geostrophic Approximation

We now examine the baroclinic instability using a geostrophic approximation. For this calculation, we analyze the oscillation modes locally, in a Cartesian box at colatitude $\theta$. We adopt a coordinate system with $\hat{x}$ in the $\phi$ direction, $\hat{y}$ in the $\theta$ direction, and $\hat{z}$ in the $r$ direction. The perturbations variables are assumed to have dependence

[TABLE]

We have the equations of motion

[TABLE]

[TABLE]

[TABLE]

where $f=2\Omega\cos\theta$. Similar to the traditional approximation, the geostrophic approximation uses $v_{z}\ll v_{x},v_{y}$ and $\Omega\ll N$ to drop Coriolis terms in equations 105 and 107. The continuity equation again uses an incompressible approximation to become

[TABLE]

We take the time derivative of the energy equation (66), applying incompressible approximation, i.e. $c_{\mathrm{s}}\rightarrow\infty$ to find

[TABLE]

Equations 105-109 can be combined to yield a geostrophic dispersion relation. Combining equations 107 and 108,

[TABLE]

Applying the WKB approximation, i.e. $k_{x},k_{y},k_{z}\gg r,H$, we take the $z$ derivative of equation 110:

[TABLE]

hence

[TABLE]

Note that in our thin box centered at $(r_{0},\theta_{0}$), $y=r_{0}(\theta-\theta_{0})$, so that

[TABLE]

and

[TABLE]

Therefore

[TABLE]

Substituting the expression for $f$, and assuming $N_{\theta}^{2}r_{0}\ll g$, we have

[TABLE]

This is the local dispersion relation for geostrophic modes in the presence of large buoyancy and baroclinicity. Note that only retrograde Rossby modes exist, as normal gravity modes have been filtered out by the geostrophic approximation. Despite prior assertions to the contrary (Spruit et al., 1983; Zahn, 1992), the baroclinic term in equation 116 can cause local growth of baroclinic modes in the geostrophic approximation. Just as in the traditional approximation, however, only modes propagating in one direction (one sign of $k_{z}$) will be able to grow. Standing modes composed of waves traveling in both directions will therefore not be able to grow, as discussed in the main text.

Considering now the case with small thermal diffusion and baroclinicity, to zeroth order the wave frequency is

[TABLE]

Expanding equation 116 in terms of thermal diffusivity,

[TABLE]

The growth rate is thus

[TABLE]

Just as our result for the traditional approximation, baroclinic modes can only grow for sufficiently small wavenumbers such that the baroclinic term outweighs the thermal damping term. We also note that Rossby modes in the traditional approximation have $F\sim\lambda$. For these modes, the traditional approximation growth rate equation 100 can be rearranged to be identical to equation 119, modulo factors of order unity.

When thermal diffusion is large and there is a finite composition gradient, we find

[TABLE]

where $a=-k_{x}\sin\theta_{0}/(2r_{0}k_{z}^{2}\Omega\cos^{2}\theta_{0})$. Equation 120 can be rearranged to show that growth requires $q\Omega\gtrsim k^{2}\kappa(k_{z}/k_{x})N_{T}^{2}/N_{\mu}^{2}$. In real stars, growth often requires wavenumbers small enough that the large thermal diffusion approximation is invalidated, i.e., there are often no growing modes when thermal diffusion is large.

Appendix C Local Analysis of Tayler Instability

Here we derive the dispersion relation and growth rate of unstable modes of Tayler instability.

C.1 Dispersion relation

Here we give a detailed discussion about the dispersion relation of Tayler instability in spherical coordinates. This generalizes the results by Tayler (1973); Spruit (1999); Zahn et al. (2007) in which only the dispersion relation near the pole is considered. Similar to preceding sections, we assume all perturbation variables $\delta Q$ have short radial wavelength due to strong buoyancy, with time/spatial dependence

[TABLE]

i.e., we make a WKB approximation. The eigenvalue of Laplacian can be simplified to

[TABLE]

As above, we make an incompressible approximation, $c_{\mathrm{s}}\rightarrow\infty$. Here, however, we ignore the baroclinic term because we shall see that thermal diffusion is generally large at the short wavelengths characteristic of the Tayler instability. The energy equation 66 then reduces to

[TABLE]

Including magnetic forces in the MHD limit, the perturbed induction equation is

[TABLE]

In differentially rotating stars, the azimuthal magnetic field is generated by winding up a weak poloidal magnetic field. The instability requires the azimuthal field to become much larger than the poloidal seed field, so we assume a background (unperturbed) field of

[TABLE]

Here, we have assumed latitudinal dependence proportional to $\sin\theta$ for the background field, which arises from our assumption of rotation constant on spherical shells, i.e., $\partial v_{\rm rot}\partial r=q\Omega\sin\theta$. Substituting equation 126 into the induction equation, we get

[TABLE]

where we applied incompressible approximation, i.e. $\nabla\cdot\vec{\xi}=0$. Due to the strong buoyancy force in radial direction, $\xi_{r}$ is small thus the second term can be neglected, then

[TABLE]

Now we calculate the perturbed Lorentz force, to first order

[TABLE]

By (126) and (128), we have

[TABLE]

and so

[TABLE]

and

[TABLE]

where we applied the WKB approximation, hence the Lorentz force is

[TABLE]

where $\omega_{\mathrm{A}}=B_{0}/\sqrt{4\pi\rho}r$ is the Alfvén frequency.

We also have the perturbed equation of motion

[TABLE]

Substituting the expressions for Lorentz force and (124) into the equation of motion gives

[TABLE]

[TABLE]

[TABLE]

We now define

[TABLE]

The continuity equation and equations (137)-(139) reduce to

[TABLE]

The dispersion relation is given by the determinant of the linear equations, i.e.,

[TABLE]

For Tayler modes, we shall see that $A\sim B\sim C\ll D$, so we further define

[TABLE]

Considering that $D\gg A$

[TABLE]

From our WKB approximation in the radial direction, $k_{r}r\gg l$, thus we can neglect the last term, and the dispersion relation reduces to

[TABLE]

Next, we define the following dimensionless variables (the same definitions used in Zahn et al. 2007) as

[TABLE]

Substituting the expressions for $A,B$ and $\tilde{N}$, and after some algebra, we arrive at the dispersion relation

[TABLE]

which is exactly a generation of the dispersion relation near the pole given by Zahn et al. (2007), whose result is only valid near the pole while our result works at arbitary latitude.

C.2 Unstable modes

Equation 154 is a little lengthy thus difficult to analyze, so we start with the more compact equation 152. Using the same notation as Spruit (1999), after a little algebra, equation 152 reduces to

[TABLE]

This is our version of equation (A8) from Spruit (1999). Note that we have set $p=1$ here, which is appropriate for a background magnetic field with $B\propto\sin\theta$ as expected for a field generated by winding a radial seed field.

For $\Omega\gg\omega_{A}$, we expect $\omega\sim\omega^{2}_{A}/\Omega$ for both real and imaginary parts, so we define

[TABLE]

and

[TABLE]

When the thermal diffusion is large, which is often the case in stars, we can take the limit $K\rightarrow\infty$. Note that the limit $K\rightarrow\infty$ is the same as the limit $k\rightarrow 0$, $N=N_{\mu}$ considered by Spruit (1999). We further limit our analysis to slow modes, i.e., we take the limit $\omega\ll\Omega$, which eliminates $r$ modes and gravity modes. The dispersion relation then reduces to

[TABLE]

This is our version of equation (A18) from Spruit (1999) for the case $p=1$. The latitudinal dependence arises only through the factor of $\cos^{2}\theta$ in the last term.

A growing mode requires equation (158) to have a solution with a positive imaginary part, which straightforwardly excludes $m=0$ modes from our interest. Though the cubic complex equation (158) is in general difficult to solve, it is easy to see that in the case with large magnetic diffusivity, i.e. $Hk_{r}^{2}\gg 1$, the solutions of $\alpha$ must have a large negative imaginary part of $\sim-Hk_{r}^{2}$ to make the LHS vanish. Hence, large magnetic diffusion or very large wavenumber $k_{r}^{2}$ will inevitably lead to damped modes.

Let us fix the azimuthal wave number $m$ and the polar angle $\theta$ for now. We start with some $n^{2},Hk_{r}^{2}$ so that the solution of (158) is a damped mode, which is always possible by setting $Hk_{r}^{2}$ large enough. By continuously changing the parameters $n^{2},Hk_{r}^{2}$, the coefficients and solutions of the cubic equation vary continuously. The boundary between damped and growing modes is defined by some set of $n^{2},Hk_{r}^{2}$ that makes $\alpha$ purely real. Hence, requiring a real solution of (158) defines an instability criterion, with mode growth only possible on one side of this boundary, and no growing modes if there is no solution for a purely real $\alpha$. By setting $\alpha$ real, we can separate (158) into its real and imaginary parts:

[TABLE]

Eliminating $\alpha$ from the equations above, we have

[TABLE]

Equation (C.2) must be satisfied for (158) to have a real solution such that growing modes exist. Since $\cos^{2}\theta H^{2}k_{r}^{4}>0$, we must have

[TABLE]

Defining $t^{2}=n^{2}/\cos^{2}\theta$, we have

[TABLE]

Since $t^{2}>0$, the maximum of the right-hand side occurs when the positive sign is chosen with $t\rightarrow 0$, where the maximum value $\rightarrow 2$. Thus, the instability can only grow if

[TABLE]

This is only possible for $m=1$ modes, with

[TABLE]

Hence the instability can only occur at high latitudes where $\theta<\pi/3$ or $\theta>2\pi/3$.

For $m=1$ and $\cos\theta=1$, instability requires

[TABLE]

which requires

[TABLE]

and it is easy to show this requires

[TABLE]

So we see that growth requires wavelengths not much larger than the Tayler instability length scale $\sim(\omega_{\rm A}/N_{\mu})r$. Looking back at equation C.2, we see that growing modes also require

[TABLE]

From above, we found that only modes $m=1$, with $n$ less than a few, and with $\cos\theta>1/2$ can grow. Then the left-hand side of the above equation is of order unity for growing modes, and growth additionally requires $Hs^{2}\lesssim 1$. Growing modes must therefore have wavelengths $\gtrsim\sqrt{\eta\Omega}/\omega_{\rm A}$. Combining this criterion with the maximum wavelength criterion from above yields the criterion of Spruit (2002) for Tayler instability that

[TABLE]

We do not attempt to solve the cubic equation for mode growth rates as was done by Zahn et al. (2007). However, we note that when growth occurs, the fastest growing modes are those with $Hs^{2}$ of order unity. The wavelength of the fastest growing modes is similar to (but slightly larger than) those on the verge of stability, with scales of $\sim\sqrt{\eta\Omega}/\omega_{\rm A}$. When $\omega_{\rm A}\sim\omega_{c}=(N_{\rm eff}^{2}\Omega\eta/r^{2})^{1/4}$, as often happens in our models, the fastest growing modes have $rk_{r}\sim N_{\rm eff}/\omega_{\rm A}$.