The Effect of Supernova Convection On Neutron Star and Black Hole Masses

Early calculations of compact remnant masses were based solely on our understanding of stellar structure, guided by nucleosynthetic yields and the use of a piston explosion model (Timmes et al. 1996). Many of the recent remnant mass distribution calculations have focused on improving one-dimensional explosion calculations to capture the engine more accurately than the piston model. However, the convective engine behind supernova explosions is an intrinsically multidimensional process and incorporating this physics is difficult in a one-dimensional simulation. The explosion is typically driven by artificially altering the energy deposition in a one-dimensional code. A large fraction of these models altered either the neutrino luminosity or absorption to increase the energy deposited (Fröhlich et al. 2006; Fischer et al. 2010; Ugliano et al. 2012; Perego et al. 2015; Ertl et al. 2016; Prša et al. 2016). These models mimic the enhanced neutrino luminosity and heating from convection within the proto-neutron star (Keil et al. 1996). Other models have implemented one-dimensional mixing above the proto-neutron star (which is believed to be much more important than the proto-stellar mixing) to drive explosions (Fryer et al. 2018; Couch et al. 2020). All of these models are simply recipes to mimic in one dimension the intrinsically multidimensional supernova engine.

To understand the difficulties in producing an accurate explosion in one dimension, we must better understand the current paradigm for the engine behind core-collapse supernovae. Herant et al. (1994) found that convection above the proto-neutron star can both enhance the conversion of internal energy produced in the bounce and reduce the pressure of the infalling material. Since these results, an increasing number of simulations have produced explosions under this paradigm, and it has slowly gained traction to become the leading theory model for the supernova engine (Fryer & Young 2007; Takiwaki et al. 2014; Lentz et al. 2015; Melson et al. 2015; Burrows et al. 2018; Fields & Couch 2021). Given that convection is key to the supernova engine, it would appear that the models including mixing would be more predictive than those that simply focus on the neutrino physics. But, as we shall discover in this paper, different one-dimensional prescriptions for mixing produce different results. In addition, these prescriptions do not capture the large-scale flows seen multidimensional models that will affect the final remnant mass distribution.

Unfortunately, multidimensional models also do not provide a final answer. The growth of convective instabilities in multidimensional models tends to be slower than the growth times predicted by turbulence models (Fryer & Young 2007). The complexity (multi-physics) and scale of the core-collapse supernova engine prevents scientists from conducting fully resolved simulations of this convection. The resulting numerical viscosity in most calculations damp out the growth of convection, delaying the development of strong convection. Because of the transient nature of supernovae, the timing of this convective growth can have a large quantitative impact on the strength of the supernova engine (Fryer & Young 2007).

On top of the numerical constraints in the engine itself, convective instabilities in the stellar progenitor also play a role in the growth of convection in the supernova engine. Early calculations justified perturbations in their numerical setup by arguing that these perturbations were on par with the convective instabilities seen in stellar convection (Herant et al. 1994; Fryer & Warren 2004; Fryer & Young 2007). More extreme initial asymmetries caused by vigorous silicon convection were proposed as the source for strong neutron star kicks (Burrows & Hayes 1996; Fryer 2004). But it was not until modelers coupled multidimensional simulations of the pre-collapse star with multidimensional engine models that the community has begun to accept the importance of stellar convection in accelerating the engine (Couch & Ott 2013; Fields & Couch 2021). Unfortunately, these pre-collapse models are also under-resolved and poorly understood, further preventing a quantitative, first-principles model for the supernova engine.

Calculating a remnant mass distribution of more massive stars depends on a different set of collapse properties, in particular, the angular momentum. If the collapsing star has considerable angular momentum and the proto-neutron star develops strong magnetic fields, a magnetar engine can drive an explosion. These explosions do not depend upon convection and will produce a different remnant mass than normal supernovae. Further, if a rapidly rotating star collapses to form a black hole, the resultant black hole accretion disk can drive a jet that is both believed to produce long-duration gamma-ray bursts (Woosley 1993; MacFadyen & Woosley 1999; MacFadyen et al. 2001) and eject much of the stellar envelope. Like our standard supernova models, simulations of magnetic field and jet production are under-resolved and are not yet predictive. Also paralleling the convective-engine paradigm, the angular momentum of the collapsing core depends upon models of stellar convection and their uncertainties. As such, it is currently beyond our computational power to model this more exotic engine and its effect on the remnant population. However, by making assumptions on the nature of this rotation, we can estimate the role rotating engines play on the compact remnant distribution.

All of these obstacles in developing a quantitative core-collapse engine model has led some groups to instead follow an approach that leverages knowledge gained from multidimensional models, one-dimensional solutions, and an understanding of turbulence growth to develop an analytic model for the remnant mass distribution (Fryer & Kalogera 2001; Fryer et al. 2012; Gogilashvili et al. 2021). In this paper, we apply a range of analytic studies of turbulence growth to core-collapse conditions to develop our intuition for the convective engine (Section 2). With this understanding of the physics behind this convective engine and its resultant uncertainties, we develop a range of mass distributions for neutron stars and low-mass black holes from supernovae (Section 3). The effect of rotation on the remnant mass is discussed in Section 4. With these models implemented into binary population synthesis models, we can then predict mass distributions in observed binary systems, including the merging binaries observed through gravitational waves (Section 5). We conclude with a summary of the primary findings of this paper.

Driven both by theoretical studies and observational insight, the basic paradigm behind the core-collapse engine has evolved with time. For type Ib/c and II supernovae, the energy source for the explosion is derived from the potential energy released when the core of a massive star collapses to nuclear densities. Although roughly 10⁵³ erg is released in the collapse, tapping just 1% of this energy has proven difficult. A clue into the nature of this explosion lies in the extensive mixing observed in the innermost ejecta (iron peak elements) of SN 1987A. Convection above the newly formed proto-neutron star both provides a means to extract the potential energy released and a way to prevent a pile-up of infalling material that would make it increasingly difficult to drive an explosion (for reviews, see Herant et al. 1994; Fryer & Young 2007; Fryer et al. 2021).

