A Theory for Neutron Star and Black Hole Kicks and Induced Spins

Radio pulsars (Manchester et al. 2005; Kaspi & Kramer 2016), magnetars (Kouveliotou et al. 1998; Kaspi & Beloborodov 2017), and low-mass and high-mass X-ray binaries (LMXBs and HMXBs; Shakura & Sunyaev 1973; Remillard & McClintock 2006) are or contain neutron stars or black holes (MacLeod & Grindlay 2023) created at the death of a massive star, often in a core-collapse supernova explosion. Excitingly, the mergers of tight binaries of stellar-mass black holes and neutron stars have recently been captured as gravitational wave sources (Abbott et al. 2016). Interestingly, a broad spectrum of observations indicates that the formation of neutron stars often leaves them with significant kick speeds that range up to ∼1500 km s⁻¹, with an average kick speed of ∼350–400 km s⁻¹ (Faucher-Giguère & Kaspi 2006; Chatterjee et al. 2009) and possible spin–kick correlations (Ng & Romani 2007; Holland-Ashford et al. 2017; Katsuda et al. 2018). However, the bound neutron star population in globular clusters suggests that some neutron stars are born with kick speeds below ∼50 km s⁻¹ (Lyne & Lorimer 1994; Arzoumanian et al. 2002), consistent with the suggestion that the kick distribution could be bimodal (Cordes & Chernoff 1998; Arzoumanian et al. 2002; Verbunt et al. 2017). In contradistinction, studies of black holes in X-ray binaries imply average kick speeds for them of less than ∼50 km s⁻¹ (Fragos et al. 2009; MacLeod & Grindlay 2023; Vigna-Gómez et al. 2023), but with a few possible exceptions (Fragos et al. 2009; Mandel 2016; Repetto et al. 2017; Atri et al. 2019; Dashwood Brown et al. 2024).

These data and inferences require explanation in the context of compact object birth, i.e., in the context of core-collapse supernova explosions (CCSNe). The Blaauw mechanism due to the unbinding of progenitor binaries as a result of explosive mass loss (Blaauw 1961; Hills 1970; Hoogerwerf et al. 2001) is quantitatively inadequate, though often subdominantly operative. Models that invoke asymmetrical radio pulsar winds (Harrison & Tademaru 1975) are not compelling (Lai et al. 2001). Models that rely on neutrino jets in the context of rapid rotation at birth (Spruit & Phinney 1998) are not consistent with modern CCSN simulations (Janka 2017; Burrows et al. 2020; Stockinger et al. 2020; Burrows & Vartanyan 2021; Coleman & Burrows 2022). The most comprehensive and well-developed theory for pulsar and black hole kicks at birth involves recoil due to the generically asymmetrical supernova explosion (Burrows & Hayes 1996; Scheck et al. 2006; Burrows et al. 2007; Nordhaus et al. 2010, 2012; Wongwathanarat et al. 2013; Janka 2017; Gessner & Janka 2018; Müller et al. 2019; Nakamura et al. 2019; Coleman & Burrows 2022) and to the accompanying aspherical neutrino emissions (Fryer & Kusenko 2006; Nagakura et al. 2019b; Rahman et al. 2022), the latter generally being subdominant for neutron star birth but important for black hole birth (Coleman & Burrows 2022).

Asphericities arise in various regions of the supernova progenitor star, and the timescales of their effects can vary from a few seconds to hours or days. Simultaneous accretion and explosion leads to interactions between ejecta and accreta, and this phenomenon is the first major physical process that influences the kick velocity experienced by the central object. When accretion terminates and the explosive shock wave has swept away most of the outer envelope a few seconds later, this process ceases, and the core has achieved its final asymptotic kick state. This marks the end of kick evolution if the explosion is sufficiently energetic, i.e., if the total energy of the outgoing matter is higher than the binding energy of the outer envelope. For less energetic explosions, part of the ejected matter may slow down due to interactions with stellar envelopes and may fall back onto the central object (Schrøder et al. 2018; Chan et al. 2020). This second process may last hours to days. Even if the progenitor fails to explode and there is no asymmetric ejecta, the asymmetric neutrino emission can still provide a modest kick (Fryer & Kusenko 2006; Nagakura et al. 2019b; Rahman et al. 2022). In this case, the accretion of the outer envelope can carry angular momenta to the compact core (Antoni & Quataert 2023).

