COMPACT REMNANT MASS FUNCTION: DEPENDENCE ON THE EXPLOSION MECHANISM AND METALLICITY

To determine the energy of our convection-enhanced, neutrino-driven engines, we assume that the explosion energy is equal to the energy stored in the convective region at the time of collapse. Fryer (2006) has already estimated this energy and we review this derivation here. Colgate et al. (1993) found that they could approximate the convective region as a roughly constant entropy atmosphere bounded by the proto-NS and the shock region of the infalling stellar envelope. The structure of this atmosphere is then well defined, with the pressure given by

$$\begin{eqnarray} P(r) &=& \big[0.25 S_0 M_{\rm NS} G (S_{\rm rad})^{-1}(1/r - 1/r_{\rm shock})\nonumber\\ &&+ P^{1/4}_{\rm shock}\big]^4\,{\rm erg \, g^{-1}}, \end{eqnarray} \tag{ 1 }$$

where M_NS is the proto-NS mass, G is the gravitational constant, S_rad is the entropy in Boltzmann's constant per nucleon, S₀ = 1.5 × 10⁻¹¹, and r_shock and P_shock are the radius and pressure of the accretion shock forming the outer bound of the convective region. P_shock is set to the ram pressure of the infalling material. If we assume that the infalling star is accreting at the free-fall rate, mass continuity gives (Fryer 2006)

$$\begin{eqnarray} &&P_{\rm shock}(r) = 1/2 \rho _{\rm shock} v^2_{\rm free\hbox{-}fall} = (2 G M_{\rm NS})^{0.5} \dot{M}_{\rm acc}\big/ \big(8 \pi r^{2.5}_{\rm shock}\big),\nonumber\\ \end{eqnarray} \tag{ 2 }$$

where the free-fall velocity v_free-fall is determined by the proto-NS mass and the accretion rate $\dot{M}_{\rm acc}$ is determined by the structure of the progenitor star assuming a pressureless collapse (which is a good approximation for the infall). For a radiation-dominated gas, the internal energy density is 3 × P(r) (Fryer 2006)

$$\begin{eqnarray} u_{\rm convection}(r) &=& 3 \bigg[4.7 \times 10^8 \frac{M_{\rm NS}}{M_\odot } \frac{10 k_{\rm B}\,{\rm nucleon}^{-1}}{S_{\rm rad}} \nonumber\\ &&\times\bigg(\frac{10^6\,{\rm cm}}{r} - \frac{10^6\,{\rm cm}}{r_{\rm shock}} \bigg) + 1.2 \times 10^6 \bigg(\frac{M_{\rm NS}}{M_\odot } \bigg)^{1/8} \nonumber \\ &&\times\bigg(\frac{\dot{M}_{\rm acc}}{M_\odot \,\rm {s}^{-1}} \bigg)^{1/4} \bigg(\frac{2\times 10^7\,{\rm cm}}{r_{\rm shock}} \bigg)^{5/8}\bigg]^4 {\rm erg \, cm^{-3}}. \nonumber\\ \end{eqnarray} \tag{ 3 }$$

Integrating over the entire atmosphere, we can derive the maximum energy stored in the convective region. If the energy in the convective region is above this maximum value, the pressure in the atmosphere will be larger than the infall pressure, causing the shock radius to expand. By the time the shock expands to 1000 km, the expansion is really an explosion. After the explosion is launched, the atmosphere is too thin to absorb much neutrino energy. Neutrinos leaking from the core can no longer add much energy to the now exploding atmosphere. The energy of the explosion from the neutrino-driven mechanism is then limited to roughly the energy when the explosion is launched.

For a typical atmospheric entropy of 10 k_B nucleon⁻¹ and shock radius of 1000 km, energies above 5 × 10⁵⁰ erg only occur for accretion rates above 1 M_☉ s⁻¹. See Fryer (2006) for a study of the dependence of the explosion energies on the atmosphere parameters. The energy that can be stored in the convective region decreases as the accretion rate decreases. The accretion rate of the infalling stellar material decreases with time, causing the total available explosion energy to decrease with time (Figure 1). If the delay in the explosion is long, the explosion will be weak.

Zoom In Zoom Out Reset image size

This model assumes the explosion energy is stored in the convective region and provides a natural explanation for why the observed explosion energy for most supernovae are found in the range ∼0.5–2 × 10⁵¹ erg even though the potential energy released in stellar collapse is ∼10⁵³ erg. The peak energy stored in the convective region is roughly a few times 10⁵¹ erg, and this engine cannot produce a stronger explosion. If the explosion occurs less than 250 ms after bounce, the energy is above 10⁵¹ erg for most progenitors. Most observed supernovae have explosion energies within these limits. Magnetohydrodynamic engines produce a much broader range of supernova energies and are thus not the likely engines behind "standard" core-collapse supernovae.

Our basic models consider only the convection-enhanced neutrino-driven supernova engine, since this is able to account for the majority of supernovae. Within this basic set, we still vary the length of the delay time prior to the explosion. One model allows for considerable contribution from the SASI engine, a current focus of many supernova groups (Blondin et al. 2003; Burrows et al. 2007; Bruenn et al. 2009; Scheck et al. 2008; Marek & Janka 2009). The SASI has a growth time that is typically over an order of magnitude longer than the Rayleigh–Taylor instabilities, and hence is unable to drive an explosion at early times. For our standard model, we allow an explosion only after the accretion rate drops below 1 M_☉ s⁻¹. At this point, we calculate the energy in the convective region. Comparing this explosion energy with the binding energy of the star as a function of radius, we can calculate the remnant mass. Here we assume that the outer layers of the star are ejected first and eject stellar material until the explosion energy and binding energy are equal. In this formalism, stars below about 15 M_☉ produce explosions enhanced by convection dominated by the instabilities developing early, whereas above this mass the explosions are enhanced by convection dominated by late-time convection with large-scale modes characterized by the SASI. For many progenitors, the latter occurs at later times. If the proto-NS exceeds 3 M_☉, we assume the star collapses promptly to a BH (a.k.a. the "direct" BH formation scenario from Fryer & Kalogera 2001).