In this paper, we focus on the uncertainties in this convection and its growth, ultimately studying its effect on compact remnant masses. Herant and collaborators initially proposed Rayleigh–Taylor instabilities as the dominant driver behind convection (Herant et al. 1994) and multiple studies have confirmed that the growth time of this convection is very rapid and it can dominant the matter motion above the proto-neutron star (Fryer & Warren 2004; Fryer & Young 2007; Couch & Ott 2013; Vartanyan et al. 2018, 2019; Fields & Couch 2021). The entropy structure set by the stalled shock and infalling stellar material has some of the same properties of accreting neutron stars (Houck & Chevalier 1991; Fryer et al. 1996). The entropy in the shock decreases with decreasing shock velocity and lower accretion rates. As the bounce shock propagates outward, the shock velocity (and, hence, the entropy at the shock) decreases. In both supernovae and neutron star accretion, the entropy profile set by the outward propagation of the shock is susceptible to Rayleigh–Taylor instabilities. Growth timescales of the Rayleigh–Taylor instability is extremely rapid (Fryer & Young 2007) and, as we shall show here, strong turbulent convection can develop in the first 3–5 ms after bounce. Although other forms of turbulence exist, e.g., standing accretion shock instabilities (Blondin et al. 2003; Abdikamalov et al. 2021), the fast growth of the Rayleigh–Taylor instability suggests it will dominate the turbulent convection in many supernova scenarios. We will focus on Rayleigh–Taylor convection in this paper.

Rayleigh–Taylor instabilities occur when there is an inversion in material properties under an acceleration (such as gravity), i.e., high-entropy material lies beneath low-entropy material (e.g., low-density fluid beneath higher-density material). Such a density/entropy gradient occurs in many simulations (the extent depends upon the equation of state) of the collapse bounce of a massive core. But for convection to occur, initial perturbations in the entropy or "seeds" to this convection must exist (presumably from asymmetries in the collapsing core). In addition, the growth time depends on the viscosity of the material (Dimonte et al. 2007; Banerjee et al. 2011; Ramaprabhu et al. 2012). The numerical viscosity in the under-resolved core-collapse calculations is many orders of magnitude higher than Spitzer viscosity from electron/ion transport (Spitzer & Härm 1953), neutrino viscosity (Fryer & Young 2007), or effective viscosities from magnetic fields (Nakano 1976; Meintjes & Breedt 2015; Parker & Constantinou 2019).

Unfortunately, the complexity of the supernova convection problem prevents state-of-the-art calculations from fully resolving the viscous scale. Most of the supernova calculations using large eddy simulations which do not capture the small scale turbulence. The numerical viscosity exceeds any physical viscosity, delaying the growth of turbulence. It has been argued that most of the energy critical to driving an explosion will ultimately reside in the large-scale flows in an attempt to justify these under-resolved calculations. While true, this resolution is often critical in determining the growth time of these large-scale flows. The dynamic nature of core-collapse supernovae (short evolutionary timescales) requires accurate calculations of the growth time of this convection to develop a quantitative picture of these explosions. Most calculations in the literature do not include any representation of the sub-grid turbulence. Sub-grid models could be used to understand the turbulent evolution below the resolution scale of the simulations. Reynolds-averaged Navier–Stokes sub-grid models are commonplace in many engineering applications but such models also have their limitations and are typically tuned to match data for specific applications. In this section, we will study these sub-grid models in more detail.

Although further studies of the supernova engine in three dimensions are critical to increasing our understanding of supernovae, if you combine these uncertainties in mixing with the vast array of additional physical uncertainties, it becomes clear that a quantitative, first-principles model of the core-collapse engine remains an ultimate goal and not a near-term goal to understanding the compact remnant population. Understanding this convective supernova engine requires a two-prong approach: both an ever-improving set of two- and three-dimensional core-collapse models and a set of one-dimensional models that are guided by these studies to extract new insight from neutron star remnant populations.

One-dimensional models of supernovae have been used for many decades, allowing scientists to study the fate of collapsing cores for a broad range of supernova progenitors (Wilson et al. 1986; Mayle & Wilson 1988; Wilson & Mayle 1988; Fryer et al. 1999, 2018; Fryer & Kalogera 2001; Fröhlich et al. 2006; Fischer et al. 2010; Ugliano et al. 2012; Perego et al. 2015; Ertl et al. 2016; Prša et al. 2016; Couch et al. 2020). Initially, explosions were initiated either through piston or thermal bombs, driving strong explosions that matched observed supernova energies. With these studies, scientists began studying additional supernova observables like nucleosynthetic yields, using these yields to probe characteristics of the progenitor. Newer studies have instead focused on varying the efficiency of neutrino energy deposition to drive an explosion (Fröhlich et al. 2006; Fischer et al. 2010; Ugliano et al. 2012; Perego et al. 2015; Ertl et al. 2016; Prša et al. 2016). These later studies could potentially be used to probe aspects of the neutrino physics. Finally, another set of calculations have focused on the convective engine and parameterized the convection to drive the explosion (Fryer & Kalogera 2001; Fryer et al. 2018). The most recent of these have used physics-based sub-grid mixing models for their explosions (Mabanta et al. 2019; Couch et al. 2020).

