Slowing the Spins of Stellar Cores

TL;DR

This paper investigates the magnetic Tayler instability as a mechanism for angular momentum transport in stellar interiors, providing improved models that match observed core rotation rates across different stellar evolutionary stages.

Contribution

It introduces a new saturation prescription for the Tayler instability, leading to more effective angular momentum transport in stellar models, aligning theory with observations.

Findings

Models reproduce nearly rigid solar and main sequence star rotation.

Core rotation rates of red giants during different burning phases are matched.

White dwarf rotation rates are accurately predicted.

Abstract

The angular momentum (AM) evolution of stellar interiors, along with the resulting rotation rates of stellar remnants, remains poorly understood. Asteroseismic measurements of red giant stars reveal that their cores rotate much faster than their surfaces, but much slower than theoretically predicted, indicating an unidentified source of AM transport operates in their radiative cores. Motivated by this, we investigate the magnetic Tayler instability and argue that it saturates when turbulent dissipation of the perturbed magnetic field energy is equal to magnetic energy generation via winding. This leads to larger magnetic field amplitudes, more efficient AM transport, and smaller shears than predicted by the classic Tayler-Spruit dynamo. We provide prescriptions for the effective AM diffusivity and incorporate them into numerical stellar models, finding they largely reproduce (1) the nearly rigid rotation of the Sun and main sequence stars, (2) the core rotation rates of low-mass red giants during hydrogen shell and helium burning, and (3) the rotation rates of white dwarfs. We discuss implications for stellar rotational evolution, internal rotation profiles, rotational mixing, and the spins of compact objects.

Equations

$$\omega_{\rm A}>\omega_{c}\sim\Omega\bigg{(}\frac{N}{\Omega}\bigg{)}^{\!1/2}\bigg{(}\frac{\eta}{r^{2}\Omega}\bigg{)}^{\!1/4}\,,$$

$${\omega }_{\mathrm{grow}}\sim \frac{{\omega }_{\mathrm{A}}^2}{\Omega }\mathrm{for}{\omega }_{\mathrm{A}}\lesssim 2\Omega .$$

$$l_r\sim \frac{1}{k_r}\lesssim l_{\perp }\frac{{\omega }_{\mathrm{A}}}{N},$$

$$N_{\mathrm{eff}}^2\simeq \frac{\eta }{K}{N}_T^2+N_{\mu }^2,$$

$${\gamma }_{\mathrm{turb}}\sim k_r^2{\eta }_{\mathrm{eff}}\sim \frac{{\omega }_{\mathrm{A}}^2}{\Omega },$$

$$\frac{\partial }{\partial t}{B}_{\varphi }=q\Omega B_r$$

$${\gamma }_{\mathrm{amp}}\sim q\Omega \frac{{\omega }_{\mathrm{A}}}{N}.$$

$$T=B_{r}B_{\phi}\sim 4\pi\rho r^{2}\Omega^{2}q^{3}\bigg{(}\frac{\Omega}{N}\bigg{)}^{4}\,.$$

$$\delta \vec{B}=(\vec{k}\cdot {\vec{B}}_{\varphi })\vec{\xi },$$

$$\delta \vec{v}\sim \frac{{\omega }_{\mathrm{A}}}{\Omega }\delta {\vec{v}}_{\mathrm{A}},$$

$$\frac{\partial}{\partial t}\vec{v}+\big{(}\vec{v}\cdot{\mbox{\boldmath$\nabla$}}\big{)}\vec{v}=2(\vec{\Omega}\times\vec{v})-\frac{{\mbox{\boldmath$\nabla$}}P}{\rho}+\frac{({\mbox{\boldmath$\nabla$}}\times\vec{B})\times\vec{B}}{4\pi\rho}-\vec{g}\,.$$

$$\delta v_{\perp }\sim \frac{{\omega }_{\mathrm{A}}}{\Omega }\delta v_{\mathrm{A},\perp }.$$

$${\gamma }_{\mathrm{cas}}\sim \frac{\delta v_{\mathrm{A}}}{r}.$$

$${\dot{E}}_{\mathrm{damp}}\sim \frac{\delta v_{\mathrm{A}}}{r}|\delta B_{\perp }{|}^2.$$

$$\frac{{\omega }_{\mathrm{A}}^2}{\Omega }\sim \frac{\delta v_{\mathrm{A}}}{r}.$$

$$\delta B_{\perp }\sim \frac{{\omega }_{\mathrm{A}}}{\Omega }{B}_{\varphi },$$

$${\dot{E}}_{\mathrm{damp}}\sim \frac{{\omega }_{\mathrm{A}}^4}{{\Omega }^3}|B_{\varphi }{|}^2.$$

$${\dot{E}}_{\mathrm{damp},\mathrm{pole}}\lesssim \frac{{\omega }_{\mathrm{A}}^4}{{\Omega }^3}{B}_{\varphi }^2,$$

$$\frac{\partial }{\partial t}{E}_{\mathrm{mag}}\sim B_{\varphi }\frac{\partial }{\partial t}{B}_{\varphi }\sim q\Omega B_{\varphi }{B}_r.$$

$$q\Omega B_{\varphi }{B}_r\sim \frac{{\omega }_{\mathrm{A}}^4}{{\Omega }^3}|B_{\varphi }{|}^2.$$

$$\frac{B_r}{B_{\varphi }}\sim \frac{{\omega }_{\mathrm{A}}}{N_{\mathrm{eff}}}.$$

$$\frac{B_{\phi}}{\sqrt{4\pi\rho r^{2}}}=\omega_{\rm A}\sim\Omega\bigg{(}\frac{q\Omega}{N_{\rm eff}}\bigg{)}^{\!1/3}\,,$$

$$\frac{B_{r}}{\sqrt{4\pi\rho r^{2}}}\sim\Omega\bigg{(}\frac{q^{2}\Omega^{5}}{N_{\rm eff}^{5}}\bigg{)}^{\!1/3}\,,$$

$$\frac{\delta B_{\perp}}{\sqrt{4\pi\rho r^{2}}}\sim\frac{\delta v_{\rm A}}{r}\sim\Omega\bigg{(}\frac{q\Omega}{N_{\rm eff}}\bigg{)}^{\!2/3}\,,$$

$$\frac{\delta v_{\perp }}{r}\sim \Omega \frac{q\Omega }{N_{\mathrm{eff}}},$$

$$\frac{\partial}{\partial t}\delta\vec{B}=(\delta\vec{B}\cdot{\mbox{\boldmath$\nabla$}})\vec{v}+(\vec{B}_{\phi}\cdot{\mbox{\boldmath$\nabla$}})\delta\vec{v}\,.$$

$$\frac{\partial }{\partial t}⟨\delta \vec{B}⟩=0.$$

$$\frac{\partial}{\partial t}\langle\delta\vec{B}\rangle=\langle(\delta\vec{B}\cdot{\mbox{\boldmath$\nabla$}})\delta\vec{v}\rangle\,.$$

$$\frac{\partial }{\partial t}⟨\delta B_r⟩\sim \frac{⟨\delta v_{\perp }\delta B_{\perp }⟩}{r}.$$

$$r^2{k}_r^2{B}_r^2\gtrsim B_{\varphi }^2.$$

Full text

1 Introduction

One of the longstanding problems in stellar astrophysics is the nature of angular momentum (AM) transport within evolving stars. After the main sequence, the stellar core contracts and spins up, while the envelope expands and spins down. The differential rotation may source various (magneto)-hydrodynamical instabilities that can transport AM outwards to slow the rotation of the stellar core, with crucial consequences for the spins of white dwarfs (WDs), neutron stars, and black holes. However, the AM transport mechanisms at work remain controversial and enigmatic.

Asteroseismic observations have revolutionized this field by measuring internal stellar rotation rates for stars at various stages of evolution. Helioseismic inversions reveal nearly rigid rotation in the Sun’s radiative zone (Howe, 2009; Gough, 2015). For low-mass ($M\!\lesssim\!3\,M_{\odot}$) stars, internal rotation rates have been measured on the main sequence (Kurtz et al., 2014; Saio et al., 2015; Benomar et al., 2015; Van Reeth et al., 2018), sub-giant/red giant branch (RGB) (Beck et al., 2012; Mosser et al., 2012; Deheuvels et al., 2014; Triana et al., 2017; Gehan et al., 2018), red clump (Mosser et al., 2012; Deheuvels et al., 2015), and finally in WD remnants (Hermes et al., 2017). The conclusion drawn from these measurements is unambiguous: core rotation rates are relatively slow, and the vast majority of AM is extracted from stellar cores as they evolve. An efficient AM transport mechanism must be at work, causing cores and compact remnants to spin orders of magnitude slower than they would in the absence of AM transport.

