{'fraction_non_alphanumeric': 0.029924477117068842, 'fraction_numerical': 0.01594185910209823, 'mean_word_length': 4.607850502394591, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Slowing the Spins of Stellar Cores

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.


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 3 M ) 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 Email: jfuller@caltech.edu 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 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 . 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.

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.

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 φ with corresponding Alfvèn frequency ωA = B φ / 4πρr 2 , (2) the angular rotation frequency Ω 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 ωA Ω N , where N is the Brunt-Väisälä frequency. As shear winds the magnetic field, B φ grows and becomes Tayler unstable when it reaches a critical strength of (Spruit 2002;Zahn et al. 2007) ωA > ωc ∼ Ω N Ω 1/2 where η is the magnetic diffusivity. The corresponding growth rate of this instability is largest for m = 1 perturbations and is approximately Due to the strong stratification in these stars, the radial length scale of the instability is limited to while the maximum horizontal length scale of the instability is l ⊥ ∼ r. 1 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 where K is the thermal diffusivity, N 2 T is the thermal component of the stratification, and N 2 µ is the compositional component. Red giant cores have large composition gradients, so N eff in much of the core (and especially at the hydrogenburning shell, the bottleneck for AM transport) is dominated by its compositional component, and thus N eff Nµ. Appendix C discusses the appropriate value of N eff when thermal diffusion is moderately important.

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 γ turb with the linear growth rate results in where the wavenumber kr is the minimum required for instability, kr ∼ ωA/(N r), and η eff is an effective turbulent diffusivity. Next, since the azimuthal field grows via winding by shear as where Br is the radial field, then the amplification rate is γamp = qΩBr/B φ . The incompressible nature of the instability implies that krBr ∼ k ⊥ B φ and thus Br ∼ (ωA/N )B φ ,   Spruit 2002 (top), and our model (bottom). Black arrows represent magnetic field lines, while red arrows represent fluid motions. Orange text describes processes that we argue operate differently than proposed by Spruit 2002. such that If the azimuthal field B φ is turbulently damped at the rate of equation 5, setting equations (5) and (7) equal determines the azimuthal field strength at saturation B φ / 4πρr 2 ∼ ωA ∼ qΩ 2 /N . The radial field strength is then Br/ 4πρr 2 ∼ q 2 Ω 4 /N 3 , so that 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 φ . Thus the axisymmetric component of Br is not necessarily related to the axisymmetric component of B φ via Br/B φ ∼ ωA/N . A second potential issue is that equation (5) may not predict the correct damping rate for a large-scale background field B φ that varies on lengthscales much larger than 1/kr. If the background field B φ 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 φ , and the saturated values of Br and B φ 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.

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 δ B and δ v. Tayler instability transfers energy from large scale magnetic fields to perturbed fields δ B that vary on the short length scale ∼ 1/kr. In the linear regime, δ B and δ v are related to each to each other via where ξ is the Lagrangian displacement associated with the instability. Using δ v −iωR ξ, where ωR ∼ ω 2 A /Ω is the real part of the perturbation frequency, and k Here we have used m = 1 and ignore geometric terms of order unity, and δ vA = δ B/ √ 4πρ is the perturbed Alfvén velocity. A similar answer can be obtained by analyzing the momentum equation 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 where δv ⊥ and δv A,⊥ are the horizontal components of the perturbed velocity and Alfvén velocity. Since we shall find ωA Ω 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 γcas ∼ δvA r .
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 γcas effectively represents a turbulent damping rate of the perturbed magnetic energy. The non-linear energy dissipation rate is thenĖ We do not expect energy in the background field B φ to be damped by a turbulent cascade to small scales. This is a key difference from Spruit (2002), who uses a damping rateĖ damp ∼ γ turb |B 2 φ | ∼ (ω 2 A /Ω)|B 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 Note that equation 15 implies that upon saturation, the perturbed and background field are related by so that the energy damping rate iṡ 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Ė damp,pole 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.
(i) Rotational shear energy is converted to magnetic energy by winding a radial field into a toroidal field.
(ii) Toroidal field energy is converted by the Tayler instability into magnetic/kinetic energy associated with the perturbed magnetic/velocity field.
(iii) 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 Equating the energy input rate to energy damping rate, we have To solve our system, we need an estimate of Br/B φ .
In Section 2.3.1, we argue that Br can grow until Lorentz forces nearly stabilize the plasma against the growth of the Tayler instability, such that 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 Br These fields can then drive AM transport and chemical mixing as further described in Section 2.4.

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, Br 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 Br can then continue to be wound by the shear, and continue the flow of energy as outlined above. The question is how this new Br is generated and how strong can it grow. As argued in Section 2.2, closing the loop between Br and B φ is difficult if only linear effects are considered. This can be seen by starting with the linearized induction equation, and taking the azimuthal average, which yields 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 To order of magnitude, we thus expect Hence, we expect some growth of an axisymmetric radial field due to non-linear induction, i.e., an α-dynamo effect. However, as shown by Braithwaite (2009), Tayler instability cannot operate if Br 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 If Br grows until the instability is quenched, then the maximum possible value of Br (corresponding to the longest length scale unstable disturbance rkr ∼ N eff /ωA) is This is identical to the condition arising from incompressibility used by Spruit (2002), but it relates the axisymmetric component of Br and B φ , whereas Spruit's relation is only valid for the non-axisymmetric component of Br.