These studies focused on only one or two instantiations of these recipes to produce a suite of models (often calibrating to SN 1987A and an assumption for this supernova's progenitor). Such studies are ideally suited to produce results for population studies, but we must bear in mind that, without a complete understanding of the physics, there are broad uncertainties in the theoretical predictions. Focusing on just the convection, we can begin to understand (and hopefully characterize) this range simply by looking at the different predictions from different sub-grid mixing algorithms, applying them to the profiles of post-bounce. For this study, we use a mixing-length theory (MLT) approach commonly used in astrophysics and a Reynolds-averaged Navier–Stokes (RANS) approach based on models developed in the engineering community Livescu et al. (2009). In a sense, MLT is a very particular RANS model, but it differs enough that it should be considered separately.

The MLT solution we apply is identical to the one described in detail by Couch et al. (2020). If we focus on the source terms (neglect the flux terms), the turbulent velocity growth is described by:

$$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho {v}_{\mathrm{turb}}^{2})}{\partial t}+\,\mathrm{Flux}\ \mathrm{Terms}=\rho {v}_{\mathrm{turb}}{\omega }_{\mathrm{BV}}^{2}{{\rm{\Lambda }}}_{\mathrm{mix}}-\rho {v}_{\mathrm{turb}}^{3}/{{\rm{\Lambda }}}_{\mathrm{mix}}.\end{eqnarray} \tag{ 1 }$$

where ρ is the density and v_turb is the turbulent velocity. In this calculation, ω_BV is the Brunt–Väisälä frequency

$$\begin{eqnarray}&&{\omega }_{\mathrm{BV}}^{2}=-{g}_{\mathrm{eff}}(1/\rho \partial \rho /\partial r-{(\rho {c}_{{\rm{s}}}^{2})}^{-1}\partial P/\partial r)\end{eqnarray} \tag{ 2 }$$

where c_s is the sound speed, P is the pressure and g_eff = −∂Φ/∂r + v_r ∂v_r /∂r is the effective gravity with Φ is the gravitational potential and v_r is the radial velocity. Finally, Λ_mix is the mixing length set to a fraction α_Λ of the pressure scale height.

Typical RANS solutions are very similar, but the notation, and driving terms are different. For this solution, we use the turbulent kinetic energy ($\widetilde{k}$). A simple driving term is in the RANS formulism (taking only the leading-order terms) is Livescu et al. (2009):

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho \widetilde{k}}{\partial t}=-{(\rho {\widetilde{U}}_{3}\widetilde{k})}_{,3}+{a}_{3}\displaystyle \frac{\partial P}{\partial r}-{R}_{33}{\widetilde{U}}_{\mathrm{3,3}}-1/2{R}_{{ii}\mathrm{3,3}}\end{eqnarray} \tag{ 3 }$$

where a₃ is:

$$\begin{eqnarray}&&{a}_{3}={{hC}}_{{\rm{t}}}{\tilde{k}}^{1/2}\end{eqnarray} \tag{ 4 }$$

where, following Livescu et al. (2009), we set C_t = 0.39 (constant defining the turbulent viscosity) and h is the pressure scale height. In these equations, ρ is the density, ${\widetilde{U}}_{3}$ is the average radial velocity from our one-dimensional calculations, R_ij is the Favre Reynolds stress described in Equation (10) of Livescu et al. (2009), $\widetilde{k}$ is the turbulent kinetic energy described in Equation (11) of Livescu et al. (2009). The 3 in these equations refers to the radial direction and, because we only have radial terms in our one-dimensional calculation, these are the only terms considered.

To compare the growth of these different instabilities, we take massive star progenitors and collapse and bounce the stars to determine post-bounce conditions. For our study, we use both 15 and 25 M_⊙ progenitor stars modeled with the Kepler code (Weaver et al. 2017). These stars are then collapsed using our one-dimensional core-collapse code (Fryer et al. 1999). This code includes a simple neutrino transport (gray flux-limited diffusion) for each of the three neutrino species, a dense nuclear equation of state, and general relativity assuming spherical symmetry. With this code, we can model the collapse and bounce of the massive star.

We use the post-bounce structure from these simulations to test the growth times of the instabilities. Figures 1 and 2 show the growth of turbulence at four different times post-bounce for 15 and 25 M_⊙ progenitors, respectively. The top panel shows the turbulent velocity using our MLT approximation using a relatively low value for the mixing length (α_Λ = 0.1). The bottom panels show the turbulent energy in our RANS formalism. Here we have taken the structure of the collapsed star at a 3, 6, 9, and 12 ms period post-bounce and, keeping these conditions fixed, determined the growth of instabilities over a 3 ms period. Both of our sub-grid models agree that Rayleigh–Taylor instabilities are strong and grow quickly in these post-bounce conditions (something not seen in the under-resolved multidimensional calculations). The models also agree on the regions where the instability grows. But they do not agree on which regions dominate this growth. Although analytical derivations show that our multidimensional models are not capturing this physics entirely, we still can not predict exactly the growth of turbulence. For this paper, we allow the growth time to vary to determine its effect on the compact remnant mass distribution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At this time, we can not predict the growth of convection exactly with either resolved multidimensional models or reduced-term analytic models. But we can understand the trends in this convective growth and use it to guide our remnant-mass models. In the next section, we show how measurements of the remnant mass distribution can help validate our convection models.

To determine how the growth of convection alters the remnant mass, we have altered the prescription used in the entropy-driven calculations of remnant mass distributions (Fryer & Kalogera 2001; Fryer et al. 2012). These calculations assume a convective region whose entropy evolves with time (mimicking the convection). In Fryer et al. (2012), the fast- and slow-convection models were developed to demonstrate the extremes allowed under the convective model. But the true answer likely lies in between these values. In this section, we produce a simple prescription to vary the remnant mass due to differences in the convective growth time.

This prescription does not include the possibility that massive star evolution may produce some stochasticity into the supernova explosion mechanism. Many stellar evolution calculations have produced wild variations in the core structure at collapse that then produce large variations in the final remnant mass (Laplace et al. 2021; Schneider et al. 2021). Whether or not these wild variations in stellar structure are real (versus numerical artifacts), the fact that the nature of the asymmetries in the collapsing core can alter the growth of convection suggests that at least some stochasticity exists in the models. We include remnant-mass prescriptions including this stochasticity at varying levels to determine whether individual stellar systems or population studies can be used to place limits on the level of this stochasticity.