In fact, the spin rates red giant cores and WDs are slower than theoretically predicted by nearly all AM transport mechanisms (Cantiello et al., 2014; Fuller et al., 2015; Belkacem et al., 2015; Spada et al., 2016; Eggenberger et al., 2017; Ouazzani et al., 2018). The MHD instability known as the Tayler-Spruit dynamo (Spruit, 2002) can provide more efficient AM transport than most other mechanisms, but prior implementations still predict spin rates roughly an order of magnitude too large because they struggle to overcome the steep composition gradient in red giants that suppresses AM mixing (Cantiello et al., 2014). Magnetorotational instability (Balbus & Hawley, 1994) may operate in some stars (e.g., Kagan & Wheeler 2014; Wheeler et al. 2015; Rüdiger et al. 2015) but it is also inhibited by composition gradients and thus has difficulty operating in red giants. Another possibility is that magnetic fields enforce rigid rotation in radiative regions of stars (Mestel, 1953), but that differential rotation develops in deep convective envelopes (Kissin & Thompson, 2015a), discussed further in Section 4. See Aerts et al. (2018) for a review of asteroseismic rotation rates and angular momentum transport mechanisms.

In this paper, we re-investigate the physics of the Tayler instability and its resulting saturation, as described in the seminal paper by Spruit 2002 (see also references therein Acheson & Gibbons 1978; Pitts & Tayler 1985; Spruit 1999; Braithwaite 2006; Denissenkov & Pinsonneault 2007; Zahn et al. 2007). We show that the instability can persist in RGB stars despite the existence of strong composition gradients, and we argue that its growth will saturate in a different manner than proposed by Spruit (2002). In our formulation, the instability can grow to larger amplitudes and produce stronger magnetic torques. We develop a convenient prescription for the effective AM/chemical diffusivity created by the instability and implement it into stellar evolution models. The core rotation rates of these models roughly match those observed in main sequence stars, red giant cores, and WDs. Hence, if the Tayler instability operates as we propose, it may largely solve the AM transport problem in stellar interiors.

2 Tayler Instability

Here we analyze the onset, growth, and saturation of the the Tayler instability. We follow the heuristic description of Spruit (2002) and use the same notation, but we address subsequent criticism by Denissenkov & Pinsonneault (2007) and Zahn et al. (2007). We begin by describing the main, generally agreed upon features of the Tayler instability and summarize how this instability is typically argued to saturate via the Tayler-Spruit dynamo. We then present an alternative way of thinking about the saturation, for which we explore the corresponding AM transport.

2.1 Instability Basics

The instability is analyzed in a rotating frame such that the local velocity is zero, but the local shear is finite. We make a number of standard assumptions that are appropriate in the context of stellar interiors, including (1) the main background component of the magnetic field is $B_{\phi}$ with corresponding Alfvèn frequency $\omega_{\rm A}=B_{\phi}/\sqrt{4\pi\rho r^{2}}$, (2) the angular rotation frequency $\Omega$ is roughly constant on spherical shells since horizontal turbulence can redistribute angular momentum latitudinally much faster than it can radially, and (3) the key frequencies are ordered such that $\omega_{\rm A}\ll\Omega\ll N$, where $N$ is the Brunt-Väisälä frequency.

As shear winds the magnetic field, $B_{\phi}$ grows and becomes Tayler unstable when it reaches a critical strength of (Spruit, 2002; Zahn et al., 2007)

[TABLE]

where $\eta$ is the magnetic diffusivity. The corresponding growth rate of this instability is largest for $m=1$ perturbations and is approximately

[TABLE]

Due to the strong stratification in these stars, the radial length scale of the instability is limited to

[TABLE]

while the maximum horizontal length scale of the instability is $l_{\perp}\sim r$.111Although Denissenkov & Pinsonneault (2007) argue the instability operates on shorter length scales, we demonstrate in a forthcoming paper (Ma & Fuller 2019, in preparation) $l_{\perp}\sim r$ is generally appropriate by deriving the dispersion relation at non-polar latitudes. The critique by Denissenkov & Pinsonneault (2007) is incorrect because it confuses the instability length scale $l_{\perp}\sim 1/k_{\perp}$ with the displacement amplitude $\xi_{\perp}$.

At the short radial length scales characteristic of Tayler instability in red giants, thermal diffusion is efficient so that the thermal stratification is largely mitigated (see Section 3.1). The main effect of this can be replicated by replacing $N$ in the above expressions with an “effective” Brunt-Väisälä frequency

[TABLE]

where $K$ is the thermal diffusivity, $N_{T}^{2}$ is the thermal component of the stratification, and $N_{\mu}^{2}$ is the compositional component. Red giant cores have large composition gradients, so $N_{\rm eff}$ in much of the core (and especially at the hydrogen-burning shell, the bottleneck for AM transport) is dominated by its compositional component, and thus $N_{\rm eff}\simeq N_{\mu}$. Appendix C discusses the appropriate value of $N_{\rm eff}$ when thermal diffusion is moderately important.

2.2 Saturation via the Tayler-Spruit Dynamo

The saturation of the Tayler instability is crucial for understanding the strength of the AM transport and chemical mixing it generates.The linear instability calculation allows us to determine the rate at which energy is transferred from background fields to perturbed fields, but energy dissipation only results from non-linear effects. This non-linear energy dissipation rate is necessary for calculating the mean amplitudes of the background and perturbed fields. One possibility is that the instability grows until it reaches a statistically stationary state in which the turbulent velocity field produces an effective viscosity or magnetic diffusivity large enough to balance the linear growth rate of the instability (Spruit, 2002). Equating the turbulent damping rate $\gamma_{\rm turb}$ with the linear growth rate results in

[TABLE]

where the wavenumber $k_{r}$ is the minimum required for instability, $k_{r}\sim\omega_{\rm A}/(Nr)$, and $\eta_{\rm eff}$ is an effective turbulent diffusivity. Next, since the azimuthal field grows via winding by shear as

[TABLE]

where $B_{r}$ is the radial field, then the amplification rate is $\gamma_{\rm amp}=q\Omega B_{r}/B_{\phi}$. The incompressible nature of the instability implies that $k_{r}B_{r}\sim k_{\perp}B_{\phi}$ and thus $B_{r}\sim(\omega_{\rm A}/N)B_{\phi}$, such that

[TABLE]

If the azimuthal field $B_{\phi}$ is turbulently damped at the rate of equation 5, setting equations (5) and (7) equal determines the azimuthal field strength at saturation $B_{\phi}/\sqrt{4\pi\rho r^{2}}\sim\omega_{\rm A}\sim q\Omega^{2}/N$. The radial field strength is then $B_{r}/\sqrt{4\pi\rho r^{2}}\sim q^{2}\Omega^{4}/N^{3}$, so that

[TABLE]

is the resulting Maxwell stress.

An important issue regarding this picture (as pointed out by Zahn et al., 2007) is that to linear order the Tayler instability grows fastest in the non-axisymmetric $m=1$ mode. Therefore the radial field generated by the instability is non-axisymmetric, and winding of this field produces no net increase in the axisymmetric component of $B_{\phi}$. Thus the axisymmetric component of $B_{r}$ is not necessarily related to the axisymmetric component of $B_{\phi}$ via $B_{r}/B_{\phi}\sim\omega_{\rm A}/{N}$.

A second potential issue is that equation (5) may not predict the correct damping rate for a large-scale background field $B_{\phi}$ that varies on lengthscales much larger than $1/k_{r}$. If the background field $B_{\phi}$ is essentially constant on this lengthscale, displacements do not mix background field lines of opposite polarity such that reconnection or dissipation occurs. Loops of background field can dissipate via reconnection if they migrate to a pole of a star where the loop has a small spatial scale, but we show in Appendix B that this mechanism produces a damping rate much smaller than equation 5. Hence, we believe equation 5 overestimates the decay rate of any large scale component of $B_{\phi}$, and the saturated values of $B_{r}$ and $B_{\phi}$ can be larger than those above.

Figure 1 presents a schematic for understanding the saturation of the Tayler instability as envisaged by Spruit (2002), and our proposed modifications discussed below.

2.3 Saturation via Magnetic Cascade

Motivated by these difficulties of calculating a turbulent/non-linear energy damping rate, we explore how turbulent cascades damp energy from the fluctuating fields $\delta\vec{B}$ and $\delta\vec{v}$. Tayler instability transfers energy from large scale magnetic fields to perturbed fields $\delta\vec{B}$ that vary on the short length scale $\sim 1/k_{r}$. In the linear regime, $\delta\vec{B}$ and $\delta\vec{v}$ are related to each to each other via

[TABLE]

where $\vec{\xi}$ is the Lagrangian displacement associated with the instability. Using $\delta\vec{v}\simeq-i\omega_{\rm R}\vec{\xi}$, where $\omega_{\rm R}\sim\omega_{\rm A}^{2}/\Omega$ is the real part of the perturbation frequency, and $\vec{k}\cdot\vec{B}_{\phi}\simeq k_{\phi}B_{\phi}\simeq imB_{\phi}/r\sin\theta$, we have

[TABLE]

Here we have used $m=1$ and ignore geometric terms of order unity, and $\delta\vec{v}_{\rm A}=\delta\vec{B}/\sqrt{4\pi\rho}$ is the perturbed Alfvén velocity. A similar answer can be obtained by analyzing the momentum equation

[TABLE]

The dominant forces in the horizontal direction of equation (11) are the Coriolis and Lorentz terms. Therefore we expect quasi-magnetogeostrophic balance, as found in rapidly rotating convective simulations by Augustson et al. (2016), such that

[TABLE]

where $\delta v_{\perp}$ and $\delta v_{\rm A,\perp}$ are the horizontal components of the perturbed velocity and Alfvén velocity. Since we shall find $\omega_{\rm A}\ll\Omega$ in most stellar applications, the perturbation energy is dominated by magnetic rather than kinetic energy.

