© The Authors 2023

1. Introduction

The greater the number of detections of gravitational wave (GW) sources, the better we know and understand the distribution of the masses of the detected remnants (e.g., Abbott et al. 2021a). These remnants are often black holes (BHs) in a binary system, which merge with each other after their binary evolution, through the emission of GWs (e.g., Abbott et al. 2016). This evolution includes the collapse into BHs, potentially accompanied by the supernovae (SNe) of both stars, which means the SN mechanism can be important in determining the final remnant mass. We are interested in the distribution of these final masses.

The top panel of Fig. 1 shows the distribution of the individual masses of the detected GW sources, from GWTC 1, 2, and 3 (Abbott et al. 2021b). The binary black holes (BBHs) seem to show a bimodal distribution, with a gap between 14 M_⊙ and 22 M_⊙, along with an additional peak at 35 M_⊙, although this peak seems to disappear when taking into account the uncertainty of the data. The bottom panel shows the distribution of the corresponding chirp masses of the black hole binaries (BHBs), represented by asymmetric Gaussians (as described in Appendix A). This distribution appears to be trimodal (Tiwari & Fairhurst 2021; Abbott et al. 2023), with one mode below 10 M_⊙, one between 10 M_⊙ and 20 M_⊙, and one above 20 M_⊙. Here, there is a gap as well, at 11 M_⊙.

Fig. 1. GW data from Abbott et al. (2021b). Top panel: mass distribution of the individual BBHs (represented by a histogram of the mean values, with a bin-width of 2 M⊙). We can also represent the datapoints as asymmetric Gaussians, as described in Appendix A. The sum of these Gaussians, ∑idNi/dm, is given as well (solid line), together with a cumulative distribution function (CDF, dotted line on a different axis). Bottom panel: chirp mass distribution (grey lines), which is made to resemble the distribution from Abbott et al. (2023). In order to do this, we use asymmetric Gaussians as well, taking into account the 90% credible interval. We also show an adjustable-width kernel density estimation (AWKDE, dashed line) from Shimazaki & Shinomoto (2010), Shimazaki et al. (2018). The Gaussians are plotted on a logarithmic axis, while the sum and the AWKDE use a linear axis. Similarly to Abbott et al. (2023), we omitted GW190814 from our analysis.

One hypothetical explanation for this is given by Broadhurst et al. (2018, 2022), who theorize that the gap (specifically the gap in the individual mass distribution) is real and caused by gravitational lensing, as a result of which the distance to some of the GW sources is underestimated and therefore the chirp mass overestimated. The gap is then located between the non-lensed and lensed sources. In order to produce a distribution similar to the top panel of Fig. 1, Broadhurst et al. (2018, 2022) have posed a merger rate that has high values at high redshift and low values at low redshift, since this produces the desired number of lensed and unlensed sources. We estimate, however, that the merger rate value they pose at high redshift implies a merger fraction that is approximately 50 times higher than the value Mandel & Farmer (2022) estimated through their Drake-like approach. Even if we make optimistic assumptions about some of the factors in this Drake-like equation, we still find a factor of 14 unaccounted for (as shown in Appendix B). Although this factor is not large enough to completely dismiss the lensing hypothesis of Broadhurst et al. (2018, 2022), we find it sufficient to deem their hypothesis improbable.

We aim to explore whether the supernova (SN) mechanism could provide an alternative and more natural explanation for the observed gap. Literature suggests that the shock-wave in a core-collapse SN (CCSN) can be stalled and not cause a successful explosion (Mazurek 1982). The star can then collapse in its totality and form a stellar remnant in a direct collapse (or “failed”) SN. Failed SNe are difficult to detect, since they do not cause a bright explosion but instead make a massive star seemingly disappear (Kochanek et al. 2008). Because of these difficulties, the properties of failed SNe are not well-known, even though there are several potential candidates (e.g. Gerke et al. 2015; Adams et al. 2017a,b; Basinger et al. 2021; Neustadt et al. 2021). In addition to the search for failed SNe, others have also looked at the implications for remnant mass distributions, for instance, Kochanek (2014a,b), who attempts to estimate the low-mass BH distribution based on the failed SN mechanism. Overall, we state that the difference between failed SNe and successful ones could potentially create a gap in the mass distribution: some sources have a successful SN and lose mass before forming a remnant and other sources have a failed SN, conserve mass, and form a more massive remnant. This means that the parameter which is important for the remnant mass distribution is the mass limit above which SNe fail.

In this work, we investigate the SN mechanism in close binaries and how it affects the BBH mass distribution (when we mention SNe we refer to CCSNe, unless specified otherwise). We start by investigating this SN mechanism in Sect. 2, wherein we describe the transitions between successful and failed SNe through an initial-final mass relation based on the works of Schneider et al. (2021) and Marchant et al. (2019; see Sect. 2.1). We also attempt to compare this relation to observed SN rates (Sect. 2.2). Then, in Sect. 3, we show that this model naturally produces a bimodal mass distribution and a trimodal chirp mass distribution, similar to Fig. 1. We note that while our paper was under review, Schneider et al. (2023) published a paper with a similar approach. In the discussion (Sect. 4) we make the relevant comparisons among the results. Finally, in Sect. 5, we summarize our conclusions.

2. Failed supernovae

2.1. Remnant function

We are interested in the differences between successful and failed SNe, as well as the relation between zero-age main-sequence (ZAMS) mass and the masses of the remnants they form. The SN explodability simulations indicate that there may not be only one mass range for successful SNe and one range for failed SNe (Ertl et al. 2016; Müller et al. 2016; Kresse et al. 2021). Instead, they find that there may be a failed SN range somewhere between 20 M_⊙ and 24 M_⊙, above which stars can have a successful SN again up until about 27 M_⊙, where the failed SNe take over again, supported by other studies of pre-SN compactness (e.g., Sukhbold & Woosley 2014; Sukhbold et al. 2016). We combine this with the results of Schneider et al. (2021), who performed a simulation and investigate stellar evolution and remnant masses. They also find that there are two ranges for failed SNe: one small range for lower masses, which they call the “BH island”, and a larger failed SN range for higher masses. They note that the two points of transition from successful to failed SNe coincide with the masses at which the stellar core change from convective to radiative carbon and neon burning, respectively. This could explain why there is a BH island in the first place. Schneider et al. (2021) also note that their relation between ZAMS mass and remnant mass gives rise to a bimodal mass distribution.