In this paper, we focus on the first process mentioned above, i.e., kicks during the simultaneous accretion and explosion phase. This is the dominant process in energetic CCSN explosions (Chan et al. 2020; Coleman & Burrows 2022). Many relevant simulations have been done in the past (Burrows & Hayes 1996; Scheck et al. 2006; Burrows et al. 2007; Nordhaus et al. 2010, 2012; Wongwathanarat et al. 2013; Janka 2017; Gessner & Janka 2018; Müller et al. 2019; Nakamura et al. 2019; Coleman & Burrows 2022), but the required sophisticated and expensive 3D simulations have not been continued to late enough times after bounce to achieve asymptotic kick speeds, with very few exceptions (Stockinger et al. 2020; Coleman & Burrows 2022). Most 3D simulations are stopped within one second or less after the bounce, far too soon to witness the final kick speed imparted. Moreover, even the exceptional longer-term 3D studies did not explore the systematics with the broad range of progenitor masses to determine the correlations between kick speed and progenitor mass and/or initial core structure, the latter frequently associated with the compactness parameter. ⁵ They focused on the few models whose recoil kick had quickly settled to nearly a final state. In contrast to papers incorporating detailed numerical simulation, papers such as Janka (2017), Mandel & Müller (2020), Bray & Eldridge (2016), and Bray & Eldridge (2018) have derived semianalytical and/or empirical relations between kick velocities and supernova/progenitor properties. Although such simple recipes are intuitively appealing, they need to be tested thoroughly using 3D numerical simulations.

Hence, a comprehensive theoretical study of the broad spectrum of systematic behaviors in the context of modern 3D core-collapse supernova theory has not yet been performed. With this paper, we attempt to remedy this situation by providing and analyzing 20 of our recent long-term state-of-the-art 3D Fornax (Skinner et al. 2019) simulations for initially nonrotating zero-age main-sequence star (ZAMS) solar-metallicity progenitors from 9 to 60 M_⊙ taken from Sukhbold et al. (2016) and Sukhbold et al. (2018). Our goal is to build a more solid foundation for the current CCSN kick theory and to provide detailed 3D models that can be used to test, calibrate, and update widely used recipes.

We emphasize that we are addressing in this paper only solar-metallicity progenitors and those most likely to be the origin of most galactic neutron stars and black holes and leave subsolar metallicity studies to future work. We also are not here to address the issue of the potential role of very massive massive stars and pulsational pair instabilities. ⁶ Importantly, many of our models have been carried in 3D beyond 4 s after the bounce and show signs of asymptoting to their final recoil kick speeds. A few of these models have or will form black holes.

In this paper, we complement our kick study with a study of the associated induced spins (Blondin & Mezzacappa 2007; Rantsiou et al. 2011; Coleman & Burrows 2022). Though a compelling explanation for pulsar and neutron star spins is spun up from an initially spinning Chandrasekhar core upon collapse, modern supernova simulations leave protoneutron star (PNS) cores spinning due to random and stochastic postexplosion accretion of angular momentum (L), even when the core is initially nonrotating. In fact, the magnitude of such induced spins (Coleman & Burrows 2022) can be comparable to measured pulsar spins (Chevalier & Emmering 1986; Lyne & Lorimer 1994; Manchester et al. 2005; Faucher-Giguère & Kaspi 2006; Popov & Turolla 2012; Igoshev & Popov 2013; Noutsos et al. 2013). Moreover, there are hints that radio pulsar spins and kicks might be correlated at birth (Johnston et al. 2005; Wang et al. 2006; Johnston et al. 2007; Ng & Romani 2007; Wang et al. 2007; Noutsos et al. 2013). Finally, angular momentum transport from the stellar core to the mantle up to the time of core collapse might perhaps be more efficient than hitherto thought (Cantiello et al. 2014), leaving the Chandrasekhar core slowly rotating just before implosion. Therefore, we continue exploring the possibility that pulsar spins might be partially reflective of induced spins during the supernova. Similarly, black holes created in supernova explosions (Burrows et al. 2023) will have both large kicks and significant induced spin parameters (a = Lc/GM²), though black holes created in a nonexploding context should have birth spins reflective of the initial progenitor spin (though perhaps in a complicated fashion). However, similar to the case with kicks, the induced spins of the central object can be influenced by different very long-term processes. In this paper, however, we focus on the crucial first few seconds postbounce.

We used the code Fornax (Burrows et al. 2019, 2020, 2023; Nagakura et al. 2019a, 2019b, 2020; Radice et al. 2019; Skinner et al. 2019; Vartanyan et al. 2019, 2022, 2023; Vartanyan & Burrows 2020) to generate the 20 3D models to late times after bounce used here to study recoil kicks and induced spins in the context both of supernova explosions and quiescent black hole formation. A more comprehensive set of studies using these new 3D simulations is in preparation, which will address other aspects and outcomes of core collapse and their dependence upon progenitor. This is the largest set of long-term (many seconds after the bounce) 3D state-of-the-art core-collapse simulations ever created, and we hope to use it to gain new insights into the core-collapse supernova phenomenon. To create this set, we employed the CPU machines TACC/Frontera (Stanzione et al. 2020) and ALCF/Theta and the GPU machines ALCF/Polaris and NERSC/Perlmutter. The specific setups and run parameters are described in Burrows et al. (2023) and Vartanyan et al. (2023). The techniques employed to derive the kicks and induced spins are described in detail in Coleman & Burrows (2022) and follow the standard approaches in the literature (Scheck et al. 2006; Wongwathanarat et al. 2013; Stockinger et al. 2020). However, for the calculation of the neutrino component, unlike in the Coleman & Burrows (2022) paper, we calculate its contribution here at the PNS radius (defined as the 10¹¹ g cm⁻³ isosurface), and not at the outer mass surface of the final protoneutron star (PNS). The result for the final total vector kick is the same, but there is a different partitioning between the neutrino and matter components. The momentum transfer between the matter and the radiation in the region between these two definitions accounts for the shift, and the total momentum is conserved under both approaches. We employ a different method here to comport with the standard practice in the literature (Stockinger et al. 2020).