Our alternate convection-enhanced supernova engine model focuses on strong convective models developing rapidly and driving an explosion within the first 250 ms (e.g., Herant et al. 1994; Fryer & Warren 2002). In this engine, we assume an explosion occurs when the accretion rate drops below 3 M_☉ s⁻¹. If the proto-NS exceeds 3 M_☉ or the delay time exceeds 250 ms, we assume the star collapses to a BH.

Figure 2 shows the supernova explosion energies for these two models as a function of progenitor mass using the Woosley et al. (2002) progenitors at solar and zero metallicity. The delayed explosions tend to be weaker than our rapid explosion model, producing many more explosions with energies below 10⁵¹ erg. The explosion energies depend upon the prescriptions for mass loss (or, more accurately, the coefficients used in the mass-loss prescription) and convection, and we discuss the differences in these models in Section 3.5. Above 30–35 M_☉, the mass loss plays the deciding role in determining the remnant mass, as we discuss this in detail in Section 3.3.

Zoom In Zoom Out Reset image size

After the initial launch of the convection-enhanced explosion, the region above the proto-NS is evacuated. At this point, the energy deposited by neutrinos in this region is drastically reduced. For our basic engine, we assume the energy in the supernova available to eject the star is limited to that stored in the convective region at the launch of the explosion.

There exist other mechanisms to extract the enormous energy reserved in the collapsed stellar core. For example, if the stellar core is rotating rapidly, strong magnetic fields may develop. As the proto-NS contracts, an enormous amount of energy can be stored in the rotation, driving a second outburst from the NS. This energy source has been invoked to explain a variety of energetic supernova explosions (e.g., Maeda et al. 2007; Kasen & Bildsten 2010). Although quantitative calculations of rotation powered explosions do not exist, we can place an upper limit based on the predicted spin-energy in the NS (∼3 × 10⁵²  erg for a millisecond NS). If such an explosion occurs, the fallback onto the proto-NS will be minimal. In our third ("magnetar") model for remnant distributions, we assume this additional explosive engine occurs in all systems.

We note that fallback can also cause a secondary outburst as the material accretes onto the proto-NS (Fryer et al. 2006). This is less explosive than a magnetar and, to show the extremes, we focus only on the magnetar explosion.

The collapsar engine behind hypernovae posits that a rapidly spinning star collapses down to a BH. The jet produced by the accretion disk that forms around the BH "drills" through the stellar envelope and ultimately ejects the star, preventing further accretion (MacFadyen & Woosley 1999). In our final scenario (collapsar), we include the mass loss from this powerful engine.

There exist many permutations of our basic convection engines combined with post-explosion engines (e.g., collapsar, magnetar). However, it is anticipated that just a small fraction of systems are affected by the post-explosion engines. To encapsulate the effect of the post-explosion models on the final mass distributions, we assume (an unreasonably high) 100% efficiency in these engines, and study only the four cases described in Table 1.

For ZAMS progenitor stars below ∼11 M_☉, the fate of collapse is the same for all of our explosion models. The fallback in the supernova explosion from these progenitors is negligible. As such, the uncertainties in the shock propagation are minimal. In addition, as we shall see below, the uncertainties in the explosion mechanism only alter the final remnant mass by, at most, 0.1 M_☉. The primary uncertainty in the remnant mass distribution from these stars arises from determining the number of these objects which, in turn, is caused by uncertainties in stellar evolution. In particular, stellar evolution models do not accurately constrain the transition mass between white dwarf (WD) and NS formation (i.e., the lower mass limit for NS formation).

The exact value of the lower limit (M^lower_NS) for core collapse has been a matter of debate for many decades (see Iben & Renzini 1983 for a review). The metallicity dependence for this limit is even more controversial. Heger et al. (2003) argued that the metallicity dependence was negligible, and that the lower limit was roughly at 9 M_☉ for all metallicities. Poelarends et al. (2008) have studied this matter in more detail, and A. J. T. Poelarends et al. (2012, in preparation) have a range of solutions for the metallicity dependence.¹⁰ Their "preferred" model predicts that M^lower_NS drops from 9 M_☉ at solar metallicity down to below 6.3 M_☉ at Z_metal = 10⁻³ Z_☉. We use the following fit to their solution:

$$\begin{eqnarray} M_{\rm NS}^{\rm lower} &=& 9.0\,M_\odot + 0.9 {\rm log_{10}} (Z_{\rm metal}/Z_{\odot }) \,M_\odot \; \nonumber \\ &&{\rm if} \, {\rm log_{10}} (Z_{\rm metal}/Z_{\odot }) > -36.3 \,M_\odot \; {\rm otherwise}. \end{eqnarray} \tag{ 4 }$$

Here we use Z_metal to represent the absolute metallicity. The fraction of NS/BH remnants (assuming a Salpeter initial mass function, IMF) produced by stars below 11 M_☉, as a function of metallicity, is shown in Figure 3. However, we note that Poelarends et al. (2008) found that this result was highly sensitive to their prescriptions for mass loss and dredge-up. Whether the minimum mass varies with metallicity from 6.3 M_☉ to 9 M_☉ (our assumed 2.7 M_☉ variation) or from 8 M_☉ to 9 M_☉ (a 1 M_☉ variation) depends upon these prescriptions.