We constructed a possible remnant function, starting by describing the BH island, based on Schneider et al. (2021), both for their models related to case A/B mass transfer as well as their case C models. They describe the BH island as a function that produces remnants between approximately 8 M_⊙ and 10 M_⊙. In order to reproduce the observed BBH distribution, however, we changed this to a linear function that produces remnants between 8 M_⊙ and 14 M_⊙. We give the exact definition of our remnant function, including these mass ranges, in Appendix C. For the other failed SNe, we used the fits given by Schneider et al. (2021), but because these fits start at a remnant mass of about 16 M_⊙, we shifted them upwards to 22 M_⊙ instead. We motivate this based on the general uncertainty in the models, coupled with the fact that rather than complete collapse or complete ejection, there could be partial fallback. Thus, some stars that are in the successful SN domain will in fact be able to form BHs. Finally, at the highest masses, we include the effect of pair-instability SNe (PISN) that leave behind no stellar remnant (Fraley 1968) and pulsational pair-instability supernova (PPISN; Woosley et al. 2007; Woosley 2017), at slightly lower masses that can form the transition between the remnants of failed SNe and the remnantless PISNe. For this, we make a polynomial fit to the PPISN results from Marchant et al. (2019). In order to create one coherent remnant function, we increase the slope of the case A/B function, connecting it to this PPISN fit, and also shift the PPISN fit for case C to lower ZAMS masses, connecting it to the case C function. The first adjustment does not influence our results significantly and the second adjustment is justified because we expect case C stars to have a more massive helium core at the end of their evolution than case A/B stars have. This is due to the fact that case C stars are able to grow more massive cores, meaning that we would expect case C stars to have a lower threshold for PPISNe. We show our complete remnant functions, with one function for case A/B and one for case C, in Fig. 2. Here, we neglect the influences of aspects such as metallicity, which we discuss in Sect. 4. Nevertheless, this remnant function suffices for our purpose: to show that the bimodal BBH distribution is to be expected based on stellar evolution (described in Sect. 3).

Fig. 2. Remnant masses (mRN), as function of ZAMS mass (mZAMS). The blue curve shows the remnant function for case A and case B mass transfer, based on the simulations by Marchant et al. (2019) and Schneider et al. (2021). The Marchant et al. (2019) results are shown as well (black dots). The triangle shows the transition point between the part of our function which is based on Schneider et al. (2021) and the part which is based on Marchant et al. (2019). The red curve shows the corresponding remnant function for case C mass transfer, where the part after the transition point is the PPISN fit translated to lower masses by 42 M⊙. The dashed lines show mRN/mZAMS = 1, 0.65 and 0.35, since the latter two are used as approximation in estimating the amount of mass transfer in our estimation. The shaded area shows the values of the remnant mass gap. Furthermore, we note that this remnant function is limited to BHs, as for our purposes, we are not interested in the neutron star (NS) masses. The remnant functions are explicitly given in Appendix C.

2.2. Supernova rates

Before we turn to the mass distribution resulting from the remnant function, we describe how we tried to detect failed SNe by investigating how they reduce SN rates. After all, since failed SNe do not produce the bright signal successful SNe do, an increase in amount of failed SNe would reduce the total SN rate. In order to make a rough estimate of the effect of failed SNe on SN rates, we assume one mass range of stars that exhibit a successful SN, from m_l to m_u. Above m_u, stars are too heavy to explode and have a failed SN, meaning they do not contribute to the SN rate. This model differs from our remnant function, since here we assume only one mass range of failed SNe, but for our purposes this suffices. In this estimation, we neglect the influences of binary interaction, effectively limiting our scope to type II SNe.

We started by describing the SN rates. Using the initial mass function (IMF), N(m)dm = κMϕ(m)dm, where κ⁻¹ = ∫mϕ(m)dm and ϕ(m) the Chabrier (2003) IMF shape, we express the total mass, M, in terms of the star formation rate (SFR) ψ(t). This comes down to the total mass of the stars with a mass equal to m, which are created a lifetime τ ago and equals M = ψ(t − τ(m))dt. Therefore, the SN rate is simply the amount of stars within a certain mass range, which reach the end of their stellar lifetime at time t and is defined as:

$$\begin{array}{c}\hfill R_{\mathrm{SN}}(t,m_l,m_u)=\kappa {\displaystyle \int }\frac{M}{\mathrm{d}t}\varphi (m)\mathrm{d}m=\kappa {\displaystyle {\int }_{m_l}^{m_u}}\psi (t-\tau )\varphi (m)\mathrm{d}m,\end{array}$$(1)

where we use τ = τ(m) = 10¹⁰ yr(m/ M_⊙)^−2.5. Since the precise shape of the SFR is often not determined accurately, it is difficult to determine R_SN using this equation. Because of this, we look at the ratio between the SN rate and the SFR. The integrand of this ratio has a factor ψ(t − τ)/ψ(t), which is why we make the approximation ψ(t − τ)/ψ(t)≈1. This approximation holds for constant SFRs and decreases in accuracy for SFRs which evolve on a short timescale. Also, we set m_l = 8 M_⊙, which is non-controversial (e.g., Kochanek et al. 2008; Kochanek 2014b). The ratio of SN rate to SFR then becomes:

$$\begin{array}{c}\hfill \frac{R_{\mathrm{SN}}(t,8M_{\odot },m_u)}{\psi (t)}\approx \frac{{\int }_{8M_{\odot }}^{m_u}\varphi (m)\mathrm{d}m}{\int m\varphi (m)\mathrm{d}m}\equiv \mathrm{\Lambda }(m_u),\end{array}$$(2)

where we used the definition of κ and introduce Λ(m_u) as shorthand for this function. The κ integral is taken from m_min = 10⁻¹ M_⊙ up to m_max = 10² M_⊙, which means the value which does not take failed SNe into account is simply Λ(m_max).

In order to determine m_u, we consider the SN simulations of Ertl et al. (2016), Müller et al. (2016), and Kresse et al. (2021), together with our remnant function based on Schneider et al. (2021). These simulations use different models of neutrino engines, which represent the collapsed core, and determine the mass dependent explodability per model. Müller et al. (2016) find a mass dependent probability distribution, where approximately m < 20.5 M_⊙ and 23.5 M_⊙ < m < 27 M_⊙ have a high probability for a successful SN. The simulations of Ertl et al. (2016) and Kresse et al. (2021) find a similar distribution, with approximately m < 21.5 M_⊙ and 25 M_⊙ < m < 27.5 M_⊙ which go SN. These are in good agreement with Schneider et al. (2021), who find that their case C results (as shown in Fig. 2) are similar to the results for single stars. Keeping in mind the IMF, we make a crude approximation and estimate, based on these simulations, that m_u ≈ 22.5 M_⊙. This is comparable to the estimation of 20 M_⊙ used by Mashian & Loeb (2017), for instance.

The value of m_u gives $\mathrm{\Lambda }(22.5M_{\odot })=9.06\times {10}^{-3}{M}_{\odot }^{-1}$. In contrast, if we neglect the influence of failed SNe, this value becomes $\mathrm{\Lambda }(m_{max})=1.18\times {10}^{-2}{M}_{\odot }^{-1}$. This means we have found a failed SNe correction of 1 − Λ(22.5)/Λ(m_max) = 0.232 ≈ 25%.

