Rotation and Magnetism of Massive Stellar Cores

TL;DR

This paper investigates the internal rotation and magnetic field evolution of massive stars, linking core dynamics to neutron star and black hole properties, with implications for magnetic field strength and binary interactions.

Contribution

It presents new models of stellar core rotation and magnetism, highlighting the impact of convective layers and binary interactions on remnant spin and magnetic fields.

Findings

Neutron star spin periods vary from 0.1 to 2 ms depending on stellar mass.

Stronger inward angular momentum pumping occurs in more massive stars.

Binary interactions can significantly spin up black hole remnants.

Abstract

The internal rotation and magnetism of evolved massive stars are considered in response to i) the inward pumping of angular momentum through deep and slowly rotating convective layers; and ii) the winding up of a helical magnetic field in radiative layers. Field winding can transport angular momentum effectively even when the toroidal field is limited by kinking. Magnetic helicity is pumped into a growing radiative layer from an adjacent convective envelope (or core). The receding convective envelope that forms during the early accretion phase of a massive star is the dominant source of helicity in its core, yielding a $\sim 10^{13}$ G polar magnetic field in a collapsed neutron star (NS) remnant. Using MESA models of various masses, we find that the NS rotation varies significantly, from $P_{\rm NS} \sim 0.1-1$ s in a 13$\,M_\odot$ model to $P_{\rm NS} \sim 2$ ms in a $25\,M_\odot$ model with an extended core. Stronger inward pumping of angular momentum is found in more massive stars, due to the growing thickness of the convective shells that form during the later stages of thermonuclear burning. On the other hand, stars that lose enough mass to form blue supergiants in isolation end up as very slow rotators. The tidal spin-up of a 40$\,M_\odot$ star by a massive binary companion is found to dramatically increase the spin of the remnant black hole, allowing a rotationally supported torus to form during the collapse. The implications for post-collapse decay or amplification of the magnetic field are also considered.

Figures

Equations

$$\xi_{M}\equiv\frac{M/M_{\odot}}{r(M)/1000~{}{\rm km}}\bigg{|}_{t=t_{\rm bounce}}.$$

$$\dot{M}=\frac{2(M_{\mathrm{ZAMS}}-3M_{\odot })}{T^2}\cdot t{M}_{\odot }{\mathrm{yr}}^{-1}.$$

$${\tau }_{\mathrm{ev}}=\min \left[\frac{l_P}{|v_r|},\frac{t_{\mathrm{cc}}-t}{3}\right].$$

$$\Omega(r)\propto\left\{\begin{array}[]{lr}r^{-2}&\quad[{\rm Co}(r)\leq 1]\\ r^{-(1+\beta)/2}&\quad[{\rm Co}(r)>1].\end{array}\right.$$

$$v_{{\rm A},r}\biggr{|}_{\rm min}={\rm max}\left[{l_{P}\over\tau_{\rm ev}}\left({r\over 2l_{P}}\Omega\tau_{\rm ev}\right)^{1/4},\;{(rl_{P})^{1/2}\over\tau_{\rm ev}}\right].$$

$${\tau }_{\mathrm{shear}}=\frac{r^2}{D_{\mathrm{shear}}}=\frac{45}{2}{N}^2{\kappa }_T^{-1}{\left(\frac{\partial \Omega }{\partial r}\right)}^{-2}.$$

$${\tau }_{\mathrm{MC}}\sim {\tau }_{\mathrm{KH}}{\left(\frac{{\Omega }_{\mathrm{K}}}{\Omega }\right)}^2,$$

$$\delta \mathcal{H}\sim \pi B_r{B}_{\varphi }{R}_{\mathrm{b}}^3{v}_{\mathrm{b}}\delta t,$$

$$\frac{B_r{B}_{\varphi }}{4\pi }\sim \frac{\rho l_P{R}_{\mathrm{b}}(\Delta \Omega {)}^2}{2\pi N_{\mathrm{dyn}}}={\epsilon }_B\rho {\Omega }^2.$$

$$\delta {\mathcal{H}}_{\mathrm{dyn}}\sim \frac{{\epsilon }_B}{N_{\mathrm{dyn}}}\cdot \pi \delta M_{\mathrm{dyn}}{\Omega }^2(R_{\mathrm{b}})R_{\mathrm{b}}^3.$$

$$\delta ({\mathcal{H}}^2)\sim \frac{\delta t}{P_{\mathrm{dyn}}}(\delta {\mathcal{H}}_{\mathrm{dyn}}{)}^2.$$

$$\delta ({\mathcal{H}}^2)$$

$$\text{mbox}\delta {\mathcal{H}}_{\mathrm{dyn}}$$

$$\frac{\delta ({\mathcal{H}}^2{)}_{l_P}}{(\delta {\mathcal{H}}_{\mathrm{dyn}}{)}^2}=\frac{\delta M_{l_P,i}}{M_{\mathrm{rad},i}-M_{\mathrm{rad},i-1}}\frac{t_i-t_{i-1}}{P_{\mathrm{dyn},i}}.$$

$$\delta ({\mathcal{H}}^2)=\sum _j\delta ({\mathcal{H}}^2{)}_{l_P,j}\frac{r_j-r_{j-1}}{{\ell }_{P,j}}.$$

$${\Phi }_r\sim {\left[\delta ({\mathcal{H}}^2{)}_{l_P}\right]}^{1/4}.$$

$${\Phi }_r\to {\mathcal{H}}_{\mathrm{con}}^{1/2}\sim N_P^{1/4}{\left[\delta ({\mathcal{H}}^2{)}_{l_P}\right]}^{1/4}.$$

$$\frac{l(r)}{\sqrt{G{M}_{\mathrm{BH}}r}}=\frac{1-2j(R_g/r{)}^{3/2}+j^2(R_g/r{)}^2}{\sqrt{1-3(R_g/r)+2j(R_g/r{)}^{3/2}}}.$$

$$\frac{R_{\mathrm{isco}}}{G{M}_{\mathrm{BH}}/c^2}=3+Z-[(3-W)(3+W+2Z){]}^{1/2},$$

$$W$$

$$\text{mbox}Z$$

$$\frac{da}{dt}=-6k_2\frac{{\tau }_{\mathrm{fric}}}{{\tau }_{\mathrm{dyn}}^2}q(1+q){\left(\frac{a}{R_{\star ,1}}\right)}^{-8}a\left(1-\frac{{\Omega }_{\star }}{{\Omega }_{\mathrm{orb}}}\right)$$

$$\frac{d{L}_{\mathrm{orb}}}{dt}=\frac{1}{2}\frac{M_{\star ,1}{M}_{\star ,2}}{M_{\star ,1}+M_{\star ,2}}{\Omega }_{\mathrm{orb}}a\frac{da}{dt}.$$

$${\tau }_{\mathrm{fric}}={\left(\frac{M_{\mathrm{env},1}{R}_{\star ,1}^2}{L_{\star ,1}}\right)}^{1/3},$$

$$k_2=20.5{\alpha }^{4/3}\frac{R_{\star ,1}^{1/3}{\tau }_{\mathrm{fric}}}{M_{\mathrm{env},1}}{\int }_{x_{\mathrm{benv}}}^1{x}^{22/3}{L}_r^{1/3}{\rho }^{2/3}{l}_{P}dx.$$

$$\frac{B_r{B}_h}{4\pi }\sim {\Sigma }_{\mathrm{rad}}{v}_h{\Omega }_{\mathrm{NS}}.$$

$$\frac{R_{\mathrm{NS}}}{v}$$

$$v_{\mathrm{kink}}\sim \frac{r}{{\tau }_{\mathrm{kink}}}\sim \frac{B_{\varphi }^2}{8\pi \rho \Omega r}.$$

$$\frac{\partial B_{\varphi }}{\partial t}\sim -r{B}_r\frac{\partial \Omega }{\partial r}-\frac{B_{\varphi }}{{\tau }_{\mathrm{kink}}}.$$

$$v_{A,\varphi }^3\sim 2\left|\frac{d\ln \Omega }{d\ln r}\right|(\Omega r{)}^2{v}_{A,r}.$$

$$\frac{B_{\varphi }{B}_r}{4\pi }r\sim \rho l_P\frac{r^2\Omega }{{\tau }_J}.$$

Full text

1. Introduction

This paper addresses some gaps in the current understanding of the rotation and magnetism of massive stars as they evolve toward core collapse. Our focus is on the transport of angular momentum by convective and magnetic stresses, and on the genesis of stable magnetic fields in radiative layers of a star. We note in particular the following:

1. Existing stellar evolution codes do not allow for the advective transport of angular momentum by extended convective plumes. This may not result in significant inaccuracies for solar-type stars with shallow convective envelopes, or for rapidly rotating stars that support vigorous magnetic dynamos (see, e.g., Augustson et al. 2016). But we suggest that the situation is different for supergiants with deep and slowly rotating convective layers: in this case, the inward advection of a small fraction of the stellar angular momentum toward the core will have a dramatic effect on the inner rotation. The presence of deep convective plumes is supported by fits of mixing-length models to giant stars (Sackmann & Boothroyd, 1991), and by anelastic calculations of non-magnetized convection (Brun & Palacios, 2009). The consequences for the rotation of stars of solar and intermediate mass were previously considered by Kissin & Thompson (2015a, b), where the largest effect was seen near the tips of the red giant branch (RGB) and asymptotic giant branch (AGB).

Consider, for example, stars more massive than about $20\,M_{\odot}$ which also retain a hydrogen-rich envelope at the moment of core collapse. The envelope and burning shells of these stars form a broad sequence of convective layers. A key result of this paper is that strong angular momentum pumping within deep convective layers can significantly compensate the rotational braking of a contracting stellar core by internal stresses, thereby allowing for more uniform rotation across radiative-convective boundaries than is seen in preceding rotational models. A similar, but overall weaker, effect will be encountered in massive helium stars.

2. The radiative layers of an evolving star tend toward a state of differential rotation, which is resisted by the winding up of an embedded magnetic field. One popular approach to magnetic field growth in radiative layers starts with a weak seed polar magnetic field, and then relies on a kink instability of the wound-up toroidal field (Tayler, 1973) to feed back on the poloidal field (Spruit, 2002). Heger et al. (2005) incorporated this formalism into evolutionary models of massive stars. On the other hand, the magnetic fields of radio pulsars may indicate the presence of a stronger embedded magnetic flux which could significantly redistribute angular momentum within the contracting stellar core (e.g. Spruit & Phinney 1998; Maeder & Meynet 2014). We re-examine this process, taking into account the kinking of the wound-up poloidal flux, and argue that kinking does not much reduce its effectiveness. The competing effect of the magnetorotational instability (MRI) is much more localized near radiative-convective boundaries, and involves small-scale fluid motions that induce strong compositional mixing and significant structural changes in the star (Wheeler et al., 2015). These effects are greatly suppressed in our model by the smoothing action of large-scale Maxwell stresses.

3. Stability of the polar magnetic field threading a radiative layer of the star depends on the presence of a buried toroidal flux that carries a net twist (Braithwaite & Spruit, 2004). The accumulation of magnetic helicity within a radiative layer cannot simply result from the freezing of a large-scale fossil magnetic field, as is sometimes invoked to explain the magnetic moments of white dwarfs and neutron stars (NSs). Instead it depends on the transport of magnetic helicity across the evolving convective-radiative boundary that, in most cases, preceded the formation of the radiative layer (Kissin & Thompson, 2015b). There is a simple connection between this helicity flux and the large-scale Maxwell stress that is needed to compensate a persistent, inhomogeneous Reynolds stress that is imposed at a convective boundary. The required Maxwell stress is a natural outcome of even a mild dynamo instability.

The method just outlined allows us to identify the evolutionary stage that contributes most to the seed magnetic field in a collapsing stellar core, and to investigate how the spin rate of the collapsed core depends on the progenitor mass and spectral type. We find that the pre-MS accretion phase is the dominant source of magnetic helicity in the core of a massive star. The outer part of a growing massive star transitions from a convective to a radiative state while the mass still below $\sim 10\,M_{\odot}$. A significant merger could create a brief convective state that re-ignites this process, but would not necessarily generate much stronger helicity in the inner part of a massive star that had already reached the main sequence. Existing treatments of magnetic field amplification during a binary stellar merger (e.g. Wickramasinghe et al. 2014) have not so far accounted for the evolution of the magnetic helicity.

The rotation of radio pulsars and stellar-mass BHs is much faster than would be expected if their progenitors were able to relax continuously to a state of solid rotation. This is commonly taken to imply that the transport of angular momentum must freeze out in some parts of the progenitor during its evolution toward core collapse.

For example, Heger et al. (2005) obtain pulsar rotation periods $P_{\rm NS}\sim 3$-15 ms from models of initial (zero-age main sequence) mass $M_{\rm ZAMS}=12$-35$\,M_{\odot}$, and solar metallicity. These results are consistent with the fastest pulsar spins, but not with the substantial fraction of pulsars born spinning in the range $P_{\rm NS}\sim 0.1$-1 s (e.g. Popov & Turolla 2012). Although the multipolar structure of the poloidal magnetic field is not clear in the formalism of Spruit (2002), making an optimistic assumption of a dipolar radial field still implies a NS polar field $\sim 10^{10}$ G, substantially weaker than is typical of radio pulsars.

In comparison, the calculations presented here imply the presence of a stronger poloidal magnetic flux within the cores of massive stars, whose large scale nature is directly tied to the presence of net magnetic helicity. We find only modest variations between progenitor models in the magnetic field entrained in the collapsing core, corresponding to $B_{p}\sim 10^{12}$-$10^{13}$ G in the NS remnant. This demonstrates that the internal transport of magnetic helicity is an essential ingredient in a model of stellar rotation; we incorporate it here for the first time in the case of massive stars.

Because the internal magnetic torques are stronger than in the calculations of Heger et al. (2005), we obtain a wider range of spin periods. We find a significant dependence of $P_{\rm NS}$ on progenitor mass and spectral type. In part, this is because more massive stars develop thicker convective burning shells and experience stronger angular momentum pumping; and also because the angular moment lost by a supergiant to a wind depends strongly on its effective temperature. In combination with the relatively uniform pre-collapse magnetic field, this suggests that the rotation of the post-collapse NS is an important variable determining how the seed poloidal magnetic field evolves after the collapse, a possibility that we discuss briefly.

1.1. Additional Transport Processes and Simplified Approach

Most stars that experience core collapse (those born with a mass $M_{\rm ZAMS}\gtrsim 8\,M_{\odot}$) rotate rapidly during the hydrogen burning phase (Huang et al. 2010; Ramírez-Agudelo et al. 2015). Early numerical studies of the rotation of post-MS stars, which either neglected angular momentum transport or included it in an elementary way, found that the core angular velocity eventually surpassed the Keplerian rate (Kippenhahn et al. 1970; Endal & Sofia 1976). This is obviously problematic and we now know unrealistic.

A suite of mixing processes have been studied and implemented in one-dimensional evolutionary calculations over the past few decades: see for example Endal & Sofia (1978), Zahn (1992), Maeder & Zahn (1998), Meynet & Maeder (2000), Heger et al. (2000), and Wheeler et al. (2015). Mechanisms of angular momentum transport can be divided into two main categories: dynamic and secular. Dynamic processes operate on the convective or rotational timescale, a familiar example being the MRI (Balbus & Hawley, 1994), which can enforce nearly solid rotational motion within spherical shells over a modest number of rotations. Convection instead generates differential rotation in latitude (see Brun & Palacios 2009 for simulations of deep envelopes). Because our leading concern is the redistribution of angular momentum in radius, we adopt a simplified description of the rotation profile, taking the angular velocity to be a function only of spherical radius, $\Omega=\Omega(r)$. Such a rotation profile is sometimes described as ‘shellular’.

Secular processes include i) radial mixing that is driven by an angular velocity gradient and facilitated by thermal diffusion (Townsend, 1958; Goldreich & Schubert, 1967; Fricke, 1968; Zahn, 1992; Maeder & Zahn, 1998; Meynet & Maeder, 2000); and ii) the transport of angular momentum across convective-radiative boundaries by gravity waves, which are sourced by convective motions and then damp within radiative layers (Goldreich & Kumar, 1990; Kumar & Quataert, 1997; Zahn et al., 1997; Talon et al., 2002). We find that turbulent hydrodynamic stresses are subdominant to the Maxwell stresses that emerge from our model for magnetic helicity growth; and, furthermore, the mean polar flux produced is consistent with the post-collapse magnetic fields measured in radio pulsars. The estimates of internal Maxwell stresses made by Maeder & Meynet (2014), which are based on outward extrapolations of pulsar magnetic fields (and may therefore underestimate the large-scale poloidal field in the outer stellar core) are consisent with this conclusion. The MRI feeds off negative $d\Omega/dr$ in thin and weakly stratified layers near radiative-convective boundaries (Menou et al., 2004), but combining its effects with those of large-scale Maxwell stresses has not yet been attempted in stellar dynamo models and is beyond the scope of this work. Finally, we note that processes that operate on the Kelvin-Helmholtz timescale $\sim GM^{2}/RL$, such as meridional circulation in a massive and rapidly rotating star, are not generally competitive with those considered here.