Zoom In Zoom Out Reset image size

Understanding stellar evolution will allow us to determine how many NSs are formed from stars with initial masses below 11 M_☉. However, the mass of the remnant that is formed in the collapse of these stars is determined by the supernova engine. For these low-mass stars (both electron-capture/ONe supernovae and iron collapse supernovae) the envelope of the star is fairly tenuous (and relatively easily ejected), and the collapsing core is close to the Chandrasekhar limit: ∼1.38 M_☉ for OMgNe stars. The remnant mass is the Chandrasekhar mass minus any ejecta mass. Simulations of the collapse of these objects have found that ECS can eject from 0.01–0.2 M_☉ of stellar material (see Fryer et al. 1999a for a review, or more recent work by Kitaura et al. 2006; Dessart et al. 2007). The favored model by Fryer et al. (1999a) predicted an ejecta mass of roughly 0.1 M_☉, leading to a baryonic mass of the collapsed NS of ∼1.28 M_☉. The recent results by Kitaura et al. (2006) suggest that the ejecta mass is closer to 0.02 M_☉, predicting a final remnant mass of ∼1.36 M_☉. Low-mass iron core collapses are likely to be similar in mass. Although we adopt the former as our default mass distribution, we discuss both in our mass distribution calculations. We expect the remnant mass to be relatively flat for all stars below 11 M_☉.

The remnant masses of stars more massive than ∼11 M_☉ are primarily determined by the amount of material that falls back onto the proto-NS after the launch of the explosion. As the ejecta push out against the surrounding stellar envelope it loses energy, and some of this material falls back onto the compact remnant. The previous piston-driven explosions both delayed and underestimated this fallback for a given supernova explosion (Young & Fryer 2007). Energy-driven explosions more accurately mimic the likely convection-enhanced explosion mechanism, and lead to a solution where most of this fallback occurs in the first 10–20 s (see Fryer 2009 for a review). For a given energy drive, the amount of fallback can be calculated fairly accurately.

The amount of fallback depends upon the energy in the explosion. As mentioned before, here we study four types of engines, producing a wide range of fallback yields. For the initial explosion energies we use either the delayed or rapid explosions discussed in Section 2.1. Figure 4 shows the remnant masses for these two engines using the Woosley et al. (2002) progenitors. The broad range of energies in the delayed explosion produces a continuous range in remnant masses, while the rapid explosion mechanism leads to a sharp transition between NSs and more massive BHs.

Zoom In Zoom Out Reset image size

For rapidly spinning progenitors, the initial fallback can be altered by a post-explosion engine. In the case of BH formation, the rapidly spinning systems will drive further outflow via the collapsar engine. This engine ultimately disrupts the entire star, typically limiting the BH mass (MacFadyen & Woosley 1999) to a maximum of 5 M_☉. Figure 4 shows a mass distribution of the rapid explosion initial model coupled to a collapsar engine.

Alternatively, a rapidly spinning magnetar model can also halt fallback in all progenitors that initially form a NS. Here we assume that this magnetohydrodynamic engine prevents further fallback. As an example, Figure 4 shows a mass distribution of the delayed explosion initial model coupled to a rapidly spinning magnetar engine.

Below ∼25–30 M_☉, the results are fairly insensitive to the metallicity. Figure 5 shows our delayed model results for a range of metallicities, with our best fit to the remnant masses for these metallicities:

$$\begin{eqnarray} M_{\rm remnant,delay} &=& 1.1 + 0.2 e^{(M_{\rm star}-11.0)/4.} \nonumber\\ &&- (2.0+Z_{\rm metal}) e^{0.4(M_{\rm star}-26.0)}, \end{eqnarray} \tag{ 5 }$$

where M_star is the ZAMS mass of the star in solar masses and Z_metal is the metallicity with respect to solar. The corresponding fit for our rapid explosion is (Figure 6)

$$\begin{eqnarray} M_{\rm remnant,rapid} &=& 1.1 + 0.2 e^{(M_{\rm star}-11.0)/7.5} + 10.0 (1.0+Z_{\rm metal}) \nonumber \\ &&\times e^{-(M_{\rm star}-23.5)^2/(1.0+Z_{\rm metal})^2} {\rm if} \,M_{\rm star} < 22\,M_\odot \nonumber \\ &=&M_{\rm remnant,delay} - 1.85 + 0.25 Z_{\rm metal} + 10.0 (1.0 \nonumber \\ &&+Z_{\rm metal}) e^{-(M_{\rm star}-23.5)^2/(1.0+Z_{\rm metal})^2} {\rm otherwise}.\nonumber\\ \end{eqnarray} \tag{ 6 }$$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

There is considerable scatter in the stellar evolution models. Although this may be real, and not numerical artifacts of explosive mixing, we have smoothed out many of these features in our fits. In our fits, we have focused on obtaining the general trends of the models.

Above 30 M_☉, the treatment of mass loss from winds and the wind metallicity dependence become the dominant uncertainty, and differences between the different explosion models fade. Figure 4 shows results for the Woosley et al. (2002) models at both solar and zero metallicity. These models use the Nieuwenhuijzen & de Jager (1990) mass-loss prescription ($\dot{M} \propto L^{1.42} M^{0.16} R^{0.81}$, where L is the luminosity, M is the stellar mass, and R is the stellar radius) that is used by many stellar evolution calculations. A metallicity dependence ($\dot{M} \propto Z_{\rm metal}^{0.5}$) from Kudritzki et al. (1989) is added to this prescription, a common approach in stellar evolution codes (for stars not in the Wolf–Rayet phase). The primary differences among stellar models is the proportionality constant in front of these quantities. This difference translates into a difference in the effective metallicity used by different stellar theorists, i.e., the Woosley et al. (2002) models at 1/3 solar are very similar to the solar metallicity Limongi & Chieffi (2006) models. For stars less than 50 M_☉, the remnant mass is fairly well fit by