3.1. Growth Time Variations

To study the effects of variations in mixing, we slowly altered the mixing algorithm in Fryer et al. (2012). By varying the growth time for the convective instabilities, we alter the energy injection and energy region. Studying variations in these properties, a wide range of explosion energies and properties can be produced for the same basic progenitor (Fryer et al. 2018). This study found that a wide range of explosion properties and remnant masses could be produced by varying the energy injection rate and region.

The final fate of these massive stars depends both on the growth time of the convection and the nature of the stellar core. Couch et al. (2020) found that, for the progenitor models in their study, they could produce an explosion with a sufficiently high mixing-length value. Increasing the mixing length accelerates the convective growth and, as this convection grows, the more likely it is to explode. This result demonstrates the same trends seen in other studies altering convection (e.g., Fryer et al. 2018) which form the basis for our models. In this section, we use these trends from one-dimensional explosion models, integrating them into prescriptions that incorporate the stellar structure, to determine the remnant mass for a given progenitor.

The stellar structure at collapse can play an even more critical role than convective growth time in determining the stellar fate. Although it is well known that a single compactness parameter is insufficient to determine the final fate of a model (e.g., Burrows et al. 2020), our models are based on a prescription that follows the time-evolution of the infalling stellar material (Fryer & Kalogera 2001; Fryer et al. 2012). To compare this method to approaches solely using compactness parameter, we plot in Figure 3 both the compactness parameter and our final remnant masses (using two different convective growth times). Although there are trends between remnant masses from our prescription and the compactness parameter, as with the models of Couch et al. (2020), Fryer et al. (2018), the ultimately fate is more complex than this.

Zoom In Zoom Out Reset image size

The final structure of the star, which dictates the fate of the collapsing star, whether using a simple compactness parameter or our more complex prescription, can vary dramatically with stellar model (Sukhbold et al. 2018). The multiple carbon and oxygen-shell-burning phases at the end of a star's life can evolve very differently depending on small amounts of mixing. The largest variation in this evolution occurs in stars with zero-age main-sequence mass lying between 20 and 30 M_⊙. These variations can alter the fate, allowing supernova explosions in these more massive stars. Fryer et al. (2012) used a single set of stellar models (Woosley et al. 2002) to guide their stellar fates. Understanding this convection is critical to understanding the fate of stars in this 20–30 M_⊙. We will discuss this further below in our parameterized models.

Final explosion energies are more difficult to calculate because one-dimensional models typically do not include additional energy from further fallback (although see Fryer et al. 2018). However, our models do calculate the energy in the convective region when it is able to blow off the infalling material and this provides a rough estimate of the explosion energy (explosion energy is lost in the ejection of the star but energy is gained through fallback, so this value can be higher or lower than the final explosion energy). The corresponding explosion energy for these zero-age main-sequence progenitors is shown in Figure 4.

Zoom In Zoom Out Reset image size

If we assume that there is a smooth trend in core properties with respect to CO core mass, we can develop a simple smooth fit by varying the mixing efficiency set by the growth time. A set of models to fill in the solutions from Fryer et al. (2012) for the final baryonic remnant mass (M_remnant) is given by:

$$\begin{eqnarray}\begin{array}{rcl}{M}_{\mathrm{remnant}} & = & 1.2+0.05{f}_{\mathrm{mix}}+0.01{({M}_{\mathrm{CO}}/{f}_{\mathrm{mix}})}^{2}\\ & & +\,{e}^{{f}_{\mathrm{mix}}({M}_{\mathrm{CO}}-{M}_{\mathrm{crit}})}\end{array}\end{eqnarray} \tag{ 5 }$$

where M_CO is the carbon/oxygen core mass, M_crit is a critical mass for black hole formation (for our models, we typically use M_crit = 5.75 M_⊙), and f_mix is a parameter describing the mixing growth time. This mass can not exceed the total stellar mass at collapse so the final remnant remnant mass (baryonic) is:

$$\begin{eqnarray}&&{M}_{\mathrm{remnant}}=\min ({M}^{\mathrm{eq}.\,5},{M}_{\mathrm{collapse}})\end{eqnarray} \tag{ 6 }$$

where M^{eq. 5} is the remnant mass from Equation (5), M_collapse is mass of the star at collapse. The resultant remnant masses using this prescription are overlaid upon our calculated remnant masses in Figure 5 where we have varied f_mix from 0.5 to 4.0 where f_mix = 4.0 corresponds to a fast growth time (rapid model in Fryer et al. 2012; blue line in Figure 5) and f_mix = 0.5 corresponds to a slow growth time (delayed model in Fryer et al. 2012; red line in Figure 5).

Zoom In Zoom Out Reset image size

We use CO core mass as the principle stellar-structure quantity for the fate of these stars because most population synthesis models include its evolution with binary mass transfer and mass loss. However, the CO core mass does not capture the full structure of the core that scientists characterize by either the compactness parameter or the more detailed infall prescription used in our models. In addition, the differences between the CO mass and this structure also depend on metallicity. Barring population synthesis models that include the evolution of these more accurate prescriptions, we are currently limited to simplified remnant descriptions depending on properties like the CO core mass. An alternate approach could use the ZAMS mass alone, but this would not include evolution caused by binary mass transfer and loss.