Our rotational models combine two novel effects: large-scale helical Maxwell stresses, which dominate the transport of angular momentum in radiative layers of the star, and the inward advection of angular momentum in deep convective layers, which enforces strong radial differential rotational. For this reason, we take an otherwise simplified approach, turning off transport by other processes, especially rotationally induced mixing and internal gravity waves. In the latter case, predictions of the sign of the angular momentum transport have varied between different authors. However, it is possible that pumping of angular momentum into radiative layers by gravity waves sometimes overwhelms the smoothing effect of the large-scale Maxwell stress, especially if the exciting convective motions have a high Mach number. The rotation periods of collapsed stellar cores estimated by Fuller et al. (2015) are competitive with those obtained here for the lowest mass stellar models, but not for higher masses. Although the neglect of gravity wave transport may introduce the strongest systematic bias in our results, the influence of hydrodynamic instabilities on core properties can be more straightfowardly quantified, and is examined briefly.

A recent numerical simulation of core collapse starting from a dynamic state of oxygen burning in a $18\,M_{\odot}$ progenitor (Müller et al., 2017) finds that a fairly rapid ($P_{\rm NS}\sim 20$ ms) neutron star spin results from an asymmetric collapse and explosion. Velocity perturbations seeded by convection in the oxygen shell have been argued to facilitate an explosion (Couch & Ott, 2013; Müller & Janka, 2015). This result appears to depend on a combination of physical processes: first, the growth of non-spherical velocity perturbations in the collapsing material (due essentially to conservation of angular momentum: Lai & Goldreich 2000), which then seeds a global buoyancy instability of the hot shocked material (Thompson, 2000). Such dynamic effects fall beyond the scope of the approach advanced here.

Finally, rotational support during the collapse provides additional channels for angular momentum transport. The formation of a quasi-Keplerian torus is a key ingredient in the collapsar model of GRBs (Woosley, 1993). In the most massive models that form BHs, we are able to calculate the mass that is directly incorporated into the BH, and the mass which may be expelled in an outflow. More generally, we divide our stellar models into those which form a NS or a BH based on estimates of the compactness of the progenitor core that emerges from the MESA calculation. Previous calculations of massive stellar evolution with different codes, which were used to estimate threshold conditions for an explosion, imply slightly different relations between ZAMS mass and the compactness and iron mass of the evolved core (O’Connor & Ott, 2011; Ertl et al., 2016; Müller et al., 2016a). This implies an intrinsic theoretical ‘fuzziness’ in the progenitor mass that will evolve to a given set of pre-collapse core conditions. For this reason, implications for the post-collapse rotation of the most massive models are considered with both outcomes in mind.

1.2. Plan of the Paper

The set-up of the MESA models used in this paper is described in Section 2. Then in Section 3 we explain our prescription for angular momentum transport in convective and radiative layers, and how the rotation profiles of successive MESA snapshots are connected to each other. The role of the Coriolis force in limiting angular momentum pumping by convection is discussed in Section 3.2, and of kinking in limiting transport by the Maxwell stress in Section 3.3. We outline how the different parts of the model stars divide into transporting and non-transporting layers, and demonstrate the effectiveness of the adopted angular momentum transport mechanisms in comparison with more popular prescriptions for rotationally driven mixing. The accumulation of magnetic helicity in growing radiative layers is described in Section 4. The results for the overall evolution of the rotation and magnetic helicity in isolated model stars are presented in Section 5, with a focus on the rotation and polar magnetic flux of the remnant NS or BH. Spin-up of the stellar core by tidal angular momentum exchange with a massive stellar companion is explored in Section 6. Our main results and conclusions are summarized in Section 7. Appendix A shows our MESA model parameter lists, and Appendix B describes how kinking limits the growth of the Maxwell stress in radiative layers of a star containing a helical magnetic field.

2. Stellar Models