$$\begin{eqnarray} M_{\rm remnant, delay} &=& {\rm min}(33.35+(4.75+1.25Z_{\rm metal}) (M_{\rm star}-34),\nonumber\\ &&M_{\rm star} - Z_{\rm metal}^{1/2} (1.3\,M_{\rm star} - 18.35)). \end{eqnarray} \tag{ 7 }$$

For our rapid explosions, the remnant is essentially identical:

$$\begin{equation} M_{\rm remnant, rapid} = M_{\rm remnant, delay} - 1.85 \,+\, Z_{\rm metal} (75-M_{\rm star})/20. \end{equation} \tag{ 8 }$$

At solar metallicity, the fate of stars more massive than 50 M_☉ depends sensitively on the prescription for mass loss. Stellar theorists run very few models at these high masses, and their results vary widely. Fitting the Woosley et al. (2002) and Heger et al. (2003) results at solar metallicities, we get identical results for both our delayed and rapid explosions:

$$\begin{eqnarray} M_{\rm remnant} &=& 1.8 + 0.04\times (90-M_{\rm star}) \quad {\rm if} \, M_{\rm star} < 90\,M_\odot \nonumber \\ &=& 1.8 + {\rm log_{10}}(M_{\rm star}-89) \; {\rm otherwise}. \end{eqnarray} \tag{ 9 }$$

For lower metallicity stars, a reasonable fit is the maximum between Equations (7) and (9). These equations do not include the effect of pulsational instabilities. Heger et al. (2003) argued that below a metallicity of ∼10⁻³ to 10⁻⁴ solar, pulsations eject mass for progenitors lying in the mass range 100–140 M_☉. At this same metallicity, stars with masses lying between 140–260 M_☉ are disrupted entirely in a pair-instability supernova, leaving behind no compact remnant whatsoever. It may be that pair-instability outbursts occur at even higher metallicity. We do not include these pair-instability explosions in this analysis.

These mass distributions do not include any additional explosions produced by magnetars or collapsars. Magnetar outbursts following a supernova explosion will prevent fallback, turning some of our BH remnants into NSs. Collapsars will limit the amount of BH accretion, turning massive BH systems into low-mass systems.

Our mass distribution, combining our analysis of stars above and below 11 M_☉ (described by Equations (1)–(4)), is shown in Figures 3, 5, and 6. Note that the primary effect of metallicity on the NS mass distribution is to alter the number of low-mass NSs. This reflects the changing lower limit for NS formation. BH mass distributions, arising from more massive stars with stronger winds, have a more dramatic dependence on metallicity. These masses are baryonic masses; the gravitational masses are likely to be ∼10% lower (e.g., Strobel & Weigel 2001).

With our initial/remnant mass relation, and assuming a value for the IMF, we can now derive mass distributions for both NSs and BHs as a function of metallicity. These mass distributions for a Salpeter (α = 2.35) IMF are shown in Figure 7. At lower metallicities, the number of low-mass NSs increases. This is because the lower limit for stellar collapse and NS formation decreases with decreasing metallicity. The maximum BH mass increases with decreasing metallicity, as mass loss becomes weaker, leading to more massive stars, and hence more massive remnants.

Zoom In Zoom Out Reset image size

The corresponding distribution versus metallicity for our rapid explosion models is shown in Figure 8. The rapid explosion produces very few compact remnants in the gap region: for the solar metallicity models of Woosley et al. (2002), ≲18% of all remnants above 2 M_☉ lie in the 2–5 M_☉ gap region. At 0.1 solar metallicity, this fraction decreases to less than 1.5%. The Limongi & Chieffi (2006) models, with their lower mass loss, predict this small fraction of remnants at solar metallicity. The BH mass distribution peaks at roughly 7–8 M_☉. On the other hand, the delayed explosion model produces far more remnants in the gap region: 33% at solar metallicity and 17% at 0.1 solar metallicity. If the observed gap represents the true compact object formation mass distribution, the rapid explosion engine is a better match to the observed gap in remnant masses between 2–5 M_☉ (Bailyn et al. 1998; Ozel et al. 2010). To directly compare with these observations, we must focus on X-ray binary systems (see Belczynski et al. 2012).

Zoom In Zoom Out Reset image size

For remnant masses below 11 M_☉, we have two sources of error. First, uncertainties in modeling mass loss and convection (dredge-up) can alter the lower limit for NS formation. We have assumed that at zero metallicity this lower limit drops to 6.3 M_☉, but a limit of 8 M_☉ also fits within the errors of the calculations (A. J. T. Poelarends et al. 2012, in preparation). Second, differences in the results of core-collapse models lead to an error in the mass of the NSs, allowing a range from 1.28–1.36 M_☉. Observations placing constraints on the lowest NS mass could reduce this error.

For remnant masses produced by stars above 11 M_☉, the uncertainties are significantly larger. Our models are based on the assumption that the standard convective engine is the correct engine for most supernovae. We assume the maximum efficiency for ejecting stellar material (i.e., the entire explosion energy is used to unbind the outer layers of the star). This assumption naturally leads to an underestimate of the final remnant mass. However, we also assume that the proto-NS lacks an energy reservoir, despite the fact that there are many mechanisms by which it could inject energy and affect the fallback. For example, proto-NS winds could generate additional energy in the explosion, thereby reducing fallback and decreasing the remnant mass. Even more important, Fryer et al. (1996, 2006) and Fryer (2009) found that fallback material will drive outflows, self-limiting the total amount of fallback. As we discussed, magnetar-strength magnetic fields can also drive outflows, preventing fallback (e.g., Maeda et al. 2007; Kasen & Bildsten 2010). All of these effects would drive our final remnant mass lower, and ignoring them overestimates the remnant mass.