As we discussed above, late-stage carbon- and oxygen-shell-burning stages depend sensitively on the progenitors and progenitor models. In general, the remnant mass increases (as in our Equation (5)) with CO mass. But this trend can turn around briefly for stars with final CO cores in the 6–9 M_⊙ mass range. In Fryer et al. (2012), final CO cores above 6 M_⊙ formed black holes and, although the details of the convection allowed some stars with CO core masses above 6 M_⊙ to produce supernova explosions, the explosions were too weak to produce neutron stars. In their models, the final black hole masses were below the mass expected from a direct collapse (no explosion), limiting the final black hole mass in a region roughly between 6.5 and 10 M_⊙ to roughly 75% of the mass of the star at collapse. One possible parameterization of this effect is to limit the mass of the black hole to a fraction (f_{CO burning}) of the collapse mass between a range M_low and M_high where, for Fryer et al. (2012), these values are roughly f_{CO burning} = 0.75, M_low = 6.5, and M_high = 10.0. We do not include this effect in the populations discussed in the rest of this paper.

If we assume stars follow the Salpeter (α_IMF = 2.35) initial mass function (IMF) and the CO and He cores from Kepler models, the prevalence of a gap in the compact remnant masses depends upon this mixing factor. Figure 6 shows the distribution using seven values of f_mix: 0.5, 0.7, 1.0, 1.4, 2.0, 2.8, and 4.0. The top panel uses M_crit = 5.75 M_⊙, the bottom panel uses M_crit = 4.75 M_⊙. For the rest of the figures in this paper, we use M_crit = 5.75 M_⊙ as this value provides a closer match to past studies (Fryer et al. 2012). However, this value can change depending on physical factors in the supernova engine (e.g., neutrino physics, equation of state). The larger mixing factors produce larger and deeper gaps in the remnant mass distributions, which can be compared to observations to constrain the growth of convection in the engine. Although all of these models produce compact remnant masses, the number of systems in the gap can vary by a factor of 10.

Zoom In Zoom Out Reset image size

This simple model produces a smooth trend in the remnant masses with progenitor mass. Variations in the remnant mass distributions arise from variations in the CO core mass in our coarse-sampled (in mass) stellar models. But the variations in stellar mixing produce even stronger variations in the CO core mass and, combined with variations in seeds of convection in the supernova engine, much larger variations could exist in the mass distribution. To study the full distribution, we need to include this stochastic variation.

3.2. Stochasticity from Convection in the Star and in the Supernova Engine

The simple smooth remnant mass/CO core mass distribution does not account for variation in the structure and uniformity of the star at collapse caused by stellar mixing. Here we introduce two modifications to our smooth profile to capture the stochasticity that is currently known to be present in supernova explosion calculations. The first captures the growth time variation that can be caused by asphericities in the collapsing core that then seed the convection. The latter captures the dramatic structural differences seen in some stars where, for example, the core of a 14.9 M_⊙ star follows our basic trend producing a NS, but the core of a 15.0 M_⊙ collapses to form a black hole.

Asymmetries in the collapsing core provide the seeds for convection that then change the growth time of the convection. To capture this effect, we sample f_mix from a Gaussian distribution:

$$\begin{eqnarray}&&F\left({f}_{\mathrm{mix}}\right)={f}_{\mathrm{norm}}{e}^{-{f}_{\mathrm{stoch}}{\left({f}_{0}-{f}_{\mathrm{mix}}\right)}^{2}}\end{eqnarray} \tag{ 7 }$$

where F(f_mix) is the fraction of stars with a value of the mixing parameter equal to f_mix, f₀ is the average value of f_mix, f_stoch denotes the variation in the growth times, and f_norm is a normalization factor. Using α_IMF = 2.35 and our Kepler models, the remnant mass distribution for f₀ = 1, f_stoch = 1, 2, 4, 8, 16 produces remnant mass distributions in Figure 7. For this distribution, we use a Monte Carlo sampling of this Gaussian distribution using 10 million sampling points. By using a low normalization factor, we have skewed these distributions toward low values of mixing.

Zoom In Zoom Out Reset image size

Broader changes are required to capture the more dramatic variation seen in calculations where slight changes during a shell-burning episode dramatically alter the size of a burning layer, ultimately altering the structure of the core at collapse and its final fate. One way to implement these changes is to use the detailed grids with stochastic results to determine CO core masses instead of using a formula fit to the CO core mass that is smooth. The evolutionary changes that alter these fates are typically reflected in the CO core mass. However, because different codes (and indeed different initial grid resolutions) can produce very different results (the nature of a stochastic process whether real or numerical), no grid of models represents the ultimate answer. It is not even obvious how to interpolate between masses in such stochastic simulations.

We instead propose a more parameterized approach where we assume some percentage f_conv of stars have a very different stellar burning phase where the fate ranges from the proto-compact object and the mass of the star at collapse. The proto-compact object mass is M_proto from Equation (18) of Fryer et al. (2012). By setting f_conv to 1%–10%, we alter the remnant mass of 1%–10% of our CO cores to a mass sampled from our entire set of systems focusing on stars with 15 < M_ZAMS/M_⊙ < 50 (ZAMS corresponds to the Zero-Age Main-Sequence mass). For an individual system, this means that some subset of stars can produce a black hole instead of a neutron star (or vice versa). Using α_IMF = 2.35 and our Kepler models, the remnant mass distribution for f_conv = 0%, 1%, 10% (again using 10 million sampling points) and f_mix = 1, 8 produces remnant mass distributions in Figure 8. Here we have used f_mix = 8 to study the extreme solutions. In general, it is likely that f_mix lies between 0.5 and 4. Although this stochasticity can dramatically change the remnant mass of a given CO core, this variation does not affect the integrated remnant mass distribution if the sample size is large.

Zoom In Zoom Out Reset image size

As we show in Figure 8, stochasticity in the progenitor is unlikely to affect the remnant mass distribution if the sample size of remnants is large (e.g., for population synthesis calculations studying a large number of binaries) because such large sampling smooths out these stochastic variations. However, the current observed distribution of compact remnants is small. Stochastic variations can explain single observations of peculiar observations: e.g., an observed system with a black hole forming from a 15 M_⊙ progenitor or a neutron star produced by a 30 M_⊙ progenitor. Population studies using small sampling sizes could consider including stochasticity. Figure 9 shows the remnant distribution with stochastic model using a sampling of only 1000 points and combining both progenitor and engine stochasticities discussed in this section. Although varying the amount of stochasticity in the progenitor (f_conv = 0%, 1%, 10%), the effect of this stochasticity remains minimal on the remnant distribution. From this analysis, unless a study is trying to explain the masses of a small number of systems or population studies are using samplings below 1000 systems, stochasticity is not an important effect.

