Bimodal Black Hole Mass Distribution and Chirp Masses of Binary Black Hole Mergers

Gravitational-wave (GW) astronomy has opened up a new window to the universe and led to unprecedented insights (Abbott et al. 2016a, 2016b, 2017a, 2017b). With almost 100 GW detections of binary neutron star (BNS), binary black hole (BBH) and neutron star–black hole (NS–BH) mergers (Abbott et al. 2019, 2021a; The LIGO Scientific Collaboration et al. 2021a), the mass distribution of stellar-mass BHs is being revealed across cosmic time. This will help to better understand many aspects relevant to the formation of BHs such as supernova (SN) explosion physics, and pre-SN evolution of massive single and binary stars. Previously, only ${ \mathcal O }(10)$ stellar-mass BHs were known, mostly residing in Galactic X-ray binaries (Casares et al. 2017) with Cygnus X-1 hosting the most massive BH of ≈ 21 M_⊙ (Miller-Jones et al. 2021).

The chirp mass of GW merger events,

$$\begin{eqnarray}&&{ \mathcal M }=\displaystyle \frac{{({m}_{1}{m}_{2})}^{3/5}}{{({m}_{1}+{m}_{2})}^{1/5}},\end{eqnarray} \tag{ 1 }$$

with component masses m₁ and m₂, is one of the most accurately known GW observables (Tiwari & Fairhurst 2021). Since the observing run O3a of the Advanced LIGO and Advanced Virgo instruments, hints of a gap in the chirp-mass distribution at 10–12 M_⊙ and peaks at 8, 14, 27, and 45 M_⊙ have been reported (Tiwari & Fairhurst 2021). They are corroborated by measurements in O3b (The LIGO Scientific Collaboration et al. 2023; Tiwari 2022) and have complementary features in the inferred BH-mass distribution at about 9, 16, 30, and 57 M_⊙ (The LIGO Scientific Collaboration et al. 2023; Tiwari & Fairhurst 2021; Tiwari 2022; Edelman et al. 2023). Given the current number of BBH merger detections, the peaks at 9 and 30 M_⊙ seem robust whereas the significance of features around ≈16 M_⊙ is less clear (Talbot & Thrane 2018; Sadiq et al. 2022; Wong & Cranmer 2022; Edelman et al. 2023; Farah et al. 2023).

BBH mergers are theoretically expected to form through isolated binary evolution and dynamically in dense stellar environments (Mapelli 2021; Kocsis 2022; Mandel & Broekgaarden 2022). In the former channel ⁴ (Figure 1), the first- and second-born BHs both originate from stars that have been stripped of their envelopes by stable binary mass transfer or a common-envelope phase (called binary-stripped stars, BSSs, hereafter). Dynamically formed BBH mergers additionally invoke BHs from single stars and repeated BH mergers.

Zoom In Zoom Out Reset image size

Envelope removal in BSSs is known to affect the core structures of stars such that they form NSs and BHs of different masses compared to single stars (Timmes et al. 1996; Wellstein & Langer 1999; Brown et al. 2001; Podsiadlowski et al. 2004; Woosley 2019; Ertl et al. 2020; Schneider et al. 2021; Laplace et al. 2021). In this paper, we compute the birth distributions of BH masses from BSSs and single stars at different metallicities. We show that the peaks around the chirp-mass gap and the corresponding peaks in the BH-mass distribution are a natural outcome of envelope stripping in binary stars.

We use published models of single stars and BSSs at solar metallicity (Schneider et al. 2021), Z = Z_⊙ = 0.0142, and newly computed models at Z_⊙/10 following the methodology of Schneider et al. (2021). We employ the stellar evolution code MESA (revision 10398; Paxton et al. 2011, 2013, 2015, 2018, 2019) and evolve models from the beginning of core-hydrogen burning until an iron core forms that collapses at velocities of >950 km s⁻¹. The pre-SN stellar structures are then used as input to the neutrino-driven, semi-analytic SN code of Müller et al. (2016) to determine whether stars explode in SNe or collapse to BHs. We apply calibrations as in Schneider et al. (2021), assume that the entire stellar mass falls into the BH, and use a maximum NS mass of 2 M_⊙. The initial masses of stars are chosen to avoid (pulsational) pair-instability SNe (for further details and all stellar model data see Appendix A).

The pre-SN core structures of BSSs depend on the timing of mass transfer (Schneider et al. 2021). We consider the removal of hydrogen-rich stellar envelopes during core-hydrogen burning (Case A), between the end of core-hydrogen and core-helium burning (Case B), and after core-helium exhaustion (Case C). Case B is further divided into early and late for donors with radiative and convective envelopes, respectively. Because Case A/B BSSs, and single stars/Case-C BSSs result in similar pre-SN core structures (Schneider et al. 2021), we simplify our figures and only show results for single stars and Case-B BSSs.