Using the one-dimensional stellar evolution code MESA (Modules for Experiments in Stellar Astrophysics: Paxton et al. 2011, version 8118), we created evolved models of several masses, beginning with pre-MS accretion onto a low-mass core and extending to core-collapse. We aimed to include models that i) produce both NS and BH remnants; ii) have a range of effective temperatures and convective penetration during the supergiant phase iii) have a wide range of peak luminosities, implying a range of convective depths within burning shells; and iv) experience a range of peak mass loss rates. We also allow for tidal interaction with a binary stellar companion by adding a source term to the stellar angular momentum.

Models of mass $M_{\rm ZAMS}$ = 13$\,M_{\odot}$ and 40$\,M_{\odot}$ were chosen as progenitors of NSs and BHs respectively, with 25$\,M_{\odot}$ as an intermediate case. Each of the stellar models has solar metallicity, with the 40$\,M_{\odot}$ model duplicated at 30% solar metallicity. These choices were guided by qualitative predictors of the outcome of core collapse. The first of these is the compactness parameter at core bounce,

[TABLE]

O’Connor & Ott (2011) found that $\xi_{M}$ evaluated at enclosed baryonic mass $M=2.5\,M_{\odot}$ implies an explosion for $\xi_{2.5}<0.4$.

More recently, Ertl et al. (2016) found a tighter relation between pre-collapse mass profile and ‘explodability’, which is the one we used. The division between NS and BH remnant corresponds to a linear curve in the plane defined by the enclosed mass and the radial mass derivative, both evaluated where the dimensionless entropy per baryon $s=4$. Figure 1 shows the four MESA models in this plane. The 13$\,M_{\odot}$ model appears to explode easily, and the solar-metallicity 40$\,M_{\odot}$ model appears to fail.

The 25$\,M_{\odot}$ model was drawn from a narrow mass range over which the core compactness drops significantly, suggesting the possibility of a successful explosion. This model lies close to the diagonal line in Figure 1 separating failed from successful explosions. The same is true of the low-metallicity 40$\,M_{\odot}$ model. This means that additional physical processes than neutrino heating, such as magnetorotational feedback, could play a significant role in determining the outcome of the core collapse. Interestingly, we find that one of these models does produce a rapidly rotating neutron star, whereas the other does not.

The 13$\,M_{\odot}$ and 25$\,M_{\odot}$ models expand to become red supergiants, whereas the two 40$\,M_{\odot}$ models remain blue supergiants. We find that varying the effective temperature during peak mass loss has a strong influence on the angular momentum carried off by the stellar wind, because convective envelopes of varying depths pump different amounts of angular momentum into the stellar interior. The same effect is encountered within the inner convective shells. At the moment of core collapse, the inner core stores more angular momentum in more massive progenitors, because they have deeper convective shells: compare the later stages of the 13$\,M_{\odot}$ model of Figure 2 with the 25$\,M_{\odot}$ model of Figure 3.

2.1. Building the Star

Each model star is assembled by gradual accretion, starting from a 3$\,M_{\odot}$ core. Growth by mergers with other stars is harder to implement in a one-dimensional evolution code. The accretion rate is taken to grow with time, $\dot{M}\propto t$, as suggested by numerical simulations and analytic calculations (Peters et al., 2011; Murray & Chang, 2015). The accretion rate is normalized so that the mass is fully assembled in $T=10^{5}$ yr,

[TABLE]

Thereafter accretion is shut off, and the only further change in mass is due to the stellar wind. Since a massive star can reach the main sequence while still accreting, we define the zero-age MS to coincide with the end of accretion.

The seed 3$\,M_{\odot}$ core is fully convective. We follow changes in the early convective structure because of its importance in magnetic helicity generation (Section 4). This structure differs little between the models. When the accretion rate reaches $\dot{M}\sim 0.5$-$1\times 10^{-4}\,M_{\odot}$ yr*-1*, the initial convection zone begins to retreat to the surface, and disappears when the stellar mass has grown to $M\sim 7.7$-$9\,M_{\odot}$. At the same time the convective core emerges. See Figure 4.

2.2. Prescription for Numerical Evolution

The MESA inlist file we used to generate the stellar models can be found in Appendix A. The free parameters with the greatest impact on the model behavior are those describing convective overshoot and mass loss. We now discuss each of these briefly.

2.2.1 Overshoot

The convective structure of massive model stars is sensitive to the prescription for convective overshoot. Strengthening the overshoot causes increased mixing, and significant changes in core mass. The amplitude of the overshoot is still poorly constrained, with some evidence for a dependence on stellar mass and metallicity: see Herwig (2000) as well as the discussion in Appendix B.7.2 of Paxton et al. (2013). We use an exponential parametrization of overshoot, meaning that the convective mixing coefficient drops on the radiative side of a radiative-convective boundary over $f_{\rm ov}=0.01$ times the pressure scaleheight $l_{P}$: $D_{\rm ov}=D_{0}e^{-|\Delta r|/f_{\rm ov}l_{P}}$, where $\Delta r$ the distance from the boundary. The amplitude $D_{0}$ is calculated a distance $f_{\rm ov,0}=5\times 10^{-4}l_{P}$ into the convection zone.

The amplitude of overshoot may well depend on a second parameter such as the magnetization of the convective material (Nordhaus et al., 2008), especially given that the buoyancy of magnetized material can be enhanced by a high radiative energy flux (Kissin & Thompson, 2015b). The latter effect is especially relevant during later stages of stellar evolution, or in massive stars.

2.2.2 Mass and Angular Momentum Loss

Mass loss plays an important role in the evolution of massive stars: our 40$\,M_{\odot}$ model loses about half its ZAMS mass by the moment of core collapse. In the absence of a surface magnetic field, as assumed here, the angular momentum lost is proportional to the ejected mass. We use the ‘Dutch’ mass loss prescription in MESA, which combines fitting formulae appropriate to different ranges of effective temperature and surface hydrogen abundance $X_{s}$, as described by Glebbeek et al. (2009).

2.2.3 Convective Stability and Semi-convection

The Ledoux criterion for convective instability is enforced. Semi-convective zones also develop within each model; although they influence the degree of mixing and the energy transport, they have a negligible effect on the angular momentum transport that we calculate in post-processing. That is because the fluid motions are relatively slow in semi-convective regions, so that angular momentum transport is dominated by the large-scale Maxwell stress.

2.3. Angular Momentum Transport and

Magnetic Field Evolution

We disable the MESA modules that handle internal rotation, as well as angular momentum transport by rotationally induced instabilities and the ‘Spruit-Tayler’ magnetic feedback process. We find that a large enough magnetic helicity is deposited in the radiative zones of our model stars to invalidate the starting assumptions of the ‘Spruit-Tayler’ process. In addition, radial angular momentum transport by the winding up of the pinned poloidal magnetic field is rapid enough to overwhelm rotationally driven mixing processes.

Instead the rotation profile is handled in post-processing, working with a spherical stellar model. One limitation of this approach is that the growth of the magnetic helicity in radiative layers of the star (Section 4), and the transport of angular momentum by convection and magnetic torques (Section 3), are decoupled from the mixing processes that influence the growth of the stellar core.

3. Rotational Evolution

We now describe the initialization of the stellar rotation, and explain our handling of angular momentum transport by convective and magnetic stresses. These processes operate in combination with each other: the effectiveness of magnetic stresses at enforcing nearly solid rotation in radiative layers of a star depends on the previous convective history.

3.1. Initialization of the Rotation

At the end of the accretion phase, we set the equatorial rotation speed to $v_{\rm rot,eq}=200$ km s*-1*, consistent with the measured average projected velocity $\langle v_{\rm rot,eq}\sin i\rangle={1\over 2}v_{\rm rot,eq}\sim 100$ km s*-1* in O and B stars (Huang et al., 2010; Ramírez-Agudelo et al., 2015). The rotation speed during the accretion phase is set to the same (fixed) fraction of the break-up speed, as determined by the evolving mass and radius.

3.2. Inward Pumping of Angular Momentum by Convection

We adopt the approach to convective angular momentum transport described in Kissin & Thompson (2015a). We make the plausible assumption that convective transport is rapid, with the rotation approaching its equilibrium profile over a modest multiple of the convective period $\tau_{\rm con}\sim l_{P}/v_{\rm con}$. We keep track of how $\tau_{\rm con}$ compares with the evolutionary time of the stellar mass profile, as defined by the Eulerian speed $v_{r}$ of the stellar material,

[TABLE]

Here $l_{P}=P/\rho g$ is the pressure scaleheight. The second term in (3) accounts for the rapid evolution of the inner burning shells just before core collapse (which happens at stellar age $t_{\rm cc}$). The last stages of nuclear burning are rapid enough that the convection in the intermediate and outer parts of the star effectively freezes out.

We restrict consideration of convective angular momentum pumping to thick shells and envelopes. These we define as having an aspect ratio $R_{\rm con+}/R_{\rm con-}\geq 2$, where $R_{\rm con+/-}$ is the radius of the top/bottom of the convective layer. This choice is motivated mainly by the observation of a small radial angular velocity gradient throughout most of the solar envelope ($R_{\rm con+}/R_{\rm con-}\sim 1.4$).

When the rotation is slow (as measured by the Coriolis parameter ${\rm Co}\equiv\Omega\tau_{con}\lesssim 1$), we allow convective plumes to conserve angular momentum, resulting in a rotation profile $\Omega\propto r^{-2}$. In faster rotating layers (${\rm Co}>1$) we account for the back reaction of the Coriolis force by considering the vorticity equation in the ‘thermal wind’ approximation (see Section 2 of Kissin & Thompson 2015a). The rotation profile in convection zones is therefore described by

[TABLE]

Here $\beta$ is the local radial power law dependence of gravitational acceleration ($g(r)\propto r^{-\beta}$). These scalings are connected to each other at radial shells where ${\rm Co}=1$.

Recently Klion & Quataert (2017) investigated core rotation in sub-giant and giant stars (radius $\sim(4-10)\,R_{\odot}$) during the early part of the first dredge up, where $\beta\sim 1$ and $\Omega\sim r^{-1}$ is predicted by Equation (4). They considered the relative splittings of pressure- and gravity-dominated modes as a probe of the relative rotation rates in the stellar core and envelope. Typical splittings were found to be marginally consistent with such an envelope rotation profile, and also to be consistent with most of the angular velocity offset concentrated in the layer between the hydrogen burning shell and the convective envelope. In the case of Kepler-56 ($R\sim 4\,R_{\odot}$, $M\sim 1.3\,M_{\odot}$ and a convective envelope extending a factor $\sim 3$ in radius), a $\sim r^{-1}$ convective angular velocity profile is not precluded; but this gradient must be extended all the way down to the burning shell, contradicting the assumption of solid rotation in all radiative layers by Kissin & Thompson (2015a). A recent application of the Klion & Quataert (2017) method to other Kepler giants by Triana et al. (2017) did not yield clear results. In this paper, we take a more general approach to rotation in radiative layers by calculating the limiting magnetic torque, as enforced by kinking of the wound up magnetic field (Section 3.3).