We now turn to SN surveys and use the works of Graur et al. (2015, 2017) and Botticella et al. (2017) in order to detect the predicted effect of failed SNe. As shown in the top panel of Fig. 3: these data consist of the type II SN rate per unit mass (R_SNuM) versus the specific SFR (sSFR), that is, the SFR per unit mass. The figure includes two fits, shaped as R_SNuM = a·sSFR^b, where the first one includes all the data-points and the second one neglects two outliers. The fitted values of b, 0.68 and 1.04 confirm that our assumption ψ(t − τ)/ψ(t)≈1 (i.e., R_SN ∝ ψ and b = 1) is appropriate. If we we assume b ≈ 1, the two points which are neglected for the second fit indeed appear to be outliers. The consequence of this assumption is that a = R_SNuM/sSFR = Λ. Our estimate $\mathrm{\Lambda }(22.5M_{\odot })=9.06\times {10}^{-3}{M}_{\odot }^{-1}$ differs significantly from the fitted values of a: $4.0\times {10}^{-6}{M}_{\odot }^{-1}$ and $1.9\times {10}^{-2}{M}_{\odot }^{-1}$. If we set b = 1, the fitted values of a become $(7.2\pm 1.7)\times {10}^{-3}{M}_{\odot }^{-1}$ and $(7.3\pm 0.9)\times {10}^{-3}{M}_{\odot }^{-1}$, respectively, which is consistent with Λ. The bottom panel of Fig. 3 shows the distribution of the individual values of R_SNuM/sSFR, including a kernel density estimation. The distribution shows an outlier at $3\times {10}^{-2}{M}_{\odot }^{-1}$, which corresponds to the outlier at 1.5 × 10⁻¹² yr⁻¹ in the top panel, which was neglected in the second fit. Also, a peak is situated at approximately $0.8\times {10}^{-2}{M}_{\odot }^{-1}$, which is approximately the value of the fitted a with b = 1.

Fig. 3. Data from SN surveys (Graur et al. 2015, 2017; Botticella et al. 2017). Top panel: Type II SN rates per unit mass (RSNuM) versus the SFR per unit mass (sSFR), together with a fit for all the data-points (dashed line) and a fit which excludes the two outliers at sSFR ≈ 1.5 × 10−12 yr−1 and 1.5 × 10−9 yr−1 (dotted line). The fits use the least-squares method on a first degree polynomial in logarithmic space and the error bars are one standard deviation. These data are corrected for the Chabrier (2003) IMF, since the original data follows the Salpeter (1955) IMF. For a constant SFR, L = κψ∫τϕL(m)dm, where L(m)∝m3.5 is the luminosity of a star with mass m. However, since τ(m)⋅L(m)∝m, it cancels out with κ and makes the SFR independent of the IMF. Since the SN rate is proportional to κ, it is multiplied by a factor κChabrier/κSalpeter ≈ 1.38, this correction therefore also multiplies RSNuM/sSFR by a factor of 1.38. The bottom panel shows a histogram of the RSNuM/sSFR ratio (with a bin-width of $2\times {10}^{-3}{M}_{\odot }^{-1}$), and a kernel density estimation (KDE).

The question arises whether uncertainty in other aspects of our model can account for a similar deviation from the standard Λ(m_max). This deviation can, for instance, be accomplished by choosing m_l to be 9.72 M_⊙, which seems high compared to the literature value of 8 M_⊙, albeit not unreasonably. Furthermore, uncertainties in the IMF could also account for such a deviation. Not only does Chabrier (2003) give different parameter values for binary systems (which would cause a deviating value of Λ), but there is also uncertainty in the high-mass slope of the IMF. Parravano et al. (2018) state that this slope has a value of $\alpha =2.{35}_{-0.15}^{+0.35}$. A value of 2.45 is enough to account for the 25% deviation, which is well within this range. This means that even if Λ(22.5 M_⊙) fits the data better than Λ(m_max), which the data may suggest, it would be difficult to attribute this to failed SNe. Because of this, we cannot conclude that we find the effects of failed SNe in the SN survey data and therefore cannot constrain the value of m_u based on this analysis. Despite not being able to confirm the failed SNe observationally, we are interested in how the remnant function (described in Sect. 2.1) shapes the BBH mass distribution.

3. Binary black hole mass distribution

We go on to show that our remnant function can indeed produce a BBH distribution similar to observation. To do this, we start by making a simple estimation of binary evolution and determine the resulting remnant mass distribution. In our estimation, we started with approximately 1.2 × 10⁶ binaries, in which the primary masses are distributed according to the Chabrier (2003) IMF and the secondary masses are distributed uniformly. This means that we define the primary star here to have the highest initial mass.

In our model, we also include mass transfer from the primary to the secondary in a Roche-lobe overflow (RLO) phase, and neglect mass transfer from the secondary to the primary. We state that the secondary accretes a fraction, η, of the mass expelled from the primary during the RLO phase. The precise value of η is uncertain (Dorozsmai & Toonen 2022), so we simply set η = 0.5, since it does not affect the general shape of the remnant distribution. In order to determine the amount of mass which is expelled from the primary, we cannot simply take the difference between ZAMS mass and remnant mass, since in successful SNe there is mass loss outside of the RLO phase. We therefore approximate the expelled mass of the primary in the RLO phase by 0.35 ⋅ m_ZAMS for case A/B and 0.65 ⋅ m_ZAMS for case C, since these curves approximate the failed SN curves, as shown in Fig. 2. The secondary, then, accretes a fraction, η, of this expelled mass and is rejuvenated, which means that we use this new mass as if it were the ZAMS mass. Besides η we also define the parameter ζ, which equals the fraction of binaries which follow the case C curve. Since ζ also does not influence the general shape of the distribution significantly, we set ζ = 0.5 as well.

The estimated mass distribution is shown in Fig. 4. The left panels show the exact results from our estimation and the right panels show the results with an added Gaussian uncertainty, as described in Appendix D. Unsurprisingly, the remnant mass distribution is indeed bimodal: it shows one peak at lower masses, caused by the BH island, and one peak at higher masses. The exact results of our estimation, without the Gaussian uncertainty, also show a PPISN peak around 44 M_⊙, caused by the flat PPISN distribution in Fig. 2. The general shape of the estimated distribution is similar to the observed distribution, although it is difficult to compare the two in detail because the relative heights of the peaks are influenced by multiple aspects. Not only is there a bias towards high-mass mergers in observation, but the parameters η and ζ also influence the heights. The effects of these parameters, including the fact that varying them does not influence our conclusions about the general distribution, as shown in Appendix E.

Fig. 4. BBH mass distribution resulting from our estimation, where we simulate approximately N = 1.2 × 106 binaries and set η = 0.5 and ζ = 0.5. The latter means that half of the binaries follow the case A/B remnant function, and the other half follow the case C remnant curve (as shown in Fig. 2). The top row shows the resulting mass distribution of the individual primary (red) and secondary (grey) masses, in a stacked histogram with a bin-width of 1 M⊙. The bottom row shows the chirp mass distribution in a similar histogram, where we distinguish between four possible configurations of the primary and secondary mass with respect to the central value of the mass gap (mG): either both BHs are below the gap (dark blue), the secondary is below the gap while the primary is above it (blue) or vice versa (beige), or both are above the gap (light blue). The left column shows the exact results of our estimation, while the right column adds a Gaussian uncertainty to the results in order to make the simulated distribution look more similar to the observed one (as shown in Fig. 1). We based the Gaussian uncertainty on the average standard deviations of the GW data, as described in Appendix D.