Understanding how energy cascades to small (or large) scales in MHD turbulence is tricky business. We look to Goldreich & Sridhar (1995); Lithwick & Goldreich (2003); Chandran (2004); Lithwick et al. (2007) for guidance, though these studies did not include the effects of stratification and rotation. In Appendix A, we attempt to account for Coriolis and buoyancy forces on the Alfvénic cascade rate to smaller spatial scales, finding

[TABLE]

Similar to the weak Alfvénic turbulence described by Lithwick & Goldreich (2003), equation 13 is determined by the rate at which energy is transferred to smaller scales when Tayler modes scatter off one another. We assume magnetic energy cascades from the large scales of the instability to small scales where it is damped, such that the cascade rate $\gamma_{\rm cas}$ effectively represents a turbulent damping rate of the perturbed magnetic energy. The non-linear energy dissipation rate is then

[TABLE]

We do not expect energy in the background field $B_{\phi}$ to be damped by a turbulent cascade to small scales. This is a key difference from Spruit (2002), who uses a damping rate $\dot{E}_{\rm damp}\sim\gamma_{\rm turb}|B_{\phi}^{2}|\sim(\omega_{\rm A}^{2}/\Omega)|B_{\phi}^{2}|$. We believe this is unphysical because the Alfvénic turbulence does not cause magnetic energy in the background field to cascade to small scales, it is only the Alfvén waves traveling along the background field (i.e., Tayler modes) that cascade to small scales where they can be damped.

As in Spruit (2002), we assume the instability saturates and reaches a statistically stationary state when the instability growth rate is matched by the turbulent damping rate such that

[TABLE]

Note that equation 15 implies that upon saturation, the perturbed and background field are related by

[TABLE]

so that the energy damping rate is

[TABLE]

Energy in the background field *can *be damped if field loops can reconnect with loops of opposite polarity, which can occur sufficiently close to the pole of the star where the loops have a small spatial scale. In Appendix B, we show that the maximum possible energy damping rate due to this effect is

[TABLE]

which is less than or equal to the energy damping rate of equation 17. Hence, both mechanisms may contribute to saturation of the instability, but equation 17 is always a good estimate of the total energy damping rate.

Next, it is useful to consider the flow of energy in this system, which is as follows.

Rotational shear energy is converted to magnetic energy by winding a radial field into a toroidal field. 2. 2.

Toroidal field energy is converted by the Tayler instability into magnetic/kinetic energy associated with the perturbed magnetic/velocity field. 3. 3.

These perturbations are damped into heat by a turbulent cascade.

At equilibrium, the energy input by winding must equal the turbulent energy dissipation rate. The energy input by winding is

[TABLE]

Equating the energy input rate to energy damping rate, we have

[TABLE]

To solve our system, we need an estimate of $B_{r}/B_{\phi}$. In Section 2.3.1, we argue that $B_{r}$ can grow until Lorentz forces nearly stabilize the plasma against the growth of the Tayler instability, such that

[TABLE]

This is the same ratio used by Spruit (2002), but it arises for different reasons. Combining this with equations (15) and (20), we expect the turbulent damping to saturate the Tayler instability at

[TABLE]

[TABLE]

[TABLE]

[TABLE]

These fields can then drive AM transport and chemical mixing as further described in Section 2.4.

2.3.1 The Importance of Non-linear Induction

Before providing prescriptions that can be used for stellar evolution calculations, it is helpful to address some conceptual challenges associated with this new approach to the saturation of the Tayler instability. Initially, $B_{r}$ can be due to a small seed field, but as this field is converted to a toroidal field and dissipation occurs, it must be replenished. This new $B_{r}$ can then continue to be wound by the shear, and continue the flow of energy as outlined above. The question is how this new $B_{r}$ is generated and how strong can it grow.

As argued in Section 2.2, closing the loop between $B_{r}$ and $B_{\phi}$ is difficult if only linear effects are considered. This can be seen by starting with the linearized induction equation,

[TABLE]

and taking the azimuthal average, which yields

[TABLE]

This is because the perturbed magnetic/velocity field is non-axisymmetric ($m=1$) to linear order, while the background magnetic field is axisymmetric. Equation (27) conveys the argument by Zahn et al. (2007) that winding of the non-axisymmetric field cannot regenerate the axisymmetric toroidal field.

However, there can be growth of the axisymmetric radial field if we include non-linear terms in the induction equation. Perturbing the induction equation to second order and taking the azimuthal average yields

[TABLE]

To order of magnitude, we thus expect

[TABLE]

Hence, we expect some growth of an axisymmetric radial field due to non-linear induction, i.e., an $\alpha$-dynamo effect.

However, as shown by Braithwaite (2009), Tayler instability cannot operate if $B_{r}$ rises above a threshold value. This occurs if magnetic tension forces due to perturbation of the radial field are larger than magnetic pressure forces driving the instability, which can be expressed as

[TABLE]

If $B_{r}$ grows until the instability is quenched, then the maximum possible value of $B_{r}$ (corresponding to the longest length scale unstable disturbance $rk_{r}\sim N_{\rm eff}/\omega_{\rm A}$) is

[TABLE]

This is identical to the condition arising from incompressibility used by Spruit (2002), but it relates the axisymmetric component of $B_{r}$ and $B_{\phi}$, whereas Spruit’s relation is only valid for the non-axisymmetric component of $B_{r}$.

2.4 Angular Momentum Transport

The torque via Maxwell stresses in the saturated state is found from combining equations (22) and (23),

[TABLE]

corresponding to an effective AM diffusivity

[TABLE]

Although these scalings apply in the case of magnetic energy dissipation balance, it is difficult to predict the exact prefactors using these analytic arguments. We therefore parameterize our result via the saturated Alfvén frequency, using a dimensionless parameter $\alpha$ such that

[TABLE]

The parameterized AM diffusivity is then

[TABLE]

We expect $\alpha$ of order unity, and indeed in Section 3 we find $\alpha\approx 1$ fits the observational data.

Combining the instability criterion given by equation (1) with the value of $\omega_{\rm A}$ in the saturated state implies a minimum shear in order for the instability to occur and saturate as outlined above. Equating (1) and (34), we find

[TABLE]

We show in Section 3 that this minimum shear appears to frequently be realized in red giant stars, such that $\omega_{\rm A}\sim\omega_{c}$ in most of the core. In this case, the core rotation rates are set mostly by the structure of the star (i.e., profiles of $N_{\rm eff}$ and $\eta$) and are very insensitive to the initial rotation rate or prior evolution of the star.

2.5 Energetics and Mixing

Note that our relations at saturation imply a hierarchy of rotational, background magnetic, perturbed magnetic, and kinetic energy densities:

[TABLE]

These hierarchies are true as long as $q\lesssim N_{\rm eff}/\Omega$, which is true in our models in Section 3 where $q\sim 1$ and $N_{\rm eff}/\Omega\sim 10^{4}$. However, in cases where $q$ is much larger the hierarchy will be altered, and the instability could saturate in a different manner.

From the divergence-free conditions on the perturbed magnetic and velocity fields, we can also calculate their radial components:

[TABLE]

[TABLE]

The radial components of the fields are typically orders of magnitude smaller than the horizontal components due to the tiny value of $\Omega/N_{\rm eff}$ in most stars.

For this reason, chemical mixing induced by the Tayler instability will likely be less important than AM transport in most stars. The effective chemical mixing diffusivity is

[TABLE]

Using the relations above, we have

[TABLE]

so that

[TABLE]

for the ratio of chemical mixing to AM transport.

In red giants, we find $\nu_{\rm mix}/\nu_{\rm AM}\sim 10^{-6}$, such that chemical mixing caused by the Tayler instability is minuscule. The chemical mixing timescale across the star is longer the Ohmic diffusion time scale, which is longer than the stellar evolution timescale (Cantiello et al., 2016), so the chemical mixing is likely negligible. The scaling of 42 is stronger than that of Eddington-Sweet circulation, so we expect chemical mixing from Tayler instabilities to be unimportant unless $q\gg 1$ or $N_{\rm eff}\ll N$.

3 Stellar Models

With our prediction for AM transport due to Tayler torques, we implement this prescription into stellar evolutionary models to predict their internal rotation rates. We then compare with asteroseismic measurements of internal rotation rates, finding generally good agreement.

3.1 Properties of Red Giant Cores

An important feature of post-main sequence stars is their steep composition gradient in and above their hydrogen burning shells. Figure 2 shows a $M=1.2\,M_{\odot}$ model on the lower RGB at $R=4.1\,R_{\odot}$ and log(g)$=3.3$. At the hydrogen burning shell, the stabilization is primarily due to the hydrogen-helium composition gradient such that $N\simeq N_{\mu}$, but even above the burning shell, we often find $N_{\mu}\sim N/5$ due to the hydrogen gradient left behind by partial pp-chain burning during the main sequence. Hence, the compositional part of the stratification is very important, even above the burning shell.

An important consideration for the operation of the Tayler instability is whether thermal diffusion will undermine the thermal component of $N^{2}$. It is useful to compare the growth rate of the instability with the thermal diffusion time scale at the instability lengthscale,

[TABLE]