The feedback of magnetic fields on the rotation of a convective layer is an important outstanding issue. The MRI must be activated at some level where ${\rm Co}>1$ (e.g. Balbus et al. 2009), which is the case during the subgiant phase probed by the Kepler asteroseismological data. By contrast, the Coriolis parameter never becomes very large in our isolated star models during a supergiant phase (or, indeed, in the fully expanded RGB/AGB models of Kissin & Thompson 2015a). Since the convective Mach number is relatively high in the envelope near full expansion, the Coriolis parameter also remains modest within interior convective shells. We expect that a reduction in peak ${\rm Co}$ due to the activation of the MRI would, therefore, have a modest effect on the rotation rate of the inner core in a red supergiant. The effect may be larger in more compact stars, especially those which gain angular momentum from a binary companion.

3.3. Magnetic Angular Momentum Transport Limited by Kink Instability

The growing concentration of mass toward the center of an evolving star generates negative $\partial\Omega/\partial r$, but in the presence of strong stable stratification this does not trigger the MRI. In radiative zones we focus on the linear winding of the embedded radial magnetic field, which generates a growing toroidal magnetic field and $r\phi$ Maxwell stress. Because the magnetic field quickly becomes tightly wound, it is susceptible to a kink instability (Tayler, 1973). The growth of the kink is, however, impeded by the Coriolis force (Pitts & Tayler 1985; Appendix B).

We find that the net effect of the kink instability is a modest increase in the radial magnetic field needed to erase most of the radial angular velocity gradient. The threshold radial Alfvén speed $v_{{\rm A},r}\equiv B_{r}/(4\pi\rho)^{1/2}$ rises from $\sim|v_{r}|$ to

[TABLE]

We take solid rotation to be established in a radiative mass shell if $v_{{\rm A},r}>v_{{\rm A},r}|_{\rm min}$ within that shell. The second term in Equation (5) allows for a situation where $\tau_{\rm ev}$ becomes so small that the stellar material is not able to execute a full rotation, $\Omega\tau_{\rm ev}\lesssim 1$.

3.4. Rotational Evolution between Model Snapshots

The rotation profile $\Omega(r)$ is evolved as follows from one MESA model snapshot to the next:

1. In each radial zone we ask if $\tau_{\rm con}<\tau_{\rm ev}$, or $v_{{\rm A},r}>v_{{\rm A},r}|_{\rm min}$. If either inequality is satisfied then angular momentum is transported between the interior and exterior zones; otherwise the angular momentum of the given zone is frozen. In this way, we divide each model snapshot into transporting and non-transporting layers. Figure 5 shows the result over the full history of the 13$\,M_{\odot}$ model. Most of the star ceases to transport angular momentum effectively in the final $\sim 0.1$ yr before core collapse. This conclusion also applies to the outer convective envelope, where $\tau_{\rm con}$ becomes longer than the time to collapse.
2. The update to the rotation rate in a non-transporting zone is obtained by matching the specific angular momentum $\frac{2}{3}\Omega r^{2}$ between corresponding mass shells in the successive snapshots. The new rotation profile in each connected transporting region (comprising multiple MESA mass shells) is obtained by i) conserving the total angular momentum in that region; and ii) fitting the profile (4) in convective parts and solid rotation in radiative parts. This fit is often non-linear, in the sense that both the normalization and the position of breaks in the slope of the convective rotation profile (surfaces where ${\rm Co}=1$) will depend on the net angular momentum of the zone.
3. The complex convective structure of massive stars generates circumstances in which thin non-transporting layers are sandwiched between two transporting zones. In such a situation we envision that the two transporting regions would successfully communicate changes in angular velocity, if and only if the non-transporting layer is thinner than a pressure scaleheight $l_{P}$.

3.5. Some Basic Results

We now highlight some consequences of the rotational model just described, before presenting more detailed results in Sections 5 and 6. Figure 6 summarizes the rotational evolution of the 25$\,M_{\odot}$ model. Plotted is the rotation period $P_{\rm NS}$ that would result if the central 1.4, 1.6, 1.8, or 2$\,M_{\odot}$ of stellar material were instantly incorporated into a neutron star, while conserving angular momentum.

The result is shown as a function of time, to provide a measure of the changing angular momentum profile. The first rise in (the effective value of) $P_{\rm NS}$ is driven by the post-MS expansion of the star, and is followed by a rapid drop as the outer convective envelope forms and deepens. The pumping effect of the envelope is supplemented by the deep convective shells that form during the later stages of thermonuclear burning, resulting in a further decrease in $P_{\rm NS}$.

Figure 7 shows for comparison the effect of turning off angular momentum transport by magnetic stresses in radiative layers, while maintaining the pumping effect of the convection. Now the core does not couple to the outer layers of the star. Conserving its angular momentum from the ZAMS, one sees that the assumption of collapse to a hydrostatic NS leads to an inconsistency: the implied $P_{\rm NS}$ is shorter than a millisecond.

3.6. Comparison with Rotationally Induced Mixing

The efficiency of angular momentum transport by convective and Maxwell stresses, as implemented here, can be compared with the rotationally induced mixing (RIM) as formulated by Zahn (1992). The latter operates at low Prandtl number $\nu/\kappa_{T}$, where $\nu$ the kinematic viscosity and $\kappa_{T}$ the thermal diffusivity, as is appropriate to the radiative zone of a star. In this situation, strong $\partial\Omega/\partial r$ combines with radiative diffusion to trigger a linear axisymmetric instability (Goldreich & Schubert, 1967; Fricke, 1968), as well as higher wavenumber turbulent mixing (Townsend, 1958; Zahn, 1992). The corresponding timescale for angular momentum transport is

[TABLE]

Here $D_{\rm shear}$ is the diffusion coefficient, and $N$ is the Brunt-Väisälä frequency. We also show for comparison the timescale for meridional circulation,

[TABLE]

where $\Omega_{\rm K}$ is the Keplerian angular frequency at the given radius and enclosed mass.

These RIM and circulation timescales are compared in Figure 8 with i) the evolutionary time (3); ii) the convective time; iii) the radial Alfvén time $\tau_{{\rm A},r}=l_{P}/v_{{\rm A},r}$ as estimated from the hemispheric magnetic flux; and iv) the (maximum) Aflvén time $l_{P}/v_{{\rm A},r}|_{\rm min}$ that allows efficient angular momentum transport after allowing for kinking of the wound-up field. We choose a snapshot of the core of the 13$\,M_{\odot}$ model around the time (about 0.1 yr before core collapse) when the core and envelope begin to decouple. The RIM timescale is consistently 3-4 orders of magnitude longer than the Alfvén and convective timescales, and 2-3 orders larger than $\tau_{\rm ev}$; the circulation time even longer.

Where the core material is radiative, the inequality $\tau_{\rm A,r}<l_{P}/v_{{\rm A},r}|_{\rm min}$ holds and the Maxwell stress is large enough to enforce solid rotation. One sees from Figure 8 that this inequality is violated outside $\sim 2\,M_{\odot}$ enclosed mass, but here the material is convective. The convective timescale remains small enough to enforce the equilibrium rotation profile, given here by Equation (4), out to $\sim 3.5\,M_{\odot}$ enclosed mass. Beyond that point, the convection has essentially frozen out at this brief interval before core collapse.

We conclude that the addition of the shear instability would not directly change the rotation profile significantly in our calculation. The secondary effect of the influence of rotationally-induced mixing on core mass and composition is addressed in the next section.

3.7. Effect of RIM on Core Mass and Compactness

The introduction of rotationally induced mixing into a one-dimensional stellar model, using a simplified prescription (Zahn, 1992; Maeder & Zahn, 1998) is known to modify the iron core mass and the inner core compactness in models with rapid initial rotation. For example, Heger et al. (2000) found that when the ZAMS model rotates at about one-half of the breakup rate (as is typical of massive stars), more than $\sim 10$% changes in iron mass are possible in comparison with non-rotating models, even after allowing for the inhibition of mixing by mean molecular weight gradients. Here we choose a similar normalization of the mixing coefficients, with the parameters $f_{c}=0.03$ and $f_{\mu}=0.05$ in Equation (53) of Heger et al. (2000). In this case, in agreement with these previous results, the total helium core mass is less sensitive to the initial rotation rate: changes in $M_{\rm He}$ are limited to a few percent unless the rotation approaches breakup, or mean molecular weight gradients are (unrealistically) ignored ($f_{\mu}=0$). We focus in this paper on massive stars that do not start as close binaries and develop extreme rotation from mergers.

The ZAMS rotation rate is too blunt a measure of the influence of rotation on core properties. We have checked systematically how the pre-collapse core properties correlate with the strength of rotation at different stellar radii and at successive evolutionary stages. By far the strongest correlation turns out to be with the maximum rotation rate within the star at the first appearance of the helium core.This quantity is expressed as ${\rm max}(\Omega/\Omega_{K})$, where $\Omega_{K}=[GM(<r)/r^{3}]^{1/2}$ and $M(<r)$ is the enclosed mass, in which case the maximum is situated typically around the base of the convective envelope. We find ${\rm max}(\Omega/\Omega_{K})\sim 0.08$ for an initial fiducial rotation rate $\Omega_{\rm fid}$ equal to ${1\over 2}$ the breakup rate, using the prescription for angular momentum transport given in Section 3. Performing the same test during later burning stages (e.g. the appearance of the carbon core) shows a weaker dependence on the local rotation rate, and so we focus on the helium core rotation rate as the second variable modulating the core properties near collapse.

To explore the effect of variable core rotation, we construct a series of rotating MESA models with the built-in prescription for RIM and the Spruit-Tayler dynamo turned on, and the convective overshoot parameters normalized to the same values as in our non-rotating MESA models. This allows us to indirectly measure the error in pre-collapse core properties introduced by our neglect of RIM. Figure 9 shows how the pre-collapse core compactness $\xi_{2.5}$ and the total mass enclosed by the Si shell depend on ${\rm max}(\Omega/\Omega_{K})$ at the first appearance of the helium core, for $M_{\rm ZAMS}=25\,M_{\odot}$. The rotating MESA model with initial (solid-body) rotation rate $\Omega_{i}=\Omega_{\rm fid}/2.5$ gives the closest correspondence with the post-processed, non-rotating model. One finds that $\xi_{2.5}$ in our post-processed model is a factor $\sim 1.5$ smaller than is implied by the trend of $\xi_{2.5}$ with ${\rm max}(\Omega/\Omega_{K})$.