The resulting chirp mass distribution also shows interesting features. Since the individual mass distribution is bimodal, with one peak at either side of the mass gap at m_G = 17 M_⊙, the chirp mass distribution shows four possible configurations, with the primary and secondary mass at the same and opposite sides of the mass gap. These configurations can clearly be found in the chirp mass distribution: there is one peak centered just below 10 M_⊙ for binaries with both BHs below the mass gap, one peak centered around 28 M_⊙ for binaries with both BHs above the mass gap, and also two peaks centered around 15 M_⊙, representing a mixed population with BHs at either side of the gap. The latter mostly consists of binaries where the primary is above the gap and the secondary below, although, depending on η, a small population of binaries where only the secondary is above the gap can be found as well. These features can be linked to the observed chirp mass distribution, because the observed distribution (Fig. 1) shows a trimodal distribution very similar to our results: one peak just below 10 M_⊙, one around 28 M_⊙, a small group of binaries around 15 M_⊙, and a gap just above 10 M_⊙. We also predict a gap just below 20 M_⊙, which is somewhat visible in our AWKDE but even more visible in the estimation by Abbott et al. (2023). Our results are not identical to the observed chirp mass distribution; it is especially above 20 M_⊙ that they seem to differ, but the general shapes of the distributions do seem to agree.

In comparing our simple estimation and the GW data, we examine the primary versus secondary mass distribution. Figure 5 shows this distribution, together with the posterior distributions of the GW observations. Both the simulated and the observed distributions clearly show the BH island and the other failed SNe. Also, there seem to be some observed mergers which can be linked to the mixed population. However, there are also differences between the two distributions. For example, there is a local maximum at m_P′ ≈ m_S′ ≈ 34 M_⊙, while the corresponding maximum in our distribution is situated around m_P′ ≈ m_S′ ≈ 28 M_⊙. Also, the observed distribution has relatively more BHs in the BH island, even though there is an observational bias towards high-mass mergers. However, we are interested in the general shape of the distribution and not normalization and details. After all, consideration of partial fallback as well as the parameters η and ζ influence the relative amount of BHs in the BH island. Another interesting feature seen in Fig. 5 is the fact that the observed mass gap seems to be situated at lower masses than the simulated mass gap; while we set the simulated mass gap at m_P′ ≈ m_S′ ≈ 17 M_⊙, the observed gap appears to be around m_P′ ≈ m_S′ ≈ 14 M_⊙. This could mean that we did not have to shift the results of Schneider et al. (2021) upwards in our remnant function at all. The difference can be explained by the fact that applying the Gaussian uncertainty to the individual masses results in a two-dimensional (2D) distribution which is more symmetric than the observed posteriors. Finally, it is interesting to compare Fig. 5 with the rate densities determined by Abbott et al. (2023). While we find three populations, namely, a BH island, other failed SNe, and a mixed population, they have found a similar distribution but for NS-NS, BH-BH, and BH-NS mergers, respectively.

Fig. 5. Primary versus secondary mass distributions, for the results of our estimation (left) and the GW data (right), both in 2D histograms with a bin-width of 1 M⊙. Left panel shows the same as the upper right panel from Fig. 4. We use the notation P′ and S′ here to denote the stars with the largest and smallest final masses, respectively, in contrast to P and S from Fig. 4, which concern the initial masses. The right panel shows the cosmologically reweighted posterior distributions of the GW data (LIGO Scientific Collaboration & Virgo Collaboration 2021; LIGO Scientific Collaboration, Virgo Collaboration, & KAGRACollaboration 2021). Since the posteriors are not normalized, we rescaled them to ${\int }_{0M_{\odot }}^{100M_{\odot }}\frac{\mathrm{d}N}{\mathrm{d}m}\mathrm{d}m\approx 1.4\times {10}^5$. Then, after adding all the posteriors, we normalized the resulting distribution in between 5 M⊙ and 65 M⊙. We omit GW200322 here, since no valid posterior for this source has been given. The grey lines are lines of equal chirp mass, with intervals of 5 M⊙.

4. Discussion

We have proposed failed SNe as a natural explanation for the shape of the BBH mass distribution, particularly the difference between failed and successful SNe. In order to do this, we posed a remnant function (Sect. 2.1), investigated type II SN rate data (Sect. 2.2), and examined the BBH distribution caused by this function (Sect. 3). However, we note that our model is rather simplistic and we do not claim that our results describe the actual binary mass distribution in detail. We are therefore careful not to draw too strong conclusions in comparing our results from Figs. 4 and 5 to the data as shown in Fig. 1. Although, the simplicity of our model does indicate that our conclusions are a robust feature of BBH distributions caused by remnant functions similar to ours.

Furthermore, we argue that our remnant function is a more natural hypothetical explanation for the BBH mass distribution than the gravitational lensing hypothesis mentioned in the introduction, for two reasons. First, although they are difficult to compare, we estimate the difference of a factor 50 between the Broadhurst et al. (2018, 2022) merger rate and an optimistic BBH fraction estimate (discussed in Appendix B) to be greater than the difference between our remnant function and the results of Marchant et al. (2019) and Schneider et al. (2021). Second, the chirp mass distribution is difficult to explain using the lensing hypothesis. After all, the two stars in the binary are either both lensed or both non-lensed. This means that there should not have been any binary where the primary is above the gap and the secondary below. Therefore, the lensing hypothesis has trouble explaining the trimodal chirp mass distribution, while our remnant function reproduces it.

As mentioned in the introduction: Schneider et al. (2023) published a paper which uses a similar method and reaches comparable conclusions, while our paper was in review. They use simulations equivalent to those of Schneider et al. (2021), showing a BH island, which resulted in a bimodal distribution of the individual BBHs. Their BBH distribution looks similar to the one we have found, except for the PPISN peak since they do not incorporate these. Moreover, they argue that the bimodal mass distribution causes a trimodal chirp mass distribution with the same argument we used: the mixing between BHs from the BH island and the higher mass BHs. They also note that fallback after the SN can extend the BH island to lower masses, due to the same reason we extended our BH island in Fig. 2 to include lower masses. Furthermore, according to Schneider et al. (2023), metallicity mainly affects the case A/B curve in Fig. 2, but because of the mass-loss through wind associated with metallicity, they state that an increase in metallicity can shift the resulting mass distribution to lower masses. Overall, we find that their results are in good agreement with our findings, which strengthens our argument.

5. Conclusions

After our analysis of the BBH mass distribution and chirp mass distribution (Sect. 3) resulting from our remnant function (Sect. 2), we draw the following conclusions:

• The mass distribution of BBHs appears to show a gap for approximately 14 M_⊙ < m < 22 M_⊙ (Fig. 1). Gravitational lensing as an explanation for this gap (Broadhurst et al. 2018, 2022) seems, initially, to assume an unreasonably high merger rate for z > 2 (Fig. 2). We find that this merger rate implies a BBH fraction which is approximately 50 times larger than the one described by Mandel & Farmer (2022). This is improbably large, but not large enough to completely dismiss this possibility.

• We investigate failed SNe as a more natural explanation for the gap. Firstly, we approximate the results of the simulations by Marchant et al. (2019) and Schneider et al. (2021) with a comprehensive remnant function which describes the relation between ZAMS mass and remnant mass, for case A/B and case C binaries (Fig. 2). This function includes a BH island, whose existence is supported by SN explodability simulations and is the cause of the bimodality in the BBH mass distribution.

• We investigate whether we can confirm the mass range for failed SNe observationally. We assume a single mass range for successful SNe, which turn into failed SNe above a certain mass limit. We take an optimistic value of 22.5 M_⊙ for this mass limit and find that this implies a R_SN/SFR reduction of approximately 25%, since failed SNe are not detected and included in the type II SN rate. However, even though this reduction is compatible with SN survey data to some degree, the data is quite uncertain and small changes in other model parameters could also account for a similar reduction in SN rate. Therefore, we cannot state that the failed SNe mass range can be confirmed by the SN survey data.