The delayed, neutrino-driven SN mechanism is the leading idea on how the collapse of the core of a massive star is reversed and leads to an SN explosion (e.g., Colgate & White 1966; May & White 1966; Arnett 1967; Fryer & Warren 2002; Mezzacappa 2005; Janka et al. 2007, 2016; Iwakami et al. 2008; Takiwaki et al. 2014; Müller 2016; O'Connor & Couch 2018; O'Connor et al. 2018; Burrows et al. 2020; Burrows & Vartanyan 2021). For the explodability of stars within this framework, an often-used proxy is the compactness of their cores at core collapse, ⁵

$$\begin{eqnarray}&&{\xi }_{M}=\displaystyle \frac{M/{M}_{\odot }}{R(M)/1000\,\mathrm{km}},\end{eqnarray} \tag{ 2 }$$

where R(M) is the radius at mass coordinate M, usually taken as M = 2.5 M_⊙ (O'Connor & Ott 2011). Low compactness ξ_2.5 suggests that stars explode, while high values indicate unsuccessful explosions and collapse to BHs (O'Connor & Ott 2011; Ugliano et al. 2012; Sukhbold & Woosley 2014; Ertl et al. 2016; Sukhbold et al. 2016; Müller et al. 2016; Chieffi & Limongi 2020; Patton & Sukhbold 2020; Schneider et al. 2021).

The compactness of massive stars and hence their final fates are mainly set by the carbon-oxygen (CO) core mass M_CO and the central carbon mass fraction X_C at the end of core-helium burning (Chieffi & Limongi 2020; Patton & Sukhbold 2020; Schneider et al. 2021). Stellar evolution beyond helium burning is driven by thermal neutrino losses that depend mostly on the density and temperature in the stellar cores, which are set by M_CO. Nuclear energy generation of the subsequent carbon burning additionally depends on the amount of available fuel, i.e., X_C .

In helium burning, three α particles fuse into carbon, and once carbon is available, α-capture via the ¹²C(α, γ)¹⁶O nuclear reaction produces oxygen and consumes carbon. The number of available α particles further depends on the growth of the convective core during helium burning, which is closely linked to shell hydrogen burning. In stars with more massive cores, X_C is smaller such that carbon burning eventually becomes neutrino dominated and turns radiative at a certain M_CO. Enhanced core contraction then leads to high compactness (Sukhbold & Woosley 2014; Chieffi & Limongi 2020; Patton & Sukhbold 2020; Schneider et al. 2021). At somewhat higher M_CO, neon and oxygen burning occur earlier, decreasing the compactness such that a compactness peak forms (see below; E. Laplace et al. 2023, in preparation). For even higher M_CO, also neon burning becomes neutrino dominated, again sharply increasing the compactness (Sukhbold & Woosley 2014; Chieffi & Limongi 2020; Patton & Sukhbold 2020; Schneider et al. 2021).

At Z_⊙, this leads to compactness peaks in single stars and BSSs at M_CO ≈ 7 M_⊙ and 8 M_⊙, respectively, where carbon burning becomes neutrino dominated, and an increase of ξ_2.5 at M_CO ≈ 13 M_⊙ and 15 M_⊙, respectively, where also neon burning becomes neutrino dominated (Figure 2(a); see also Schneider et al. 2021). The M_CO-shifts are caused by a larger X_C in BSSs compared to single stars because the convective helium-burning cores of BSSs grow less in mass as a consequence of having lost the hydrogen-rich envelopes by binary mass transfer (Brown et al. 1996, 2001; Podsiadlowski et al. 2004; Woosley 2019; Chieffi & Limongi 2020; Patton & Sukhbold 2020; Sukhbold & Adams 2020; Laplace et al. 2021; Schneider et al. 2021, E. Laplace et al. 2023, in preparation). At Z_⊙/10, we find the same qualitative ξ_2.5–M_CO relations (Figure 2(b)).

Zoom In Zoom Out Reset image size

Despite the difference of a factor of 10 in metallicity, we find that the CO-core masses at the compactness peaks are similar (M_CO ≈ 7.5 M_⊙); the peak of single stars is shifted by 0.5 M_⊙ to lower M_CO compared to the Z_⊙ models and by 0.8 M_⊙ in BSSs (Figure 2). At lower metallicity with weaker winds, the convective helium-burning core grows relatively more in mass than at higher metallicity such that more α particles are mixed into it and more carbon is burnt into oxygen. With lower X_C , the compactness peaks shift to smaller M_CO (Figure 2).