This rotation-dependent shift in core properties is combined with a strong sensitivity to progenitor mass. The latter effect has mainly been explored so far by neglecting the effects of RIM (Ertl et al., 2016; Müller et al., 2016a), and is also seen in our MESA models. The implication here is that RIM will shift the ZAMS mass that produces a given value of $\xi_{2.5}$. To gauge the magnitude of this shift, we plot in Figure 10 the pre-collapse core compactness for a range of $M_{\rm ZAMS}$, and both $\Omega_{i}=0$ and $\Omega_{i}=\Omega_{\rm fid}/2.5$. One sees that models producing the same $\xi_{2.5}$ are shifted in initial mass by 1-2$\,M_{\odot}$. The rapid variations seen in pre-collapse compactness as a function of $M_{\rm ZAMS}$, when combined with the effects of RIM, imply an intrinsic fuzziness in the relation between $M_{\rm ZAMS}$ and post-collapse core properties of this order.

4. Magnetic Helicity Accumulation

The inner parts of all our stellar models pass through multiple convective phases, which leave behind radiative material as they contract. There are two main contributions to ${\cal H}$ in these growing radiative layers (Kissin & Thompson, 2015b). The flux of magnetic helicity across the convective-radiative boundary is proportional to the magnetic torque which acts against growing differential rotation. A minimal source of differential rotation comes from the gradually changing mass profile of the star. But a stronger source is provided by inhomogeneous latitudinal rotation within the convective layer. The latitudinal gradient in rotation sources a radial gradient on the opposing (radiative) side of the boundary, where the magnetic field tends to enforce solid rotation. This process generates a toroidal magnetic field, and is a component of the hydromagnetic dynamo operating near the convective-radiative boundary.

The helicity is a conserved topological charge that, absent a boundary flux of magnetic twist, can only decay on a long resistive timescale. Net helicity is needed to stabilize the magnetic field in a radiative layer. The generic helical field configuration involves a finite open poloidal flux $\Phi_{r}$ that is surrounded by a twisted toroidal loop carrying flux $\Phi_{\phi}\sim{\cal H}/\Phi_{r}$ (Braithwaite & Spruit, 2004). The toroidal flux is confined by loops of poloidal field that close within the radiative material. We assume, following Kissin & Thompson (2015b), that the twisted magnetic field isotropizes as the helicity accumulates, so that $\Phi_{r}\sim\Phi_{\phi}\sim{\cal H}^{1/2}$.

The first convective phase is encountered during the pre-MS evolution, as the accreting massive star transitions from a fully convective state to the MS configuration of a radiative envelope and a convective core (Figure 4). Net magnetic helicity is left behind in the inner $\sim 8$-9$\,M_{\odot}$ of material because the convective envelope offers an escape route for the compensating helicity, e.g. via flaring activity. The remainder of the star is assembled after it has developed a radiative envelope, and its magnetization is determined by the interaction with the accretion disk. We (somewhat arbitrarily) set the magnetic field threading this outer material to zero.

We now describe the process of magnetic helicity pumping in a general way, and then turn to describe specific episodes of convective retreat in more detail. We first need an estimate of the helicity that would be stored in a shell adjacent to the convective-radiative boundary (at radius $R_{b}$). The shell thickness is determined by the Lagrangian speed $v_{\rm b}$ of the boundary through the stellar material. The helicity accumulated in a single hemisphere over a time $\delta t$ is then

[TABLE]

assuming that the magnetic field maintains a constant shape and strength over this interval. See Equation (27) of Kissin & Thompson (2015b).

The strength of the magnetic field in the dynamo layer is estimated by i) relating the toroidal field to the poloidal field through the linear winding term in the induction equation; and ii) assuming that the poloidal field is just strong enough to transfer angular momentum across a distance $\sim l_{P}$ over the dynamo period. Then

[TABLE]

Here $\Delta\Omega$ is the mismatch between the (nearly uniform) angular velocity of the radiative layer, and the angular velocity on the opposing convective side of the boundary. We take the dynamo period to be $N_{\rm dyn}\cdot\Delta\Omega^{-1}$, with $N_{\rm dyn}=10^{2}$. Further assuming that $\Delta\Omega\sim\Omega$, as appropriate for a slowly rotating and deeply convective layer (Brun & Palacios, 2009), one finds $\varepsilon_{B}\sim 10^{-3}$.

The integral of the helicity flux over latitude does not generally vanish within a single hemisphere, and maintains a uniform sign as long as the pole-equator angular velocity difference also maintains a constant sign (Kissin & Thompson, 2015b). There is, however, a cancellation between hemispheres if the magnetic field is reflection symmetric about the rotational equator. This cancellation can only be approximate if the rotation near the convective boundary is sustained by a modest number of deeply penetrating plumes. The strength of the cancellation is normalized here to $N_{\rm dyn}^{-1}\sim 0.01$ over a dynamo period, with a sign that varies stochastically over multiple cycles.

Then the net helicity deposited over the dynamo period $P_{\rm dyn}$ is

[TABLE]

The mass of radiative material added over the dynamo period is $\delta M_{\rm dyn}=4\pi R_{\rm b}^{2}v_{\rm b}P_{\rm dyn}\rho(R_{b})$. Given that the sign of the imbalance between hemispheres fluctuates randomly over multiple dynamo periods, the net helicity accumulated over an interval $\delta t\gg P_{\rm dyn}$ is

[TABLE]

This expression is easily implemented in a sequence of MESA model snapshots (labelled $i$), each of age $t_{i}$:

[TABLE]

This expression is independent of the time spacing between snapshots, as long and this spacing is fine enough to resolve the large-scale changes in convective structures.

The normalization of the helicity flux given by Equation (4) is supplemented by a threshold condition for the Coriolis parameter, measured close to the radiative-convective boundary. When ${\rm Co}<{\rm Co}_{\rm crit}$, the dynamo is shut off and the helicity flux vanishes. In this paper we take ${\rm Co}_{\rm crit}=0.1$.

We will need to ascribe a radial magnetic flux $\Phi_{r}$ to each radial zone in each snapshot. As the radiative material grows, we do this by summing $\delta({\cal H}^{2})$ over the last shell of thickness $l_{P}$ added to the zone. Thus $\delta({\cal H}^{2})\sim\Phi_{r}^{4}$ is the Lagrangian variable which follows each radiative mass shell, and which is modified by successive convective structures. The mass in this shell is $\delta M_{l_{P}}=4\pi R_{\rm b}^{2}\rho(R_{\rm b})l_{P}$, meaning that

[TABLE]

The sum (4) is recovered from the Lagrangian variable $\delta({\cal H}^{2})_{l_{P}}$ by summing over the mass shells (labelled $j$) of the radiative material:

[TABLE]

The hemispheric flux threading a single mass shell is estimated as

[TABLE]

Several pressure scale heights of material (number $N_{P}$) are mixed together when convection is excited in the core during hydrogen and helium burning. The total helicity contained in the core is ${\cal H}_{\rm con}\sim N_{P}^{1/2}[{\cal H}^{2}_{l_{P}}]^{1/2}$, and the flux increases after mixing to

[TABLE]

4.1. Successive Convective Structures

We calculate the contribution to ${\cal H}$ from three successive convective episodes: i) the transition from a convective to a radiative envelope during the pre-MS phase; ii) the contraction of the He-rich convective core at the end of the MS; and iii) the contraction of the C-rich convective core at the end of core He burning. The first, pre-MS, contribution dominates in all the models we consider. Figures 11-13 show the effective hemispheric magnetic flux, as corrected by convective mixing where appropriate (Equation 16), and the flux distribution at the moment of core collapse.

The relative importance of the pre-MS contribution to ${\cal H}$ arises from the fast recession of the envelope: in the 25$\,M_{\odot}$ this takes place in $\Delta t\sim 3\times 10^{4}$ yr, as compared with $\sim 7\times 10^{6}$ yr for the recession of the convective core at the end of the MS. The cancellation factor $\sim P_{\rm dyn}/\Delta t$ in ${\cal H}^{2}$ is much stronger in the later phase.

We now summarize some of the finer details of this process.

End of Core He Burning. Even though the recession of convective He burning core has the potential to add helicity to the inner 2-12$\,M_{\odot}$ (depending on the ZAMS mass), in all of our models the Coriolis parameter during this phase is $<0.1$. We therefore shut off this contribution.

Multiple convective shells. The complex convective behavior of massive stars creates numerous instances of receding convection zones, which we cannot track individually. To account for this contribution, we arbitrary add a term ($10^{25}\,{\rm Mx})^{4}$ to ${\cal H}^{2}$. This turns out to have a negligible effect on the total helicity and angular momentum transport.

Expulsion of helicity from a slowly rotating convective envelope. We sometimes find that the entire supergiant envelope is slowly rotating, meaning that it cannot sustain an active dynamo but still can lose magnetic helicity through the stellar surface, especially during strong mass loss. We allow this ejection process to occur if the radiative buffer at the surface is thinner than five pressure scaleheights. The effects of surface helicity ejection can be seen in Figures 11-13, which show a radial cutoff in helicity marking the deepest penetration of the slowly rotating convective envelope.

5. Combined Magnetic and Rotational Evolution to Collapse

The mechanisms of angular momentum transport described in Section 3 work in concert with the mechanism of magnetic helicity deposition outlined in Section 4. In this section, we contrast the behavior of the various stellar models, as manifested especially by the magnetization and rotation of the post-collapse remnants.

The rate of helicity deposition in growing radiative layers depends on the star’s angular velocity profile and history of angular momentum loss. Conversely, the embedding of magnetic helicity facilitates a rotational coupling between core and envelope, through intermediate radiative and semi-convective zones.

5.1. Dependence on Effective Temperature

Both 40$\,M_{\odot}$ models become blue supergiants ($T_{\rm eff}\sim 8$-9000 ${\rm K}$), and we find that they end up with far less angular momentum than the lower-mass models, which expand to become red supergiants ($T_{\rm eff}\sim 4000\,{\rm K}$). Figure 14 shows the relation between the angular momentum and mass carried away by the stellar wind. During the supergiant phase the rotation is slow enough that $\Omega(r)\propto r^{-2}$ throughout most of the envelope. This means that the depth of the envelope has a strong influence on the angular momentum stored close to the surface. The result is specific angular momentum several times higher near the surfaces of blue supergiants.