Zoom In Zoom Out Reset image size

Two primary black hole formation scenarios have been identified (Fryer 1999): fallback and direct collapse formation. Most studies of compact remnant formation focus on determining the dividing lines between different compact remnant fates (NS, fallback BH, direct BH) and the amount of mass accreted in the fallback BH formation. It is generally assumed that in the direct-collapse BH formation that all of the mass collapses to form a black hole. This is not the case for rapidly rotating stars. Since the first models of collapsar jets through a massive star, it has been known that the jet drives a shock that propagates around the star, ejecting much of the stellar mass (MacFadyen & Woosley 1999; MacFadyen et al. 2001). This jet-driven mass ejection limits the final mass of the black hole. For progenitor stars with sufficient angular momentum to form a disk that drives a jet, this mass ejection process can eject much of the star even if the star collapses directly to a black hole. To determine the mass ejection from disk- and jet-driven outflows, we must understand disk formation and the angular momentum in stars.

To form a disk, the star must have sufficient angular momentum such that infalling material of the imploding star is centrifugally supported in a disk. In general, the specific angular momentum of a star increases with radius. For example, if the star were strictly coupled with a constant angular velocity (ω), the angular momentum of the star increases with the square of the radial position in the star (j = I ω ∝ ω r²). Many stellar evolution models assume only mild coupling in the stars, allowing the inner core to spin more rapidly than the outer layers. But even for these stars, the angular momentum of a star increases as we move outward in mass coordinate (see, for example Heger et al. 2000). The imploding star continues to accrete onto its core until the angular momentum of infalling stellar material is sufficiently high to produce a disk. Because the angular momentum increases with radius, even if the innermost material does not have enough angular momentum to form a disk, the outer layers of the star may. The innermost material will accrete onto the black hole, increasing its mass until the higher-angular momentum material forms a disk. ⁷

We can use the angular momentum profile to determine the mass of the core when a disk forms. This angular momentum also sets the spin of the black hole produced when the disk forms. Figure 10 shows the mass of the remnant at the time of the formation of a 300 km disk as a function of the final black hole mass spin period: ${a}_{\mathrm{BH}}=({{cJ}}_{\mathrm{BH}})/({{GM}}_{\mathrm{BH}}^{2})$ where c is the speed of light, J_BH is the black hole angular momentum, G is the gravitational constant and M_BH is the black hole mass. This figure uses two stars from the GENEC stellar evolution code (Maeder & Meynet 2000) with weak coupling where we artificially raise and lower the spin period to produce a range of angular momentum profiles. The faster the spin rate, the earlier a disk forms and the higher the value of the final black hole spin. After the disk forms, the core will continue to accrete until a disk-driven jet is able to eject the star. In this manner, the mass of the remnant mass at disk formation in Figure 10 shows the minimum remnant mass of a collapsing core as a function of its angular momentum. If the coupling is strong, as suggested by (e.g., Fuller & Ma 2019; Fuller et al. 2019), the angular momentum in the helium core is insufficient to form a disk.

Zoom In Zoom Out Reset image size

This black hole accretion disk is the power source invoked in the collapsar gamma-ray burst engine. We will assume that the efficiency at which the disk ejects mass is proportional to the energy expected in this black hole accretion disk engine. There have been many studies of these black hole accretion disks and the energy of their outflows. For our study, we use a simple formula derived by Heger et al. (2003) based on the disk models from Popham et al. (1999) using the Blandford & Znajek (1977) energy mechanism:

$$\begin{eqnarray}\begin{array}{rcl}{L}_{\mathrm{BZ}} & = & {10}^{50}{a}_{\mathrm{spin}}^{2}{10}^{0.1/(1-{a}_{\mathrm{spin}})-0.1}\\ & & \times {\left(\displaystyle \frac{{M}_{\mathrm{BH}}}{3\ {M}_{\odot }}\right)}^{-3}\displaystyle \frac{{dM}/{dt}}{0.1{M}_{\odot }\,{{\rm{s}}}^{-1}}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\end{array}\end{eqnarray} \tag{ 8 }$$

where dM/dt is the accretion rate. ⁸ Figure 11 shows this energy as a function of the final black hole spin using the same models as Figure 10.

Zoom In Zoom Out Reset image size

Even with these approximations, it is still not straightforward to determine the final black hole mass because we must estimate the amount of disk material that accretes while its outflow and jet ejects the rest of the star. For powerful jet-driven explosions (high accretion rates), the outflow can eject much of the infalling star. In such an extreme, we assume that the accretion of 1 M_⊙ of material is sufficient to eject the entire star. If this occurs, instead of our standard black hole remnant mass formula (for direct collapse, this value is close to the mass at collapse), the remnant mass for a rotating system ${M}_{\mathrm{BH}}^{\mathrm{rot}}$ can be as low as:

$$\begin{eqnarray}&&{M}_{\mathrm{BH}}^{\mathrm{rot}}=\max (3.0,7.5-15.0({a}_{\mathrm{BH}}-{a}_{\mathrm{crit}})){M}_{\odot }\end{eqnarray} \tag{ 9 }$$

for a_BH > a_crit where a_crit is set by the models (∼0.7, 0.45 for the 32, 40 M_⊙ stars in our study) and a_BH is the final black hole spin. If a_BH < a_crit, we assume the angular momentum in the star is insufficient to produce a disk that produces a sufficiently energetic outflow to alter our standard remnant mass.