• Using this remnant function, a bimodal BBH distribution can be estimated which look similar to observation. Also, our simplistic estimation produces a trimodal chirp mass distribution, which corresponds to observation. The primary versus secondary mass distribution resulting from our estimation corresponds to some degree with the posterior samples of the observed GW sources. This distribution clearly shows a mass gap, and a comparison of estimation and observation indicates that our remnant function does not need to deviate from literature values in the degree that it does, strengthening our conclusions.

Based on our results, we conclude that failed SNe, and in particular the relation between ZAMS mass and remnant mass which includes a BH island, can provide a natural explanation for the apparent bimodality in the observed BBH mass distribution. We therefore state that, when trying to explain the general shape of the observed BBH distribution, a consideration of stellar evolution is shown to be fruitful.

Future research could expand or improve multiple aspects of our work. Reproducing the observed BBH distribution in detail, for instance, would be an interesting topic of research. Not only could the remnant function perhaps stay closer to literature values and still produce the desired mass gap, as implied by Fig. 5, but other aspects such as a consideration of partial fallback, PPISNe, and the values of η and ζ also influence the relative heights of the different peaks. It would be an interesting endeavor to create a model which includes these aspects in such a way that it reproduces the observed distribution in more detail. A relevant question may be related to the chirp mass distribution and whether or not it corresponds to such a model in detail. Finally, more (precise) SN survey and GW data would improve the observational verification, which could help in identifying the precise effects of failed SNe.

Acknowledgments