Until recently, 3D models experiencing kicks had been calculated to postbounce times insufficient to actually witness their final values. Moreover, due to computational limitations, the systematics with progenitor core structure and total ZAMS mass were obscured, since only one or a few models early during the accumulating recoil could be simulated. This is akin to viewing the entire progenitor panorama through a straw. The stochasticity of the final values of any explosion observable inherent in the chaotic turbulence of the CCSN phenomenon also served to compromise extrapolation and interpretation. However, with our new set of 20 long-term 3D full-physics simulations, we are now able to extract what we suggest are real trends across the progenitor continuum and find correlations among observables hitherto only imperfectly glimpsed. Hence, though for all core-collapse supernova observables, there is expected to be a distribution of values, even for the same progenitor, the availability of a raft of such models across the progenitor continuum mitigates in part against this modest impediment to determining system-wide trends. This is the philosophy of our multimodel study, enabled by the efficiency of Fornax and the availability of high-performance computing resources.

Table 1 provides the summary numbers relevant to this study. The top rank displays the data for the exploding models that leave neutron stars (NS), while the bottom rank displays the corresponding data for the models that leave black holes (BH; see Section 2.1), either quiescently or explosively. Most of the models have either asymptoted to the final state (in particular the lowest-mass, lowest-compactness progenitors) or nearly so. The left-hand side of Figure 1 depicts the temporal evolution of the absolute value of the vector kick, including both the matter and neutrino recoils and all momentum flux, pressure, and gravitational forces (Coleman & Burrows 2022). The lowest-mass progenitors generally flatten/asymptote earliest, while the models with the largest compactness (frequently more massive) take longer. Such long-duration 3D simulations have not been available in the past.

Zoom In Zoom Out Reset image size

Table 1. Results for the Simulations Presented in This Work

MZAMS	ξ1.75	${t}_{\max }$	Type	Mbaryonic	Mgrav	pkick/1040	${v}_{\mathrm{kick}}^{\mathrm{total}}$	${v}_{\mathrm{kick}}^{\nu }$	L/1046	Pcold	$\hat{L}\cdot \hat{v}$
(M⊙)		(s)		(M⊙)	(M⊙)	(g cm s−1)	(km s−1)	(km s−1)	(g cm2 s−1)	(ms)
9(a)	6.7 × 10−5	1.775	NS	1.347	1.237	3.235	120.7	87.4	1.046	749.1	0.402
9(b)	6.7 × 10−5	2.139	NS	1.348	1.238	2.106	78.6	59.1	0.190	4132	−0.564
9.25	2.5 × 10−3	3.532	NS	1.378	1.263	3.842	140.1	93.2	0.344	2346	0.742
9.5	8.5 × 10−3	2.375	NS	1.397	1.278	5.796	208.6	77.5	1.636	501.5	−0.756
11	0.12	4.492	NS	1.497	1.361	20.83	699.4	90.5	23.28	38.54	0.266
15.01	0.29	4.384	NS	1.638	1.474	5.668	173.9	10.6	6.995	144.2	−0.743
16	0.35	4.184	NS	1.678	1.505	15.62	468.0	63.5	20.23	51.45	0.566
17	0.74	4.473	NS	2.044	1.785	30.90	759.8	61.6	91.52	14.9	0.914
18	0.37	6.778	NS	1.684	1.510	22.97	686.0	81.4	10.23	102.2	−0.221
18.5	0.80	5.250	NS	2.139	1.854	32.43	762.4	57.7	16.88	86.14	−0.184
19	0.48	4.075	NS	1.803	1.603	18.89	526.7	30.6	68.05	16.85	0.360
20	0.79	5.503	NS	2.143	1.857	24.21	567.9	6.08	16.17	90.22	0.568
23	0.74	6.228	NS	1.959	1.722	10.94	280.8	75.9	9.294	138.2	0.825
24	0.77	3.919	NS	2.028	1.773	36.79	912.1	34.9	37.58	35.88	−0.990
25	0.80	5.611	NS	2.106	1.830	26.52	633.0	55.2	34.38	41.37	0.979
60	0.44	6.386	NS	1.791	1.594	30.81	864.8	46.0	119.1	9.537	−0.380
MZAMS	ξ1.75	${t}_{\max }$	Type	Mbaryonic	Mfinal	pkick/1040	${v}_{\mathrm{kick}}^{\mathrm{total}}$	${v}_{\mathrm{kick}}^{\nu }$	L/1046	a	$\hat{L}\cdot \hat{v}$
(M⊙)		(s)		(M⊙)	(M⊙)	(g cm s−1)	(km s−1)	(km s−1)	(g cm2 s−1)
12.25	0.34	2.090	BH (F)	1.815	∼11.1	1.441	39.9 (10.7)	39.9 (6.90)	3.621	0.002 (0.0004)	−0.677
14	0.48	2.824	BH (F)	1.990	∼12.1	1.684	42.5 (9.59)	42.5 (7.9)	4.638	0.003 (0.0004)	−0.920
19.56	0.85	3.890	BH (S)	2.278	?	82.76	1826 (?)	28.1 (?)	12.83	0.006 (?)	0.966
40	0.87	1.76 (21)	BH (S)	2.434	∼3.5	72.85	1504(1034)	79.5 (83.0)	84.32	0.039 (0.6)	−0.883

