Modeling the outcome of supernova explosions in binary population synthesis using the stellar compactness

Following the collapse of their cores, some of the massive binary stars that populate our Universe are expected to form merging binaries composed of black holes and neutron stars. Gravitational-wave observations of the resulting compact binaries can reveal precious details on the inner workings of the supernova mechanism and the subsequent formation of compact objects. Within the framework of the population-synthesis code MOBSE, we present the implementation of a new supernova model that relies on the compactness of the collapsing star. The model has two free parameters, namely the compactness threshold that separates the formation of black holes and that of neutron stars, and the fraction of the envelope that falls back onto the newly formed black holes. We compare this model extensively against other prescriptions that are commonly used in binary population synthesis. We find that the cleanest signatures of the role of the pre-supernova stellar compactness are (i) the relative formation rates of the different kinds of compact binaries, which mainly depend on the compactness threshold parameter, and (ii) the location of the upper edge of the mass gap between the lightest black holes and the heaviest neutron stars, which mainly depends on the fallback fraction.


Introduction
Gravitational wave (GW) observations of merging compact binaries offer unprecedented insights into the life of massive stars. The black holes (BHs) and neutron stars (NSs) observed by LIGO and Virgo (Abbott et al., 2019(Abbott et al., , 2021 constitute the end product of stellar collapse -the same cosmic events that are well understood to be behind supernova (SN) explosions. The direct observation of compact binaries thus provides a novel opportunity to probe the inner workings of the SN explosion mechanism.
Stars with zero-age-main-sequence (ZAMS) masses greater than ∼ 10 M are robustly predicted to undergo core-collapse SN once their cores reach the Chandrasekhar limit (for a recent review on SN theory see Burrows & Vartanyan 2021). Lighter stars might instead explode as electron-capture SN. In those cases, the degenerate oxygen-neon core reaches the critical mass of ∼ 1.38 M and electron-capture reactions destabilize the inner region (Miyaji et al., 1980;Nomoto, 1984). While NSs are the only possible outcome of electron-capture SNe, core-collapse SNe can also produce BHs, with the mass of the resulting compact remnant increasing with the ZAMS mass overall. For progenitors heavier than ∼ 70 M , however, the collapsing core becomes unstable to pair production. This removes radiation pressure support from the star which, in turn, ignites explosive carbon-oxygen (CO) burning. The core is either partially (pulsational pair-instability SN; Woosley et al. 2007) or entirely (pair-instability SN;Heger et al. 2003) disrupted, thus introducing a characteristic upper limit of ∼ 50M (Farmer et al., 2019;Woosley & Heger, 2021) to the masses of BH remnants that can be produced by conventional stellar evolution.
The key GW observables to probe the physics of SN are the so-called gaps in the BH mass spectrum. The existence or absence 1 School of Physics and Astronomy & Institute for Gravitational Wave Astronomy, University of Birmingham, Birmingham, B15 2TT, UK. * mxd773@star.sr.bham.ac.uk of compact objects with masses between ∼ 3M and ∼ 5M ("lower mass gap") separating BHs and NSs has been long investigated using X-ray binary data (Bailyn et al., 1998;Farr et al., 2011;Özel et al., 2010) and is now under active scrutiny after some of the most recent LIGO/Virgo events, notably GW190814 (Abbott et al., 2021). At the high-end of the BH mass spectrum ("upper mass gap"), current GW measurements point to a significant drop in the merger rate which is consistent with prediction from pair-instability SN theory. At the same time, some of the observed BH masses, notably those of GW190521 (Abbott et al., 2021), are well within this forbidden region, perhaps hinting at a different origin for this event (for a review see Gerosa & Fishbach 2021). Besides the masses, the orientations of the BH spins can also provide insights on SN physics. In particular, they are sensitive to the asymmetric emission of mass and neutrinos occurring during the SN and the subsequent kick imparted to the newly formed compact object (e.g. Kalogera 2000;Gerosa et al. 2013;O'Shaughnessy et al. 2017).
In this paper, we consider the standard formation scenario where GW sources originate from massive binary stars (e.g. Postnov & Yungelson 2014) and explore the impact of a new SN model on the resulting compact binaries. The most common prescriptions used to predict the SN outcome implemented in state-of-the-art population-synthesis codes are based on the CO mass after the carbon burning stage and the pre-supernova total mass of the progenitor star (e.g. SEBA, Portegies Zwart & Verbunt 1996;STARTRACK, Belczynski et al. 2020;COSMIC, Breivik et al. 2020;MOBSE, Giacobbo et al. 2018;COMPAS, Stevenson et al. 2017).
In particular, the two leading SN models are the so-called rapid and delayed prescriptions first proposed by Fryer et al. (2012). The main difference between these two models is the time after core bounce at which the explosion is launched, which is set to < 250 ms and > 500 ms in the rapid and delayed case, respectively. In terms of the properties of the resulting BH bina-ries, the rapid (delayed) model does (does not) predict the lower mass gap between BHs and NSs (e.g. Belczynski et al. 2012). Additional recipes to model the SN outcome include the effect of envelope stripping in binaries (Schneider et al., 2021) and probabilistic couplings between natal kicks and remnant masses (Mandel & Müller, 2020).
Hydrodynamical simulations (e.g. O'Connor & Ott 2011) suggest that, for a given equation of state, the most critical parameter to estimate the outcome of a SN is given by the compactness of the stellar core at the onset of the explosion, where R(M) is the radius that encloses a mass M at core bounce. We explore the impact of such "compactness model" of compact-object formation on the resulting single and binary mass spectrum by means of population-synthesis simulations. In Sec. 2, we describe an updated version of the MOBSE code where we have implemented a treatment of SNe based on the stellar compactness. In Sec. 3, we discuss the effect of the compactness model on the formation and evolution of merging compact binaries and compare our results against the more common rapid and delayed prescriptions. In Sec. 4, we summarize our findings and illustrate prospects for future work.

Methods
After a brief introduction to the MOBSE code (Sec. 2.1), we describe the implementation of our compactness model (Sec. 2.2) and the simulation setup used to explore its relevance (Sec. 2.3).

Massive Objects in Binary Stellar Evolution
Both single and binary stars are evolved with the rapid population-synthesis code MOBSE ("Massive Objects in Binary Stellar Evolution", . Built on top of the older BSE (Hurley et al., 2002), MOBSE's key updates include up-todate prescriptions for stellar winds in massive stars , natal kicks , and modeling of pulsational pair-instability and pair-instability SNe (Spera & Mapelli, 2017). The version of MOBSE presented in this paper has been assigned number v1.1. The code is publicly available at mobse-webpage.netlify.app.

Compactness model
While the prescriptions for the rapid and delayed models are straightforward to implement in a population-synthesis code, the compactness parameter ζ M of Eq. (1) depends on the internal structure of the star before the SN. This information is typically not available with a population-synthesis approach, which does not track the details of stellar structure but only some of the star's global properties.
As in O'Connor & Ott (2011), we will refer to the reference value of ζ at M = 2.5M , resulting in the parameter ζ 2.5 . Limongi & Chieffi (2018) demonstrated that there is a strong correlation between ζ 2.5 and the carbon-oxygen mass M CO of the pre-SN, which is a readily accessible quantity. More recently, Mapelli et al. (2020) showed that this relation can be well represented by the simple fit (2) We make use of ζ 2.5 to distinguish between the formation of a NS and that of a BH. The precise threshold value is still debated, with current estimates ranging from ∼ 0.2 (Horiuchi et al., 2014) to ∼ 0.45 (O'Connor & Ott, 2011). In our fiducial model, we assume a treshold value of ζ 2.5 = 0.365 which was chosen to match the NS-to-BH transition predicted by the rapid model (cf. Sec. 3.1). Hence, stellar progenitors with ζ 2.5 > 0.365 will form BHs while progenitors with ζ 2.5 ≤ 0.365 will form NSs.
The compactness itself, however, does not predict the mass of the resulting compact remnant. We thus adopt the same strategy of Mapelli et al. (2020). NS masses are drawn from a Gaussian distribution with mean µ = 1.33 M and standard deviation σ = 0.09 M , which agrees with the observed masses of NSs in binary systems (Özel & Freire, 2016). BH masses are instead given by where M He is the mass of the helium core and M fin is the total stellar mass at the onset of collapse (both of which are tracked in standard population synthesis). The quantity f H ∈ [0, 1] is a free parameter that describes the fraction of the hydrogen envelope that is accreted by the BH. We adopt a fiducial value f H = 0.9, i.e. assuming that 90% of the star's hydrogen envelope falls back onto the remnant after core collapse. This fiducial value is motivated by recent studies (e.g. Fernández et al., 2018;Lovegrove & Woosley, 2013;Sukhbold et al., 2016) that show that not all of the hydrogen envelope is accreted onto the newborn BH, even in the case of direct collapse. This is because the outermost layers of the envelope are only weakly bound to the core. Figure 1 illustrates the relationship between the initial mass of the star M ZAMS and the mass of the compact remnant M rem obtained with our fiducial compactness model (ζ 2.5 = 0.365 and f H = 0.9). Much like in the rapid case, the compactness model also produces a gap in the remnant masses which separates BHs and NSs. As already explored at length (Belczynski et al., 2010;Neijssel et al., 2019;Spera et al., 2015), both the maximum BH mass and the upper edge of the mass gap strongly depends on the progenitor's metallicity Z. This is because higher metallicities drive larger mass loss via stellar winds, thus leading to smaller remnant masses. The peak at M ZAMS ∼ 70M for Z 0.004 is due to pulsational pairinstability SNe while the sharp decrease at M ZAMS ∼ 140M marks the onset of the proper pair-instability regime.

Simulation setup
Starting from the fiducial model we just described, we further explore the parameter space spanned by ζ 2.5 and f H . We consider four different values for the threshold parameter ζ 2.5 = 0.2, 0.3, 0.365 (fiducial), and 0.4. For each of these, we vary f H = 0.1, 0.3, 0.5, 0.7, and 0.9 (fiducial). For all the resulting combinations of ζ 2.5 and f H we evolve 10 6 massive binaries varying their metallicity Z = 0.0002, 0.0004, 0.0008, 0.0012, 0.0016, 0.002, 0.004, 0.006, 0.008, 0.012, 0.016, and 0.02. We also run control cases using the rapid and delayed models. The details of their precise implementation in MOBSE have been presented by Santoliquido et al. (2021) and include some minor modifications with the respect to the original prescriptions by Fryer et al. (2012).
The initial condition of the simulated binaries are generated as follows : primary ZAMS masses are extracted from an initial mass function (IMF)

Results
We now present the results of our simulations. We first discuss the impact of the compactness SN model on the mass spectrum of compact object (Sec. 3.1). We then illustrate the formation of compact object binaries (Sec. 3.2) and focus on the sub-sample of merging systems (Sec. 3.3).

Mass spectrum
In Fig. 2, we compare the mass spectrum for our fiducial compactness model (ζ 2.5 = 0.365, f H = 0.9) against those obtained with the more standard rapid and delayed models. Both the rapid and the compactness models produce a gap in the remnant masses. However, while in the rapid model the mass gap ranges from ∼ 2M to ∼ 5M independently of the metallicity, in the compactness model the size of the mass gap strongly depends on both Z and f H . In particular, the shaded areas in Fig. 2 indicate the range of predicted masses by compactness models with 0.1 ≤ f H ≤ 0.9. The impact of f H is largest at lower metallicities. This is because stars at higher metallicities tend to end their life as Wolf-Rayet stars that have light hydrogen envelopes, so the mass of the resulting BH is almost independent of f H . On the other hand, at lower metallicities stellar winds are less efficient and stars approach the SN phase with heavier envelopes.
Furthermore, the compactness model does not predict a sharp decrease in M rem at 20M M ZAMS 30M which instead characterizes the rapid model. With some dependency on metallicity and fallback, the compactness model tends to, overall, predict heavier remnants in that intermediate region (while for the rapid model those stars do not undergo direct collapse, Fryer et al. 2012). At the high-end of the mass spectrum M zams 50M (but the precise value depends on Z), the rapid, delayed and fiducial compactness model all predict the forma-tion of BHs via direct collapse and thus return similar remnant masses.
The four panels of Fig. 3 illustrate the impact of ζ 2.5 on the mass spectrum. As ζ 2.5 is reduced, the NS-to-BH transition moves toward lower ZAMS masses, i.e. BHs are formed from lighter progenitors. This has important consequences for the relative abundances of NSs and BHs (cf. Sec. 3.2). This is also sensitive to the metallicity. At Z = 0.02 (Z = 0.0002) the ZAMS mass transition value changes from ∼ 14M (∼ 13M ) for ζ 2.5 = 0.2 to ∼ 37M (∼ 23M ) for ζ 2.5 = 0.4.
Finally, let us notice that f H has a strong impact on the maximum BH mass. Smaller values of f H lead to lighter BHs by construction, because the parameter f H describes the amount of fallback material. More interestingly, for massive progenitors ( 60 M ) and lower values of f H ∼ 0.1 , the dependence of the maximum BH mass on the metallicity is not monotonic. In particular, the heaviest BHs are formed at intermediate metallicities Z ∼ 0.006.

Compact binaries
From our samples of simulated massive stars (see Sec. 2.3), we select the systems that form binaries composed of compact objects and classify them as binary black holes (BBHs), black-hole neutron-star binaries (BHNSs), and binary neutron stars (BNSs). The solid lines in Fig. 4 shows the resulting distribution of their chirp masses M chirp = (M 1 M 2 ) 3/5 /(M 1 + M 2 ) 1/5 (where here M 1,2 are the masses of the two compact objects). This is the mass combination parameter that is measured best in GW observations (e.g. Thorne 1987). Larger values of ζ 2.5 result in a narrower mass range for both BBHs and BHNSs. For BBHs, the minimum chirp mass ranges from ∼ 2.5 M for ζ 2.5 = 0.2 to ∼ 7 M for ζ 2.5 = 0.4. In particular, for ζ 2.5 = 0.365 (ζ 2.5 = 0.2) our compactness model matches the range predicted by the rapid (delayed) prescription. The observationally based NS mass prescription we implemented (cf. Sec. 2.2) does not depend on ζ 2.5 . Consequently, the effect of ζ 2.5 on the BHNS chirp masses is mitigated compared to BBHs, while the masses of BNSs are entirely unaffected.  Figure 4 shows results for our two most extreme cases of fallback retention, f H = 0.1 and 0.9, for each type of compact binaries and values of ζ 2.5 . As already highlighted in Sec. 2.2, the impact of f H on the remnant masses depends on both Z and ζ 2.5 . In particular, the mass spectrum of lighter stars (M ZAMS 40) presents a shallower slope for higher f H . This translates ito well defined peaks in the chirp masses of BBHs for f H = 0.9. At low metallicities, the chirp-mass distribution of both BBHs and BHNSs shows a shortage of the most massive systems for f H = 0.1 compared to f H = 0.9. This is a direct consequence of the fact that, for larger ZAMS masses, the parameter f H only affects BH masses at low metallicities.

Merging systems
Only a fraction of the compact binaries that form will merge within a Hubble time (here taken to be 13.6 Gyr). In Fig. 4, we compare the chirp masses of such merging systems against those of the entire simulated samples. In agreement with  ble time and that the fraction of merging systems decreases for increasing metallicities. Figure 5 further illustrates the total number of merging systems as a function of metallicity, separating the different kinds of compact binaries (BBHs, BHNSs, and BNSs). The total number of merging BBHs and merging BHNSs strongly depends on the metallicity of the progenitors. In particular, most models show a monotonic trend that predicts more mergers at lower metallicities. The BHNSs formed in models with ζ 2.5 0.3 are an exception and tend to have a mild peak at intermediate metallicities. In contrast, the metallicity has a weaker impact on the BNSs rate, with fewer mergers predicted at intermediate metal- The number of merging BNSs tends to slightly decrease as ζ 2.5 increases, while that of both BBHs and BHNSs increases substantially, especially at high metallicities. This happens because the NS-to-BH transition moves toward higher masses for higher values of ζ 2.5 and this effect is more prominent for higher metallicities (cf. Fig. 3). Because the IMF is bottom-heavy, the location of the NS-to-BH transition has a strong impact on the relative formation rates when outcomes are classified in terms of BBHs, BHNSs, and BNSs.
The fraction of merging system is largely independent of the fallback parameter f H . Indeed, results from simulations with different values of f H overlap almost perfectly in Fig. 5 (at least within the resolution of the figure). In order to merge within a Hubble time, most systems need to evolve through episode(s) of ). The bottom row shows control results from the standard rapid and delayed SN models. In all cases, colors indicate the metallicity, with Z = 0.0002, 0.006, 0.012, 0.02 from darker to lighter. We further separate between the full sample of compact binaries produced in our evolutions (empty histograms) and the subset of those that merge within a Hubble time (filled histograms). mass transfer and/or common envelope, with a consequent loss of their hydrogen envelope. The progenitors of merging binaries are most likely Wolf-Rayet stars with light envelopes and, consequently, the resulting BHs are relatively unaffected by f H .
Finally, it is interesting to note that compactness model predicts a larger fraction of merging BBHs compared to either the rapid or the delayed models for all values of ζ 2.5 and f H . This is related to the IMF, which in our case is a power law with negative spectral index. Moving the NS-to-BH transition at smaller ZAMS masses favors the formation of a larger number of BHs, which in turn translates into a larger number of (merging) BBHs. The merger time due to GW emission also depends on the component masses, with massive systems merging in a shorter time (Peters, 1964). This effect is also present in the case of merging BNSs: the rapid model tends to produce lighter NSs and thus has the lowest number of merging BNSs for all metallicities. This also explains why all runs performed with the compactness model have similar distributions of merging BNSs. We find that the effect of ζ 2.5 is strongest for BHNSs. This because shifting the formation threshold between NSs and BHs toward higher masses affects the kicks imparted to the newly formed compact objects. In models with larger ζ 2.5 , relatively massive stars that would have formed BHs now collapse to NSs, with a consequent large mass ejection during the explosion. Conservation of linear momentum then translates the large mass of the ejecta into strong natal kicks with a consequent suppression of the merger rate.

Conclusions
SNe are a key process to understand the formation of compact binaries, with different explosion mechanics producing different observable features (e.g. the presence of mass gaps). In this work, we have investigated the impact of the pre-SN stellar compactness (O'Connor & Ott, 2011) on the mass spectrum of compact objects and the resulting population of GW sources. In particular, we presented a new version of the MOBSE populationsynthesis code that includes a new model of BH and NS formation. Our "compactness model" has two free parameters : • The quantity f H describes the fraction of the hydrogen envelope that falls back onto BHs after their formation. BHs are lighter for smaller f H while, in our simple model at least, the masses of NSs are unaffected. This might have important consequences when combining GW sources formed from isolated binaries to those formed in dynamical formation channels (e.g. Bouffanais et al. 2019;Wong et al. 2021;Zevin et al. 2021). For compactness models with small f H , the heaviest BHs in the sample become exclusive to dynamical channels where, for example, more massive objects can be assembled via repeated mergers (Gerosa & Fishbach, 2021).
• The parameter ζ 2.5 instead marks the stellar compactness where stars transition from forming NSs to forming BHs. We find that this parameter mostly affects the relative abundance of BHs and NSs that are expected to form in realistic populations of massive stars.
We compared our models against two of the most commonly used prescriptions for SN explosions in binary population synthesis: the rapid and delayed models first introduced by Fryer et al. (2012). In particular, we find that distinguishing between these models and our compactness model with GW data will be facilitated by the detection of BHs which are relatively light (M 30M ). For heavier BHs, predictions from most of the models we tested tend to overlap. The cleanest signature of the compactness model is the location of the upper edge of the lower mass gap between the lightest black holes and the heaviest neutron stars, which depends mainly on f H , together with the metallicity Z. The merger rates might also help us discriminate between the different models, although this claim needs to be verified with an analysis which includes the metallicity-dependent star-formation history and the LIGO/Virgo selection effects.
The compactness model we implemented in MOBSE suffers from several simplifications. First, we use the stellar compactness only to discriminate between NSs and BHs, but then rely on ad-hoc prescriptions for their masses. Furthermore, we use a monotonic expression for the compactness as a function of the CO mass, see Eq.
(2), although recent studies have shown that their interplay is more complex (Chieffi & Limongi, 2020;Chieffi et al., 2021;Sukhbold et al., 2018). Both these aspects will be refined in future work together with estimates of the resulting merger and detection rates. The consequences of the SN models presented in this paper in terms of GW emission and heavyelement nucleosynthesis also need to be further explored.