Here, we have used the maximum radial lengthscale for Tayler instability $l_{r}\sim r(\omega_{\rm A}/N_{\rm eff})$, and the thermal diffusivity $\chi=16\sigma_{\rm SB}T^{3}/(3\rho^{2}c_{v}\kappa)$. Thermal diffusion strongly reduces the effective thermal stratification when $\gamma\gtrsim\omega$, where $\omega$ is the real part of the frequency of the overstable oscillations. Zahn et al. (2007) show that $\omega\sim\omega_{\rm A}^{2}/\Omega$. Using our saturated field strength (equation 34), we find thermal diffusion is very important when

[TABLE]

A comparison of $\omega_{t}$ and $\Omega$ in Figure 2 shows that the former is larger throughout the radiative core, such that thermal diffusion is very important. This is almost always the case in our post-main sequence models. In this case, as discussed by Zahn et al. (2007), the effective stabilization is given by equation 4. In most regions of our models, we find $N_{\rm eff}^{2}\simeq N_{\mu}^{2}$.

3.2 Comparison with Measurements

We expect the AM diffusivity of equation (35) to capture the scaling of magnetic torques in differentially rotating stars, but we must still evaluate the appropriate value of $\alpha$. To that end, we construct rotating stellar models using the MESA stellar evolution code (Paxton et al., 2011; Paxton et al., 2013, 2015, 2018). We assume rotation constant on spherical shells and near rigid rotation in convective zones. These models include AM transport via the diffusivity of equation (35) applied to gradients in angular rotation frequency, which is included if equation (36) is satisfied. Our models also include hydrodynamic AM transport mechanisms (which are usually negligible compared to our revised TS torques), but we do not use MESA’s default prescription for TS torques. A full inlist can be found in Appendix D.

To calculate core rotation rates $P_{\rm core}$ for comparison with asteroseismology, we compute the average core rotation period $P_{\rm core}=2\pi/\Omega_{\rm core}$ as sensed by a gravity wave in the WKB limit,

[TABLE]

The bounds of the integral in equation 46 correspond to the boundaries of the core gravity mode cavity where $\omega_{\rm g}\!<\!N$, and we consider gravity waves with frequency $\omega_{\rm g}=2\pi\nu_{\rm max}$, where the frequency of maximum power is calculated via classical scaling relations $\nu_{\rm max}=3090\,\mu{\rm Hz}\,(M/M_{\odot})(R/R_{\odot})^{-2}(T/T_{\odot})^{-1/2}$. For our WD models, we simply set $P_{\rm core}$ equal to the central rotation rate.

We evolve solar metallicity models ranging from $1.2-6\,M_{\odot}$ from the zero-age main sequence (ZAMS) to the WD phase. We initiated each model with a spin rate $P_{i}=2$ days, except for the $1.2\,M_{\odot}$ model for which we used $P_{i}=20$ days to account for main sequence magnetic braking. Figure 3 shows evolution of the core and surface rotation rates of a $1.6\,M_{\odot}$ model with a ZAMS rotation period $P_{\rm i}=2$ days. We also denote typical measured rotation rates of cores of stars ascending the RGB (Mosser et al., 2012; Gehan et al., 2018), stars on the red clump (Mosser et al., 2012), and WDs (Hermes et al., 2017), all of which descended predominantly from main sequence stars in the range $1\,M_{\odot}\lesssim M\lesssim 3\,M_{\odot}$. Typical core rotation rates are $P_{\rm core}\sim 10-20$ days on the lower RGB, $P_{\rm core}\sim 50-200$ days on the red clump, and $P\sim 0.5-4$ days for WDs.

Our models generally exhibit very similar rotation rates to observations for $\alpha\approx 1$, a reasonable value since we expect $\alpha\sim 1$. The agreement is very good along the RGB, red clump, and in the WD phase. We also plot the model’s surface rotation rate, which shows that nearly rigid rotation is maintained beyond the end of the main sequence. The rigid rotation is maintained until the ratio of $\Omega/N_{\rm eff}$ becomes sufficiently small that the Tayler instability cannot fully prevent the spin-up of the core. After this point, differential rotation develops during the late sub-giant/early red giant phase as the core contracts and tries to spin up, while the envelope expands and spins down. However, the AM transport is strong enough that the core rotation period actually increases between the main sequence and the tip of the RGB. In contrast, a model including the default prescription for TS torques exhibits core spin-up along the RGB, spinning an order of magnitude too fast compared to observations, in agreement with the results of Cantiello et al. (2014). Models with only hydrodynamic prescriptions for AM transport have even faster core rotation and are totally incompatible with observations.

Our models diverge from those with different AM transport prescriptions along the sub-giant branch and lower RGB because most of the core AM extraction in our models occurs during these phases. At later stages of evolution of low-mass stars ($M\lesssim 2\,M_{\odot}$), the stabilizing composition gradients are so large in comparison to the local rotation rates that very little AM transport occurs after the RGB bump. This result agrees with Cantiello et al. (2014), who find that red clump and WD rotation rates require approximate conservation of core AM after the RGB bump.

To understand the dependence of our results on the parameter $\alpha$ and the star’s initial spin rate, Figure 4 shows core rotation rates for models with different values of $\alpha$ and $P_{i}$. Remarkably, the core spin rate on the RGB and red clump is nearly independent of the initial spin. The reason is the strong dependence of AM transport on the local spin rate, with $\nu_{\rm AM}\propto\Omega^{3}$. Rapidly rotating cores experience a stronger spin-down torque while slowly rotating cores feel a weaker spin-down torque, causing convergent migration in the post-main sequence core spin rate. Note that the WD spin rate does exhibit some dependence on initial spin rate, largely because this determines the AM of the accreted material on the clump and asymptotic giant branch (AGB).

Figure 4 shows that the post-main sequence spin period is roughly proportional to $\alpha$. Nearly rigid rotation is maintained along the main sequence, regardless of $\alpha$, except for very slow rotators in which the value of $\nu_{\rm AM}$ is much smaller (see discussion in Section 4). The models slightly diverge from one another on the lower RGB, with smaller values of $\alpha$ allowing faster core rotation. We find the main effect of the value of $\alpha$ in our models is not the prefactor in equation 35, but rather in determining the minimum shear $q_{\rm min}$ (equation 36) required for Tayler instability to saturate as we have outlined. When $q>q_{\rm min}$, efficient AM transport generally decreases the core rotation and shear, thereby reducing $q$ until $q\sim q_{\rm min}$, as shown in Appendix D. In this limit, equation 36 predicts that the core rotation scales as $\Omega\propto\alpha^{-12/13}$, in line with our numerical finding that the core spin period is approximately proportional to $\alpha$.

Figure 5 shows the rotation profile of a $1.2\,M_{\odot}$ model at the base of the RGB for $\alpha=1$. The shear is strongest at the hydrogen burning shell where $N_{\mu}$ is largest. There is very little shear in the helium core where almost no compositional stratification exists. Significant shear also exists in the radiative region above the burning shell due to the composition gradient left over from incomplete hydrogen burning outside the central core during the main sequence. Encouragingly, this rotational profile is very similar to that inferred for the RGB star KIC 4448777 (Di Mauro et al. 2016, Figure 11) at nearly the same phase of evolution, though we caution that the actual rotation profile is poorly constrained by asteroseismic data.

Our model makes important predictions for core rotation rates as a function of progenitor mass, as shown in Figure 6. We find core rotation rates on the RGB in the range $10\,{\rm days}\!\lesssim P_{\rm core}\!\lesssim\!30\,{\rm days}$ regardless of mass. On the clump, the differences between stars of different masses are slightly larger: our models predict slower core rotation rates for low-mass clump stars, with $P_{\rm core}\sim 100-200\,{\rm days}$ for $\approx 1.2\,M_{\odot}$ stars. We predict faster rotation for secondary clump stars, with $P_{\rm core}\sim 50$ days for $M\approx 2.2\,M_{\odot}$. The trend of faster rotation for higher-mass clump/secondary clump stars indeed appears to be present in the results of Mosser et al. (2012). We also predict very mild core spin-down during helium-burning for lower-mass stars, whereas we predict significant core spin-down during the helium-burning phase of secondary clump stars (by a factor of $\sim\!2$).

We also make predictions for WD spin rates as a function of WD mass. Figure 7 shows asteroseismicly measured WD spin periods from Hermes et al. (2017) as a function of WD mass, along with predictions from our stellar models evolved down the WD cooling track to the ZZ-ceti instability strip. The ZAMS models have masses $M_{\rm ZAMS}=(1.2,1.6,2.2,3.0,4.0,5.0,6.0)\,M_{\odot}$ and produce carbon-oxygen WDs with masses $M_{\rm WD}=(0.54,0.56,0.58,0.64,0.81,0.87,0.95)\,M_{\odot}$. Figure 7 demonstrates that our predicted WD spin rates are very similar to those observed, with $P_{\rm WD}\sim 1-3$ days for WDs with $M\approx 0.6\,M_{\odot}$. W predict that more massive WDs rotate faster, a trend indeed observed in Hermes et al. (2017), but our highest mass models rotate much slower than the observed high-mass ($M\gtrsim 0.7\,M_{\odot}$) WDs.