At Z = Z_⊙/10, the second major increase in compactness of single stars is at almost the same M_CO of ≈13 M_⊙ as in the Z_⊙ models. In BSSs at Z_⊙/10, the transition is also at M_CO ≈ 13 M_⊙ but is smaller by 2 M_⊙ at Z_⊙. At the higher metallicity, wind mass loss during core-helium burning in BSSs with M_CO ≳ 10 M_⊙ is so strong that the entire helium layer and parts of the CO-rich layer are lost, delaying the compactness increase to higher M_CO and reducing the total stellar mass.

Applying our SN model, we find that stars in the compactness peak at M_CO ≈ 7–8 M_⊙ and those with high compactness at M_CO > 13–15 M_⊙ produce BHs (Figure 3). The same conclusion is reached by applying other explodability criteria based on neutrino-driven SNe: a compactness threshold of 0.44 predicts the same explodability as in our analysis in >95% of cases, and the two-parameter explodability criterion of Ertl et al. (2016) agrees in >92% of cases (Appendix A). NS masses reach up to 2 M_⊙, where the lowest masses are from stars with the smallest compactness. In single stars, NS and BH formation occur at almost the same initial masses for both metallicities (at M_ini ≲ 35 M_⊙ and outside the compactness peak at M_ini ≈ 22 M_⊙; Figure 3). In BSSs, we find NS formation at Z_⊙ for initial masses of up to about 70 M_⊙, while this is up to about 40 M_⊙ at Z_⊙/10.

Zoom In Zoom Out Reset image size

In general, at higher Z, winds are stronger such that the BH masses are smaller. Because BSSs additionally lose their hydrogen-rich envelopes, their final and thus BH masses are smaller than those of single stars (Figure 3).

The BH masses of BSSs in the compactness peaks at Z_⊙ and Z_⊙/10 are 8.7–9.8 M_⊙, with the BHs at Z_⊙/10 being more massive by about 0.6 M_⊙ (Figure 3). As shown above, the compactness peaks occur at a certain M_CO, which is set by the size of the convective helium-burning core. This in turn is given by the total mass of the BSS. Because the CO-core masses M_CO at the end of core helium of BSSs in the compactness peak at Z_⊙ and Z_⊙/10 are similar, this implies that also the total helium-star mass and thus the BH mass must be similar, regardless of the exact mass-loss history. Compactness-peak BHs from BSSs thus have a characteristic mass across metallicity. The above argument no longer holds if mass loss during core-helium burning reduces the helium-core mass. The latter can happen in stars with Z > Z_⊙, and beyond some Z, BH formation may no longer be possible (Heger et al. 2003).

The same argument holds for the BH masses of BSSs of ≈ 16 M_⊙ at the second major compactness increase (Figure 3). We have thus identified two characteristic BH masses from BSSs of ≈ 9 M_⊙ and ≈ 16 M_⊙ that we call "low-mass" (BH_L ) and "high-mass" (BH_H ) BHs hereafter.

In about 10%–20% of our Z_⊙/10 models beyond the compactness peak, the SN is not energetic enough to unbind the entire stellar envelope, and SN fallback is expected. At Z_⊙, this fraction is 20%–40%. It remains unknown how much mass falls back and hence what is the resulting BH mass. Furthermore, some BH-forming stars may experience mass loss at core collapse because of neutrino losses leading to a weak SN-like transient (Nadezhin 1980; Lovegrove & Woosley 2013), and some BHs may accrete mass from their companion star during later binary-star evolution (see Figure 1). Such mechanisms are not accounted for in our models and can affect the BH masses.

Assuming a Salpeter initial mass function (Salpeter 1955) for single stars and primary stars in binaries and assuming that initial binary mass ratios and logarithmic orbital separations follow uniform distributions (Öpik 1924), we compute the mass spectrum of BHs (Figure 4). We mark single stars that spent more than 5 × 10⁴ yr in regions of the Hertzsprung–Russell diagram where S-Doradus luminous blue variables (LBVs) are found (Smith et al. 2004; Davies et al. 2018) because they likely experience unaccounted-for mass loss. They may lose their entire envelope and behave like BSSs.

Zoom In Zoom Out Reset image size