It is worth emphasizing that this result does not depend on the detailed reasons why both 40$\,M_{\odot}$ models fail to become a red supergiant. (See Woosley & Weaver 1995 and Ekström et al. 2012 for a summary of current thinking on this issue.) The stronger wind mass loss may be responsible. Metallicity does not appear to be the key variable, and we are not including rotationally driven mixing effects, or (in this section of the paper) allowing for angular momentum injection from a binary companion. The handling of convective overshoot may also be relevant, but a full exploration of its effects is beyond the scope of this paper.

5.2. Dependence on Stellar Mass

A significant difference in the rotation rate of the inner core is observed between the 13$\,M_{\odot}$ and 25$\,M_{\odot}$ models: the center of the lower-mass model rotates much more slowly. To explain this result, we show in Figure 15 the angular momentum profiles at different times preceding core collapse. The burning shells develop much thicker convective layers in the 25$\,M_{\odot}$ model. They are slowly rotating enough to be effective inward pumps of angular momentum, with a large ratio of inner to outer angular velocities. These convective layers correspond to the flat parts of the specific angular velocity curves in Figure 15. As a guide, we mark off the radii of several enclosed mass coordinates ($m=1.4$, 1.7 and 2$\,M_{\odot}$). Close to core collapse, the specific angular momentum is about two orders of magnitude larger at $m=1.4\,M_{\odot}$ in the 25$\,M_{\odot}$ model than in the 13$\,M_{\odot}$ model.

5.3. Neutron Star Remnants: Rotation

and Dipole Magnetic Flux

Consider now the spin period $P_{\rm NS}$ of a NS that forms from a collapsed baryonic mass $M_{\rm col}$. This depends on $M_{\rm col}$ as well as on the angular momentum profile of the pre-collapse core. A first estimate of $M_{\rm col}$ is obtained from the observation that the infall of the oxygen shell causes an outward expansion of the post-collapse standing shock, which when combined with neutrino heating can drive an explosion (e.g. Müller et al. 2016b). Taking $M_{\rm col}$ to be the mass enclosed by the silicon shell $+\,0.1\,M_{\odot}$, we obtain $\sim 1.7\,M_{\odot}$ for the 13$\,M_{\odot}$ model and $\sim 2.0\,M_{\odot}$ for the 25$\,M_{\odot}$ model. Figure 15 shows that the lever arm is raised by a factor $\sim 2$, and the specific angular momentum by a factor $\sim 4$, when moving outward from $m=1.4$ to 1.7$\,M_{\odot}$ in the 13$\,M_{\odot}$ core at a short interval ($\sim 10^{-3}$ yr) before core collapse. Increasing the enclosed mass from 1.4 to 2$\,M_{\odot}$ raises the specific angular momentum by approximately the same factor in the 25$\,M_{\odot}$ pre-collapse core.

The dependence of $P_{\rm NS}$ on $M_{\rm col}$ is shown in Figures 16 and 17. Here the evolving angular momentum profile is probed by assuming contraction to a cold neutron star over a range of times in advance of the actual core collapse. The angular momentum of the collapsed mass is equated with $2\pi P_{\rm NS}^{-1}I_{\rm NS}$. The cold NS moment of inertia111The energy lost to neutrinos is neglected here, meaning that the gravitational mass of the NS is slightly overestimated; but in compensation our choice of $R_{\rm NS}$ is slightly smaller than most modern equations of state would suggest. is taken to be $I_{\rm NS}=0.35M_{\rm col}R_{\rm NS}^{2}$, and the radius $R_{\rm NS}=10$ km.

The rotation of the 13$\,M_{\odot}$ remnant ($P_{\rm NS}\sim 0.2$ s for $M_{\rm col}=M_{\rm Si}+0.1\,M_{\odot}$) is not atypical of pulsars: a significant fraction of NSs appear to be born with spin periods in the range 0.1-1 s (Popov & Turolla, 2012). On the other hand, the 25$\,M_{\odot}$ model has a predicted spin period around 1.5 ms, making it a strong candidate for further dynamo amplification of the magnetic field post-collapse (Thompson & Duncan, 1993). We also include, for completeness, the $P_{\rm NS}$ that would be obtained from the 40 $M_{\odot}$ models if they were able to explode and leave behind a stable NS below the maximum mass. As expected from the greatly increased angular momentum loss in these models, some 2-4 orders of magnitude less angular momentum is stored in the core than in the lower-mass models (Figure 17). The remnant spin is calculated in Section 5.4 in the more likely scenario of BH formation.

The blue lines in Figures 11 and 12 show the magnetic flux threading the hydrogen-depleted core just before collapse. The thin red line represents the convective core on the ZAMS, which has nearly the same amplitude but covers a wider range of mass; this demonstrates that the contribution from later convective stages is small. Here we assign a single number $\Phi_{r}={\cal H}_{\rm con}^{1/2}$ for the polar flux, where ${\cal H}_{\rm con}$ is the magnetic helicity stored in the convective core. This must overestimate the NS poloidal flux: only the inner $\sim 15$-$30\%$ of the core material is incorporated into the NS. For example, if the core had a constant ratio of flux to mass before the collapse, then $\Phi_{r}$ should be corrected downward by a factor $\sim 0.15^{2/3}$-$0.3^{2/3}\sim 0.3$-0.4.

Applying this correction to the polar magnetic field $B_{p}\sim\Phi_{r}/\pi R_{\rm NS}^{2}\sim{\cal H}_{\rm con}^{1/2}/\pi R_{\rm NS}^{2}$, we obtain $B_{p}\sim 2\times 10^{13}$ G in the cold NS. This is in the upper part of the pulsar dipole field distribution: the dipole fields of isolated pulsars are typically $B={1\over 2}B_{p}\sim(3\times 10^{11}$ - $3\times 10^{13})R_{\rm NS,6}^{-3}$ G (Bhattacharya & van den Heuvel 1991, with $R_{\rm NS,6}=R_{\rm NS}/10^{6}$ cm). This translates into a hemispheric flux $\Phi_{p}\sim\pi B_{p}R_{\rm NS}^{2}\sim(2\times 10^{24}$ - $2\times 10^{26})R_{\rm NS,6}^{-1}$ Mx. For comparison, Heger et al. (2005) used the dynamo model of Spruit (2002) and obtained a much lower flux, $\Phi_{p}\sim 1\times 10^{22}$ Mx. The implications for the origin of pulsar magnetism, including the growth or decay of the dipole field post-collapse, are discussed further in Section 7.3.

5.4. Black Hole Remnants: Mass and Spin

Black holes are expected to form during the collapse of the most massive stars, as represented by our solar-metallicity 40$\,M_{\odot}$ model, and possibly one or other of the 25$\,M_{\odot}$ model and the $0.3\,Z_{\odot}$, 40$\,M_{\odot}$ model (see Figure 1). The remnant of such a collapse may engulf a large part of the progenitor. The fraction of the progenitor that is accreted is regulated by the spin up of the infalling material: conserving angular momentum, it may become rotationally supported. An energetic outflow from the black hole ergosphere and the inner parts of an orbiting disk would expel outer mass shells which have not yet become strongly bound to the hole, thereby limiting its mass.

As the black hole builds up by accretion, its mass $M_{\rm BH}$ and angular momentum $J_{\rm BH}$ are nearly equal to those of the accreted precollapse core material (enclosed mass $m$). Some loss of energy and angular momentum to neutrinos occurs in the first part of the collapse, before the event horizon forms, but this represents a small correction after the hole has grown to several solar masses. We therefore neglect neutrino losses here. The first mass shell that forms a centrifugally supported disk can be estimated by comparing the specific angular momentum $j_{\star}$ of the infalling material with that of the innermost stable circular orbit (ISCO). The specific angular momentum in an orbit of radius $r$ is a function of $M_{\rm BH}=m$ and $J_{\rm BH}=J(m)$ (Frolov & Novikov, 1998),

[TABLE]

Here the spin parameter $j\equiv J_{\rm BH}c/GM_{\rm BH}^{2}$ and $R_{g}\equiv GM_{\rm BH}/c^{2}$. We evaluate $l_{\rm isco}\equiv l(R_{\rm isco})$ at the radius of the ISCO,

[TABLE]

where

[TABLE]

Figure 18 compares $j_{\star}$ and $l_{\rm isco}$ as functions of $m$ in the 40$\,M_{\odot}$ models. In both cases, so much angular momentum has been lost to the wind that each star should almost entirely collapse into a BH (of mass 18.3$\,M_{\odot}$ in the $Z=Z_{\odot}$ model, versus 21.6$\,M_{\odot}$ for $Z=0.3Z_{\odot}$).

6. Angular Momentum Injection from a

Binary Companion

Most massive stars form in multiple systems: for example, Kobulnicky & Fryer (2007) find that the multiple fraction is greater than 80%. Interaction with a binary companion is therefore the rule rather than the exception. Strong interactions involving conservative mass transfer, a merger, or a common envelope phase are beyond the scope of this paper.

There is, however, a range of binary separations over which there is a gradual transfer of angular momentum by tides to the spin of the more massive (primary) star during its supergiant phase. This effect is implemented in our model stars by applying a positive torque to the hydrogen envelope. We assume a companion mass equal to ${1\over 2}$ the primary ZAMS mass, and consider a range of binary separations with the orbital plane aligned with the initial rotation of the primary.

The orbital separation $a$ and angular momentum $L_{\rm orb}$ evolve according to (Hut, 1981)

[TABLE]

and

[TABLE]

Here the primary is labeled ‘1’ and the secondary ‘2’, the mass ratio is $q\equiv M_{\star,2}/M_{\star,1}$, and $\Omega_{\rm orb}$ is the circular orbital frequency. The tidal friction time is

[TABLE]

where $M_{\rm env,1}$ is the envelope mass of the primary, and we define $\tau_{\rm dyn}=(R_{\star,1}^{3}/GM_{\star,1})^{1/2}$. The Love number is obtained from an integral over the convective envelope,

[TABLE]

Here $\alpha=2$ is the convective mixing length parameter, $x=r/R_{\star,1}$, $x_{\rm benv}$ marks the base of the convective envelope, and $L_{r}$ is the luminosity at a given radius. As we integrate expressions (20) and (21) for $a$ and $L_{\rm orb}$, we linearly interpolate the various stellar parameters between MESA model snapshots.

The tidal interaction is concentrated near the maximum expansion of the primary, given the strong dependence of the torque on the aspect ratio $R_{\star,1}/a$. This simplifies the calculation, because the hydrogen envelope is mostly radiative during the early expansion, complicating the handling of the tidal dissipation. On the other hand, a deep outer convective envelope redistributes angular momentum rapidly, allowing us to apply the rotation model described in Section 3.2.