Unfortunately, we do not know the angular momenta of massive stars prior to collapse. Stars are born rotating and, if the different burning layers do not couple, the core could be spinning rapidly at collapse. A weak coupling (typically assumed to be produced by magnetic fields), has been argued as an explanation of the spin periods behind pulsars. If we assume a set fraction of our black-hole-forming stars have a spin value that is sufficiently high to cause the collapsed core to eject half of its final mass, we can use our IMF to predict BH mass distributions for our two different mass-limit extremes (Figure 12). This figure varies the number of high-rotating systems from 0 to 20% of the black-hole-forming stars.

Zoom In Zoom Out Reset image size

For stars that lie above the pair-instability mass gap corresponding to a ZAMS mass above ∼260 M_⊙, it is possible that, unlike BH-forming scenarios where first a Proto-NS is formed that accretes until it collapses to a BH (first forming a low-mass BH), the star can first form a proto-BH (Fryer et al. 2001). In such a case, the remnant mass will be closer to the proto-black hole mass, producing gravitational masses above 60–70 M_⊙. We add the distribution from a model assuming proto-BH formation in our plot of black hole masses (Figure 12).

Although measurements of X-ray binaries suggests that many collapsing stars have high angular momenta, these measurements are both difficult to make and limited to a select set of formation scenarios requiring short-period, accreting binaries (Miller & Miller 2015; Belczynski et al. 2021). Gravitational wave measurements of black hole binaries suggest that many black hole systems have low spin values (Abbott et al. 2019), suggesting that angular momentum coupling in massive stars is strong (e.g., Fuller & Ma 2019; Fuller et al. 2019). If this coupling is strong, forming fast-spinning black holes and collapsed cores with disks will be limited to systems that can be spun up prior to collapse. One method to spin up these stars is through tidal locking (Kushnir et al. 2016, 2017; Hotokezaka & Piran 2017; Qin et al. 2018; Zaldarriaga et al. 2018). Incorporating a solution only using this tidal-locking spin-up, Belczynski et al. (2020) and Olejak & Belczynski (2021) were able to use binary population models to match the observed χ_eff values observed by gravitational wave observatories (Abbott et al. 2019).

In this section, we study the role of this tidal-locked angular momentum model using the results of binary population synthesis models. We use the population synthesis code StarTrack (Belczynski et al. 2008) to calculate evolution of four populations of massive binary stars. Each population consists of 10,000 massive binary systems, in which initial primary mass is drawn from power-law IMF with exponent α = −2.3 and in the ZAMS mass range M_ZAMSA = 5–200 M_⊙, and secondary initial mass is in sampled from a ZAMS range M_ZAMSB = 3–200 M_⊙ and is obtained from uniform mass ratio (q = M_ZAMSB/M_ZAMSB) distribution. For each population we study two metallicities, Z = 0.1 Z_⊙ and 0.01 Z_⊙ (with Z_⊙ =0.014), and two mixing parameters: f_mix = 0.5, 4.0. These massive binaries are potential progenitors of double compact objects (NS-NS, BH-NS, and BH-BH systems). Note that formation of double compact objects is expected only from a very small fraction of massive binaries.

In our calculations, we assume standard wind losses for massive stars: O/B star winds (Vink et al. 2001), LBV mass loss, and WR star winds (specific prescriptions for these winds are listed in Section 2.2 of Belczynski et al. 2010). We treat the accretion onto compact objects during the Roche lobe overflow and from stellar winds using the analytic approximations presented by King et al. (2001) and by Mondal et al. (2020), and limit accretion during the common envelope (CE) phase to 5% of the Bondi rate (MacLeod et al. 2017). In calculations of remnant mass, we employ the following scheme: for NSs and low-mass BHs we employ the prescriptions presented in this study where we calculate the BH/NS mass with f_mix = 0.5, 4.0 but do not account for effects of stochasticity. We allow for formation of BHs up to 90 M_⊙ (as detailed in Belczynski 2020) taking into account recent advancements in understanding effects of pair-instabilities on evolution of massive stars (Farmer et al. 2020; Costa et al. 2021; Farrell et al. 2021). We do not produce black holes from systems above the pair-instability limit (stars above ∼260 M_⊙). We employ the fallback decreased NS/BH natal kicks with σ = 265 km s⁻¹, we do not allow CE survival for Hertzsprung gap donors (submodels B in our past calculations). Note that we allow for standard CE development criteria and most double compact objects to form through the CE channel (Belczynski et al. 2016), and not through the stable Roche lobe overflow channel (van den Heuvel et al. 2017; Olejak et al. 2021). The most updated description of StarTrack is given by Belczynski et al. (2020) and the model M30 in this study describes our standard choices of input physics.

From the progenitors of BH-BH binaries from these calculations, we obtain orbital properties prior to the collapse of each star. These systems undergo binary interactions that tighten the orbits and eject mass. To determine the spin of these stars, we assume direct tidal locking of the star at collapse for systems with periods less than 1.3 days (assuming the timescale for tidal locking is longer than the evolution timescale for binaries with greater orbital periods as did Belczynski et al. 2020). For each of the binaries in our population synthesis model, we calculate the orbital period prior to the collapse to a black hole of both the primary and the secondary star. For these calculations, we refer to the primary star in the system as the star that has a larger ZAMS mass in the binary (the secondary star has the lower ZAMS mass in the system). The specific angular momenta at a radius of 10¹⁰ cm in each of these stars is shown in Figure 13. In this figure, we compare results from binary population synthesis models using two metallicities and two mixing parameters. We use 10¹⁰ cm because, for the massive stars in this study, this lies well within the helium layer and provides a maximum extent (highest angular momentum) for systems whose compact remnant is not so massive that the formation of the disk will still drive powerful jets. For the material to hang up in a disk, the centrifugal acceleration must exceed that of gravity above the innermost stable circular orbit (given for a non-rotating black hole):