Angular Momentum Transport
The torque via Maxwell stresses in the saturated state is found from combining equations (22) and (23), corresponding to an effective AM diffusivity 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 α such that The parameterized AM diffusivity is then We expect α of order unity, and indeed in Section 3 we find α ≈ 1 fits the observational data.
Combining the instability criterion given by equation (1) with the value of ω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 We show in Section 3 that this minimum shear appears to frequently be realized in red giant stars, such that ωA ∼ ω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 eff and η) and are very insensitive to the initial rotation rate or prior evolution of the star.

Energetics and Mixing
Note that our relations at saturation imply a hierarchy of rotational, background magnetic, perturbed magnetic, and kinetic energy densities: These hierarchies are true as long as q N eff /Ω, which is true in our models in Section 3 where q ∼ 1 and N eff /Ω ∼ 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: The radial components of the fields are typically orders of magnitude smaller than the horizontal components due to the tiny value of Ω/N 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 Using the relations above, we have so that for the ratio of chemical mixing to AM transport. In red giants, we find νmix/νAM ∼ 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 , 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 1 or N eff N .

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.  Figure 2. Important frequencies as a function of radius in a 1.2 M , 4 R model at the base of the RGB. We show the Brunt-Väisälä frequency N , its compositional component Nµ, and its effective value N eff (equation 4) when thermal diffusion is important. We also plot the angular rotation frequency Ω, the saturated Alfvén frequency ω A (equation 34), the minimum Alfvén frequency required for Tayler instability ωc (equation 1), and the thermal diffusion frequency at the instability length scale ωt (equation 43). Note that ω A Ω throughout the interior such that Tayler instability occurs in the rapidly rotating regime. Because ωt ω A , the instability occurs in the limit where thermal diffusion is important.

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 model on the lower RGB at R = 4.1 R and log(g)= 3.3. At the hydrogen burning shell, the stabilization is primarily due to the hydrogen-helium composition gradient such that N Nµ, but even above the burning shell, we often find Nµ ∼ 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, Here, we have used the maximum radial lengthscale for Tayler instability lr ∼ r(ωA/N eff ), and the thermal diffusivity χ = 16σSBT 3 /(3ρ 2 cvκ). Thermal diffusion strongly reduces the effective thermal stratification when γ ω, where ω is the real part of the frequency of the overstable oscillations. Zahn et al. (2007) show that ω ∼ ω 2 A /Ω. Using our saturated field strength (equation 34), we find thermal dif- White dwarfs P surf P core , P i = 2 d, α = 1 P core , P i = 2 d, TS Only Figure 3. Post-main sequence rotational evolution of a 1.6 M star, with a ZAMS rotation rate P i = 2 days and including our updated prescription for AM transport with α = 1. We plot the surface rotation rate (blue line), and core rotation rate as sensed by mixed modes (red line). We also include the core rotation rate using a prior prescription for the TS dynamo (black line). Our AM transport scheme closely matches observations along the red giant branch (Mosser et al. 2012;Gehan et al. 2018), red clump (Mosser et al. 2012;Deheuvels et al. 2015), and white dwarf phase (Hermes et al. 2017).
fusion is very important when A comparison of ωt and Ω 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 2 eff N 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 α. To that end, we construct rotating stellar models using the MESA stellar evolution code (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 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 Pcore for comparison with asteroseismology, we compute the average core rotation period Pcore = 2π/Ωcore as sensed by a gravity wave in the WKB limit, The bounds of the integral in equation 46 correspond to the boundaries of the core gravity mode cavity where ωg < N , and we consider gravity waves with frequency ωg = 2πνmax, where the frequency of maximum power is calculated via classical scaling relations νmax = 3090 µHz (M/M )(R/R ) −2 (T /T ) −1/2 . For our WD models, we simply set Pcore equal to the central rotation rate. We evolve solar metallicity models ranging from 1.2 − 6 M from the zero-age main sequence (ZAMS) to the WD phase. We initiated each model with a spin rate Pi = 2 days, except for the 1.2 M model for which we used Pi = 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 model with a ZAMS rotation period Pi = 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 M 3 M . Typical core rotation rates are Pcore ∼ 10 − 20 days on the lower RGB, Pcore ∼ 50 − 200 days on the red clump, and P ∼ 0.5 − 4 days for WDs.
Our models generally exhibit very similar rotation rates to observations for α ≈ 1, a reasonable value since we expect α ∼ 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 Ω/N 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 2 M ), 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 α and the star's initial spin rate, Figure 4 shows core rotation rates for models with different values of α and Pi. 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 νAM ∝ Ω 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-  Figure 4. Core rotation rate for the same star as Figure 3, but varying the initial rotation rate and parameter α. The post-main sequence rotation period scales approximately as α, but is relatively insensitive to initial rotation rate.
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 α. Nearly rigid rotation is maintained along the main sequence, regardless of α, except for very slow rotators in which the value of νAM is much smaller (see discussion in Section 4). The models slightly diverge from one another on the lower RGB, with smaller values of α allowing faster core rotation. We find the main effect of the value of α in our models is not the prefactor in equation 35, but rather in determining the minimum shear qmin (equation 36) required for Tayler instability to saturate as we have outlined. When q > qmin, efficient AM transport generally decreases the core rotation and shear, thereby reducing q until q ∼ qmin, as shown in Appendix D. In this limit, equation 36 predicts that the core rotation scales as Ω ∝ α −12/13 , in line with our numerical finding that the core spin period is approximately proportional to α. Figure 5 shows the rotation profile of a 1.2 M model at the base of the RGB for α = 1. The shear is strongest at the hydrogen burning shell where Nµ 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 days Pcore 30 days regardless of mass. On the clump, the differences between stars of different masses are slightly  Figure 6. Core rotational evolution (as sensed by mixed modes) for models of several masses. Each model uses α = 1 and is initiated with P i = 2 days, except the 1.2 M model which has P i = 20 days. The highest plateaus of each model correspond to the core helium-burning phase, and each model ends as a carbonoxygen white dwarf. larger: our models predict slower core rotation rates for lowmass clump stars, with Pcore ∼ 100 − 200 days for ≈ 1.2 M stars. We predict faster rotation for secondary clump stars, with Pcore ∼ 50 days for M ≈ 2.2 M . 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 ∼ 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 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 ∼ 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.