In general, the observed population on WDs appears to exhibit more scatter than our model predictions, some of which may be inherited from the scatter in progenitor rotation rate as shown in Figure 4. Additionally, our models do not take into account binary effects such as mergers (either during stellar evolution or WD mergers) that may produce faster rotating stellar cores and WDs. Kilic et al. (2018) suggest that $\sim\!10\%$ of WDs, especially higher mass WDs, are likely to be merger products. We speculate some of the faster rotating WDs shown in Figure 7 resulted from stellar mergers during post-main sequence evolution. Finally, our models do not make reliable predictions for descendants of magnetic Ap/Bp stars, whose strong internal fields likely increase AM transport and may keep their cores more slowly rotating than our predictions.

4 Discussion

Our AM transport prescription predicts extremely short AM transport times $t_{\rm AM}$ within radiative zones of main sequence stars,

[TABLE]

For a fast rotating young Sun, we find that the instability occurs in the nearly adiabatic limit, such that $N_{\rm eff}\approx N$. Evaluating equation (47) in the radiative zone of a young solar model rotating at $P=3$ days yields a typical AM transport timescale $t_{\rm AM}\sim 10$ years. In the current Sun, the rotation rate is much slower and $t_{\rm AM}\sim 10^{4}\,{\rm years}$, but this still enforces nearly rigid rotation in agreement with helioseismic measurements (Howe, 2009; Gough, 2015). In most cases, we predict nearly rigid rotation for main sequence stars, although modest differential rotation may exist in very slowly rotating stars. For rotation rates of about 100 days, equation (47) predicts $t_{\rm AM}\sim 10^{6}\,{\rm yr}$, which may be longer than the timescale for shear to develop due to other effects such as internal gravity waves (Rogers et al., 2013; Fuller et al., 2014; Townsend et al., 2018). Hence, differential rotation can persist in slowly rotating stars, and this could explain why some very slowly rotating stars (see Kurtz et al. 2014; Saio et al. 2015; Triana et al. 2015; Kallinger et al. 2017; Sowicka et al. 2017) appear to exhibit some degree of differential rotation, while more rapidly rotating stars main sequence stars appear to be nearly rigidly rotating (Aerts et al., 2018).

The short AM transport time for main sequence stars may seemingly contradict observations of rotational evolution of young $\approx\!1.0\,M_{\odot}$ stars, for which several works (e.g., Denissenkov et al. 2010; Gallet & Bouvier 2015; Lanzafame & Spada 2015) find evidence for core-envelope coupling times in the range $t_{\rm AM}\!\sim\!10-100$ Myr. These coupling times are deduced by fitting models including magnetic braking and core-envelope decoupling to the surface rotational evolution of cluster stars at a variety of ages. The wide distribution of surface rotation rates extending to ages older than $100$ Myr can only be fit using a value of $t_{\rm AM}\sim 10^{7}-10^{8}\,{\rm yr}$ that varies with mass and rotation rate. These models all utilize relatively simple and deterministic magnetic braking laws, but the bimodal rotation rates of low-mass cluster stars (e.g., Rebull et al. 2016, 2017; Rebull et al. 2018) cannot be explained by such models. Instead, it appears that magnetic braking is strongly influenced by surface magnetic field morphology, such that rapidly rotating stars with more complex fields can spin down more gradually than slower rotating stars with mostly dipolar fields, which can explain the bimodality and rotational evolution of low-mass cluster stars (Brown, 2014; Garraffo et al., 2015, 2016, 2018), even assuming rigid internal rotation. In light of the relatively short AM transport time scales needed to explain the slow rotation of red giant cores, we find it most likely that $t_{\rm AM}$ is indeed very small for main sequence stars such that they rotate nearly rigidly, but that magnetic braking can be a more complex process than previously assumed, especially for young stars.

Our results have important consequences for mixing processes that depend on stellar rotation, such as meridional circulation and various shear instabilities. Mixing resulting directly from Tayler instability (equation 41) is typically quite small as long as there is a composition gradient such that $\Omega/N_{\rm eff}\ll 1$, which is often the case in radiative regions that have undergone any nuclear processing. An exception to this is horizontal mixing. Unless $q\ll 1$, the horizontal circulation given by equation 25 is more rapid than the Eddington-Sweet circulation. This helps to justify the assumption of Zahn (1992) that horizontal mixing is much faster than vertical circulation currents, and so supports the conclusion that vertical chemical mixing is slow relative to the vertical advection rate. More importantly, our models exhibit slower core rotation and smaller shears than previous predictions, resulting in less mixing via shear instabilities. We thus suspect that rotational mixing has been overestimated in many previous works. To quantify this statement, more thorough calculations must be performed, incorporating AM/chemical transport via Tayler instability, meridional circulation, shear instabilities, convective overshoot, etc. The coupled effects of AM tranpsort and mixing can then be compared with abundance/rotation measurements (see e.g., Somers & Pinsonneault 2016) for stars in clusters. We hope to explore mixing effects and make detailed predictions for surface abundances in future work.

Our models assume that convection zones are nearly rigidly rotating, which may not be true for deep convective zones where asymmetric convective energy/AM fluxes may cause deeper layers of the convective envelope to rotate faster (Brun & Palacios, 2009; Kissin & Thompson, 2015a). Indeed, some degree of envelope differential rotation may be necessary to explain rotation rates of horizontal branch stars (Sills & Pinsonneault, 2000). However, envelope differential rotation does not always change our predictions for rotation in radiative cores for two reasons. First, the strong dependence of AM transport on local rotation rate, $\nu_{\rm AM}\propto\Omega^{3}$, causes core rotation to converge to a rate only weakly dependent on surface rotation rate. Second, in our models we find the core rotation rate often converges to a state marginally unstable to Tayler instability such that the core rotation rate is set by equation (36). Preliminary tests indicate envelope differential rotation may allow for slightly faster rotation rates of WDs, similar to the effect of decreasing the initial spin period and allowing the core to accrete more AM during the AGB.

Currently the most viable alternative to our model is that of Kissin & Thompson (2015a), which posits rigid rotation in radiative zones enforced by magnetic torques and differential rotation in convective zones due to convective AM pumping. The clear prediction from our model is that differential rotation is mostly in the core, while the Kissin & Thompson (2015a) model predicts differential rotation confined to the evelope. Asteroseismic observations appear to disfavor envelope differential rotation (Di Mauro et al., 2016; Klion & Quataert, 2017; Di Mauro et al., 2018), though currently their ability to distinguish between the models is limited. Both models may have some tension with observations, as Kissin & Thompson (2015a) predicts rotation rates that are too slow for low-mass ($M\lesssim 1.2M_{\odot}$) RGB stars, and Kissin & Thompson (2015b) appear to predict anomalously slow rotation rates for some white dwarfs.

A potential problem with our mechanism is that it may underpredict the scatter in observed core rotation rates, as our models converge to a similar rotation rate regardless of initial conditions. Additionally, we predict significant spin-up of red giant cores (by a factor of $\sim$2) along the lower RGB, whereas Gehan et al. (2018) find no clear spin-up/spin-down as a function of evolutionary state. More work predicting spin rates for a population of stars (incorporating changes in initial spin-rate, metallicity, binarity, etc.), along with more asteroseismic measurements and a better understanding of measurement biases222We are concerned that measurement bias limits the number of core rotation measurements for stars with more rapidly rotating cores higher up the RGB, where the rotational frequency splitting becomes comparable to the mixed mode period spacing, and the asteroseismic power spectrum becomes difficult to interpret (Deheuvels et al., 2017). will help to distinguish between the competing models, though we note that differential rotation in both the convective envelope and the radiative core is possible.

Another obstacle for our model is that the radiative core must have very weak fossil fields in order for the Tayler instability to dominate AM transport. Figure 8 plots various components of the magnetic fields in a stellar model. Of particular importance is the radial field $B_{r}$, which we predict to have a strength of $B_{r}\sim 10^{-2}\,{\rm G}$ through much of the radiative core. Recall that for the predicted value of $B_{\phi}$, Tayler instability cannot occur if there is a fossil field with strength greater than equation 23. If there is a fossil component of $B_{r}$ larger than that shown in Figure 8, the azimuthal component $B_{\phi}$ must be amplified by shear to larger field strengths before Tayler instability kicks in. By the time this occurs, the Maxwell stress $B_{r}B_{\phi}$ will be larger than predicted by our model, bringing the radiative core closer to a state of rigid rotation. So, even relatively weak fossil fields ($B_{r}\gtrsim 10^{-2}\,{\rm G}$) may enforce nearly rigid rotation of the radiative core. While the internal field strengths of red giants are not well known333While some red giants with suppressed dipole oscillation modes may have very strong ($B_{r}\gtrsim 10^{5}G$) magnetic fields (Fuller et al., 2015; Stello et al., 2016), those whose internal rotation has been measured must have weaker fields in order for gravity waves to propagate in their cores such that the core rotation rate can be measured., a rigidly rotating core enforced by fossil fields would necessitate large differential rotation in the convective envelope, as advocated by Kissin & Thompson (2015a).

5 Conclusions