$$\begin{eqnarray}&&{r}_{\mathrm{ISCO}}=6{{GM}}_{\mathrm{BH}}/{c}^{2}\end{eqnarray} \tag{ 10 }$$

where G is the gravitational constant, c is the speed of light, and M_BH is the enclosed mass that forms a black hole. The corresponding minimum angular momentum (j_crit)is:

$$\begin{eqnarray}&&{j}_{\mathrm{crit}}=\sqrt{6}{{GM}}_{\mathrm{BH}}/c.\end{eqnarray} \tag{ 11 }$$

For a 5 M_⊙ black hole, this critical specific angular momentum is 5.4 × 10¹⁶ cm² s⁻¹. As we shall see in Figure 14, the mass enclosed within 10¹⁰ cm for BH-forming stars exceeds 5 M_⊙. For nearly all of our tidally locked stars, the angular momentum at 10¹⁰ cm is insufficient to produce a disk. The specific angular momentum increases with radius squared, so the outer layers in some cases will produce a disk. The stars with the highest angular momentum are a subset of the collapsing secondaries in tightly bound orbits. These stars are likely to form disks around low-mass BH cores and are the most likely to form strong jets (see Equation (8)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To determine the radius (and mass) at which the infalling material will form a disk, we must first study the structure of the star at collapse. Figure 14 shows the radius versus enclosed mass for different black hole forming progenitors using three different stellar evolution codes: KEPLER (Woosley et al. 2002), GENEC (Maeder & Meynet 2000), and MESA (Paxton et al. 2011). The MESA and KEPLER calculations follow the evolution through silicon core burning and to stellar collapse. The GENEC code does not model the star through silicon burning and, not surprisingly, is slightly more extended than the KEPLER and MESA simulations in the inner core (within 4 M_⊙). Beyond the inner core, the results vary between different stellar mass, rotation, and code used. Given the uncertainties in the models, we use an approximate fit for the enclosed mass above 2 M_⊙ (M_enc) to capture the features of this structure to calculate the angular momentum of the tidally locked star:

$$\begin{eqnarray}&&{M}_{\mathrm{enc}}=7({\mathrm{log}}_{10}(r+{10}^{9})-9){M}_{\odot }\end{eqnarray} \tag{ 12 }$$

where r is the stellar radius in cm.

With this fit to the stellar model, we can use the spins and mass/radius relation of the stars to determine the position where the angular momentum exceeds that needed to form a disk for the enclosed mass. Many of our systems are not spinning rapidly enough to form a disk. Figure 15 shows the enclosed mass when the angular momentum from tidal locking causes the star to form a disk as a function of collapse mass for those systems that do form a disk (7%–10% of the binaries in our four population synthesis models). The strength of the jet produced in these disks decreases rapidly with the mass of the interior compact remnant (see Equation (8)). Only 1% of our binary systems produce black holes less massive than 10 M_⊙ when the disk forms. Because the power of the jet is proportional to the energy in the disk and the disk energy drops as the event horizon increases (with remnant mass), it is difficult to make a strong gamma-ray burst jet with more massive black holes. As such, only 1% of our models would produce strong collapsar Gamma-Ray Bursts. None of our systems had enough angular momentum to produce a disk around the core when the core was less than 3 M_⊙ and these systems would not produce rotationally powered supernovae.

Zoom In Zoom Out Reset image size

In this paper, we build on our latest understanding of stellar evolution and the core-collapse supernova engine to study the compact remnant distribution. We summarize the primary results here:

• 1.  
  State-of-the-art simulations can not resolve the growth of convection. Although common RANS models used to capture this sub-grid convection agree on the short growth times and the dominate regions of convective instabilities, we demonstrate that the relative strengths depend on the leading terms in these models.
• 2.  
  We demonstrate that the convective growth time can produce a range of solutions for the remnant versus CO core mass. The number of remnants within the NS/BH mass gap region depends on the growth timescales. We provide a simple parameterized formula for these remnant masses.
• 3.  
  We discussed two forms of stochasticity in the remnant mass distribution. The first captures the growth time variation that can be caused by asphericities in the collapsing core that then seed the convection. The latter captures the dramatic structural differences seen in some stars (e.g., small differences in shell burning lead to large changes in pre-collapse star structure). Although stochasticity can alter the fate of an individual system, its effect on large populations is much less important.
• 4.  
  Although evidence continues to grow that the standard engine behind supernovae is from the convective model studied here, rapidly rotating models can explode via alternative mechanisms. We studied the effect of rapid rotation on the black hole remnant mass distribution and its potential to fill in the pair-instability mass gap.
• 5.  
  We studied the effect of rotation on the black hole mass distribution if rapid rotation in stars is only produced through tidal locking using recent population synthesis models. In such a spin-up scenario, roughly 10% of stars in close binaries achieve sufficiently high angular momenta to form a disk. Only 1% of these stars form disks around remnants with masses less than 10 M_⊙. If this is the only way to produce fast-spinning stars, rotationally driven explosions are common in only a small fraction of BH-forming stars and will not produce NSs.

The remnant mass distribution can be used to help us better constrain the growth timescales of convection in supernovae, reducing this key uncertainty in our understanding of the supernova engine. Upcoming measurements both through lensing (Mroz et al. 2021) and gravitational wave observations of mergers with next-generation detectors (Ng et al. 2021) are set to grow the number (and fidelity) of remnant distributions. Coupled with better supernova calculations and binary studies (e.g., Olejak et al. 2022), we can use these observations to drive breakthroughs in our understanding of the supernova engine and its implications on all of the astrophysics affected by these explosions.

We would like to thank Eric Burns, Mathieu Renzo, Lieke van Son, and an anonymous referee for useful comments improving this paper. The work by C.L.F. was supported by the US Department of Energy through the Los Alamos National Laboratory. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218CNA000001). K.B. and A.O. acknowledge support from the Polish National Science Center (NCN) grant Maestro (2018/30/A/ST9/00050).