This research has made use of data or software obtained from the Gravitational Wave Open Science Center (https://gwosc.org/), a service of LIGO Laboratory, the LIGO Scientific Collaboration, the Virgo Collaboration, and KAGRA. LIGO Laboratory and Advanced LIGO are funded by the United States National Science Foundation (NSF) as well as the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO600 detector. Additional support for Advanced LIGO was provided by the Australian Research Council. Virgo is funded, through the European Gravitational Observatory (EGO), by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale di Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by institutions from Belgium, Germany, Greece, Hungary, Ireland, Japan, Monaco, Poland, Portugal, Spain. The construction and operation of KAGRA are funded by Ministry of Education, Culture, Sports, Science and Technology (MEXT), and Japan Society for the Promotion of Science (JSPS), National Research Foundation (NRF) and Ministry of Science and ICT (MSIT) in Korea, Academia Sinica (AS) and the Ministry of Science and Technology (MoST) in Taiwan. This research is supported by the Netherlands Organisation for Scientific Research (NWO). We thank both Onno Pols for useful discussions and the anonymous referees who provided comments which helped to improve this paper.

References

Appendix A: Asymmetric Gaussians

In order to approximate the posterior chirp mass distribution from Abbott et al. (2023) in Fig. 1, we use Gaussian probability density functions with two different standard deviations. This is equivalent to using two halves of different Gaussians which meet at the mean. These functions are described by the standard Gaussian formula:

$$\begin{array}{c}\hfill f^{\pm }(x)=\frac{1}{{\sigma }_{\pm }\sqrt{2\pi }}exp(-\frac{1}{2}{\left(\frac{x-\mu }{{\sigma }_{\pm }}\right)}^2),\end{array}$$(A.1)

where f⁻(x) describes the function below the mean μ with standard deviation σ₋ and f⁺(x) describes the function above the mean with standard deviation σ₊. In other words, σ_± = σ₋ if x ≤ μ and σ_± = σ₊ if x > μ. Also, $\overline{\sigma }$ is simply the average standard deviation (σ₋+σ₊)/2. The combination of these halves can then be constructed as

$$\begin{array}{c}\hfill f(x)=k_1·\{\begin{array}{cc}{f}^{-}(x)& \mathrm{for}x\le \mu \\ k_2{f}^{+}(x)& \mathrm{for}x>\mu \end{array},\end{array}$$(A.2)

where the constant k₂ makes sure f⁻(μ) = k₂f⁺(μ) and k₁ ensures normalization. From this description of k₂ and Eq. A.1 it immediately follows that k₂ = f⁻(μ)/f⁺(μ) = σ₊/σ₋. Also, the normalization requirement can be used to determine k₁: ${\int }_{-\infty }^{\infty }f(x)dx=1$ gives $k_1={\sigma }_{-}/\overline{\sigma }$. Filling in the values of k₁ and k₂ gives, then, the overall asymmetric Gaussian

$$\begin{array}{c}\hfill f(x)=\frac{1}{\overline{\sigma }\sqrt{2\pi }}exp(-\frac{1}{2}{\left(\frac{x-\mu }{{\sigma }_{\pm }}\right)}^2),\end{array}$$(A.3)

which has the following antiderivative:

$$\begin{array}{c}\hfill F(x)=\frac{{\sigma }_{\pm }}{2\overline{\sigma }}\mathrm{erf}\left(\frac{x-\mu }{\sqrt{2}{\sigma }_{\pm }}\right)·\end{array}$$(A.4)

However, the asymmetry brings with it some ambiguity. When only μ and the 90% credible deviations (L₋ and L₊) are given, it is not immediately obvious how to define a 90% interval in Eq. A.3. As a first option, it is possible to assume that the interval, μ − L₋ ≤ x ≤ μ + L₊, has borders proportional to the standard deviations with a constant a (i.e., L_± = aσ_±). This means: ${\int }_{\mu -a{\sigma }_{-}}^{\mu +a{\sigma }_{+}}f(x)dx=90\%$, or

$$\begin{array}{c}\hfill 90\%=\mathrm{erf}\left(\frac{a}{\sqrt{2}}\right)⟹a=\sqrt{2}{\mathrm{erf}}^{-1}(90\%)\approx 1.65.\end{array}$$(A.5)

Alternatively, it is possible to define the 90% interval as “cutting of” 5% on either side. The relation between the interval borders and the standard deviations then no longer has a single proportionality constant but two: b₋ and b₊ (i.e., L_± = b_±σ_±). The fact that ${\int }_{-\infty }^{\mu -L_{-}}f(x)dx={\int }_{\mu +L_{+}}^{\infty }f(x)dx=5\%$ translates into:

$$\begin{array}{c}\hfill 5\%=\frac{{\sigma }_{\pm }}{2\overline{\sigma }}-\frac{{\sigma }_{\pm }}{2\overline{\sigma }}\mathrm{erf}\left(\frac{L_{\pm }}{\sqrt{2}{\sigma }_{\pm }}\right)=\frac{1-\mathrm{erf}(b_{\pm }/\sqrt{2})}{1+{\sigma }_{\mp }/{\sigma }_{\pm }}.\end{array}$$(A.6)

This equation can be transformed, after substitution of σ_± = L_±/b_±, into the two (±) equations:

$$\begin{array}{c}\hfill b_{\pm }=\sqrt{2}{\mathrm{erf}}^{-1}(95\%-5\%\frac{b_{\pm }}{b_{\mp }}\frac{L_{\mp }}{L_{\pm }}),\end{array}$$(A.7)

which can be solved numerically to obtain b₋ and b₊. It is clear that Eq. A.7 satisfies σ₋ = σ₊ → b₋ = b₊ = a. Figure A.1 shows (similarly to Fig. 1), the chirp mass distribution represented by both methods of asymmetric Gaussians (to avoid a chaotic plot, only a subset of the data is shown). It is clear that both methods, using either a or b_±, give almost indistinguishable results for the data-set. In fact, the total distributions differ less than 5%. Because of this, we use the a-method (Eq. A.5) in this work, which is mathematically more straightforward.

Fig. A.1. Randomly selected subset of the chirp masses (Abbott et al. 2021b), represented by asymmetric Gaussians (as in Fig. 1), which are similar to the posterior distributions of Abbott et al. (2023). The 90%-interval is determined by either the a-method from Eq. A.5 (top panel) or the b±-method from Eq. A.7 (bottom panel).

Appendix B: Merger rate

According to Broadhurst et al. (2018, 2022), gravitational lensing can explain the apparent bimodality in the BBH mass distribution, as shown in the top panel of Fig. 1. They note that gravitational lensing focuses the GW signal, causing the distance of the merger to be underestimated by a factor equal to the square root of the magnification. This means that the inferred chirp mass is overestimated. After all, an increased distance means a higher redshift, or a lower frequency, which gives a higher inferred chirp mass. They state that the gap in BBH masses can be explained by the hypothesis that a significant part of the observed BBHs is gravitationally lensed. In order to explain the observed BBH mass distribution, their model needs a merger rate which has a high value at high redshift, since there has to be a significant amount of mergers to be lensed by massive objects at lower redshift. This merger rate, then, has to decline rapidly towards lower redshift. They pose this hypothetical merger rate, not as a result coming from a physical model, but because the desired bimodal BBH distribution follows from it. They explicitly assume all mergers with a chirp mass above 20M_⊙ to be gravitationally lensed (Broadhurst et al. 2022).

The merger rate they used takes the shape of R_merger(z) = a ⋅ S(z)⋅F(z), where the merger rate (R_merger) is given in terms of a constant a = 6 ⋅ 10³yr⁻¹Gpc⁻³, as well as a function S(z) that traces the Madau & Dickinson (2014) SFR and a modulation function F(z), which is equal to 1 for z > 2 and causes an exponential decline for 2 > z > 0.5. Figure B.1 shows both the merger rate and the Madau & Dickinson (2014) SFR, together with the SN rate (which we define in Section 2.2).

Fig. B.1. Madau & Dickinson (2014) star formation rate per M⊙, the corresponding SN rate with ml = 8M⊙ and mu = 30M⊙ (through Eq. 1) and the Broadhurst et al. (2018, 2022) BBH merger rate.

The merger rate seems to be extremely high for z > 2. In order to investigate if this is at all possible, we turn to the work of Mandel & Farmer (2022), who used a “Drake-like” equation to estimate the fraction of binaries that produce a merging BBH (f_BBH). In order to produce a BBH merger, the following fractions need to be multiplied: the fraction of primary stars which fall within the SN range (f_prim), the fraction of their secondary stars which fall inside this range (f_sec), the fraction of these binaries which have a separation suitable for merging (f_init sep), the fraction of these binaries which survive the primary SN (f_{survive SN1}) as well as the common envelope phase (f_CE), and the secondary SN (f_{survive SN2}), and the fraction of these binaries that successfully merge (f_merge). In total, this gives:

$$\begin{array}{cc}\hfill f_{\mathrm{BBH}}& =f_{\mathrm{prim}}·f_{\sec }·f_{\mathrm{init}\mathrm{sep}}·f_{\mathrm{survive}\mathrm{SN}1}·f_{\mathrm{CE}}·f_{\mathrm{survive}\mathrm{SN}2}·f_{\mathrm{merge}}\hfill \\ \hfill & \approx 0.001·0.5·0.5·1·0.1·1·0.2=5·{10}^{-6},\hfill \end{array}$$(B.1)

where we fill in the estimations of the fractions made by Mandel & Farmer (2022).

In order to compare the merger rate inferred by Broadhurst et al. (2018, 2022) to Eq. 1, this merger rate has to be converted to a BBH fraction (which we call $f_{\text{BBH}}^{\prime }$), this is equal to the ratio between the number of binary stars that merge and the total number of stars. The number of binary stars that merge is calculated by taking twice the merger rate (since a merger consists of two stellar remnants). The latter is determined through the SFR. However, we need to consider the fact that the SFR is determined through the observed luminosity. We define the binary to have a certain mass ratio q = m_S/m_P, which means that the total luminosity of the binary can be described by $L\propto m_P^{3.5}+m_S^{3.5}=\left(1+q^{3.5}\right)m_P^{3.5}$. In other words, the SFR implies 1 + q^3.5 stars, while there are actually 2. Therefore, the total number of stars has to be adjusted by a factor 2/(1+q^3.5). Here, we make the approximation 2/(1+q^3.5) ≈ 2, which is valid even if q takes values as high as 0.6 or 0.7.

We use the Chabrier (2003) IMF and compute the IMF integrals from m_min = 10⁻¹ M_⊙ to m_max = 10² M_⊙. The total amount of stars is, then, simply κM∫ϕ(m)dm. If we assume a Madau & Dickinson (2014) SFR, which is the function S(z) multiplied by a constant b = 1.5 ⋅ 10⁷ M_⊙yr⁻¹Gpc⁻³, we get for z > 2:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill f_{\mathrm{BBH}}^{\prime }& =\frac{N_{\mathrm{BBH}}}{N_{\mathrm{tot}}}=\frac{2R_{\mathrm{merger}}/dt}{2\kappa \psi (z)/dt\int \varphi (m)dm}=\frac{aS(z)F(z)}{bS(z)}\frac{\int m\varphi (m)dm}{\int \varphi (m)dm}\hfill \\ \hfill & \approx \frac{6·{10}^3}{1.5·{10}^7}·0.67=2.7·{10}^{-4},\hfill \end{array}\end{array}$$(B.2)

where ψ(z) is the SFR at redshift z, which means that M = ψ(z)/dt. We do not take into account the fact that there is a time difference between the stellar formation and the merger events, effectively assuming that the bulk of the mergers occur on a relatively short timescale, compared to the SFR.

In order to analyze $f_{\text{BBH}}^{\prime }$, we use fractions similar to those used in Eq. B.1. Leaving out f_{survive SN1} and f_{survive SN2}, since they are equal to 1, the remaining factors have to equal $f_{\text{BBH}}^{\prime }$ = 2.7 ⋅ 10⁻⁴ and account for the deviation from f_BBH of $f_{\text{BBH}}^{\prime }$/f_BBH = 56 ≈ 50. At this point, we must question whether this is even possible. Using IMF integrals, combined with the 2/(1+q^3.5) correction, we estimate that reasonable adjustments to f_prim and f_sec could both account for a factor of approximately 2. This means that the remaining $f_{\text{initsep}}^{\prime }$, $f_{\text{CE}}^{\prime }$ and $f_{\text{merge}}^{\prime }$ have to account for a factor of 56/4 ≈ 14, or an average factor of 14^1/3 ≈ 2.4 per fraction. This factor is large, but not large enough to completely dismiss the Broadhurst et al. (2018, 2022) merger rate at z > 2. After all, the factor of 14 can be distributed between the remaining fractions without one of them becoming greater than 1.

In contrast, if we look at the merger rate at z < 0.5, the fraction seems to become quite small. This is due to the factor F(z) in the merger rate. If we approximate the merger rate below z = 0.5 as a constant value (specifically as the value at z = 0.2), we find that this value equals F(0.2)⋅$f_{\text{BBH}}^{\prime }$(z > 2)≈1.4 ⋅ 10⁻³ ⋅ 2.7 ⋅ 10⁻⁴ = 3.8 ⋅ 10⁻⁷. This is approximately 8% of the Mandel & Farmer (2022) fraction (Eq. 1). Since f_BBH consists of 7 fractions, each has to be adjusted by a factor of 0.08^1/7 ≈ 0.7, on average, which is well within their uncertainty.

Appendix C: Remnant function

The remnant function we pose in Fig. 2 consists of the following equations. First, we follow Woosley (2019) for the relation between ZAMS mass and the final mass of the helium core:

$$\begin{array}{c}\hfill f_{\mathrm{He}}(m)=\{\begin{array}{cc}0.0385m^{1.603}\hfill & \mathrm{for}m\le 30M_{\odot }\hfill \\ 0.50m-5.87M_{\odot }\hfill & \mathrm{for}m>30M_{\odot }\hfill \end{array}.\end{array}$$(C.1)

We use this relation to rewrite the results of Marchant et al. (2019) as a function of ZAMS mass instead of helium core mass. This takes the form of two fourth degree polynomials, shaped as

$$\begin{array}{c}\hfill f_{\mathrm{PPISN}}(m)=\{\begin{array}{cc}{M}_{\odot }\sum _{i=0}^4{a}_i{\left(\frac{m}{M_{\odot }}\right)}^i& \mathrm{for}m\le 151.74M_{\odot }\\ M_{\odot }\sum _{i=0}^4{b}_i{\left(\frac{m}{M_{\odot }}\right)}^i& \mathrm{for}m>151.74M_{\odot }\end{array},\end{array}$$(C.2)

where the coefficients a_i are determined through a polynomial fit through the data points at f_He(m)/M_⊙ = 100, 115.74, 127.74 and 151.74, and b_i is given as a fit through the data points at f_He(m)/M_⊙ = 151.74, 167.74, 179.74 and 189.74, as given by Marchant et al. (2019). The two polynomials connect at $m=f_{\mathrm{He}}^{-1}(70M_{\odot })=151.74M_{\odot }$, with the restriction that $\frac{d}{\mathit{dm}}{f}_{\mathrm{PPISN}}(m)=0$ here. The values for the coefficients are then:

$$\begin{array}{c}\hfill \left(\begin{array}{c}{a}_0\\ a_1\\ a_2\\ a_3\\ a_4\end{array}\right)=\left(\begin{array}{c}-3.38751·{10}^2\\ 1.01736·{10}^1\\ -1.06666·{10}^{-1}\\ 5.18494·{10}^{-4}\\ -9.74407·{10}^{-7}\end{array}\right)\mathrm{and}\left(\begin{array}{c}{b}_0\\ b_1\\ b_2\\ b_3\\ b_4\end{array}\right)=\left(\begin{array}{c}-3.47418·{10}^4\\ 8.53251·{10}^2\\ -7.84251\\ 3.20206·{10}^{-2}\\ -4.90175·{10}^{-5}\end{array}\right),\end{array}$$(C.3)

this fit is valid for at least 100M_⊙ ≤ m ≤ 192M_⊙. Now, in order to describe the failed SN population, apart from the BH island, we use the work of Schneider et al. (2021). Their results entail both a relation between ZAMS mass and CO core mass and a fit for the failed SN population, which depends on this relation. For the case A/B mass transfer, this comes down to:

$$\begin{array}{c}\hfill f_{\mathrm{FSNe}}^{\mathrm{AB}}(m)=15.832M_{\odot }exp\left(0.058[\frac{3}{20}\frac{m}{M_{\odot }}+4]\right)-14.784M_{\odot },\end{array}$$(C.4)

where the part in square brackets is a linear fit to the CO core mass relation and the last term determines the upper border of the mass gap. Similarly, the failed SN function for case C mass transfer is equal to:

$$\begin{array}{c}\hfill f_{\mathrm{FSNe}}^{\mathrm{C}}(m)=9.727M_{\odot }exp\left(0.041[\frac{16}{31}\frac{m}{M_{\odot }}-\frac{128}{31}]\right)+4.594M_{\odot }.\end{array}$$(C.5)

Equations C.4 and C.5 determine the upper border of the mass gap, at 22M_⊙. Now, for the BH island we simply use a linear function which produces remnants ranging from 8M_⊙ up to the lower border of the mass gap, at 14M_⊙. These different functions are valid in ZAMS mass ranges which are similar to the values used by Schneider et al. (2021) and Marchant et al. (2019), although (small) deviations from these values should not influence the general shape of the remnant distribution. For case A/B mass transfer, this comes down to

$$\begin{array}{c}\hfill f_{\mathrm{RN}}^{\mathrm{AB}}(m)=\{\begin{array}{cc}0& \mathrm{for}m<30.5M_{\odot }\hfill \\ 1.2m-28.6M_{\odot }& \mathrm{for}30.5M_{\odot }\le m\le 35.5M_{\odot }\hfill \\ 0& \mathrm{for}35.5M_{\odot }<m<70M_{\odot }\hfill \\ f_{\mathrm{FSNe}}^{\mathrm{AB}}(m)& \mathrm{for}70M_{\odot }\le m\le 100M_{\odot }\hfill \\ f_{\mathrm{PPISNe}}(m)& \mathrm{for}100M_{\odot }<m\le 192M_{\odot }\hfill \\ 0& \mathrm{for}m>192M_{\odot }\hfill \end{array},\end{array}$$(C.6)

where we assume that not only PISNe but also regular successful SNe do not form remnants, since we are not interested in the distribution of neutron stars. The function for case C is similar to Eq. C.6:

$$\begin{array}{c}\hfill f_{\mathrm{RN}}^{\mathrm{C}}(m)=\{\begin{array}{cc}0& \mathrm{for}m<20.5M_{\odot }\hfill \\ 2.0m-33.0M_{\odot }& \mathrm{for}20.5M_{\odot }\le m\le 23.5M_{\odot }\hfill \\ 0& \mathrm{for}23.5M_{\odot }<m<35.5M_{\odot }\hfill \\ f_{\mathrm{FSNe}}^{\mathrm{C}}(m)& \mathrm{for}35.5M_{\odot }\le m\le 58M_{\odot }\hfill \\ f_{\mathrm{PPISNe}}(m+42M_{\odot })& \mathrm{for}58M_{\odot }<m\le 150M_{\odot }\hfill \\ 0& \mathrm{for}m>150M_{\odot }\hfill \end{array},\end{array}$$(C.7)

where the PPISN function is shifted to lower masses in order to connect with the failed SN function. Equations C.6 and C.7 give the curves that are shown in Fig. 2.

Appendix D: Model uncertainties

In order to apply a Gaussian uncertainty to the results of our estimation in Fig. 4, we make lognormal fits to the average relative standard deviations of the GW data (Abbott et al. 2021b). Figure D.1 shows the data together with our fits. Applying this to our model, we pick a standard deviation randomly from the lognormal distributions and replace the result by a randomly chosen point from a Gaussian with the result as mean and the chosen standard deviation. However, the lognormal fits have a sharp enough peak that taking a single value for the relative standard deviations, for instance, at the maximum, would not influence the results significantly.

Fig. D.1. Relative average standard deviations in the GW data (Abbott et al. 2021b) for both individual masses (left) and chirp masses (right).

Appendix E: Mass transfer efficiency and type

Figures E.1 and E.2 show the results of our estimation (Fig. 4), but for varying η and ζ. The figures show that both parameters do not greatly influence the final mass distributions. The strongest influence of these factors concerns the peak at m ≈ 22M_⊙ and the gap at ℳ ≈ 19M_⊙, but these effects are not big enough to influence our conclusions.

Fig. E.1. Results of the estimation from Fig. 4, for varying η and ζ = 0.5.

Fig. E.2. Results of the estimation from Fig. 4, for varying ζ and η = 0.5.

All Figures

Fig. 1. GW data from Abbott et al. (2021b). Top panel: mass distribution of the individual BBHs (represented by a histogram of the mean values, with a bin-width of 2 M⊙). We can also represent the datapoints as asymmetric Gaussians, as described in Appendix A. The sum of these Gaussians, ∑idNi/dm, is given as well (solid line), together with a cumulative distribution function (CDF, dotted line on a different axis). Bottom panel: chirp mass distribution (grey lines), which is made to resemble the distribution from Abbott et al. (2023). In order to do this, we use asymmetric Gaussians as well, taking into account the 90% credible interval. We also show an adjustable-width kernel density estimation (AWKDE, dashed line) from Shimazaki & Shinomoto (2010), Shimazaki et al. (2018). The Gaussians are plotted on a logarithmic axis, while the sum and the AWKDE use a linear axis. Similarly to Abbott et al. (2023), we omitted GW190814 from our analysis.

In the text

Fig. 2. Remnant masses (mRN), as function of ZAMS mass (mZAMS). The blue curve shows the remnant function for case A and case B mass transfer, based on the simulations by Marchant et al. (2019) and Schneider et al. (2021). The Marchant et al. (2019) results are shown as well (black dots). The triangle shows the transition point between the part of our function which is based on Schneider et al. (2021) and the part which is based on Marchant et al. (2019). The red curve shows the corresponding remnant function for case C mass transfer, where the part after the transition point is the PPISN fit translated to lower masses by 42 M⊙. The dashed lines show mRN/mZAMS = 1, 0.65 and 0.35, since the latter two are used as approximation in estimating the amount of mass transfer in our estimation. The shaded area shows the values of the remnant mass gap. Furthermore, we note that this remnant function is limited to BHs, as for our purposes, we are not interested in the neutron star (NS) masses. The remnant functions are explicitly given in Appendix C.

In the text

Fig. 3. Data from SN surveys (Graur et al. 2015, 2017; Botticella et al. 2017). Top panel: Type II SN rates per unit mass (RSNuM) versus the SFR per unit mass (sSFR), together with a fit for all the data-points (dashed line) and a fit which excludes the two outliers at sSFR ≈ 1.5 × 10−12 yr−1 and 1.5 × 10−9 yr−1 (dotted line). The fits use the least-squares method on a first degree polynomial in logarithmic space and the error bars are one standard deviation. These data are corrected for the Chabrier (2003) IMF, since the original data follows the Salpeter (1955) IMF. For a constant SFR, L = κψ∫τϕL(m)dm, where L(m)∝m3.5 is the luminosity of a star with mass m. However, since τ(m)⋅L(m)∝m, it cancels out with κ and makes the SFR independent of the IMF. Since the SN rate is proportional to κ, it is multiplied by a factor κChabrier/κSalpeter ≈ 1.38, this correction therefore also multiplies RSNuM/sSFR by a factor of 1.38. The bottom panel shows a histogram of the RSNuM/sSFR ratio (with a bin-width of $2\times {10}^{-3}{M}_{\odot }^{-1}$), and a kernel density estimation (KDE).

In the text

Fig. 4. BBH mass distribution resulting from our estimation, where we simulate approximately N = 1.2 × 106 binaries and set η = 0.5 and ζ = 0.5. The latter means that half of the binaries follow the case A/B remnant function, and the other half follow the case C remnant curve (as shown in Fig. 2). The top row shows the resulting mass distribution of the individual primary (red) and secondary (grey) masses, in a stacked histogram with a bin-width of 1 M⊙. The bottom row shows the chirp mass distribution in a similar histogram, where we distinguish between four possible configurations of the primary and secondary mass with respect to the central value of the mass gap (mG): either both BHs are below the gap (dark blue), the secondary is below the gap while the primary is above it (blue) or vice versa (beige), or both are above the gap (light blue). The left column shows the exact results of our estimation, while the right column adds a Gaussian uncertainty to the results in order to make the simulated distribution look more similar to the observed one (as shown in Fig. 1). We based the Gaussian uncertainty on the average standard deviations of the GW data, as described in Appendix D.

In the text

Fig. 5. Primary versus secondary mass distributions, for the results of our estimation (left) and the GW data (right), both in 2D histograms with a bin-width of 1 M⊙. Left panel shows the same as the upper right panel from Fig. 4. We use the notation P′ and S′ here to denote the stars with the largest and smallest final masses, respectively, in contrast to P and S from Fig. 4, which concern the initial masses. The right panel shows the cosmologically reweighted posterior distributions of the GW data (LIGO Scientific Collaboration & Virgo Collaboration 2021; LIGO Scientific Collaboration, Virgo Collaboration, & KAGRACollaboration 2021). Since the posteriors are not normalized, we rescaled them to ${\int }_{0M_{\odot }}^{100M_{\odot }}\frac{\mathrm{d}N}{\mathrm{d}m}\mathrm{d}m\approx 1.4\times {10}^5$. Then, after adding all the posteriors, we normalized the resulting distribution in between 5 M⊙ and 65 M⊙. We omit GW200322 here, since no valid posterior for this source has been given. The grey lines are lines of equal chirp mass, with intervals of 5 M⊙.

In the text

	Fig. A.1. Randomly selected subset of the chirp masses (Abbott et al. 2021b), represented by asymmetric Gaussians (as in Fig. 1), which are similar to the posterior distributions of Abbott et al. (2023). The 90%-interval is determined by either the a-method from Eq. A.5 (top panel) or the b±-method from Eq. A.7 (bottom panel).
In the text

	Fig. B.1. Madau & Dickinson (2014) star formation rate per M⊙, the corresponding SN rate with ml = 8M⊙ and mu = 30M⊙ (through Eq. 1) and the Broadhurst et al. (2018, 2022) BBH merger rate.
In the text

	Fig. D.1. Relative average standard deviations in the GW data (Abbott et al. 2021b) for both individual masses (left) and chirp masses (right).
In the text

	Fig. E.1. Results of the estimation from Fig. 4, for varying η and ζ = 0.5.
In the text

	Fig. E.2. Results of the estimation from Fig. 4, for varying ζ and η = 0.5.
In the text