The pioneering work of Spruit (2002) has shown that Tayler instabilities naturally occur in stellar interiors and may dominate internal angular momentum (AM) transport. However, the saturation of the instability, and the resulting AM transport, remain poorly understood. Whereas Spruit (2002) posits that energy in the background field is dissipated at the instability growth rate, we argue that Tayler instability saturates via the turbulent dissipation of unstable magnetic field perturbations. Our saturation mechanism results in a smaller energy dissipation rate, such that the magnetic fields reach larger mean amplitudes. The stronger fields produce larger Maxwell stresses and more efficient AM transport. Crucially, our proposed saturation condition does not depend on the closure of a dynamo loop and thus avoids the problems pointed out by Zahn et al. (2007). Another important difference is that the minimum shear for significant AM transport (equation 36) is smaller than that of equation (26) of Spruit (2002), and thus Tayler instability can occur at much lower shear as long as equation (1) is satisfied.

When Tayler instability operates, we find that it produces an effective AM diffusivity given by equation (35). In our models, the resultant torque often reduces the shear to a state of marginal stability given by equation (36). Implementation of our results into stellar evolution codes shows that a reasonable saturation parameter $\alpha\approx 1$ leads to core rotation rates in good agreement with asteroseismic measurements for main sequence stars, red giants, and white dwarfs across a wide range in mass. Hence, these findings may be a key step toward solving the AM transport problem within stars, and they open the door to realistic predictions of internal stellar rotation rates during phases of evolution prohibitively difficult to observe.

Our results have major implications for the core rotation rates of massive stars and their compact progeny. Prior estimates (Heger et al., 2005) based on the original TS dynamo prescription predicted neutron star rotation rates of $P_{\rm NS}\sim 10\,{\rm ms}$, somewhat faster than typical pulsar birth periods $P_{\rm NS}\gtrsim 10-50$ ms (Faucher-Giguère & Kaspi, 2006; Igoshev & Popov, 2013; Gullón et al., 2014). We expect our updated AM prescription to yield significantly slower NS rotation than prior predictions, and in future studies we will investigate the core rotation rates of evolving massive stars to predict the natal spin rates of neutron stars and black holes.

6 Acknowledgments

We thank the referee, Dr. Henk Spruit, for constructive and much-deserved criticism. We thank Peter Goldreich, Marc Pinsonneault, Matteo Cantiello, and Daniel Lecoanet for useful discussions. This research is funded in part by a Rose Hills Innovator Grant, the Gordon and Betty Moore Foundation through Grant GBMF7392, and by the National Science Foundation under Grant No. NSF PHY-1748958.

Appendix A Turbulent Magnetic Cascade

In this Appendix we summarize the nonlinear damping arguments of Lithwick & Goldreich (2003) and extend them to the instability discussed in this paper. We begin with the equations of incompressible MHD in the absence of rotation or buoyancy:

[TABLE]

The pressure $P$ is the total pressure, including both gas and magnetic contributions. If the magnetic field is composed of a large-scale slowly-varying component and a small-scale fluctuating one we may decompose it, or equivalently the Alfvén velocity, as

[TABLE]

where $\boldsymbol{v}_{{\rm A},0}$ is the slowly-varying background copmponent and $\delta\boldsymbol{v}_{{\rm A}}$ is the fluctuating one. With this decomposition we define combinations of fluid velocity and magnetic field fluctuation

[TABLE]

The equations of motion may then be cast as

[TABLE]

As usual in incompressible systems the pressure is not an independent degree of freedom, and may be used to ensure that the second of these equations is satisfied. The net result is that $P$ serves to project the non-linear term in the first equation into the subspace specified by the second equation.

Studying only linear terms, equation (54) just describes the advection of combined variations in the fluid velocity and the magnetic field. That is,

[TABLE]

These combined fluctuations evidently only propagate along the large-scale magnetic field. With the addition of the non-linear terms, we see that packets of $\boldsymbol{w}_{+}$ may scatter off of those of $\boldsymbol{w}_{-}$ and vice-versa, but because the linear evolution is constrained to be along $\boldsymbol{v}_{{\rm A},0}$ only those packets which are bound to the same field line may scatter.

When these scattering events occur, Lithwick & Goldreich (2003) showed that they result in a bending of the field lines by an angle of order $w_{\pm}/v_{{\rm A},0}$. By repeatedly bending the field-lines of a wave-packet of $\boldsymbol{w}_{\pm}$ the packet may be disrupted, such that its energy cascades non-linearly to different scales. If the packets of $\boldsymbol{w}_{\pm}$ have wavelength $\Lambda$ parallel to $\boldsymbol{v}_{{\rm A},0}$ and wavelength $\lambda$ transverse to it then this disruption occurs if the field lines are displaced transversely by an amount of order $\lambda$, or by an angle of order $\lambda/\Lambda$. If scattering events are equally likely to bend field lines in all directions the process of non-linear interactions may be described as a diffusive random walk with step size $w_{\pm}/v_{{\rm A},0}$, so that it takes of order

[TABLE]

scattering events to cause a cascade. Note that the “strong” regime of Lithwick & Goldreich (2003) just corresponds to the point where the above expression yields $N\la 1$. Because each scattering event takes time of order $\Lambda/v_{{\rm A},0}$, the time-scale over which wave-packets are disrupted is

[TABLE]

The specific energy in the system, neglecting the bulk magnetic field, is of order $w_{+}^{2}+w_{-}^{2}$. Both $w_{+}$ and $w_{-}$ at any given scale damp due to the cascade, with packets disrupting after time-scale of order $t_{\rm cas}$. It follows that the damping rate of the energy in the system is of orderr

[TABLE]

or equivalently, the non-linear loss rate is

[TABLE]

We now generalize this argument to the Taylor instability. In this case there are two additional phenomena that must be considered. First, in a rotating system the Coriolis effect adds an acceleration to equation (48), so that

[TABLE]

The new acceleration term is not fully absorbed into the pressure gradient because it is not fully directed along the wave-vector $\boldsymbol{k}$, and hence, at least to order of magnitude, it must be kept. Upon making the same change to equation (54) we find

[TABLE]

Crucially, this is no longer of the form of a simple advection equation like (56). Rather, packets of $w_{+}$ now mix with those of $w_{-}$ over time. This means that the non-linear interaction, which may only act between $w_{+}$ and $w_{-}$, may act on a single wave-packet as it transitions between $w_{\pm}$.

Assuming that $\boldsymbol{w}$ is not oriented nearly-parallel to the rotation axis, the time-scale over which wave-packets transition is

[TABLE]

Recalling that we work in the limit of magnetogeostrophic balance,

[TABLE]

so

[TABLE]

While in the mixed state, the cascade proceeds with time-scale

[TABLE]

Analogous to the strong and weak regimes of Lithwick & Goldreich (2003), we therefore find two regimes. In the first, $\omega_{\rm A}t_{\rm cas}^{*}\ll 1$, so that scattering occurs rapidly once the wave-packets mix. This results in an effective cascade time

[TABLE]

In the opposing limit, $\omega_{\rm A}t_{\rm cas}^{*}\ll 1$, scattering is slow and it makes sense to average $\boldsymbol{w}_{\mp}\cdot\nabla\boldsymbol{w}_{\pm}$ over the mixing time $\omega_{\rm A}^{-1}$. In effect $w_{+}$ and $w_{-}$ mix quickly, a small amount of scattering occurs, and then they are unmixed again. This repeats until the amount of scattering is of order unity, so

[TABLE]

Putting the two regimes together we find the overall effect of rotation is to reduce the cascade time-scale to

[TABLE]

The above argument may also be cast in terms of new linear combinations of $\boldsymbol{v}$ and $\delta\boldsymbol{v}_{\rm A}$ which do follow an advection-like equation. However these new linear combinations do not preserve the structure of the non-linear interaction term, and generically give rise to interactions of the form $\boldsymbol{w}_{+}\cdot\nabla\boldsymbol{w}_{+}$, and likewise for $\boldsymbol{w}_{-}$. These new self-interaction terms result, following the arguments above, in the same cascade time-cale $t_{\rm cas}\approx\Lambda/w_{\pm}$.

The second modification we must consider is that of buoyancy. This works in much the same way. We define

[TABLE]

as the Lagrangian displacement of a fluid element. With this, and working in the Boussinesq limit, equation (61) becomes

[TABLE]

where we have taken the entropy gradient to be in the radial direction and omitted the Coriolis effect for simplicity. Once more inserting the new acceleration into equation (54) we obtain

[TABLE]

We may approximate the displacement as

[TABLE]

where $\omega$ is the linear frequency associated with any given mode. Hence

[TABLE]

It follows that this term, like the Coriolis one, produces mixing between $w_{\pm}$. Because the real and imaginary parts of $\omega$ are both of order $\omega_{\rm A}^{2}/\Omega$, this may be written as

[TABLE]

Noting that

[TABLE]

we see that

[TABLE]

We again use magnetogeostrophic balance to obtain

[TABLE]

so

[TABLE]

Finally we must project away the component along $\boldsymbol{k}$, because this is eliminated by the pressure gradient in geostrophic balance. Because $k_{\perp}\approx k_{r}\omega_{\rm A}/N$ we find

[TABLE]

It follows that this acceleration produces the same mixing and hence the same cascade rate as the Coriolis effect.

The cascade rate we have found determines the non-linear damping of the magnetic energy, so that

[TABLE]

Or phrased in terms of the linear magnetic field,

[TABLE]

Substituting

[TABLE]

we find

[TABLE]

such that the effetive damping rate is

[TABLE]