In addition, the mass distribution is sensitive both to the prescription used for mass loss in the stellar evolution codes and the stability of the stellar evolution codes. Figure 9 compares the results from three different sets of stellar models using three different stellar evolution codes: Woosley et al. (2002) solar metallicity (solid line), Limongi & Chieffi (2006) solar metallicity (dotted), and the binary models of Young et al. (2009) at solar metallicity (open squares).

Zoom In Zoom Out Reset image size

Different stellar models utilize different prescriptions for mass loss, accounting for the huge differences between the Limongi & Chieffi (2006) and Woosley et al. (2002) results above 30 M_☉. Much of this difference is just a scaling factor: The primary difference in the mass-loss prescription is the coefficient in front of the metallicity-dependent term. For example, the Limongi & Chieffi (2006) solar metallicity models are very similar to the Woosley et al. (2002) models at 1/3 solar metallicity. To determine the remnant mass distribution of the Limongi & Chieffi (2006) models, we would use our remnant mass prescriptions (e.g., Equations (7) and (8)) at 1/3 solar metallicity.

The jagged nature of our calculated mass estimate comes from instabilities in the mixing within the stellar evolution code. Because of this, the core masses are not a smooth function of progenitor mass. The high remnant mass at 23 M_☉ in our rapid explosion models in Figure 6 is a demonstration of how important this mixing can be. It is possible that better convection models will reduce the non-smooth dependence of core mass on progenitor mass (see Young et al. 2005), but this remains to be established.

We have also neglected the effects of binaries in the discussion above. Because fallback occurs prior to the shock hitting the hydrogen envelope, and because the hydrogen envelope has very little binding energy, the primary effect of binaries on the mass distribution will be to reduce the maximum mass of the remnant to that of the helium core mass. Fryer & Kalogera (2001) have studied this effect and find that it is negligible in the case of NSs. However, the broader remnant mass distribution is sensitive to binary effects, and we study these next (Section 4).

Here we present a brief description of the computation of compact object masses in the StarTrack population synthesis code (for the full description see Belczynski et al. 2008, 2010a).

We determine the mass of a NS/BH remnant using information on the final CO mass M_CO combined with the knowledge of the pre-supernova mass of the star M. For a given initial ZAMS mass, the final CO core mass is obtained from the original Hurley et al. (2000) formula, while we use the models of Timmes et al. (1996) with the addition of Si shell mass to estimate final FeNi core mass (which we here refer to as proto-compact object mass). The proto-compact object mass is obtained from

$$\begin{eqnarray} \left\lbrace \begin{array}{@{}ll} M_{\rm proto}= 1.50\,M_{\odot }& M_{\rm CO} < 4.82\,M_{\odot }\\ M_{\rm proto}= 2.11\,M_{\odot }& 4.82 \le M_{\rm CO} < 6.31\,M_{\odot }\\ M_{\rm proto}= 0.69 M_{\rm CO} - 2.26\,M_{\odot }& 6.31 \le M_{\rm CO} < 6.75\,M_{\odot }\\ M_{\rm proto}= 0.37 M_{\rm CO} - 0.07\,M_{\odot }& M_{\rm CO} \ge 6.75\,M_{\odot }. \end{array} \right. \nonumber\\ \end{eqnarray} \tag{ 10 }$$

The fallback of material after the launch of the explosions adds mass to the remnant. To calculate the amount of fallback, M_fb, for a given core mass, M_CO, we employ

$$\begin{eqnarray} \left\lbrace \begin{array}{@{}ll} M_{\rm fb} = 0\,M_{\odot }& M_{\rm CO} < 5.0\,M_{\odot }\\ f_{\rm fb} = 0.378 M_{\rm CO}-1.889\,M_{\odot }& 5.0 \le M_{\rm CO} < 7.6\,M_{\odot }\\ f_{\rm fb} = 1.0 & M_{\rm CO} \ge 7.6\,M_{\odot }\end{array} \right. \nonumber\\ \end{eqnarray} \tag{ 11 }$$

with M_fb = f_fb(M − M_proto) in the mass range for which the fractional fallback f_fb is given.¹²

The final remnant (baryonic) mass is calculated from

$$\begin{equation} M_{\rm rem,bar} = M_{\rm proto} + M_{\rm fb}, \end{equation} \tag{ 12 }$$

and we convert baryonic to gravitational mass (M_rem) using

$$\begin{equation} M_{\rm rem,bar} - M_{\rm rem} = 0.075\ M_{\rm rem}^2 \end{equation} \tag{ 13 }$$

for NSs (Lattimer & Yahil 1989; see also Timmes et al. 1996), while for BHs we simply approximate the gravitational mass with

$$\begin{equation} M_{\rm rem} = 0.9\ M_{\rm rem,bar}. \end{equation} \tag{ 14 }$$

To calculate the final mass of a compact object we need to know the mass of a star, M, and its CO core mass, M_CO, at the time of core-collapse/supernova explosion. For the rapid explosion mechanism, the explosion either occurs quickly or not at all. We set the proto-compact object mass

$$\begin{equation} M_{\rm proto} = 1.0\,M_{\odot }, \end{equation} \tag{ 15 }$$

independent of exploding star mass according to hydrodynamical simulations of supernova explosions (Woosley et al. 2002). Depending on the amount of mass above the proto-compact object (M − M_proto) and the strength of the explosion, potential fallback may increase the mass of the compact object. We calculate the amount of fallback, M_fb, for a given core mass, M_CO:

$$\begin{eqnarray} \left\lbrace \begin{array}{@{}ll} M_{\rm fb} = 0.2\,M_{\odot }& M_{\rm CO} < 2.5\,M_{\odot }\\ M_{\rm fb} = 0.286 M_{\rm CO}-0.514\,M_{\odot }& 2.5\,M_{\odot } \le M_{\rm CO} < 6.0\,M_{\odot }\\ f_{\rm fb} = 1.0 & 6.0\,M_{\odot } \le M_{\rm CO} < 7.0\,M_{\odot }\\ f_{\rm fb} = a_1 M_{\rm CO} + b_1 & 7.0\,M_{\odot } \le M_{\rm CO} < 11.0\,M_{\odot }\\ f_{\rm fb} = 1.0 & M_{\rm CO} \ge 11.0\,M_{\odot }\end{array} \right. \nonumber\\ \end{eqnarray} \tag{ 16 }$$

with a₁ = 0.25 − (1.275/M − M_proto), b₁ = −11a₁ + 1, and M_fb = f_fb(M − M_proto) in the mass range for which f_fb is given.

The final baryonic mass of the remnant is then

$$\begin{equation} M_{\rm rem,bar} = M_{\rm proto} + M_{\rm fb} \end{equation} \tag{ 17 }$$

and the gravitational mass of the remnant is obtained from Equations (13) and (14).

The calculation of the final mass of a compact object for the delayed mechanism is similar to the rapid supernova case detailed above. The formulae are based on the delayed supernova calculations discussed in Sections 3.2 and 3.3. Again, for binary population synthesis, we use the CO core mass employed in most population synthesis calculations. First, we calculate the proto-compact object mass

$$\begin{equation} \left\lbrace \begin{array}{@{}ll} M_{\rm proto} = 1.2\,M_{\odot }& M_{\rm CO} < 3.5\,M_{\odot }\\ M_{\rm proto} = 1.3\,M_{\odot }& 3.5\,M_{\odot } \le M_{\rm CO} < 6.0\,M_{\odot }\\ M_{\rm proto} = 1.4\,M_{\odot }& 6.0\,M_{\odot } \le M_{\rm CO} < 11.0\,M_{\odot }\\ M_{\rm proto} = 1.6\,M_{\odot }& M_{\rm CO} \ge 11.0\,M_{\odot }. \end{array} \right. \end{equation} \tag{ 18 }$$

This assumes that the delay increases for more massive cores, causing more material to accrete onto the proto-NS during the explosion.

The amount of fallback is given by

$$\begin{eqnarray} \left\lbrace \begin{array}{@{}ll} M_{\rm fb} = 0.2\,M_{\odot }& M_{\rm CO} < 2.5\,M_{\odot }\\ M_{\rm fb} = 0.5 M_{\rm CO}-1.05\,M_{\odot }& 2.5\,M_{\odot } \le M_{\rm CO} < 3.5\,M_{\odot }\\ f_{\rm fb} = a_2 M_{\rm CO} + b_2 & 3.5\,M_{\odot } \le M_{\rm CO} < 11.0\,M_{\odot }\\ f_{\rm fb} = 1.0 & M_{\rm CO} \ge 11.0\,M_{\odot }\end{array} \right. \nonumber\\ \end{eqnarray} \tag{ 19 }$$

with a₂ = 0.133 − (0.093/M − M_proto), b₂ = −11a₂ + 1, and M_fb = f_fb(M − M_proto) in the mass range for which f_fb is given.

The final baryonic mass of the remnant is then

$$\begin{equation} M_{\rm rem,bar} = M_{\rm proto} + M_{\rm fb}, \end{equation} \tag{ 20 }$$

and the gravitational mass of the remnant is obtained from Equations (13) and (14).

In our calculations for the masses of compact objects, as detailed above, we allow for NS formation through ECS (e.g., Podsiadlowski et al. 2004), following the approach in Belczynski et al. (2008). From Hurley et al. (2000), we use the He core mass at the asymptotic giant branch base to set the limit for the formation of various CO cores. If the He core mass is smaller than M_cbur1, the star forms a degenerate CO core, and ends up forming a CO WD. If the core is more massive than M_cbur2 = 2.25 M_☉, the star forms a non-degenerate CO core with subsequent burning of elements until the formation of an FeNi core, which ultimately collapses to a NS or a BH. Stars with cores between M_cbur1 and M_cbur2 may form partially degenerate CO cores. If such a core reaches a critical mass (M_{co, crit} = 1.08 M_☉; Hurley et al. 2000), it ignites CO off-center and non-explosively burns CO into ONe, forming a degenerate ONe core. If in subsequent evolution the ONe core increases its mass to M_ecs = 1.38 M_☉, the core collapses due to electron capture on Mg and forms a NS. We refer to these as ECS NS, to distinguish them from NS from regular iron core collapse. The ECS NSs are assumed to have unique masses of M_{rem, bar} = M_ecs. If the ONe core mass remains below M_ecs, the star forms an ONe WD.

Hurley et al. (2000) suggested M_cbur1 = 1.66 M_☉ corresponding to M_zams = 6.5 M_☉ for Z = 0.02. Later calculations with an updated evolutionary code (Eldridge & Tout 2004a, 2004b) indicated that ECS may occur for higher initial masses (M_zams ≳ 7.5 M_☉). For our standard model we adopt M_cbur1 = 1.83 M_☉ (M_zams = 7.0 M_☉), and this results in ECS NS formation above M_zams = 7.6 M_☉ (for masses M_zams = 7.0–7.6 M_☉, the ONe core does not reach M_ecs and an ONe WD is formed). It is noted that binary evolution through Roche-lobe overflow may either decrease the initial mass of the ZAMS star required to form the core of mass M_cbur1 (due to rejuvenation) or increase it (due to mass loss). Therefore, binary evolution effectively leads to wider initial progenitor masses for ECS NS formation. Metallicity and wind mass loss may also influence the ECS NS formation range.