Not considering SN fallback, the BH-mass distributions at Z_⊙ and Z_⊙/10 are bimodal with a pronounced BH-mass gap, and the lowest-mass BHs are from stars in the compactness peak (Figure 4). In single stars, BHs are more massive than ≈12 M_⊙, and the gap is at ≈ 14–17 M_⊙ for Z_⊙ and at ≈ 20–27 M_⊙ for Z_⊙/10. The masses of the compactness-peak BHs (≈ 9 M_⊙) and the position of the BH-mass gap at ≈10–15 M_⊙ are similar in BSSs at both metallicities. Compactness-peak BHs (BH_L ) and BHs from the second major compactness increase (BH_H ) can thus explain the inferred peaks at 9 and 16 M_⊙ in the BH-mass distribution of BBH mergers (Tiwari & Fairhurst 2021; The LIGO Scientific Collaboration et al. 2023; Tiwari 2022).

The BH-mass spectrum of BSSs at ≳20 M_⊙ depends on metallicity because of Z-dependent Wolf–Rayet (WR) winds (Figure 4; Vink & de Koter 2005; Belczynski et al. 2010; Spera et al. 2015). At Z_⊙, the maximum BH mass in our models barely reaches 20 M_⊙, while it is almost 45 M_⊙ at Z_⊙/10. This limits the maximum mass of BHs in X-ray binaries formed via binary-star evolution and naturally explains why no BH with higher mass than in Cygnus X-1 has been found in the Milky Way (21.2 ± 2.2 M_⊙; Miller-Jones et al. 2021).

Assuming that 50% of a star's envelope falls back, BSS BHs in the range 5–8 M_⊙ and 6–9 M_⊙ form in our Z_⊙ and Z_⊙/10 models, respectively (Figure 4). Fallback is thus required to explain the BH masses near the observationally suggested NS–BH mass gap at 2–5 M_⊙ (Özel et al. 2010; Farr et al. 2011; Wyrzykowski & Mandel 2020; Schneider et al. 2021).

The two newly identified characteristic BH masses of BSSs of ≈9 M_⊙ (BH_L ) and ≳16 M_⊙ (BH_H ) may lead to three features in the chirp-mass distribution of BBH mergers (BH_L +BH_L , BH_L +BH_H , and BH_H +BH_H ). At Z_⊙, they are at ${ \mathcal M }\approx 7.5\mbox{--}8.1\,{M}_{\odot }$, 9.9–11.7 M_⊙, and 13.2–17.3 M_⊙, while they are at ${ \mathcal M }\approx 7.9\mbox{--}8.9\,{M}_{\odot }$, 10.7–18.1 M_⊙, and 14.6–42.2 M_⊙ at Z_⊙/10 (Figure 5). SN-fallback BHs will broaden the first contribution and extend it to the observed NS–BH mass gap.

Zoom In Zoom Out Reset image size

Isolated binary-star evolution suppresses BH_L +BH_H mergers (Schneider et al. 2021; Appendix B). When forming the less-massive BH_L from stars in the compactness peak first, the initially less-massive companion star cannot gain enough mass by mass transfer to result in a high-mass BH_H . A high-mass BH_H must form first, but then the parameter space for the companion to fall into the compactness peak and give rise to a low-mass BH_L is significantly smaller than forming another high-mass BH_H . Therefore, our models of BSSs predict a dearth of BBH mergers at chirp masses of about 9–13 M_⊙.

We predict NS+BH_L mergers with ${ \mathcal M }\approx 2.7\mbox{--}3.7\,{M}_{\odot }$ and NS+BH_H mergers with ${ \mathcal M }\approx 3.6\mbox{--}7.0\,{M}_{\odot }$ from our BSS models for a minimum NS mass of 1.3 M_⊙ and a maximum BH mass of 45 M_⊙. As explained above, forming the BHs in NS–BH mergers first is favored by binary-star evolution. Hence, NS–BH mergers with BH masses of 9 M_⊙ are predicted to occur more frequently than those with >16 M_⊙ BHs because of the stellar initial mass function.

The NS–BH mergers GW200105 with component masses ${8.9}_{-1.5}^{+1.2}\,{M}_{\odot }$ and ${1.9}_{-0.2}^{+0.3}\,{M}_{\odot }$ and GW200115 with ${5.7}_{-2.1}^{+1.8}\,{M}_{\odot }$ and ${1.5}_{-0.3}^{+0.7}\,{M}_{\odot }$ (Abbott et al. 2021b) appear consistent with a compactness-peak and an SN-fallback BH, respectively. The high mass ratio event GW190814 of a 22.2–24.3 M_⊙ BH and a 2.50–2.67 M_⊙ compact object (Abbott et al. 2020b) might be an example of a high-mass BH_H merging with a very massive NS. Double NS mergers are possible with a total mass of up to 4 M_⊙ in our models, consistent with GW190425 (total mass of ${3.4}_{-0.1}^{+0.3}\,{M}_{\odot }$; Abbott et al. 2020a).