For Tayler instability, the wavelength $\Lambda$ of the fastest growing modes is $\Lambda\approx 1/k_{\phi}\approx r/m\approx r$. In the main text, we show that $\delta v_{\rm A}/r\ll\omega_{\rm A}$, such that the effective damping rate is

[TABLE]

We note that a similar result can be obtained using the heuristic argument of Lithwick & Goldreich (2003). For Tayler instability, $\Lambda\approx\lambda\approx r$. Unlike isotropic magnetic turbulence, rotating Tayler instability is composed of magnetic perturbations that travel at group speed $v_{\rm g}\approx\omega_{\rm A}^{2}r/\Omega$. Then each scattering event occurs over time scale $t_{\rm scat}\approx r/v_{\rm g}\approx\Omega/\omega_{\rm A}^{2}$. Following the same argument used to derive equation 58, the cascade rate is then

[TABLE]

Using this result (instead of equation 86) in equation 13 yields an identical result. Additionally, we note that our saturated solution entails that

[TABLE]

where $\chi>1$ entails strong MHD turbulence and $\chi<1$ entails weak MHD turbulence, as defined by Chandran (2004). For strong MHD turblence, the cascade rate is $t_{\rm cas}\approx\delta v_{\rm A}/\lambda\approx\delta v_{\rm A}/r$, again equal to our result above.

Appendix B Energy Dissipation by Diffusively aided Field Migration

Magnetic energy in axisymmetric loops can be dissipated near the poles of the star by reconnection that connects magnetic field lines of opposing polarity. This only happens sufficiently close to the pole, where either diffusion can act across a magnetic field loop, or horizontal displacements compare to the size of the loop. Below we will show the latter length scale is larger and hence the relevant scale where dissipation occurs. The overall picture is that magnetic loops “migrate” both poleward and equatorward due to reconnection following a Tayler displacement. After one oscillation cycle, the maximum distance a loop can migrate in the latitudinal direction is $\xi_{\perp}$, the horizontal displacement of a loop caused by the Tayler instability. Its value is

[TABLE]

where we have used the fact that the Tayler instability growth rate and oscillation frequency (i.e., the imaginary and real components of the frequency) are both $\omega\sim\omega_{\rm A}^{2}/\Omega$. In what follows we assume $\omega_{\rm A}<\Omega<N$ as expected in stars.

Now, as shown in the text, the horizontal velocity is related to the perturbed magnetic field by $\delta v_{\perp}\sim(\omega_{\rm A}/\Omega)\delta v_{\rm A}$. Then

[TABLE]

We argue in the text that growth and damping of the instability are balanced when $\omega_{\rm A}^{2}/\Omega\sim\delta v_{\rm A}/r$. We show below that $\delta v_{\rm A}/r$ remains the relevant damping rate of the instability in spite of magnetic dissipation near the pole and any loop migration. Then we have

[TABLE]

These calculations are meant to be a mid latitudes where the cylindrical coordinate $R$ is comparable to the radial coordinate $r$. Some quantities will have different values very near the pole where $R\ll r$, but the migration time is dominated by mid latitudes where $R\sim r$, so magnetic energy can only be dissipated at the pole as fast as it migrates from mid latitudes.

Assuming loops of azimuthal field are totally dissipated near the pole, their effective damping rate is equal to their migration rate $\gamma_{\rm mig}$. Because the loop migration is essentially a random walk process, the migration timescale is

[TABLE]

where $N_{\rm step}$ is the migration length divided by a step length, and $t_{\rm step}$ is the time it takes to complete each step. The number of steps is $N_{\rm step}\sim r/\xi_{\perp}$. The time of each step is a magnetic diffusion time across a radial wavelength, $t_{\rm step}^{-1}\sim k_{r}^{2}\eta$. In order for the the instability to grow, this diffusion rate must be smaller than the growth rate $\omega_{\rm A}^{2}/\Omega$. So we have $t_{\rm step}\gtrsim\Omega/\omega_{\rm A}^{2}$. Then using equations 91 and 92, we have

[TABLE]

This maximum migration rate will be realized when $\omega_{\rm A}\sim\omega_{c}$, with the critical field strength $\omega_{c}$ defined by equation 1. When $\omega_{\rm A}\sim\omega_{c}$, field loops can reconnect with loops of opposite oscillation phase (i.e., those separated by radial distance $\sim 1/k_{r}$) after $\sim 1$ oscillation cycle, such that they can migrate by a distance $\sim\xi_{\perp}$ each oscillation cycle. When $\omega_{\rm A}>\omega_{c}$, reconnection requires many oscillation cycles, the migration rate will depend on the magnetic diffusivity, and it will be smaller than equation 93.

Why do we still think the instability damping rate is $\delta v_{\rm A}/r$? Let’s consider whether migration of magnetic loops toward the pole can destroy them at faster rates. Assuming a loop reconnects with its neighbor after being displaced horizontally by $\xi_{\perp}$ at a rate $\omega_{\rm A}^{2}/\Omega$, the two loops have moved apart from one another at a speed $v\sim\omega_{\rm A}^{2}\xi_{\perp}/\Omega\sim\delta v_{\perp}$. Then the maximum rate at which the instability can be damped due to loops migrating to the pole is

[TABLE]

But since $\omega_{\rm A}<\Omega$, $\gamma_{\rm diss}$ is smaller than the damping rate $\delta v_{A}/r$. The actual destruction rate is likely much slower due to the random walk process discussed above, and is given by equation 93. This means that unstable perturbations will damp faster by weak turbulence than they will by migrating toward the pole. So the instability is still limited by weak turbulence, and setting the growth rate equal to the damping rate still implies

[TABLE]

The migration rate of equation 93 implies that the background field, whose energy density is $\sim B_{\phi}^{2}$, is destroyed at this rate. The assosiated energy damping rate per unit volume is

[TABLE]

The turbulent energy damping discussed in the paper will operate regardless of the loop migration, and will be more important when the reconnection timescale is longer than an oscillation times scale. Hence, while damping from loop migration may be relevant when $\omega_{\rm A}\sim\omega_{c}$, we do not expect any of our scaling arguments or results to be altered.

Appendix C Effective Stratification

As discussed in Spruit (2002), the effective Brunt-Väisälä frequency $N_{\rm eff}$ depends on the thermal diffusion timescale across the Tayler instability lengthscale. This thermal diffusion timescale in turn depends on $N_{\rm eff}$. Spruit (2002) considered the limit of pure thermal/compositional stratification, but in red giant cores, both components are important. Here we derive an improved method for incorporating thermal diffusion in the general case.

Following the suggestion by Spruit (2002), thermal diffusion reduces the thermal component of the effective stratification $N_{\rm T}$ by roughly

[TABLE]

where $k$ is the instability wavenumber, $\chi$ is the thermal diffusivity, and $\omega$ is the instability timescale. The fastest growing modes have $rk\sim N_{\rm eff}/\omega_{\rm A}$ and $\omega=\omega_{\rm A}^{2}/\Omega$. Using the saturated value of $\omega_{\rm A}$ from equation 34, we have

[TABLE]

The effective stratification is $N_{\rm eff}^{2}=N_{\rm T,eff}^{2}+N_{\mu}^{2}$, where $N_{\mu}$ is the compositional component of the stratification. Substituting for $N_{\rm T,eff}$, we find

[TABLE]

Equation 99 can be solved for the appropriate value of $N_{\rm eff}$ given a stellar structure. Inspection reveals that it reduces in the appropriate limits. When $\chi\rightarrow 0$, we recover $N_{\rm eff}=N$. When $\chi\rightarrow\infty$, we recover $N_{\rm eff}=N_{\mu}$. And when $N_{\mu}\rightarrow 0$ and thermal diffusion is large, $N_{\rm eff}\ll N$ and we find $N_{\rm eff}=(r^{6}q^{4}\alpha^{12}\Omega^{13}N^{6}/\chi^{3})^{1/16}$. We have not yet implemented numerical solutions of equation 99 into our MESA routines, but we plan to do this in future work.

Appendix D MESA Model Inlists

We use the MESA stellar evolution code Paxton et al. (2011); Paxton et al. (2013, 2015, 2018) version 10108 to generate our stellar models. The inlist for our models is as follows:

&star_job

  pgstar_flag = .true.

  new_rotation_flag = .true.
  change_rotation_flag = .true.
  change_initial_rotation_flag = .true.

  new_omega = 3.64e-5
  set_initial_omega = .true.

/ ! end of star_job namelist