At the time of birth, NSs and BHs receive a natal kick due to asymmetries in the supernova explosions. We use the distribution inferred from observed velocities of radio pulsars by Hobbs et al. (2005): a single Maxwellian with σ = 265 km s⁻¹.

Compact objects formed without fallback receive full kicks drawn from the above distribution. This is the case for most NSs, with the exception of ECS NSs, for which we adopt either no natal kicks (standard model) or full kicks (to explore a range of parameters). The same approach is applied to NSs formed via accretion induced collapse of WD to NS in an accreting binary system.

Additionally, compact objects formed with small amounts of fallback (M_fb < M_{fb, small}) receive full kicks (this includes some low-mass BHs). In our standard approach, we adopt M_{fb, small} = 0.0 M_☉ for the StarTrack scheme and M_{fb, small} = 0.2 M_☉ for the rapid and delayed supernovae. However, to fully explore parameter space, we also consider full kicks with a broad fallback mass range (as high as M_{fb, small} ∼ 1 M_☉).

For compact objects formed with noticeable fallback (M_fb ⩾ M_{fb, small}), kicks are lowered proportional to the amount of fallback:

$$\begin{equation} V_{\rm kick} = (1-f_{\rm fb}) \sqrt{V_{\rm x}^2 + V_{\rm y}^2 + V_{\rm z}^2}, \end{equation} \tag{ 21 }$$

where V_x, V_y, V_z are the three velocity components drawn from Hobbs et al. (2005) distribution, and f_fb is the fraction (from 0 to 1) of the stellar envelope that falls back.

For the most massive BHs, which are formed silently (without a supernova explosion) via direct collapse (f_fb = 1) of a massive star, we assume that no natal kicks are imparted. In this formalism, the kick is solely a function of the fallback fraction (f_fb), and not of the mass of the compact remnant.

The mass distributions from our prescriptions for remnant masses at solar metallicity are shown in Figure 10. The dashed line shows the results from a standard StarTrack calculation (e.g., Belczynski et al. 2010a), the solid line shows the results from our delayed engine, and the dot-dashed line shows the results for the rapid engine. The rapid explosion will produce many more BHs of ∼8 M_☉ than the delayed explosion. For intermediate-mass stars that produce BHs (M_zams ≈ 40–80 M_☉), the StarTrack BHs are more massive than BHs predicted in more recent simulations (for both delayed and rapid models).

Zoom In Zoom Out Reset image size

Most interesting are the differences between these models in producing compact remnants below 5 M_☉. The rapid explosion mechanism produces a range of compact remnant masses up to 2 M_☉, but makes very few remnants of mass in the range ∼2–5 M_☉. The StarTrack models partially fill this gap, while the delayed explosion models contribute a significant population of remnants in the gap. These trends persist at lower metallicities (Figure 11).

Zoom In Zoom Out Reset image size

This striking gap in compact object mass (∼2–5 M_☉) results from the sensitivity of the rapid models to the structure of the core. In the Woosley et al. (2002) and Woosley & Heger (2007) models, enhanced burning in the late evolutionary stages leads to sharply different internal structures of stars with mass M_zams ≈ 20–25 M_☉. The internal structure in turn affects the supernova explosion and dramatically changes the mass of the remnant formed in the collapse. The changes are such that in this mass range we get an abrupt transition of compact object mass from 2 M_☉ to much higher mass, with no compact objects in between (see Figures 10 and 11). This result agrees with studies of O'Connor & Ott (2011), who find that the bounce compactness changes dramatically for these stars, altering their fate. The ultimate fate of these stars is sensitive to the late-stage burning, and if the burning changes with newer stellar evolution codes, the remnant mass distribution will change commensurately.

One of the strongest observational constraints is the gap in BH mass between 2 and 5 M_☉ (Farr et al. 2011; Ozel et al. 2010). However, these studies do not include systems which do fall in the gap, such as 4U1700-37 (Clark et al. 2002) with a measured mass of 2.44 ± 0.27 M_☉. We expect some remnants to fall in the gap, but have difficulty estimating the number. Our model predictions range from ∼15% for our rapid collapse to ∼40% for our delayed model in the 2–4 M_☉ mass range (Section 4.4). No model can predict a complete absence of BHs in the 2–4 M_☉ mass range. Observations of IGR J17091-3624 (Altamirano et al. 2011) suggest that some BHs exist in this range, but these systems appear to be rare.

The high space velocity of 4U1700-37 points out another difficulty in using X-ray binary mass observations to constrain the remnant mass of supernovae: BH and NS kicks. Many of the proposed NS kick mechanisms will also impart kicks onto BHs at formation. To impart a kick on the compact remnant from stellar collapse, momentum conservation requires asymmetries. Most compact remnant kick mechanisms either argue for asymmetries in the baryonic ejecta from the supernova or in the neutrino emission. Here we review the implications for compact remnant kicks for our two explosion models utilizing these two kick mechanisms.