Here a significant difference emerges between the two 40$\,M_{\odot}$ models (Figure 19). The solar metallicity model reaches a maximum radius and then contracts while mass loss continues. This means that most of the angular momentum deposited by tides (or by a merger with a lower-mass companion) near the maximum expansion would be subsequently lost. The situation is different with the $0.3\,Z_{\odot}$ model, which maintains a nearly uniform radius as a supergiant, and then experiences a late secondary expansion during which the tidal torque peaks but there is little time for additional mass loss.

The binary separation is started with a range of values, but the integration is stopped (and the evolutionary track discarded) if the internal rotation of the primary approaches breakup at any point within the convective envelope. Applying this cut also allows us to discard binaries which merge.

The resulting core rotation rate is shown in Figure 20 for the 13$\,M_{\odot}$ and 25$\,M_{\odot}$ models. In both cases, the maximum spin-up of the core, compared with an isolated star, is about a factor $\sim 2$.

The lower metallicity 40$\,M_{\odot}$ model experiences a much greater relative spin-up (a factor $\sim 10^{2}$-$10^{4}$: compare in Figure 21 results for $a_{i}=\infty$ to the other results), a consequence of the strong angular momentum loss to a wind. The consequences for the rotation and accretion of a BH remnant are shown in Figure 22. Now the infalling material becomes centrifugally supported outside an enclosed mass $\sim 12.6\,M_{\odot}$.

An energetic outflow from the accretion torus that forms as a result could limit much further growth of $M_{\rm BH}$. Here it is interesting to note that the outflow power is a strong increasing function of the black hole spin (Tchekhovskoy et al., 2010). The accretion of an additional $\sim 1.8\,M_{\odot}$ of high angular momentum core material (following the first appearance of the torus) is needed to bring $J_{\rm BH}$ close to the limiting value $GM_{\rm BH}^{2}/c$ (Figure 23). We also deduce that a successful magnetorotational explosion, leading to the formation of a NS, cannot be expected for these binary parameters.

The observed properties of Galactic binary stellar-mass BHs are consistent with such a scenario. Fitting of thermal X-ray spectra suggests in many cases a high spin rates for the hole, with $cJ_{\rm BH}/GM_{\rm BH}^{2}$ approaching unity (e.g. McClintock et al. 2006, 2011). Such a high spin could not have been gained by mass transfer from the companion star. A similar effect should occur in closer binaries in which the primary envelope is ejected during a common envelope phase.

7. Conclusions

We have investigated the evolving magnetism and rotation of massive stars, starting in the pre-MS accretion phase and extending to the last stages of nuclear burning and the onset of core collapse. Our focus is on two effects which have received little attention in models of stellar evolution: i) inhomogeneous rotation in deep convective layers, especially slowly rotating layers with a small Coriolis parameter $\Omega\tau_{\rm con}$; and ii) the pumping of magnetic helicity into a growing radiative layer at a moving convective-radiative boundary. The impact of these processes on the magnetism and rotation of lower-mass giant stars was considered previously by Kissin & Thompson (2015a, b).

The inward advection of even a small amount of angular momentum through a deep convective layer can have a major effect on the rotation rate of the stellar core and its collapsed remnant. This effect is more important for red than blue supergiants, and for more massive stars with thicker convective burning shells.

We account for the limiting effect of a kink instability on the transport of angular momentum by the $r\phi$ Maxwell stress through radiative layers of a star. Our approach to magnetic field growth in radiative layers differs substantially from that of Spruit (2002) in that it also incorporates a large-scale poloidal flux, which is stablized by buried magnetic twist.

We expand on previous efforts by tracking the growth of the magnetic field in the stellar core during the pre-MS accretion phase, when the star acquires most of its mass. This first phase of convective retreat turns out to be the dominant contributor to the magnetic helicity stored in the core, yielding a poloidal flux density around $10^{13}$ G in the NS remnant.

A relatively simple rotation profile is found over much of the evolutionary history of our model stars, corresponding to a significant redistribution of angular momentum between the inner and outer parts of the star. Radiative layers tend to rotate as solid bodies, and deep convective layers approach constant specific angular momentum, as motivated by anelastic calculations of deep and very slowly rotating stellar envelopes (Brun & Palacios, 2009).

To test these effects, we chose two MESA models which are likely to leave behind a NS remnant (13$\,M_{\odot}$ and possibly the 25$\,M_{\odot}$, based on the explosion criterion of Ertl et al. 2016 as implemented in the MESA output). These were supplemented by the two 40$\,M_{\odot}$ models, the more metal rich of which likely forms a BH. The rotation of the collapsed remnant (either NS or BH) varies much more between these models than does the magnetic flux. That is because the rotation is sensitive to the depth of the outer convective envelope during later stages of mass loss, and because the most massive stars form thicker convective shells. The predicted NS rotation period ranges from the millisecond range (in the 25$\,M_{\odot}$ model), up to $0.1-1$ s (in the 13$\,M_{\odot}$ model). These cover the full range of the estimated initial rotation rates of pulsars (Popov & Turolla 2012), implying that radio pulsars can form from isolated massive stars.

Much longer spin periods are found in the highest mass models, which end up as blue supergiants and collapse to black holes. This means that the entire star is likely to collapse through the event horizon of the BH, corresponding to a mass $M_{\rm BH}\sim 20\,M_{\odot}$ for our 40$\,M_{\odot}$ progenitors.

We find that tidal angular momentum exchange from a binary companion with semi-major axis 3000-20000$\,R_{\odot}$ only spins up the core by a factor $\sim 2$ in the two lower-mass models. On the other hand, enough angular momentum is deposited in the 40$\,M_{\odot}$, $0.3\,Z_{\odot}$ model star to produce a rapidly rotating BH: the outer third of the collapsing star forms a rotationally supported torus and may be expelled in a disk-driven outflow, leaving a remnant of mass $M_{\rm BH}\sim 13$-15$\,M_{\odot}$.

These calculations of the magnetism and rotation of the remnants of massive stars bear a promising resemblance to observations. It is possible that one or both of mechanisms investigated here is manifested more weakly in real stars. Lower magnetization of the radiative layers of a star could partly compensate weaker inward pumping of angular momentum through its convective layers. Nonetheless, an intermediate level of differential rotation in the convective envelope would still have a profound effect on the core rotation rate during the supergiant phase. And even a factor $\sim 10^{-2}$ reduction in the amplitude of the helicity flux (compared with our approach) would leave behind a poloidal magnetic flux well within the range observed in radio pulsars.

7.1. Comparison with Previous Models

The most direct comparison is with the models of Heger et al. (2005) and Wheeler et al. (2015). Both models include the effects of magnetic stresses as evolved according to the Spruit (2002) model, to which the second adds mixing by the MRI excited near compositional boundaries. As was explained in Section 1, the core magnetic flux preceding collapse cannot be usefully extracted from these earlier calculations, because the magnetic helicity is not computed, and the magnetic instabilities considered have a high spatial wave number. However, one does deduce that the poloidal magnetic field is systematically weaker, by at least a few orders of magnitude, than is suggested by pulsar magnetic fields, and is found in our calculations of helicity growth during pre-MS accretion.

We therefore focus on comparing core rotation rates. Both the preceding calculations combine an assumption of solid rotation in the deep convective envelope and convective hydrogen-depleted layers of an evolved star, with an incomplete rotational coupling across radiative-convective boundaries. Our model has opposing properties, meaning that any near agreement between the results for central core rotation should be viewed as fortuitous. One can certainly imagine a sequence of rotational models with decreasing angular velocity gradient in the convective envelope, which is compensated by a growing (negative) angular velocity jump(s) across internal compositional boundaries. The difficulty in constructing such intermediate models lies in developing a deterministic approach to reconnecting the large-scale poloidal magnetic field at compositional boundaries, and in prescribing the envelope angular velocity profile over a wide range of Coriolis parameter. A purely hydrodynamic approach to the latter problem was developed by Kissin & Thompson (2015a), but the result probably also depends on magnetic feedback, especially from the MRI where ${\rm Co}\gtrsim 1$. Nonetheless, it should be kept in mind that solid rotation in the envelope of a red supergiant would imply extremely low Co throughout the envelope, pushing conditions toward the regime where hydrodynamic stresses may dominate.

That being said, Wheeler et al. (2015) find a specific equatorial angular momentum $j\sim 3\times 10^{14}$ cm2 s*-1* at an enclosed mass $\sim 1.8\,M_{\odot}$ in a 15 $M_{\odot}$ model, and $j\sim 3\times 10^{13}$ cm2 s*-1* in a 20 $M_{\odot}$ model. By comparison, Figure 15 shows $j\sim 10^{14}$ cm2 s*-1* at the same enclosed mass in our 13 $M_{\odot}$ model, increasing to $3\times 10^{15}$ cm2 s*-1* in the 25 $M_{\odot}$ model. Taking into account the signficant scatter in core compactness and mass that is expected with a changing ZAMS mass, these results can at best be viewed as overlapping; and both show some promise in comparison with pulsar rotation rates. The calculations of Heger et al. (2005) do not include magnetic torques of any type at radiative-convective boundaries, and therefore typically yield stronger angular velocity jumps, and systematically faster core rotation.

7.2. Systematic Uncertainties in the Model

Much of the systematic uncertaintly in our approach is related to the still developing understanding of differential rotation in deeply convective layers. Here we explore how our predictions would be altered in the presence of a slightly different pattern of differential rotation.

First, the inward pumping of angular momentum by deep convective plumes could result in a shallower angular velocity profile than the one ($\Omega(r)\propto r^{-2}$) suggested by anelastic calculations of slowing rotating envelopes (Brun & Palacios, 2009). Supposing that the index softens to $-2+\alpha$, and taking into account that a large fraction of the stellar angular momentum is contained in the envelope, one finds that the core rotation period is increased by a factor $\simeq(1+4\alpha/3)^{-1/2}(R_{\star}/R_{\rm base})^{\alpha}\sim 8$ for $\Omega(r)\propto r^{-3/2}$ and $R_{\star}\sim 10^{2}R_{\rm base}$. (The coefficient here corresponds specifically to a density profile $\rho(r)\propto r^{-3/2}$ in the envelope.) Even with such a revision, the angular velocity at the base of the convective envelope remains dramatically faster relative to the surface than in existing evolution codes such as MESA.