Note. M_ZAMS is the initial ZAMS mass of the progenitor model. ${t}_{\max }$ is the maximum postbounce time reached by the simulation. For the black hole formers, the times are the time of black hole formation, and the second time (if shown) is the latest time after black hole formation; we have carried that model hydrodynamically. The types of central objects are summarized. NS means neutron star, BH (F) means black hole formed with a failed explosion, and BH (S) means black hole formed accompanied by a successful CCSN explosion. M_baryonic is the baryonic mass of the central object at ${t}_{\max }$, while M_grav is the gravitational mass of a cold, catalyzed NS. M_final is the estimated final black hole mass. The final masses of the nonexploding black hole formers are the total masses in the progenitor models at collapse, assuming all matter will eventually be accreted, while the final masses of the exploding black hole formers are derived from our hydro-only blast wave simulations after black hole formation. p_kick gives the linear momentum of the object at ${t}_{\max }$. ${v}_{\mathrm{kick}}^{(\mathrm{total})}$ and ${v}_{\mathrm{kick}}^{(\nu )}$ are the total kick speeds and the kick speeds induced by neutrinos only, respectively. Both are measured at ${t}_{\max }$. Numbers in brackets show the estimated final values (kick speed, spin parameter) of the black holes. Note that the 40 M_⊙ model has been simulated to 21 s after bounce, at which point the spin parameter has reached ∼0.35, and that we estimate its final gravitational mass and spin parameter (a) will be ∼3.5 M_⊙ 0.6, respectively, achieved after a further ∼25 s of evolution. We have yet to determine the corresponding numbers for the 19.56 M_⊙ model. These numbers still need to be checked. L_central is the angular momentum of the central object at ${t}_{\max }$, and P_cold is the inferred final spin period of the cold final neutron star, assuming $P=\tfrac{2\pi {L}_{\mathrm{central}}}{{I}_{\mathrm{NS}}}$. The neutron star moment of inertia is estimated using ${I}_{\mathrm{NS}}=(0.237+0.674{\xi }_{\mathrm{NS}}+4.48{\xi }_{\mathrm{NS}}^{4}){{MR}}_{\mathrm{NS}}^{2}$, where ${\xi }_{\mathrm{NS}}=\tfrac{{GM}}{{R}_{\mathrm{NS}}{c}^{2}}$ is the compactness parameter of the neutron star itself (Breu & Rezzolla 2016). We assume here that the final neutron star radius is R_NS = 12 km. For black holes, the relativistic spin parameter $a=\tfrac{{Lc}}{{{GM}}^{2}}$ is shown, and the number in brackets is our guess for its final value. The final column contains the cosine of the spin/kick relative angle at the end of each simulation.

Download table as:  ASCIITypeset image

The right-hand side of Figure 1 depicts the final kicks versus the final gravitational mass of the corresponding neutron stars. This is a plot of observables that have not been available before and speaks to the importance of long-term 3D simulations to provide testable predictions for modern CCSN theory. We note that the lowest-mass progenitors (blue) are all clustered at the low end of the kick range. This is correlated with the greater degree of sphericity of models that explode early before the postshock turbulence can grow to a significant enough strength and Mach number to manifest the growth of large bubbles that, when explosion ensues, lead it and set an explosion direction. For these initially nonrotating progenitors, the direction of the explosion is random, but once determined, a crude axis is set. This is not the case for the lowest-compactness progenitors that explode early and, hence, more spherically.