For ejecta-driven kicks, no kick is imparted on BHs formed without supernova explosions (since there is then no ejecta). Low-mass BHs, on the other hand, are formed when the supernova explosion is delayed and the explosion is weak. If the convection in the supernova engine develops low-mode instabilities, the resulting explosion can be asymmetric (Herant et al. 1992; Buras et al. 2003; Blondin et al. 2003). The SASI engine is characterized by low-mode convection and several groups have argued that these large asymmetries will drive strong kicks (Buras et al. 2003; Blondin et al. 2003; Scheck et al. 2006; Wongwathanarat et al. 2010); see Fryer & Young (2007) for a different interpretation of the low-mode convection. If the low-mode convection in the SASI engine produces strong kicks, the explosions with the longest delays will have the largest ejecta asymmetries and, hence, the strongest kick momenta. But these delayed supernova engines also produce the weakest supernova explosions, forming massive NSs and low-mass BHs. We thus expect massive NSs and low-mass BHs to have the highest momentum "kick" of any compact remnant.

For the rapid explosion model most BHs are formed without a supernova explosion, and only a small fraction of low-mass BHs will receive kicks. This is to be contrasted with the delayed model, where a much larger fraction of BHs are formed from fallback, and these low-mass BHs will have large kicks.

Likewise, for most asymmetric neutrino kick mechanisms, the longer the delay in the supernova engine is, the stronger the momentum asymmetry will be (Socrates et al. 2005; Kusenko 2004; Fryer & Kusenko 2006). Again, massive NSs and low-mass BHs receive the largest kicks of any of our compact remnants under our delayed explosion model.

Strong kicks unbind binary systems. If strong kicks occur in high-mass NSs and low-mass BHs, these systems are more likely to be disrupted in the massive star collapse. This bias will make it difficult to use remnant mass distributions derived from X-ray binary systems without detailed population synthesis models. We defer this study to a later paper (Belczynski et al. 2012).

Although our focus has been on estimating the compact remnant mass distribution, our models also make predictions for the distribution of explosion energies for core-collapse supernovae. The explosion energy distributions for our two convection-enhanced engines at both solar and zero metallicity are shown in Figure 12. We have assumed that core-collapse supernovae from low-mass stars and ECS all produce the same low-energy, 0.55 × 10⁵¹ erg explosion. Since the minimum mass decreases with metallicity, these supernovae dominate the explosions at low metallicity, even with the relatively shallow Salpeter IMF (α = 2.35).

Zoom In Zoom Out Reset image size

From Figure 12 we note that, except for a low-mass supernovae such as ECS, our rapid explosion model produces only strong explosions (above 10⁵¹ erg). The delayed model produces a broader range of explosion energies, with peak energies below that of the rapid explosion model. Even at solar metallicity, less than 50% of the explosions from the delayed engine have energies in excess of 10⁵¹ erg. In principal, the distribution of supernova explosion energies, if sufficiently accurate, could distinguish between our two engines. For example, a low-energy explosion from a massive star is almost certainly arising from a delayed explosion. Similarly, any explosion above 1.5 × 10⁵¹ erg produced by a slowly rotating star is likely to result from a rapid explosion.

It is likely that a combination of both the delayed and rapid explosion engines will be required to explain supernovae, and future observations may determine the relative fraction.

Gravitational waves are expected to provide new constraints on the formation of compact remnants. We use the StarTrack population synthesis code, incorporating the binary prescriptions discussed above, to calculate delayed (e.g., SASI) and rapid models for the initial–remnant mass relation. Additionally, we employ two models of common envelope evolution. In the first (model A), we employ the standard energy balance (Webbink 1984) to calculate the outcome of common envelope evolution with a Hertzsprung gap donor. In the second model (B) we assume that, when a Hertzsprung gap star fills its Roche lobe and starts the common envelope phase, it merges with its companion and ends its binary evolution. The same energy balance is used to compute the common envelope for any type of donor in both models. Since many BH–BH progenitors experience common envelope evolution with a Hertzsprung gap donor, the birth and merger rates of BH–BH binaries are greatly reduced in model B (for details and the physical motivation see Belczynski et al. 2007). We consider two metallicities Z_☉ = 0.02 and 0.1 Z_☉ = 0.002 (approximately covering the range of metallicity observed in local universe). We focus on binary stars that form double BH systems (BH–BH) at the end of their evolution.

We predict the distribution of chirp masses for the two supernova engine models (see Figures 13 and 14), and find a striking difference. For the delayed model, the chirp masses start from M_chirp ≈ 2.4(2.5) M_☉, whereas for the rapid model, the lowest chirp masses are ≈4.9(5.4) M_☉ for models A(B). These very different chirp mass distributions are the result of the differing initial–remnant mass relations presented in Figures 10 and 11. In the delayed model BHs start forming from M_BH = 2.5 M_☉, and their mass continuously increases as the initial mass of the progenitors increases. The chirp mass distribution starts at a low value, M_chirp ≈ 2.4 M_☉ (the chirp mass of two 2.5 M_☉ BHs is 2.2 M_☉). In contrast, in the rapid supernova engine model, there are no BHs with mass below 5 M_☉, and therefore the chirp mass distribution starts at M_chirp ≈ 5 M_☉. The origin of this remnant gap is described and referenced in Section 4.6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Table 2 we present Galactic merger rates for our models for BH–BH binaries. Rates are calculated for a Milky Way type galaxy (10 Gyr of continuous star formation at a rate of 3.5 M_☉ yr⁻¹), as described in Belczynski et al. (2007). The change of rates with metallicity, and between our two models for common envelope evolution, is fully discussed in Belczynski et al. (2010b). Here we note that both supernova models produce comparable numbers of close BH–BH binaries, hence the similar merger rates in Table 2. Therefore, observations of BH–BH binaries will illuminate the underlying supernova engine only if the physical properties of the mergers can be extracted from the gravitational waves. If the gap in the mass distribution from X-ray binaries is confirmed by gravitational wave observations of merging binaries, this gap will rule out delayed supernova mechanisms for most systems.