1 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, 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\(_\odot\) are robustly predicted to undergo core-collapse SN once their cores reach the Chandrasekhar limit (for a recent review on SN theory, see Burrows and 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_\odot\) 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_\odot\), 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 \(\,50 M_\odot\) (Farmer et al. 2019; Woosley and 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 of compact objects with masses between \(\,3 M_\odot\) and \(\,5 M_\odot\) (“lower mass gap”) separating BHs and NSs has been long investigated using X-ray binary data (Bailyn et al. 1998; Özel et al. 2010; Farr et al. 2011) 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 and 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 and 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 and 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 binaries, 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 and Müller 2020).

Hydrodynamical simulations (e.g., O’Connor and 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,

$$\begin{aligned} \zeta _{\mathrm{M}} = \frac{M/M_\odot }{R(M)/1000\,\mathrm{km}}, \end{aligned}$$
(1)

where R(M) is the radius that encloses a mass M at core bounce.

Fig. 1
figure 1

Mass spectrum of compact remnants as a function of the progenitor initial mass \(M_{\mathrm{ZAMS}}\) in our fiducial compactness SN model (\(\zeta _{2.5}=0.365\) and \(f_{\mathrm{H}}=0.9\)). Different colors represent different metallicities \(Z=0.0002, \dots , 0.02\) (dark to light)

Full size image

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 Sect. 2, we describe an updated version of the mobse code where we have implemented a treatment of SNe based on the stellar compactness. In Sect. 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 Sect. 4, we summarize our findings and illustrate prospects for future work.

2 Methods

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

2.1 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”, Giacobbo et al. 2018). Built on top of the older bse (Hurley et al. 2002), mobse’s key updates include up-to-date prescriptions for stellar winds in massive stars (Giacobbo and Mapelli 2018), natal kicks (Giacobbo and Mapelli 2019, 2020), and modeling of pulsational pair-instability and pair-instability SNe (Spera and 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.

2.2 Compactness model

While the prescriptions for the rapid and delayed models are straightforward to implement in a population-synthesis code, the compactness parameter \(\zeta _{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 and Ott (2011), we will refer to the reference value of \(\zeta\) at \(M=2.5 M_\odot\), resulting in the parameter \(\zeta _{\mathrm{2.5}}\). Limongi and Chieffi (2018) demonstrated that there is a strong correlation between \(\zeta _{2.5}\) and the carbon-oxygen mass \(M_{\mathrm{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

$$\begin{aligned} \zeta _{2.5} \simeq 0.55 - 1.1 \left( \frac{M_\mathrm{CO}}{1\, M_{\odot }}\right) ^{-1}\,. \end{aligned}$$
(2)

We make use of \(\zeta _{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 and Ott 2011). In our fiducial model, we assume a threshold value of \(\zeta _{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 \(\zeta _{2.5} > 0.365\) will form BHs, while progenitors with \(\zeta _{2.5} \le 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 \(\mu = 1.33\) M\(_\odot\) and standard deviation \(\sigma = 0.09\) M\(_\odot\), which agrees with the observed masses of NSs in binary systems (Özel and Freire 2016). BH masses are instead given by

$$\begin{aligned} M_{\mathrm{BH}} = M_{\mathrm{He}} + f_{\mathrm{H}}(M_{\mathrm{fin}} - M_{\mathrm{He}}), \end{aligned}$$
(3)

where \(M_{\mathrm{He}}\) is the mass of the helium core and \(M_{\mathrm{fin}}\) is the total stellar mass at the onset of collapse (both of which are tracked in standard population synthesis). The quantity \(f_\mathrm{H} \in [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_{\mathrm{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 and 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.

Fig. 2
figure 2

Comparison between different SN explosion models. Red dash-dotted, purple dotted and green solid lines mark the fiducial compactness model, the rapid model, and the delayed model, respectively. Each panel illustrates how the mass of the remnant \(M_{\mathrm{rem}}\) depends on the mass of the initial star \(M_{\mathrm{ZAMS}}\) for different metallicities Z. The shaded areas indicate the range of masses obtained with \(f_{\mathrm{H}}=0.9\) (fiducial, thick red lines) and \(f_{\mathrm{H}}=0.1\) (thin orange lines). In this figure, the NS/BH compactness threshold is fixed to \(\zeta _{2.5}=0.365\)

Full size image

Figure 1 illustrates the relationship between the initial mass of the star \(M_{\mathrm{ZAMS}}\) and the mass of the compact remnant \(M_{\mathrm{rem}}\) obtained with our fiducial compactness model (\(\zeta _{2.5} = 0.365\) and \(f_{\mathrm{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; Spera et al. 2015; Giacobbo et al. 2018; Neijssel et al. 2019), 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_{\mathrm{ZAMS}}\sim 70 M_\odot\) for \(Z\lesssim 0.004\) is due to pulsational pair-instability SNe while the sharp decrease at \(M_{\mathrm{ZAMS}}\sim 140 M_\odot\) marks the onset of the proper pair-instability regime.

2.3 Simulation setup

Starting from the fiducial model we just described, we further explore the parameter space spanned by \(\zeta _{2.5}\) and \(f_\mathrm{H}\). We consider four different values for the threshold parameter \(\zeta _{2.5} =\) 0.2, 0.3, 0.365 (fiducial), and 0.4. For each of these, we vary \(f_{\mathrm{H}} =\) 0.1, 0.3, 0.5, 0.7, and 0.9 (fiducial). For all the resulting combinations of \(\zeta _{2.5}\) and \(f_{\mathrm{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 (Giacobbo and Mapelli 2018): primary ZAMS masses are extracted from an initial mass function (IMF) \(p(M_1)\propto M_1^{-2.3}\) with \(M_1 \in [5,150] M_\odot\), mass ratios \(q=M_2/M_1\) are drawn according to \(p(q|M_1) \propto q^{-0.1}\) with \(q \in [0.1,1]\), periods P are drawn according to \(p(P) \propto \log (P/\mathrm{days})^{-0.55}\) with \(P\in [10^{0.15},10^{5.5}]\), and eccentricities are drawn according to \(p(e) \propto e^{-0.45}\) with \(e \in [0,1]\) (Sana et al. 2012).

Fig. 3
figure 3

BH mass spectrum as a function of the initial mass \(M_{\mathrm{ZAMS}}\) for different metallicities \(Z=\) 0.02, 0.012, 0.006, and 0.0002 (colors). Each panel is produced assuming a different value of the NS/BH compactness threshold \(\zeta _{2.5}=\) 0.2, 0.3, 0.365, 0.4 (top to bottom). The shaded areas indicate the range of masses obtained with \(f_{\mathrm{H}} = 0.9\) (fiducial, thick lines) and \(f_{\mathrm{H}} = 0.1\) (lower edge of the colored regions)

Full size image

3 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 (Sect. 3.1). We then illustrate the formation of compact object binaries (Sect. 3.2) and focus on the sub-sample of merging systems (Sect. 3.3).

Fig. 4
figure 4

Distribution of chirp masses \(M_{\mathrm{chirp}}\) of BBHs (two leftmost columns), BHNSs (two middle columns), and BNSs (two rightmost columns). The top four rows show results obtained with the compactness model introduced in this paper, with \(\zeta _{2.5}=\) 0.2, 0.3, 0.365, 0.4 increasing from top to bottom. The left (right) column of each pair is obtained with \(f_\mathrm{H}=0.1\) (\(f_{\mathrm{H}}=0.9\)). 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)

Full size image

3.1 Mass spectrum

In Fig. 2, we compare the mass spectrum for our fiducial compactness model (\(\zeta _{2.5} = 0.365\), \(f_\mathrm{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 \(\sim 2\) to \(\sim 5 M_\odot\) independent of the metallicity, in the compactness model the size of the mass gap strongly depends on both Z and \(f_{\mathrm{H}}\). In particular, the shaded areas in Fig. 2 indicate the range of predicted masses by compactness models with \(0.1\le f_{\mathrm{H}} \le 0.9\). The impact of \(f_{\mathrm{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_{\mathrm{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_{\mathrm{rem}}\) at \(20 M_\odot \gtrsim M_{\mathrm{ZAMS}} \gtrsim 30 M_\odot\) 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_{\mathrm{zams}}\gtrsim 50 M_\odot\) (but the precise value depends on Z), the rapid, delayed and fiducial compactness model all predict the formation of BHs via direct collapse and thus return similar remnant masses.

The four panels of Fig. 3 illustrate the impact of \(\zeta _{2.5}\) on the mass spectrum. As \(\zeta _{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. Sect. 3.2). This is also sensitive to the metallicity. At \(Z=0.02\) (\(Z=0.0002\)) the ZAMS mass transition value changes from \(\,14 M_\odot\) (\(\,13 M_\odot\)) for \(\zeta _{2.5}=0.2\) to \(\,37 M_\odot\) (\(\,23 M_\odot\)) for \(\zeta _{2.5}=0.4\).

Finally, let us notice that \(f_{\mathrm{H}}\) has a strong impact on the maximum BH mass. Smaller values of \(f_{\mathrm{H}}\) lead to lighter BHs by construction, because the parameter \(f_{\mathrm{H}}\) describes the amount of fallback material. More interestingly, for massive progenitors (\(\gtrsim 60\) M\(_\odot\)) and lower values of \(f_\mathrm{H}\sim 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 \sim 0.006\).

3.2 Compact binaries

From our samples of simulated massive stars (see Sect. 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_\mathrm{chirp} = {(M_1M_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 \(\zeta _{2.5}\) result in a narrower mass range for both BBHs and BHNSs. For BBHs, the minimum chirp mass ranges from \(\,2.5\) M\(_\odot\) for \(\zeta _{2.5}=0.2\) to \(\,7\) M\(_\odot\) for \(\zeta _{2.5}=0.4\). In particular, for \(\zeta _{2.5} = 0.365\) (\(\zeta _{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. Sect. 2.2) does not depend on \(\zeta _{2.5}\). Consequently, the effect of \(\zeta _{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_{\mathrm{H}}=\) 0.1 and 0.9, for each type of compact binaries and values of \(\zeta _{2.5}\). As already highlighted in Sect. 2.2, the impact of \(f_{\mathrm{H}}\) on the remnant masses depends on both Z and \(\zeta _{2.5}\). In particular, the mass spectrum of lighter stars (\(M_{\mathrm{ZAMS}} \lesssim 40\)) presents a shallower slope for higher \(f_{\mathrm{H}}\). This translates into well-defined peaks in the chirp masses of BBHs for \(f_{\mathrm{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_\mathrm{H}=0.1\) compared to \(f_{\mathrm{H}}=0.9\). This is a direct consequence of the fact that, for larger ZAMS masses, the parameter \(f_{\mathrm{H}}\) only affects BH masses at low metallicities.

Fig. 5
figure 5

Number of compact binary mergers as a function of metallicity for BBHs (top panel), BHNSs (middle panel), and BNSs (bottom panel). Results from the compactness model are shown in different shades of blue for \(\zeta _{2.5}=0.2,0.3, 0.365,0.4\) (light to dark). The fiducial model with \(\zeta _{2.5}=0.365\) is marked with circles, while the additional model variations are marked with triangles. Red pluses and yellow crosses indicate the standard rapid and delayed models, respectively. For all the compactness models, results obtained for \(f_{\mathrm{H}}=0.1, 0.3, 0.5, 0.7,\) and 0.9 are largely indistinguishable on the scale of this plot

Full size image

3.3 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 Giacobbo and Mapelli (2018); Spera et al. (2019); Wiktorowicz et al. (2019) , we find that the most massive systems do not merge within a Hubble 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 \(\zeta _{2.5} \lesssim 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 metallicities \(Z\sim 10^{-3}\).

The number of merging BNSs tends to slightly decrease as \(\zeta _{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 \(\zeta _{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_{\mathrm{H}}\). Indeed, results from simulations with different values of \(f_{\mathrm{H}}\) overlap almost perfectly in Fig. 5 (at least within the resolution of the figure). To merge within a Hubble time, most systems need to evolve through episode(s) of 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_{\mathrm{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 \(\zeta _{2.5}\) and \(f_{\mathrm{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 \(\zeta _{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 \(\zeta _{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.

4 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 and 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 population-synthesis code that includes a new model of BH and NS formation. Our “compactness model” has two free parameters (Mapelli et al. 2020):

  • The quantity \(f_{\mathrm{H}}\) describes the fraction of the hydrogen envelope that falls back onto BHs after their formation. BHs are lighter for smaller \(f_{\mathrm{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; Zevin et al. 2021; Wong et al. 2021). For compactness models with small \(f_{\mathrm{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 and Fishbach 2021).

  • The parameter \(\zeta _{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\lesssim 30 M_\odot\)). 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_{\mathrm{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 (Sukhbold et al. 2018; Chieffi and Limongi 2020; Chieffi et al. 2021). 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 heavy-element nucleosynthesis also need to be further explored.