Figure 2 depicts this dichotomy. Its left-hand side portrays the blast as traced by the 10% iso-⁵⁶Ni surface of the explosion of the 9 M_⊙ progenitor, colored by Mach number. Purple is positive (outward), and blue is inward (infall). The shock wave is the spherical bluish veil around the ⁵⁶Ni bubbles. In contrast, the right-hand side of Figure 2 depicts the corresponding plot for the exploding 23 M_⊙ model. ⁷ This model has been carried to ∼6.2 s after the bounce and is only slowly asymptoting, while the 9 M_⊙ progenitor flattened within 2 s (Coleman & Burrows 2022). The blast structure of the 23 M_⊙ model is more generic, and this is reflected in the separation of the blue dots on the right-hand side of Figure 1 from the rest. The latter are spread between ∼300 and ∼1000 km s⁻¹ and the former are found near ∼100−200 km s⁻¹; the difference reflects the dichotomy illustrated in Figure 2. ⁸

Zoom In Zoom Out Reset image size

Binary neutron stars are typically found to be, on average, on the low-mass side (Fan et al. 2023), and they seem to have a distinct population. One explanation is that the lower-mass neutron stars are birthed (on average) with lower kick speeds so that such binary systems can survive. This correlation emerges naturally, without tweaking, from our simulations and is consistent with observations (Tauris et al. 2017). Such a correlation was anticipated by Müller et al. (2019) and explored empirically by Bray & Eldridge (2016) and Bray & Eldridge (2018). The right-hand side of Figure 1 confirms this correlation, and also suggests that there is a neutron star mass threshold, roughly around 1.35 M_⊙, above which the kick velocities show much larger variations. It is also the case that even lower-mass progenitors not in our sample (such as the 8.8 M_⊙ progenitor of Nomoto (1984) and accretion-induced collapse (AIC) systems) would naturally experience small kicks, and these systems could play a role. Such progenitors are also likely to explode quickly and experience only very modest kicks (though if rapidly rotating could experience asymmetric MHD jets that could result in interesting recoils). The lowest progenitor mass range and AIC systems might be relevant as sources for the neutron stars bound in globular clusters, for which the kick speeds need to be smaller than the potential well depths (≤50 km s⁻¹). Though we believe the qualitative picture is now clear, we are still some distance from a fully quantitative explanation of the binary neutron star conundrum, as suggested in population synthesis studies. More 3D numerical simulations are required to generate a complete quantitative theory that can fruitfully be compared to the observed mass–kick distribution of neutron stars. Moreover, we emphasize that the progenitor models we use in this work are all single-star models (Sukhbold et al. 2016, 2018), and that binary effects may significantly influence progenitor structures and, hence, kick velocities (Podsiadlowski et al. 2004).

The left-hand side of Figure 3 depicts the temporal evolution of the corresponding neutrino contribution to the kick. This is its absolute value, and in the total kick provided in Figure 1 the matter and neutrino terms are vectorially added; the two are not generally aligned, and their individual instantaneous directions vary with time. We see that the lowest ZAMS mass progenitors (such as the 9 M_⊙ model) have some of the highest neutrino kicks (Coleman & Burrows 2022), but that most other models at higher compactness whose explosions are more delayed have lower neutrino kick contributions. Moreover, since the matter contributions to the kicks for the more delayed explosions are generally larger, the neutrino contributions are proportionately less for them and can be as low as ∼3%. They are, on the other hand, as high as ∼80% for the lowest-mass/compactness progenitors. We emphasize, however, that the chaos of the CCSN phenomenon can result in a range of values for all observables and that this range has not yet been well determined.

Zoom In Zoom Out Reset image size

On the right-hand side of Figure 3 we provide a plot of the total kick magnitude derived versus the mass dipole of the ejecta. The latter is a measure of the explosion asymmetry. We see that the two are roughly correlated, as expected for the general recoil model of kicks, but that there is a bit of scatter. In Figure 4, we provide a similar plot of the kick magnitude versus the explosion energy. The rough trends of kick with explosion energy and with ejecta mass dipole are clear and generally expected (Burrows et al. 2007; Janka 2017).

Zoom In Zoom Out Reset image size

2.1. Black Hole Kicks

Though the rotation profiles of the cores of massive stars at stellar collapse are not known, we know that massive star surfaces rotate. The loss of angular momentum to winds is guaranteed, but its magnitude is a function of surface B-fields, which themselves are less well known. Furthermore, the redistribution of angular momentum in the stellar interior is no doubt via magnetic torques, but the processes and couplings are not firmly constrained. The fact that red giant interiors are rotating more slowly (Cantiello et al. 2014) than inferred using a Taylor–Spruit dynamo (Spruit 2002) suggests a more efficient spindown process than employed by Heger et al. (2005) in estimating the spin at the collapse of massive-star models. The latter found using a Taylor–Spruit prescription for angular momentum redistribution that the Chandrasekhar cores might have initial periods at collapse near ∼30–60 s, which would translate by angular momentum conservation into neutron star periods of ∼30–60 ms (Ott et al. 2006), with the lowest-mass massive stars spinning the slowest. Hence, since these periods are not much different from those expected for radio pulsars (Chevalier & Emmering 1986; Lyne & Lorimer 1994; Manchester et al. 2005; Faucher-Giguère & Kaspi 2006; Popov & Turolla 2012; Igoshev & Popov 2013; Noutsos et al. 2013), it would seem that the issue of the initial spin period of the cores of massive stars is resolved, at least in broad outline.