Predicting merger rates requires population-synthesis modeling including the cosmological star formation history and thus the contributions of stellar populations of different Z. This is beyond the scope of this work. Nevertheless, because the BH_L and BH_H masses of BSSs are so similar across a large range of metallicity, we predict peaks in the observed BH and chirp-mass distribution. This is in agreement with the location of the observed chirp-mass gap at 10–12 M_⊙ and also the peaks at 8 and 14 M_⊙ in the chirp-mass distribution of 70 BBH mergers to date ⁶ (Figure 5).

The peaks and gap in the BH-mass distribution of single stars vary strongly with Z and can thus not give rise to the peaks and gap in the observed chirp-mass distribution. Moreover, BBH mergers from single stars populate the observed chirp-mass gap and rather contribute at chirp masses ≳13 M_⊙ (Figure 5). Hence, the low chirp-mass regime in our models is dominated by BHs from BSSs regardless of whether the BBH mergers are from the isolated-binary or dynamical-formation channel. This is not so surprising as BSSs naturally have lower masses than single stars and thus also form less-massive BHs. Furthermore, this interpretation is in line with models of dynamically formed BBH mergers in dense stellar environments that have their dominant contribution at higher masses and underpredict the observed rate at ≈ 10 M_⊙ (Rodriguez et al. 2019; Yang et al. 2019; Antonini et al. 2023). Hierarchical (repeated) BH mergers of BH_L and BH_H may help explain the observed peaks at 30 and 57 M_⊙ in the BH-mass distribution (Tiwari & Fairhurst 2021; Mahapatra et al. 2022; Tiwari 2022), and other mechanisms (e.g., binary evolution) may add further features (van Son et al. 2022b) but leave our predicted characteristic BH masses unaltered.

The compactness landscape found in this work with a peak at M_CO ≈ 7–8 M_⊙ and high compactness at M_CO ≳ 13 M_⊙ agrees with the results obtained by other groups using different stellar evolution codes and/or making different physics assumptions (see e.g., Sukhbold & Woosley 2014; Chieffi & Limongi 2020; Patton & Sukhbold 2020; Schneider et al. 2021; Takahashi et al. 2023). Moreover, Woosley (2019) and Aguilera-Dena et al. (2023) model helium stars at different metallicities as proxies of BSSs, and their results confirm that high compactness is found at characteristic M_CO for a large range of metallicities and that the final (i.e., BH) masses of the helium stars at these M_CO are also similar, as in our work. Ertl et al. (2020) study the collapse of the solar-metallicity helium stars of Woosley (2019) with their one-dimensional, neutrino-driven SN code and confirm our results of a bimodal BH-mass distribution with a large ≈ 9 M_⊙ BH component from compactness-peak BSSs and a gap at about 11 M_⊙ (see also Woosley et al. 2020).

The characteristic BH masses of BSSs of 9 and 16 M_⊙ and hence their bimodal BH-mass spectrum depend on stellar, SN, and nuclear physics (see also Appendix C). Besides metallicity and wind strength (see above), the masses also depend, e.g., on convective core-boundary mixing and the ¹²C(α, γ)¹⁶O nuclear reaction rate. While convective boundary mixing during core-hydrogen burning mainly affects the number of BHs as it relates helium-core masses to initial masses, the size of the convective helium-burning core directly affects the relation of CO-core to BH mass. Changing the convective core overshooting in our single-star models by a factor of 2 results in a 0.3–0.5 M_⊙ shift in M_CO of the compactness peak and thus a change in the BH masses of BSSs by a similar amount. The ¹²C(α, γ)¹⁶O reaction controls X_C and thus at which CO-core masses stars form BHs. Using the slower rate of Kunz et al. (2002) leads to an increase of the compactness-peak M_CO by ≈ 0.9 M_⊙. This reaction is also important for the pair-instability SN–BH mass gap (Farmer et al. 2020). The chirp-mass distribution will help constrain it and other essential physics determining the formation of stellar-mass BHs.

From GW detections of compact-object mergers, one can directly infer their luminosity distances and use them as "standard sirens" (Schutz 1986; Chernoff & Finn 1993; Holz & Hughes 2005). Using universal features in the compact-object mass distribution, such as those predicted here, one can then obtain the redshift to the sources and hence measure the cosmological redshift–luminosity distance relation. This allows for an independent determination of the Hubble constant (Farr et al. 2019; Ezquiaga & Holz 2022). The characteristic BH masses identified here may thus help us better understand the expansion history of the universe.