DISCUSSION
Our AM transport prescription predicts extremely short AM transport times tAM within radiative zones of main sequence stars, For a fast rotating young Sun, we find that the instability occurs in the nearly adiabatic limit, such that N eff ≈ N . Evaluating equation (47) in the radiative zone of a young solar model rotating at P = 3 days yields a typical AM transport timescale tAM ∼ 10 years. In the current Sun, the rotation rate is much slower and tAM ∼ 10 4 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 tAM ∼ 10 6 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 ≈ 1.0 M 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 tAM ∼ 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 tAM ∼ 10 7 − 10 8 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. 2016Rebull et al. , 2017Rebull 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. 2015Garraffo et al. , 2016Garraffo et al. , 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 tAM 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 Ω/N eff 1, which is often the case in radiative regions that have undergone any nuclear processing. An exception to this is horizontal mixing. Unless q 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, νAM ∝ Ω 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 1.2M ) 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 ∼2) along the lower RGB, whereas Gehan et al. (2018) find no clear spinup/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 biases 2 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 Br, which we predict to have a strength of Br ∼ 10 −2 G through much of the radiative core. Recall that for the predicted value of B φ , Tayler instability cannot occur if there is a fossil field with strength greater than equation 23. If there is a fossil component of Br larger than that shown in Figure 8, the azimuthal component B φ must be amplified by shear to larger field strengths before Tayler instability kicks in. By the time this occurs, the Maxwell stress BrB φ 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 (Br 10 −2 G) may enforce nearly rigid rotation of the radiative core. While the internal field strengths of red giants are not well known 3 , a rigidly rotating core enforced by fossil fields would necessitate large differential rotation in the 2 We 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). 3 While some red giants with suppressed dipole oscillation modes may have very strong (Br 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.

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 α ≈ 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 PNS ∼ 10 ms, somewhat faster than typical pulsar birth periods PNS 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.

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.
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: 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 where vA,0 is the slowly-varying background copmponent and δvA is the fluctuating one. With this decomposition we define combinations of fluid velocity and magnetic field fluctuation The equations of motion may then be cast as 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 (A7) just describes the advection of combined variations in the fluid velocity and the magnetic field. That is, ∂tw± ± vA,0 · ∇w± = 0. (A9) These combined fluctuations evidently only propagate along the large-scale magnetic field. With the addition of the nonlinear terms, we see that packets of w+ may scatter off of those of w− and vice-versa, but because the linear evolution is constrained to be along vA,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±/vA,0. By repeatedly bending the field-lines of a wave-packet of w± the packet may be disrupted, such that its energy cascades non-linearly to different scales. If the packets of w± have wavelength Λ parallel to vA,0 and wavelength λ transverse to it then this disruption occurs if the field lines are displaced transversely by an amount of order λ, or by an angle of order λ/Λ. 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±/vA,0, so that it takes of order 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 1. Because each scattering event takes time of order Λ/vA,0, the timescale over which wave-packets are disrupted is 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 tcas. It follows that the damping rate of the energy in the system is of orderr or equivalently, the non-linear loss rate iṡ 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 (A1), so that The new acceleration term is not fully absorbed into the pressure gradient because it is not fully directed along the wave-vector k, and hence, at least to order of magnitude, it must be kept. Upon making the same change to equation (A7) we find ∂tw± ± vA,0 · ∇w± + Ω × (w+ + w−) Crucially, this is no longer of the form of a simple advection equation like (A9). 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±.
Assuming that w is not oriented nearly-parallel to the rotation axis, the time-scale over which wave-packets transition is Recalling that we work in the limit of magnetogeostrophic balance, While in the mixed state, the cascade proceeds with timescale Analogous to the strong and weak regimes of Lithwick & Goldreich (2003), we therefore find two regimes. In the first, ωAt * cas 1, so that scattering occurs rapidly once the wavepackets mix. This results in an effective cascade time In the opposing limit, ωAt * cas 1, scattering is slow and it makes sense to average w∓ · ∇w± over the mixing time ω −1 A . 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 Putting the two regimes together we find the overall effect of rotation is to reduce the cascade time-scale to The above argument may also be cast in terms of new linear combinations of v and δvA 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 w+ · ∇w+, and likewise for w−. These new selfinteraction terms result, following the arguments above, in the same cascade time-cale tcas ≈ Λ/w±.
The second modification we must consider is that of buoyancy. This works in much the same way. We define as the Lagrangian displacement of a fluid element. With this, and working in the Boussinesq limit, equation (A14) becomes ∂tv + v · ∇v = −(1/ρ)∇P −rr · xN 2 + vA · ∇vA, (A24) 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 (A7) we obtain We may approximate the displacement as where ω is the linear frequency associated with any given mode. Hencê It follows that this term, like the Coriolis one, produces mixing between w±. Because the real and imaginary parts of ω are both of order ω 2 A /Ω, this may be written aŝ we see thatrr We again use magnetogeostrophic balance to obtain Finally we must project away the component along k, because this is eliminated by the pressure gradient in geostrophic balance. Because k ⊥ ≈ krωA/N we find 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 nonlinear damping of the magnetic energy, so that Or phrased in terms of the linear magnetic field, such that the effetive damping rate is For Tayler instability, the wavelength Λ of the fastest growing modes is Λ ≈ 1/k φ ≈ r/m ≈ r. In the main text, we show that δvA/r ωA, such that the effective damping rate is We note that a similar result can be obtained using the heuristic argument of Lithwick & Goldreich (2003). For Tayler instability, Λ ≈ λ ≈ r. Unlike isotropic magnetic turbulence, rotating Tayler instability is composed of magnetic perturbations that travel at group speed vg ≈ ω 2 A r/Ω. Then each scattering event occurs over time scale tscat ≈ r/vg ≈ Ω/ω 2 A . Following the same argument used to derive equation A11, the cascade rate is then Using this result (instead of equation A39) in equation 13 yields an identical result. Additionally, we note that our saturated solution entails that where χ > 1 entails strong MHD turbulence and χ < 1 entails weak MHD turbulence, as defined by Chandran (2004). For strong MHD turblence, the cascade rate is tcas ≈ δvA/λ ≈ δvA/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 ξ ⊥ , the horizontal displacement of a loop caused by the Tayler instability. Its value is 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 ω ∼ ω 2 A /Ω. In what follows we assume ωA < Ω < N as expected in stars.
Now, as shown in the text, the horizontal velocity is related to the perturbed magnetic field by δv ⊥ ∼ (ωA/Ω)δvA. Then We argue in the text that growth and damping of the instability are balanced when ω 2 A /Ω ∼ δvA/r. We show below that δvA/r remains the relevant damping rate of the instability in spite of magnetic dissipation near the pole and any loop migration. Then we have 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 r, but the migration time is dominated by mid latitudes where R ∼ 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 γmig. Because the loop migration is essentially a random walk process, the migration timescale is where Nstep is the migration length divided by a step length, and tstep is the time it takes to complete each step. The number of steps is Nstep ∼ r/ξ ⊥ . The time of each step is a magnetic diffusion time across a radial wavelength, t −1 step ∼ k 2 r η. In order for the the instability to grow, this diffusion rate must be smaller than the growth rate ω 2 A /Ω. So we have tstep Ω/ω 2 A . Then using equations B3 and B4, we have This maximum migration rate will be realized when ωA ∼ ωc, with the critical field strength ωc defined by equation 1. When ωA ∼ ωc, field loops can reconnect with loops of opposite oscillation phase (i.e., those separated by radial distance ∼ 1/kr) after ∼ 1 oscillation cycle, such that they can migrate by a distance ∼ ξ ⊥ each oscillation cycle. When ωA > ωc, reconnection requires many oscillation cycles, the migration rate will depend on the magnetic diffusivity, and it will be smaller than equation B5. Why do we still think the instability damping rate is δvA/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 ξ ⊥ at a rate ω 2 A /Ω, the two loops have moved apart from one another at a speed v ∼ ω 2 A ξ ⊥ /Ω ∼ δv ⊥ . Then the maximum rate at which the instability can be damped due to loops migrating to the pole is But since ωA < Ω, γ diss is smaller than the damping rate δvA/r. The actual destruction rate is likely much slower due to the random walk process discussed above, and is given by equation B5. 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 The migration rate of equation B5 implies that the background field, whose energy density is ∼ B 2 φ , is destroyed at this rate. The assosiated energy damping rate per unit volume isĖ 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 ωA ∼ ω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 eff depends on the thermal diffusion timescale across the Tayler instability lengthscale. This thermal diffusion timescale in turn depends on N 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 NT by roughly where k is the instability wavenumber, χ is the thermal diffusivity, and ω is the instability timescale. The fastest growing modes have rk ∼ N eff /ωA and ω = ω 2 A /Ω. Using the saturated value of ωA from equation 34, we have The effective stratification is N 2 eff = N 2 T,eff + N 2 µ , where Nµ is the compositional component of the stratification. Substituting for N T,eff , we find Equation C3 can be solved for the appropriate value of N eff given a stellar structure. Inspection reveals that it reduces in the appropriate limits. When χ → 0, we recover N eff = N . When χ → ∞, we recover N eff = Nµ. And when Nµ → 0 and thermal diffusion is large, N eff N and we find N eff = (r 6 q 4 α 12 Ω 13 N 6 /χ 3 ) 1/16 . We have not yet implemented numerical solutions of equation C3 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. (2011Paxton et al. ( , 2013Paxton et al. ( , 2015Paxton et al. ( , 2018 version 10108 to generate our stellar models. The inlist for our models is as follows: 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 steplike 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.