However, it is still possible that many of the cores of massive stars destined to undergo collapse may be torqued down significantly. This thought is motivated not only by the results of Cantiello et al. (2014), but by the observation that white dwarfs are born rotating very slowly (Kawaler 2015; Fuller et al. 2019). Could this be the fate of a good fraction of the white dwarf cores of massive stars when they die? In order to explore this possibility, one first needs theoretical insight into what induced spins are possible in the context of core-collapse supernovae. This is what we provide in this paper.

The left-hand side of Figure 5 renders the evolution of the angular momentum imparted to the PNS cores during their postbounce evolution. This spinup is due to the stochastic accretion of plumes of infalling matter into the PNS, with varying impact parameters. If the explosion is very asymmetrical, as is the case when the kick speeds are significant, or black hole formation is accompanied by an explosion, the consequently asymmetrical accretion during the explosion (mostly then from one side) can spin the PNS up to interesting periods. For the black holes that form in an explosion, these spins can reach high values (see Table 1), even though their cores were initially not rotating. This implies that such objects could be collapsar candidates (Heger & Woosley 2003), despite their lack of initial spin and their solar-metallicity progenitors.

Zoom In Zoom Out Reset image size

The right-hand side of Figure 5 portrays the dependence of the final induced spin period of the neutron star upon ZAMS mass. We see a trend, roughly recapitulated with the kick systematics (see Table 1), of decreasing period with increasing ZAMS mass, but with a lot of scatter. Such scatter should be expected in the context of the chaotic turbulence seen in 3D simulations of CCSNe. We emphasize that the lowest ZAMS mass and compactness progenitors, which explode early and approximately spherically, have the slowest induced spin rates, while those exploding later and at higher compactness have greater induced spin rates. In fact, our highest compactness models can leave cold neutron stars with high angular momenta and spins. However, while we are comfortable with the compactness/kick and PNS mass/kick correlations we derive, we are much less sanguine that induced spins are the whole story of compact object birth spins. Importantly, we know that the Crab explosion involved a low-mass progenitor near ∼9−10 M_⊙ (Stockinger et al. 2020), but that its birth spin was ∼15 ms (Lyne et al. 2015). This is far from the induced periods (many seconds) we find for this mass range of progenitors, not only suggesting but requiring the story to be more complicated. Spinup by binary interaction is also distinctly possible (Fuller & Lu 2022). Therefore, the relative roles of induced and initial spins in the distribution of observed birth periods for compact objects remains open.

Nevertheless, we find a weak but intriguing correlation between induced spin and recoil kicks, demonstrated in Figure 6. The greater the kick, the greater the induced spin, again with significant scatter. Is there something of importance in this apparent correlation? Is there rough support for this, even if complicated by the overlapping effect of the initial spin, in the pulsar database?

Zoom In Zoom Out Reset image size

The induced spins for the nonexploding black hole channel (see Table 1) are quite low, again resulting solely from the torques due to the radiated neutrinos. Spin parameters of close to zero are inferred from our simulations and model extrapolations. However, should there be any significant angular momentum in the accreting stellar envelope that, on long timescales, will bulk up the core (something quite likely), these numbers can be radically changed. Hence, even under the assumption that the initial core spin rate is zero, we can not say with any confidence what the final spin parameter of black holes born through this channel might be.

However, the induced spins for the exploding black hole channel can be large. This can follow from the significant asymmetry in their explosions, itself related to the large recoil kick. A spin parameter as large as ∼0.6 emerges for the 40 M_⊙ model (see Table 1), though we have yet to determine the corresponding number for the 19.56 M_⊙ model. Their large compactnesses and related large postexplosion asymmetrical mass accretion rates are direct causative factors. In addition, accretion continues to power the driving neutrino emissions, which for this black hole formation channel are quite large. Accretion settles into biased accretion onto one general side and spins up the periphery of the PNS before black hole formation. After black hole formation, continued accretion parallels continued spinup until, after many seconds, the black hole is isolated (for the 40 M_⊙ model after ∼50 s ?). This general black hole formation scenario is discussed in Burrows et al. (2023).