F.R.N.S. and E.L. acknowledge support from the Klaus Tschira Foundation. This work has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 945806). This work is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy EXC 2181/1-390900948 (the Heidelberg STRUCTURES Excellence Cluster). This research has made use of data or software obtained from the Gravitational Wave Open Science Center (gw-openscience.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.

Software: NumPy (Oliphant 2006), SciPy (Virtanen et al. 2020), Matplotlib (Hunter 2007), Jupyter Notebooks (Kluyver et al. 2016).

All models are nonrotating, and the initial helium mass fraction of the newly computed models is Y = 0.24904. This choice assumes that the helium abundance increases linearly over cosmic time, starting at Y = 0.24668 set by big bang nucleosynthesis (Planck Collaboration et al. 2016) and reaching the solar value of Y_⊙ = 0.2703 (Asplund et al. 2009). Convective boundary mixing during core-hydrogen and core-helium burning assumes step overshooting of 0.2 pressure scale heights. The semi-convection efficiency is 0.1, and the convective mixing-length parameter is 1.8. Nuclear burning follows MESA's approx21_cr60_plus_co56 reaction network, and the nuclear reaction rates are taken from the JINA REACLIB database V2.2 (Cyburt et al. 2010), i.e., the ¹²C(α, γ)¹⁶O reaction rate is from the NACRE II compilation (Xu et al. 2013). All models burn helium as red supergiants. Stellar wind mass-loss rates $\dot{M}$ are as in MESA's "Dutch" wind scheme but with a metallicity scaling of $\dot{M}\propto {Z}^{0.5}$ for stars with effective temperatures T_eff < 10⁴ K (Mauron & Josselin 2011) and following the model of Vink & de Koter (2005) for winds of WR stars. Models are computed until iron-core collapse. If stars do not directly collapse into BHs but emit neutrinos during the formation of a proto-NS, the neutrinos may carry away a significant fraction of the binding energy of a star such that some envelope mass is ejected (Nadezhin 1980; Lovegrove & Woosley 2013). This would systematically decrease our BH masses. Similarly, binary mass accretion after the formation of compact objects can in principle increase their mass. However, in the case of Eddington-limited accretion onto NSs and BHs, one may assume that the masses of NSs and BHs hardly change once they are formed.

To test the dependence of NS and BH formation in our models regarding the employed SN code, we apply two widely used explodability criteria based on (i) the compactness ξ_2.5 (Equation (2)) and (ii) the two-parameter explodability criterion of Ertl et al. (2016). The latter uses the normalized mass inside a shell of specific entropy per nucleon of s = 4,

$$\begin{eqnarray}&&{M}_{4}=\displaystyle \frac{m(s=4)}{{M}_{\odot }},\end{eqnarray} \tag{ 3 }$$

and the normalized mass derivative at this point in the star,

$$\begin{eqnarray}&&{\mu }_{4}=\displaystyle \frac{{\rm{d}}m/{M}_{\odot }}{{\rm{d}}r/1000\,\mathrm{km}}.\end{eqnarray} \tag{ 4 }$$

Both criteria (as well as our SN code) rely on the neutrino-driven SN mechanism. For example, if stars explode via some engine-driven mechanism, the explodability criteria would be different (e.g., they could be related to the rotational properties of a star at iron-core collapse). For a compactness threshold for BH formation of ξ_2.5 > 0.44, we find the same explodability as in our analysis in >95% of cases, and the two-parameter explodability criterion of Ertl et al. (2016), μ₄ < k₁ M₄ μ₄ + k₂ with constants k₁ = 0.200 and k₂ = 0.072, agrees in >92% of cases. Our results regarding successful SN explosions forming NSs and failed ones producing BHs thus seem robust against the exact explodability criteria.

Key stellar properties of the Z_⊙ and Z_⊙/10 models used throughout the paper are provided in Tables 1 and 2, respectively.

We consider isolated binary-star evolution through the common-envelope and stable mass-transfer channels as illustrated in Figure 1. The evolution starts with two massive stars of initial masses M_ini,1 and M_ini,2, where M_ini,1 > M_ini,2. The first mass-transfer phase has to be stable to form a BBH merger (Belczynski et al. 2016; Mapelli 2021). The accretor is usually a main-sequence star that rejuvenates, i.e., that adapts its interior structure to the new total mass M₂ after accretion. The first mass-transfer phase ends after the envelope mass M_env,1 of star 1 has been transferred to star 2. After envelope loss, star 1 produces a BH from a Case-AB binary-stripped star. The second mass-transfer phase is either stable or leads to a common-envelope episode. In both cases, star 2 is also stripped of its envelope and produces another BSS BH. The first-formed BH is usually assumed to not accrete significantly during the second mass-transfer phase because of Eddington-limited accretion.

Binary BH mergers invoking low-mass BH_L and high-mass BH_H are suppressed by the isolated binary-star evolution described above and shown as follows. Let M_ini,1 ≈ 35 M_⊙ (M_ini,1 ≈ 25 M_⊙) at Z = Z_⊙ (Z = Z_⊙/10) such that star 1 produces a low-mass BH_L because it falls into the compactness peak (Figure 3; masses/values in parentheses are for the lower-metallicity models). In our models, star 1 has an envelope mass of M_env,1 ≈ 17 M_⊙ (M_env,1 ≈ 13 M_⊙) such that the accretor can at most reach a mass of ${M}_{2,\max }\approx 35\,{M}_{\odot }+17\,{M}_{\odot }=52\,{M}_{\odot }$ if mass transfer is fully conservative (${M}_{2,\max }\approx 25\,{M}_{\odot }\,+13\,{M}_{\odot }=38\,{M}_{\odot }$). Such massive stars cannot experience Case C mass transfer in our models (Tables 1 and 2), and the second mass transfer thus leads to BSSs that do not produce high-mass BH_H —our models either predict NS or BH_L formation (Figure 3). Hence, first forming BH_L in BH_L +BH_H mergers is forbidden in our models at both metallicities.

Next, we consider that the more massive BH_H forms first. This requires initial masses of M_ini,1 > 70 M_⊙ (M_ini,1 > 40 M_⊙) because star 1 will otherwise produce an NS or a low-mass BH_L (Figure 3). As before, we assume that the initial masses of the secondary star, M_ini,2, are drawn from a uniform mass-ratio distribution, i.e., all M_ini,2 are equally likely for a given M_ini,1. For an arbitrary, fixed mass-transfer efficiency, there is only a very limited range of initial masses for star 2 (M_ini,2 < M_ini,1), to obtain a mass of M₂ ≈ 35 M_⊙ ± ΔM (M₂ ≈ 25 M_⊙ ± ΔM) after the first mass-transfer phase such that star 2 falls into the compactness peak of width ΔM and then leaves behind a low-mass BSS BH_L . In the majority of cases, star 2 has a different mass and does not fall into the compactness peak such that it forms a BSS NS or BH_H after the second mass-transfer phase. Let M_min and M_max be the minimum and maximum (initial) mass of a star to form a NS or BH, respectively, and ${M}_{{\mathrm{BH}}_{{\rm{H}}}}$ the (initial) mass threshold for BH_H formation. The probability of forming a compactness-peak BH_L after the second mass transfer is then ${\rm{\Delta }}M/({M}_{\max }-{M}_{\min })$, while the probabilities to form an NS and BH_H are $({M}_{{\mathrm{BH}}_{{\rm{H}}}}-{M}_{\min }-{\rm{\Delta }}M)/({M}_{\max }-{M}_{\min })$ and $({M}_{\max }-{M}_{{\mathrm{BH}}_{{\rm{H}}}})/({M}_{\max }-{M}_{\min })$, respectively. BH_L +BH_H mergers where BH_H forms first are thus suppressed compared to NS+BH_H and BH_H +BH_H mergers by factors of $({M}_{{\mathrm{BH}}_{{\rm{H}}}}-{M}_{\min }-{\rm{\Delta }}M)/{\rm{\Delta }}M$ and $({M}_{\max }-{M}_{{\mathrm{BH}}_{{\rm{H}}}})/{\rm{\Delta }}M$, respectively. Assuming ${M}_{\min }\approx 10\,{M}_{\odot }$, ${M}_{\max }\approx 100\,{M}_{\odot }$, and ΔM ≈ 5 M_⊙ and taking ${M}_{{\mathrm{BH}}_{{\rm{H}}}}\approx 70\,{M}_{\odot }$ at Z_⊙ and ≈40 M_⊙ at Z_⊙/10, we find suppression factors of about 5–10. We conclude that BH_L +BH_H mergers are not forbidden but suppressed in our models, and the models predict that some BBHs in the 10–12 M_⊙ chirp-mass gap should be detected in the future.

The exact bimodal structure of the BH-mass spectrum of single stars and BSSs depends on various physical processes and as such offers the possibility to better understand and constrain them with the help of observations. In Figure 6, we schematically show how the masses of the low-mass BH_L and high-mass BH_H may change and how the overall number of BHs in comparison to NSs may vary for changes in the wind mass-loss rates of stars, convective boundary mixing and the ¹²C(α, γ)¹⁶O nuclear reaction rate. There are further physical processes and other uncertain ingredients in stellar models that may affect the BH-mass spectrum, but we limit our discussion to the aforementioned aspects.

Zoom In Zoom Out Reset image size

Wind mass-loss rates $\dot{M}$ are still considerably uncertain and could in reality be lower or higher compared to what is assumed in our models (Smith 2014). Similarly, mass-loss rates of stars are higher at higher metallicity and lower at smaller metallicity. The indicated changes in Figure 6(a) are thus meant to cover systematic uncertainties in $\dot{M}$ and stars at different metallicities. In this work, we have shown that BH_L and the lower end of BH_H masses depend only mildly on $\dot{M}$ and that they are systematically smaller for higher $\dot{M}$ and thus larger metallicity (and vice versa). In contrast, the maximum BH mass depends strongly on the wind strength and metallicity, with BH masses being larger at lower $\dot{M}$ and smaller metallicity, and vice versa (Vink & de Koter 2005; Belczynski et al. 2010; Spera et al. 2015).

In our models, we consider extra mixing at the boundary of convective cores during core-hydrogen and core-helium burning. This extra mixing is parameterized as step overshooting by an amount of f_ov as measured in units of the local pressure scale height. In BSSs, i.e., in helium stars, the overshooting during core-helium burning sets the size of the resulting CO-core mass and thus the ratio of total helium-star to CO-core mass. Larger overshooting then implies a smaller total and hence smaller BH mass for the characteristic CO-core masses that result in low-mass BH_L and the lowest-mass BH_H (Figure 6(b)). For the highest-mass BHs, larger overshooting during both core-hydrogen and even more so during core-helium burning implies larger stellar luminosity for a given total mass, hence higher wind mass loss and smaller BH masses.

Step overshooting during core-hydrogen burning sets the helium-core masses of stars. It thus links the helium masses of BSSs to their initial mass and thus to their overall fraction as given by the stellar initial mass function. A larger overshooting value means that initially less-massive stars can produce the same helium-star mass; hence, the overall fraction of BHs increases with larger overshooting during core-hydrogen burning (Figure 6(b)). In BSSs, overshooting during core-hydrogen burning does not directly affect the resulting BH_L and BH_H masses as they are rather connected to the CO-core masses set by core-helium burning.

The ¹²C(α, γ)¹⁶O nuclear reaction rate ${r}_{{}^{12}{\rm{C}}}$ governs the abundance of ¹²C remaining at the end of core-helium burning (Langer 1989; Woosley et al. 1993; Brown et al. 2001). A faster reaction rate converts more carbon into oxygen and vice versa. Fewer carbon atoms during core carbon burning mean that neutrinos overcome energy generation from nuclear burning already at lower CO-core masses, and so the compactness peak systematically shifts to lower CO-core masses. This in turn results in less-massive compactness-peak BHs (Figure 6(c)). The abundance of neon after core carbon burning is less if there are fewer carbon atoms after helium burning. Hence, also the increase of compactness related to neutrinos overcoming the energy generation from core neon burning is found at lower CO-core masses such that the corresponding BH masses are smaller.

The maximum BH mass shown in Figure 6 is connected to wind mass loss and the occurrence of pair-instability SNe (PISNe). Assuming that it is set by PISNe that leave no BH remnants, we indicate in Figure 6 the finding of Farmer et al. (2020) that the lower end of the PISN BH mass gap is at higher masses for a slower ${r}_{{}^{12}{\rm{C}}}$ and vice versa. The masses of BHs beyond the PISN mass gap are not shown.

Quantitatively, decreasing Z by a factor of 10 as done in this work, i.e., decreasing $\dot{M}$ by about a factor of 3, decreases M_CO of the compactness peak by 0.5 M_⊙ in single stars and 0.8 M_⊙ in BSSs. In the BSSs, this translates into an increase of the mass of compactness-peak BHs by a similar amount (≈ 0.6 M_⊙).

Increasing and decreasing the overshooting parameter in our single-star models by factors of 2, we find that the compactness-peak M_CO increases by ≈0.5 M_⊙ and decreases by ≈0.3 M_⊙, respectively. As for the change in Z, we expect a similar change in the mass of compactness-peak BHs of BSSs.

When using the nuclear reaction rate ${r}_{{}^{12}{\rm{C}}}$ of Kunz et al. (2002), which is about 20% slower than the one used in our models (Xu et al. 2013), the M_CO of the compactness peak gets ≈0.9 M_⊙ larger. Similarly, boosting ${r}_{{}^{12}{\rm{C}}}$ in our models by 10% leads to ≈0.4 M_⊙ smaller compactness-peak M_CO. Again, the corresponding BH masses of BSSs are expected to change by a similar amount of mass.