Second, we chose a strong latitudinal angular velocity gradient at the base of the convective layer, again as observed in anelastic calculations of deeply convective layers. A weaker level of differential rotation would reduce the magnetic helicity flux in proportion to ${\cal H}\propto(\partial\Omega/\partial\theta)^{2}$ (following Equations (9) and (10)), and the remnant hemispheric magnetic flux as ${\cal H}^{1/2}\propto\partial\Omega/\partial\theta$. A reduction by an order of magnitude of the strength of differential rotation corresponds to a polar magnetic field $B_{p}\sim{\cal H}^{1/2}/\pi R_{\rm NS}^{2}=2\times 10^{12}$ G, which is still well within the range observed in young radio pulsars.

Third, we explored in Section 3.7 how the introduction of RIM into one-dimensional MESA calculations would modify the compactness and mass of the collapsing stellar core. We estimated a $\sim 0.1\,M_{\odot}$ error in mass enclosed by the silicon layer from the dependence of $M_{\rm Si}$ on peak rotation rate (relative to Keplerian) at the formation of the helium core (Figure 9). At a fixed progenitor mass, the simplest measure $\xi_{2.5}$ of the core compactness is estimated to rise from $\sim 0.2$ to $\sim 0.35$. However, much of these changes can be compensated by a modest 1-$2\,M_{\odot}$ adjustment of the progenitor mass, as shown in Figure 10.

7.3. Implications for Growth or Decay of the Magnetic Field Post-Collapse

The dipole magnetic fields of radio pulsars carry a minuscule fraction of the neutron star binding energy (about $10^{-12}(B_{p}/10^{12}~{}{\rm G})^{2}$, where $B_{p}$ is the polar flux density). This complicates our understanding of their origin: a variety of processes might contribute to such a relatively weak field. The $\sim 10^{13}$ G NS magnetic field that our models produce (over a wide range of progenitor masses) is moderately stronger than the spindown field of most radio pulsars, but weaker than the magnetic fields of active magnetars. Our next task is therefore to consider how the field may be modified post collapse.

Rapid neutrino-driven convection in the proto-NS, with an overturn time $\tau_{\rm con}\sim 3$ ms, can have two competing effects on the magnetic field (Thompson & Duncan, 1993). First, the entrained field lines will diffuse across the surface of the star over a reasonably short timescale, which depends on the initial rotation period. A mixing together of radial fluxes of opposing signs from the two magnetic hemispheres would reduce the external magnetic moment. We show that such a reduction is most feasible in slow rotators.

The convective motions do not extend fully to the surface of the proto-NS (Lattimer & Mazurek, 1981). A stably stratified layer of mass $M_{\rm rad}\sim 0.1\,M_{\odot}$ maintains significant inertia. The magnetic field lines which thread both the convective material and this more inert shell are stretched in the non-radial direction near the interface between the two, producing a strong horizontal field $B_{h}\gg B_{r}$.

The shuffling motion is slow enough to be further impeded by the Coriolis force. Integrating this through the outer shell (column $\Sigma_{\rm rad}=M_{\rm rad}/4\pi R_{\rm NS}^{2}$) gives an estimate of the speed $v_{h}$ of the shuffling motions,

[TABLE]

The shuffling timescale is then

[TABLE]

For the estimated post-collapse field $B_{r}\sim 2\times 10^{13}$ G, and for a spin period longer than $\sim 0.1$ s, the timescale (7.3) is comparable to the Kelvin-Helmholtz timescale of a few seconds. Slower rotators can therefore experience greater dipole cancellation.

The second effect is encountered in the regime of faster rotation. A proto-NS spinning with $P_{\rm NS}\lesssim 3$-10 ms would generate strong toroidal magnetic fields from the mean poloidal field. This happens preferentially after the bounce shock has expanded away and the high-entropy mantle surrounding the star has collapsed (over perhaps the first $\sim 0.3$ s); but before the bulk of the interior has cooled. During this phase, a positive radial angular velocity gradient develops in the outer parts of the star, driven by the collapse of the mantle. This shear is strong enough for the radial field to experience linear winding and the toroidal magnetic field to approach the dynamical limit, $B_{\phi}\sim(4\pi\rho)^{1/2}r\Delta\Omega=2\times 10^{16}\,\rho_{14}^{1/2}(\Delta\Omega/\Omega)(P_{\rm NS}/10~{}{\rm ms})^{-1}$ G, all before the NS cools. A field this strong induces a strong temperature perturbation, which is rapidly erased by charged-current neutrino reactions, thereby generating buoyant motions and feedback on the seed poloidal field (Thompson & Murray, 2001). In this way a dynamo feedback loop becomes possible.

Appendix A MESA inlist

Below we include the inlist parameters we changed from their default values in our MESA runs:

&star_job

create_pre_main_sequence_model = .true.

kappa_file_prefix = ’gs98’

change_v_flag = .true.

new_v_flag = .true.

warn_run_star_extras = .false.

/ ! end of star_job namelist

&controls

initial_mass = 3

initial_Y = 0.25

initial_Z = 2d-2 (6d-3 in the 40$\,M_{\odot}$, $0.3Z_{\odot}$ model)

velocity_logT_lower_bound = 7

max_dt_yrs_for_velocity_logT_lower_bound = 1

mesh_delta_coeff_for_highT = 2

okay_to_reduce_gradT_excess = .true.

cool_wind_RGB_scheme = ’Dutch’

cool_wind_AGB_scheme = ’Dutch’

Dutch_scaling_factor = 1

cool_wind_full_on_T = 1d8

cool_wind_full_off_T = 1.1d8

include_dmu_dt_in_eps_grav = .true.

use_Type2_opacities = .true.

Zbase = 2d-2 ! must set this in the main inlist

mixing_length_alpha = 2

MLT_option = ’Henyey’

use_Ledoux_criterion = .true.

alpha_semiconvection = 0.1

thermohaline_coeff = 2

overshoot_f_*** = 1d-2     ! This applies to all variations of ’overshoot_f’

overshoot_f0_*** = 5d-4    ! This applies to all variations of ’overshoot_f0’

min_timestep_limit = 1d-12

/ ! end of controls namelist

Appendix B Angular Momentum Transport by Winding a Poloidal Magnetic Field

Here we consider how the kinking of a wound-up magnetic field may limit the transport of angular momentum through the radiative layers of a star. The background state has finite magnetic helicity ${\cal H}$ and is threaded by large-scale toroidal and poloidal fields. Each field component carries a finite flux, which was deposited during a transition from a convective to the present radiative state (Section 4). Neither the hemispheric poloidal flux $\Phi_{r}$ nor the toroidal flux $\Phi_{\phi}$ is directly modified by differential rotation, even while the toroidal field energy may increase substantially.

Magnetic twist is stored on loops of poloidal field that close within the radiative layer (Braithwaite & Spruit, 2004). These structures act as barriers to the mixing and reconnection of open poloidal fluxes from the opposing hemispheres. In the application to evolving stellar interiors, we posit an initial relaxation to a roughly isotropic state with $\Phi_{r}\sim\Phi_{\phi}\sim{\cal H}^{1/2}$. This maximizes the poloidal flux for a given ${\cal H}$, as well as the rate of redistribution of angular momentum by poloidal torsional magnetic waves.

Winding by differential rotation leaves the magnetic field susceptible to an ideal hydromagnetic ‘kink’ instability (Tayler, 1973). In contrast with the model analysed by Spruit (2002), in which the mean poloidal field is absent, the kink is not the primary source of poloidal flux: it only induces high-wavenumber distortions of the poloidal field. We therefore simplify the problem by assuming a fixed poloidal field, and consider the action of the kink instability only on the wound-up component of the toroidal field.

Three characteristic timescales can be distinguished here: i) the evolution time of the stellar mass profile, $\tau_{\rm ev}={\rm min}[l_{P}/|v_{r}|,(t_{\rm cc}-t)/3]$, which varies significantly between the inner and outer parts of the star during post-MS evolution; ii) the growth time $\tau_{\rm kink}$ of the kink instability; and iii) the timescale $\tau_{J}$ for the redistribution of angular momentum between different shells within a slowly rotating radiative layer. (Convective angular momentum transport is generally much faster, taking a few eddy overturns.) We are interested in deducing the minimum $\Phi_{r}$ that will allow nearly solid rotation to be established, corresponding to $\tau_{J}\sim\tau_{\rm ev}$.

To do this, we first must address the growth of the kink. Two inequalities which can be justified ex post facto are i) $\tau_{\rm kink}<\tau_{\rm ev}$ unless the rotation is extremely slow ($\Omega\lesssim\tau_{\rm ev}^{-1}$), meaning that the growth of the toroidal field will be limited by kinking; and ii) $\tau_{\rm kink}>\Omega^{-1}$, so that the Coriolis force must be taken into account in evaluating $\tau_{\rm kink}$. Then the kink is associated with a hydromagnetic displacement of speed (Pitts & Tayler, 1985)

[TABLE]

To obtain a relationship between $B_{\phi}$ and $B_{r}$ we use the induction equation,

[TABLE]

Here we adopt a simplified analysis which focuses on the zone near the rotational equator. Then when $\tau_{\rm kink}<\tau_{\rm ev}$, one has

[TABLE]

Here $v_{A,\phi(r)}=B_{\phi(r)}/(4\pi\rho)^{1/2}$ are the toroidal (poloidal) Alfvén speeds.

The timescale for angular momentum transport is obtained by balancing the Maxwell torque exerted through a lever arm $\sim r$ against the change in angular momentum,

[TABLE]

Rotational equilibrium corresponds to $\tau_{J}\sim\tau_{\rm ev}$, which we substitute into Equation (B4) along with Equation (B3) to get an expression for $|d\ln\Omega/d\ln r|$. Requiring that the angular velocity gradient be weak, $|d\ln\Omega/d\ln r|\lesssim 1$, we obtain

[TABLE]

One can also substitute Equation (B3) into (B1) and show that

[TABLE]

The inequality is obtained from the threshold (B5) for nearly solid rotation.

One sees that the kink is excited self-consistently when this threshold is reached, unless the rotation of the star is extremely slow, or the evolution time is very short, corresponding to $\Omega\lesssim\tau_{\rm ev}^{-1}$. In this second regime, one can alternatively substitute $B_{\phi}=-B_{r}|d\ln\Omega/d\ln r|\Omega\tau_{\rm ev}$ into the torque formula, and obtain the threshold condition for angular momentum transport,

[TABLE]

A more general condition is obtained by taking the minimum of the right-hand sides of Equations (B5) and (B7). The main point to take away from this analysis is that the kink only causes a slight reduction in the Maxwell stress: the threshold value of $v_{A,r}$ is lengthened compared with $\sim l_{P}/\tau_{\rm ev}$ only by a factor $\sim(r\Omega\tau_{\rm ev}/2l_{P})^{1/4}$.