Figure 2 portrays the Mach number for ⁵⁶Ni surfaces and the crudely axial explosions when the compactness is not small and when there is a greater delay to the explosion. It also indicates that there is continued accretion roughly in the perpendicular plane. Hence, when the compactnesses and ZAMS masses are not small, explosion and accretion are roughly perpendicular. Accretion then spins up the periphery of the PNS, and the net spin vector is either roughly parallel or antiparallel to the kick (which is along or antiparallel to the explosion axis). Hence, we would think a spin–kick correlation (Ng & Romani 2007; Holland-Ashford et al. 2017; Katsuda et al. 2018) would naturally emerge. Figure 7 depicts the dot product of the final spin and kick unit vectors for our model set versus ZAMS mass (left) and compactness (right). We see that though there seems to be a deficit in the perpendicular direction and there is a slight preference for an (anti-)aligned orientation, the null hypothesis can not be rejected. Indeed, more rigorously, we perform a Kolmogorov–Smirnov (K-S) test with the null hypothesis that the simulated angles are drawn from a uniform distribution on the sphere. We see a D-statistic of 0.19, which means the cumulative distributions are not very different. In addition, due to the small number of models, the p-value of this K-S test is 0.51, which means that there is a chance that the deficit in the perpendicular direction is simply due to statistical fluctuations. Therefore, however intriguing, many more simulations are required to properly address a potential spin–kick correlation in the initially nonrotating case.

Zoom In Zoom Out Reset image size

Nevertheless, such a correlation may arise when the Chandrasekhar core is initially rotating, since such rotation may slightly bias (depending upon its magnitude) the axis and ejecta mass dipole vector of the subsequent neutrino-driven explosion. If such is the case, then the spin–kick correlation may emerge. Continued accretion is very roughly perpendicular to the explosion direction, itself roughly (anti-)parallel to the kick direction, resulting in a (anti-)correlation between spin and kick. Note that we have not invoked B-fields nor MHD jets and that this is an expected outcome in the neutrino-driven explosion context. However, we emphasize that we have yet to perform the requisite rotating simulations.

With a focus on the recoil kicks imparted to the residues of stellar core collapse in the supernova and black hole formation contexts, we have analyzed 20 state-of-the-art 3D long-term core-collapse simulations generated using the code Fornax. In the past, most 3D models were not of sufficient duration to witness the cessation of net acceleration and the asymptoting of the kicks to their final values. This includes our own previous work (Coleman & Burrows 2022), as well as that of other workers in the field. However, and for the first time for a large and uniform collection of 3D supernova models, we have either asymptoted the kicks or come within 20% of doing so. As a result, we obtain an integrated and wide-angle perspective of the overall dependence of the recoil kicks and induced spins upon progenitor mass and their Chandrasekhar-like core structures, the latter indexed approximately by compactness. As had been postulated in the original recoil kick model (Janka & Mueller 1994; Burrows & Hayes 1996), the recoil is always opposite in direction to the bulk of the ejecta.

We find that lower-mass and lower compactness progenitors that explode quasi-spherically due to the short delay to explosion experience, on average, smaller kicks in the ∼100−200 km s⁻¹ range, while higher mass and higher compactness progenitors that explode later and more aspherically give birth to neutron stars with kicks in the ∼300−1000 km s⁻¹ range. We also find that these two classes can be correlated with the gravitational mass of the residual neutron star. The correlation of a lower kick speed with lower neutron star gravitational mass suggests that the survival of binary neutron star systems, for which their mean mass is lower than average, may in part be due to the lower kick speeds we see for them. With some scatter, there is a correlation of the kick with both the mass dipole of the ejecta and the explosion energy. We also find, as did Coleman & Burrows (2022), that the contribution of the recoil due to anisotropic neutrino emission can be important for the lowest-mass/compactness progenitors.

However, unlike Coleman & Burrows (2022), we partition the fractional contribution of the neutrino component to the overall kick differently, which results in an apparent diminution of its contribution for higher ZAMS mass and compactness progenitors, for which the matter recoil effect then predominates. The difference with Coleman & Burrows (2022) is a consequence of our use here of the standard literature radius at which to calculate the neutrino contribution and not due to some major differences with the overall results in Coleman & Burrows (2022). Between the effective radius at which the neutrino component is calculated in Coleman & Burrows (2022) and here, there is momentum transfer between the matter and neutrino radiation that is properly taken into account in the vector sum of the two provided in both papers.

We find two channels for black hole birth. The first leaves masses of ∼10 M_⊙, is not accompanied by a neutrino-driven explosion, and experiences small kicks, perhaps below ∼10 km s⁻¹. The second (Burrows et al. 2023) is through the vigorous, neutrino-driven explosion of a high-compactness progenitor that leaves behind a black hole with a mass of perhaps ∼2.5−3.5 M_⊙ that is kicked to high speeds, perhaps above ∼1000 km s⁻¹. This second black hole formation channel is novel and unexpected, but if true, it challenges most notions of what can explode. Both channels challenge prior notions of the mapping between progenitor and outcome and must still be verified. However, in the context of our CCSN simulations, the associated mass ranges for black hole formation depend entirely upon the as-yet unconverged progenitor suite model. In particular, the ZAMS mass range for our quiescent black hole formation channel could shift significantly when binarity, overshoot, shell-merger, mixing processes, and the ¹²C(α,γ) rate are properly understood.