&controls

  !------------------------  MAIN

  initial_mass = 1.6
  initial_z = 0.02
  use_Type2_opacities = .true.
  Zbase = 2.d-2

  set_min_D_mix = .true.
  min_D_mix = 1d1

  mesh_delta_coeff = 0.7
  varcontrol_target = 0.7d-3

  predictive_mix(1) = .true.
  predictive_superad_thresh(1) = 0.005
  predictive_avoid_reversal(1) = ’he4’
  predictive_zone_type(1) = ’any’
  predictive_zone_loc(1) = ’core’
  predictive_bdy_loc(1) = ’top’

  dX_div_X_limit_min_X = 1d-4
  dX_div_X_limit = 5d-1
  dX_nuc_drop_min_X_limit = 1d-4
  dX_nuc_drop_limit = 1d-2

  !---------------------  Rotation

  am_nu_ST_factor = 0
  use_other_am_mixing = .true.

  am_time_average = .true.
  premix_omega = .true.
  recalc_mixing_info_each_substep = .true.
  am_nu_factor = 1
  am_nu_non_rotation_factor = 1d0
  am_nu_visc_factor = 0.333
  angsml = 0.0

  !-------------------------  WIND

  cool_wind_RGB_scheme = ’Reimers’
  cool_wind_AGB_scheme = ’Blocker’
  RGB_to_AGB_wind_switch = 1d-4
  Reimers_scaling_factor = 0.2
  Blocker_scaling_factor = 0.5
  use_accreted_material_j = .true.
  accreted_material_j = 0

  !-------------------  OVERSHOOTING

  overshoot_f_above_nonburn_core = 0.015
  overshoot_f0_above_nonburn_core = 0.005
  overshoot_f_above_nonburn_shell = 0.015
  overshoot_f0_above_nonburn_shell = 0.005
  overshoot_f_below_nonburn_shell = 0.015
  overshoot_f0_below_nonburn_shell = 0.005

  overshoot_f_above_burn_h_core = 0.015
  overshoot_f0_above_burn_h_core = 0.005
  overshoot_f_above_burn_h_shell = 0.015
  overshoot_f0_above_burn_h_shell = 0.005
  overshoot_f_below_burn_h_shell = 0.015
  overshoot_f0_below_burn_h_shell = 0.005

  overshoot_f_above_burn_he_core = 0.015
  overshoot_f0_above_burn_he_core = 0.005
  overshoot_f_above_burn_he_shell = 0.015
  overshoot_f0_above_burn_he_shell = 0.005
  overshoot_f_below_burn_he_shell = 0.015
  overshoot_f0_below_burn_he_shell = 0.005

/ ! end of controls namelist

Some important controls include the use of predictive mixing to help mitigate “breathing pulses” in the size of the convective helium-burning core during the clump. Additionally, the use of am_time_average, premix_omega, and recalc_mixing_info_each_substep help reduce numerical artifacts related to AM transport. Smoothing the shear and AM diffusivity (see next section) also help reduce these numerical instabilities. The artifacts arise because large MESA timesteps can cause AM transport to artificially create step-like features in the stellar rotation profile. The steps arise where AM transport in some grid cells is slightly more efficient than neighboring grid cells due to the discrete grid size and inaccurate numerical derivatives. A large time step will cause the rotation profile to flatten in grid cells with larger AM diffusivity, and steepen in neighboring grid cells with smaller AM diffusivity. The controls above help mitigate these effects, but in some cases enforcing smaller time steps may be useful.

The initial masses and rotation rates are adjusted as described in the text. In some models we adjust the wind scaling factors on the AGB in order to avoid late helium flashes. We also enable MLT++ to evolve more massive stars from the AGB to the WD cooling track:

okay_to_reduce_gradT_excess = .true.
gradT_excess_max_change = 1d-2

and in some cases we remove the last few hundredths of a solar mass of the hydrogen envelope using

remove_H_wind_mdot = 1d-4
remove_H_wind_H_mass_limit = 1d-5

We have performed some basic resolution testing to verify our results are very insensitive to the model’s grid resolution and time stepping.

D.1 Implementation of Angular Momentum Transport

Our run_star_extras.f code for implementation of AM transport in our MESA models is as follows:

subroutine TSF(id, ierr)

integer, intent(in) :: id integer, intent(out) :: ierr type (star_info), pointer :: s integer :: k,j,op_err,nsmooth,nsmootham real(dp) :: alpha,shearsmooth,nu_tsf,nu_tsf_t,omegac,omegag,omegaa,omegat real(dp) :: difft,diffm,brunts,bruntsn2,logamnuomega,alphaq

call star_ptr(id,s,ierr) if (ierr /= 0) return

alpha=1d0 nsmooth=5 nsmootham=nsmooth-3 shearsmooth=1d-30 op_err = 0

!Calculate shear at each zone, then calculate TSF torque do k=nsmooth+1,s% nz-(nsmooth+1)

nu_tsf=1d-30
nu_tsf_t=1d-30
!Calculate smoothed shear, q= dlnOmega/dlnr
shearsmooth = s% omega_shear(k)/(2.*nsmooth+1.)
do j=1,nsmooth
  shearsmooth = shearsmooth + (1./(2.*nsmooth+1.))*( s% omega_shear(k-j) + s% omega_shear(k+j) )
end do

!Magnetic diffusivity
diffm =  diffmag(s% rho(k),s% T(k),s% abar(k),s% zbar(k),op_err)
!Thermal diffusivity
difft = 16d0*5.67d-5*(s% T(k))**3/(3d0*s% opacity(k)*(s% rho(k))**2*s% Cv(k))
!Alfven frequency at saturation
omegaa = s% omega(k)*(shearsmooth*s% omega(k)/sqrt(abs(s% brunt_N2(k))))**(1./3.)
!Thermal damping rate assuming adiabatic instability
omegat = difft*pow2(sqrt(abs(s% brunt_N2(k)))/(omegaa*s% r(k)))
!Suppress thermal part of brunt
brunts = sqrt(abs( s% brunt_N2_composition_term(k) +
   (s% brunt_N2(k)-s% brunt_N2_composition_term(k))/(1d0 + omegat/omegaa) ))
!Effective brunt for isothermal instability
bruntsn2 = sqrt(abs( s% brunt_N2_composition_term(k) +
   (s% brunt_N2(k)-s% brunt_N2_composition_term(k))*min(1d0,diffm/difft) ))
!Choose max between suppressed brunt and isothermal brunt
brunts = max(brunts,bruntsn2)
!Don’t let Brunt be smaller than omega
brunts = max(s% omega(k),brunts)
!Recalculate omegaa
omegaa = s% omega(k)*abs(shearsmooth*s% omega(k)/brunts)**(1./3.)

!Calculate nu_TSF
if (s% brunt_N2(k) > 0.) then
  if (pow2(brunts) > 2.*pow2(shearsmooth)*pow2(s% omega(k))) then
    !Critical field strength
    omegac = 1d0*s% omega(k)*((brunts/s% omega(k))**0.5)*(diffm/(pow2(s% r(k))*s% omega(k)))**0.25
    !Suppress AM transport if omega_a<omega_c
    nu_tsf = 5d-1+5d-1*tanh(5d0*log(alpha*omegaa/omegac))
    !nu_omega for revised Tayler instability
    nu_tsf = nu_tsf*alpha**3*s% omega(k)*pow2(s% r(k))*(s% omega(k)/brunts)**2
  end if
  ! Add TSF enabled by thermal diffusion
  if (pow2(brunts) < 2.*pow2(shearsmooth)*pow2(s% omega(k))) then
    nu_tsf_t = alpha*abs(shearsmooth)*s% omega(k)*pow2(s% r(k))
  end if
  s% am_nu_omega(k) = s% am_nu_omega(k) + max(nu_tsf,nu_tsf_t) + 1d-1
end if

end do

Ψ !Smooth nu_omega logamnuomega=-3d1 do k=nsmootham+1,s% nz-(nsmootham+1)

!Don’t smooth convective diffusivity into non-convective zones
if (s% mixing_type(k)==1) then
   s% am_nu_omega(k) = s% am_nu_omega(k)
   !Smooth zones if not including a convective zone
else
  logamnuomega = log10(s% am_nu_omega(k))/(2.*nsmootham+1.)
end if

do j=1,nsmootham
  !Don’t smooth convective diffusivity into non-convective zones
  if (s% mixing_type(k-j)<3.5) then
    logamnuomega = log10(s% am_nu_omega(k))
    !Smooth zones if not including a convective zone
  else
    logamnuomega = logamnuomega + (1./(2.*nsmootham+1.))*log10(s% am_nu_omega(k-j))
  end if
end do

do j=1,nsmootham
  !Don’t smooth convective diffusivity into non-convective zones
  if (s% mixing_type(k+j)<3.5) then
    logamnuomega = logamnuomega
    !Smooth zones if not including a convective zone
  else
    logamnuomega = logamnuomega + (1./(2.*nsmootham+1.))*log10(s% am_nu_omega(k+j))
  end if
end do

s% am_nu_omega(k) = 10.**logamnuomega

end do

!Values near inner boundary do k=s% nz-nsmootham,s% nz s% am_nu_omega(k) = s% am_nu_omega(k-1) end do

!Values near outer boundary do k=nsmootham,1 s% am_nu_omega(k) = s% am_nu_omega(k-1) end do

end subroutine TSF

These controls work well for our models, but we caution that they may not work well in different situations. For instance, we estimate $N_{\rm eff}$ in a way which is accurate for our models but may be problematic in some stars. Magnetic diffusivity is calculated via the modules included in MESA’s default implementation for TS torques. To disable our AM transport prescription when $\omega_{\rm A}\!<\!\omega_{c}$, we use a tanh function to smoothly transition from no torque at $\omega_{\rm A}\!<\!\omega_{c}$ to full torque at $\omega_{\rm A}\!>\!\omega_{c}$. Additionally, we smooth the dimensionless shear by 5 grid cells on each side, and we smooth the AM diffusivity by 32 grid cells on each side. In our models, this level of smoothing helps suppress numerical instabilities but does not strongly affect the evolution because larger smoothing lengths deliver nearly identical results.