Associated with kicks are induced spins that accompany them in the context of the modern theory of turbulence-aided, neutrino-driven supernova explosions. We find that the induced spins of nascent neutron stars range from a few seconds at low compactness and low progenitor mass to ∼10 ms at higher compactness and higher ZAMS mass. However, this can not be the whole story of the birth spins of neutron stars, and the relative roles of initial and induced spins have yet to be determined. Nevertheless, we find that it is, in principle, possible to explain the spin rates of some pulsars with the postbounce asymmetrical accretion of net angular momentum in the context of asymmetrical supernova explosions for which spherical symmetry is generically broken. Moreover, we find that a spin/kick correlation for pulsars is not yet suggested by the theory. However, if there is an initial spin to the collapsing core, its vector direction could bias the direction of the explosion. If this is the case, a spin/kick correlation may emerge. However, this speculation has yet to be verified, but if true, a spin/kick correlation (in a statistical sense) would naturally emerge as a byproduct of the turbulence-aided neutrino mechanism of core-collapse supernova explosions.

In the nonexploding case of black hole formation, the induced spins are quite low due to anisotropic neutrino emission. Any slight neutrino torques would be due to a random stochastic residual effect. Under these circumstances, an initial spin would be determinative. However, for the exploding context of black hole formation, the induced spin is likely to be large. For our 40 M_⊙ progenitor, we find that the spin parameter could reach ∼0.6 (the maximum is 1.0). This is in the collapsar range and suggests that even an initially nonrotating, solar-metallicity, high-compactness progenitor could be a seat of long-soft gamma-ray bursts with a vigorous explosive precursor. This is a controversial claim, and one that needs much further scrutiny. However, this is what is suggested by our simulations.

It must be remembered that the chaos in the turbulence generic and fundamental to core-collapse phenomena will result in a spread of outcomes, even for a given progenitor. This translates into distribution functions for the explosion energies, nucleosynthesis, morphology, kicks, and induced spins that are not delta functions and are currently unknown. Hence, there is natural scatter in the kicks and induced spins for a given initial model that has yet to be determined, but that will spread observed values. Some of the scatter seen in this study no doubt originates in this behavior. Such are the wages of chaos. Nevertheless, though the study of a single 3D model could distort one's understanding, the study of a wide range of progenitor structures should provide a more complete picture of the theoretical trends and correlations that will minimize some of the confusion due to such scatter. In any case, this is the philosophy of this 3D multimodel, a long-term exploration of the core collapse, but much more remains to be done. In addition, we note that the mapping of observables to progenitor structure, modeled loosely as "compactness," seems more robust than the mapping with ZAMS progenitor mass, tethered as it is to the progenitor suite employed. Given this, the fact that we are in this study using the Sukhbold set of progenitors limits what one can say about the observable/ZAMS mass correlations and the guidance we can provide to population synthesis modelers. The effects of binarity (mass transfer), overshoot, shell mergers, perturbations, rotation, the ¹²C(α,γ) nuclear rate, and wind mass loss have yet to be fully understood, and their influence retired. Moreover, physical processes that might change residual black hole properties on timescales of hours, days, or years, such as "fallback" (Schrøder et al. 2018; Chan et al. 2020) or random late-time angular momentum accretion (Antoni & Quataert 2023), are not covered in this work.

Though Fornax is a sophisticated tool, it incorporates no approach to neutrino oscillations, uses a moment-closure approach to the radiative transfer, and incorporates approximate, not precise, general relativity. Nevertheless, the greatly expanded panorama the new capabilities and this new numerical database provide is qualitatively more comprehensive and predictive, due to its breadth and temporal extent, than has heretofore been possible and appears to explain many observed pulsar properties by default.

We thank Chris White for previous insights, collaboration, and conversations. We acknowledge support from the U.S. Department of Energy Office of Science and the Office of Advanced Scientific Computing Research via the Scientific Discovery through Advanced Computing (SciDAC4) program and grant DE-SC0018297 (subaward 00009650) and support from the U.S. National Science Foundation (NSF) under grants AST-1714267 and PHY-1804048 (the latter via the Max Planck/Princeton Center (MPPC) for Plasma Physics). Some of the models were simulated on the Frontera cluster (under awards AST20020 and AST21003), and this research is part of the Frontera computing project at the Texas Advanced Computing Center (Stanzione et al. 2020). Frontera is made possible by NSF award OAC-1818253. Additionally, a generous award of computer time was provided by the INCITE program, enabling this research to use resources of the Argonne Leadership Computing Facility, a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357. Finally, the authors acknowledge computational resources provided by the high-performance computer center at Princeton University, which is jointly supported by the Princeton Institute for Computational Science and Engineering (PICSciE) and the Princeton University Office of Information Technology, and our continuing allocation at the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the U. S. Department of Energy under contract DE-AC03-76SF00098.

The numerical data underlying this article will be shared upon reasonable request to the corresponding author.