© ESO 2021

1. Introduction

Stars in binary systems are common in the Universe (Sana et al. 2012), but the details of their evolution are uncertain. For massive binaries, it is difficult to observationally constrain the details of physical processes, such as mass-transfer (MT), as the lifetimes of these interacting binary phases are short and hence it is unlikely to observe many of them directly. However, with gravitational-wave (GW) observations, one can search for the imprints of these processes on the properties of their stellar remnant populations, such as the binary black hole (BBH) population.

The LIGO Scientific and Virgo Collaboration (LVC) has recently released the new GW catalogue GWTC-2 (Abbott et al. 2020a), which includes 37 new potential BBH detections¹ from the first half of the third observing run (O3a). In total, GWTC-2 contains 47 BBHs detections (Abbott et al. 2019, 2020a), and the intrinsic rate density of BBH mergers is currently estimated to be $23.9_{-8.6}^{+14.9}$ Gpc⁻³ yr⁻¹ (The LIGO Scientific Collaboration & the Virgo Collaboration 2020a). Each GW detection can constrain some combination of the astrophysical BH parameters: spin and mass. Convenient well-measured quantities are the chirp mass

$$\begin{array}{c}\hfill M_{\mathrm{chirp}}=\frac{{(m_1{m}_2)}^{3/5}}{{(m_1+m_2)}^{1/5}},\end{array}$$(1)

where m₁ and m₂ (m₁ ≥ m₂) are the BHs masses, the binary mass ratio q = m₂/m₁ ≤ 1, and the effective inspiral spin parameter

$$\begin{array}{c}\hfill {\chi }_{\mathrm{eff}}=\frac{m_1{\mathit{a}}_1+m_2{\mathit{a}}_2}{m_1+m_2}·\widehat{\mathit{L}},\end{array}$$(2)

where $\widehat{\mathit{L}}$ is the orbital angular momentum (AM) unit vector and a₁ and a₂ are the BH dimensionless spin vectors. The dimensionless spin vectors are defined as

$$\begin{array}{c}\hfill {\mathit{a}}_{\mathit{i}}=\frac{c{\mathit{J}}_i}{G{m}_i^2},i\in \{1,2\},\end{array}$$(3)

where c is the speed of light, G is the gravitational constant and J_i is the spin AM vector of the BH. There is a degeneracy between χ_eff and q which limits the accuracy to which each quantity can be measured independently (Poisson & Will 1995; Hannam et al. 2013). Nevertheless, the combination of the three observables provide a robust constraint on the properties of a BBH.

Multiple formation channels have been proposed to explain the origin of merging BBHs. They can be divided into two broad categories: (i) isolated binary evolution and (ii) dynamical assembly.

The former occurs during isolated stellar evolution in the field under some specific binary evolution interactions. Interacting binaries that after the formation of the first BH go through (A) stable mass transfer (SMT; e.g. van den Heuvel et al. 2017; Inayoshi et al. 2017; Neijssel et al. 2019) or (B) unstable mass transfer leading to a common envelope (CE) phase (e.g. Smarr & Blandford 1976; van den Heuvel 1976; Tutukov & Yungelson 1993; Kalogera et al. 2007; Postnov & Yungelson 2014; Belczynski et al. 2016) have been shown to form merging BBHs. Another possibility is the formation of BBHs from massive stars with low metallicities and orbital period less than ∼4 days, which due to their tidal interaction, can maintain an almost critical rotation and are going to evolve (C) chemically homogeneously (e.g. de Mink et al. 2009; Mandel & de Mink 2016; Marchant et al. 2016; du Buisson et al. 2020).

The second category of formation channels category occurs in dense stellar environments where stars and binaries can dynamically interact with each other and assemble new binary systems with more massive BHs and tighter orbits, that may eventually merge within the Hubble time. This formation path is present in (D) globular, open, and nuclear stellar clusters (e.g. Kulkarni et al. 1993; Sigurdsson & Hernquist 1993; Portegies Zwart & McMillan 2000; Miller & Lauburg 2009; Banerjee et al. 2010; Rodriguez et al. 2015, 2019; Antonini & Rasio 2016; Arca-Sedda & Gualandris 2018; Fragione & Kocsis 2018) and (E) active galactic nuclei disks (e.g. Bartos et al. 2017; Stone et al. 2017; McKernan et al. 2018; Tagawa et al. 2020). Finally (F) triple or higher-order stellar systems can also lead to the formation of BBHs (e.g. Silsbee & Tremaine 2017; Rodriguez & Antonini 2018; Gupta et al. 2020; Toonen et al. 2020). Within their uncertainties, almost all of these formation channels have been shown to have rate estimates consistent or marginally consistent with the empirical LVC rates.

In this study we consider the formation of BBHs in isolated binary evolution (A and B) though the SMT and CE phase. In these formation channels, two massive stars are born in a relatively wide binary (orbital separations of order ∼1000 R_⊙), where binary interactions happen after the more massive star leaves the main sequence (MS). At this stage, the star expands to become a red supergiant, and inflates its hydrogen-rich envelope beyond its Roche lobe, leading to the first MT episode. The MT stops when the entire stellar envelope is lost, leaving behind a naked helium (He)-star which eventually collapses to form a BH. When the companion reaches the end of its MS, the process repeats itself for the companion star. This MT phase can be either stable or unstable, with the latter leading to the formation of a CE of gas engulfing the binary. If the stripping of the secondary’s envelope is successful, we are left with either a tight BH–He-star system in the case of CE, or with a somewhat wider system in the case of SMT. Eventually the secondary star also collapses to form the second-born BH, and due to energy and AM loss from GW emission (Peters 1964), the BBH system coalesces to form a single, more massive BH.

In the SMT and CE formation channels, the spin of the first-born BH is determined by the AM transport efficiency during the evolution of the progenitor star. Measurements of neutron star and white dwarf spins (Heger et al. 2005; Suijs et al. 2008) and asteroseismology studies (Fuller et al. 2014; Cantiello et al. 2014) suggest that this mechanism must be efficient (Spruit 1999, 2002; Fuller & Ma 2019). Thus, upon expansion, the initial AM of the star is mostly transported to the outer layers which are subsequently lost due to MT and wind mass loss. This leads to the formation of slowly spinning BHs (a₁ ≃ 0), as initially suggested in the context of BH low-mass X-ray binary formation by Fragos & McClintock (2015) and subsequently quantitatively shown in Qin et al. (2018), Fuller & Ma (2019) and Belczynski et al. (2020). In the case of the SMT channel, during the second MT episode, the first born BH may accrete material and spin up (Thorne 1974), depending on the accretion efficiency rate. On the other hand, the spin of the second-born BH is determined by the net effect of the stellar wind and the tidal interaction of the BH–He-star binary system. Because of the efficiency of the AM transport, the He-star emerges from the second MT event with a negligible spin. If the orbital separation is small enough and stellar winds do not widen the system significantly, the He-star can be spun up by tides. These conditions are met at low metallicities for BBHs formed through the CE formation channel (Bavera et al. 2020). In contrast, in the case of SMT, the orbits shrink less efficiently leading to less tidally spun up second-born BHs compared to the CE channel.

All formation channels can be investigated through population synthesis studies which adopt stellar and binary models to rapidly evolve millions of binary stars. This approach gives us insights on the overall population observables given a set of physical assumptions. To explore a wide landscape of parameter values and generate many realisations of the studied population, we need to efficiently evolve millions of binaries. This can be achieved through parametric population synthesis codes which employ fits of single stellar evolution with analytical models to simulate the binary interactions at the expense of coarser approximations when modelling these interactions. Despite this limitation, rapid population synthesis still allows for investigation into how the observable distributions of a population change within the astrophysical model uncertainties. Bavera et al. (2020) recently showed how, given a specific theoretical framework, one can adopt detailed stellar and binary simulations in a population synthesis study to obtain new observable estimates such as the BH spin distributions. In this paper we study how these model predictions are affected by the uncertainties of MT physics such as MT stability and efficiency, CE efficiency and initial orbital distributions.

In Sect. 2 we present the framework we use to generate the population of BBHs and how we convolve the synthetic BBH population with the redshift- and metallicity-dependent star formation rate, as well as incorporate GW detector selection effects. We also summarise the key differences between this work and Bavera et al. (2020). The BBH observable distributions for different MT and CE efficiencies are presented in Sect. 3, where we also show how the χ_eff, M_chirp and q distributions change for these different physical assumptions and determinate thorough model selection which CE efficiency is supported by GWTC-2 BBHs events. The impact of MT stability and initial orbital distributions on the uncertainties of our models is discussed in Sect. 4. We conclude by summarising our findings in Sect. 5.

2. Methods

We generate our populations of BBHs by modelling isolated binary evolution with the POSYDON code.² POSYDON, among many other functionalities, can run and combine detailed stellar and binary evolution simulations performed with the MESA code (Paxton et al. 2011, 2013, 2015, 2018, 2019) to existing parametric binary population synthesis codes. This integration lets us target particular evolutionary phases with more detailed modelling. Similar hybrid approaches have been used in previous population synthesis studies (e.g. Nelson 2012; Chen et al. 2014; Fragos et al. 2015; Shao & Li 2015; Shao et al. 2019). In this work, we use COSMIC (Breivik et al. 2020) to model the evolution of binaries starting from zero age MS (ZAMS) until the formation of the BH–He-star system. We then use MESA to model in detail the subsequent evolution until the formation of the BBH which is the evolutionary phase that determines the spin of the second-born BH (Qin et al. 2018; Bavera et al. 2020).

2.1. Binary black hole population

We create synthetic BBH populations similar to Bavera et al. (2020), but with some key differences.

We assume similar initial binary properties with the exception of the initial orbital periods which here are drawn from an extended Sana et al. (2012) log-power law with coefficient π = −0.55 in the range p ∈ [10^0.15, 10^5.5] days and extrapolated down to p = 0.4 days assuming a log-flat distribution (as the power law is not defined for p < 10^0.15 day). This extension includes the portion of the parameter space leading to chemical homogeneous evolution (du Buisson et al. 2020). All initial binary property assumptions are explained in Appendix A and discussed in Sect. 4.1.

To determine the MT stability, we adopt critical mass ratios q_crit as in Neijssel et al. (2019) with one exception. For stars in the giant branch (GB) and asymptotic giant branch (AGB) we use the same q_crit fits as in Neijssel et al. (2019) but do not adopt Soberman et al. (1997) radial response to adiabatic mass loss for evolved stars beyond the Hertzsprung gap (HG) because they are not currently available in COSMIC. For our reference models, the stable mass-accretion efficiency onto degenerate objects is Eddington-limited. This leads to a highly non-conservative mass-transfer phase where the first-born BH accretes a negligible amount of matter and cannot spin up due to accretion. Unstable MT is parameterised with the standard α_CE–λ formalism (see e.g. Ivanova et al. 2013, for a review). In contrast to Bavera et al. (2020), we adopt λ fits as in Claeys et al. (2014) while we explore different α_CE efficiencies: α_CE ∈ [0.2, 0.35, 0.5, 0.75, 1.0, 2.0, 5.0]. Since, approximately, α_CE scales linearly with the orbital separation post CE, see Eq. (B.2), low CE efficiencies lead to tighter orbital separations post CE. Therefore, we expect that more BH–He-star systems will undergo tidal spin up at lower α_CE. We describe the details of our COSMIC model, MT stability and CE in Appendix B. In Sects. 4.2 and 4.3 we discuss how our BBH population distributions and rates are affected by these assumptions.

The late-end phase of the binary evolution of BBHs formed through CE and SMT channels are BH–He-star systems. We update our MESA models (Qin et al. 2018; Bavera et al. 2020) to match the stellar model assumptions of du Buisson et al. (2020). In contrast to Bavera et al. (2020), we relax the He-star models to zero age helium MS (ZAHeMS) before initiating the binary evolution. This ensures that the He-star model is in thermal and hydrostatic equilibrium when the binary interactions begin. In order to verify that the He-star will not overfill the L2 Roche volume throughout the binary evolution, we include the prescription of Misra et al. (2020). The ingredients of our MESA model are explained in Appendix C.

Once the He-star reaches carbon depletion, the MESA simulations are stopped. We then collapsed the profile of the He-star according to the procedure used in Bavera et al. (2020) which accounts for disk formation. Here we adopted a different treatment of neutrino mass loss where we assume that the innermost 3 M_⊙ forms a proto-neutron star which collapses to form a BH of 2.5 M_⊙ while 0.5 M_⊙ are converted into neutrinos and escape the system carrying away AM (cf. Zevin et al. 2020a). The complete procedure used to collapse the He-star profiles is explained in detail in Appendix D.

We use our detailed binary stellar models to cover the four-dimensional parameter space of initial metallicity Z, BH mass M_BH, He-star mass M_{He − star} and orbital period p. We run grids for 30 different metallicities ranging from log₁₀(Z) = − 4.0 to log₁₀(1.5 Z_⊙) ≃ − 1.593 in steps of log₁₀(Z) ≃ 0.083 where we adopt the solar reference Z_⊙ = 0.017 (Grevesse et al. 1996). For each metallicity we run 11 BH masses in the log-range [2.5, 54.4] M_⊙, 17 He-star masses in the log-range [8, 80] M_⊙ and 20 initial binary periods in the log-range [0.09, 8] days. In total, we calculated roughly 110, 000 new binary evolution sequences. These grids were used to determine the final outcomes and final parameters of the late-end evolution stage of the binary systems through linear interpolation for each metallicity, independently. The features of these grids and the interpolation accuracy are discussed in Appendix E.

In Fig. 1 we show an example of a two-dimensional slice from our four-dimensional grid. The parameter space is sliced at Z = 0.001 and M_BH = 29.4 M_⊙. We show the final second-born BH mass and spin, as well as the BBH merger timescale as a function of initial orbital period and He-star mass at ZAHeMS. We see that the binary interactions determining the spin of the second-born BH create a gradual transition between tidally locked systems and non-spinning systems. The complex interactions of stellar winds, tides, internal differential rotation and (in some cases) mass transfer are important in determining the spin of the second-born BH, and therefore require a detailed treatment which traditional rapid population codes cannot offer.

Fig. 1. Example of a two-dimensional slice of the four-dimensional grid showing their initial BH–He-star orbital period pi (in days) and initial He-star mass mHe − star (in M⊙) for log10(Z) ≃  − 3 and mBH1 = 29.39 M⊙. The final mass mBH2 and spin a2 of the second-born BH, as well as merger timescale Tmerger, are coloured according to the legend of each panel. All successful MESA simulation stopped at carbon depletion (square markers) while other termination flags are shown in the bottom legend. The merger timescale colour bar is capped at 14 Gyr.

Once a BBH system is formed, GW inspiral leads to the system’s eventual merger. We calculate the merger timescale T_merger as in Peters (1964) accounting for eccentricity. In our models, this timescale is anti-correlated with the observable χ_eff. This is caused by tides, as they are the only mechanism able to spin up the progenitor of the second-born BH. Their efficiency is highly dependent on the orbital separation, see Eq. (C.1). For tidally locked systems one recovers $T_{\text{merger}}\propto a_2^{-8/3}\propto {\chi }_{\text{eff}}^{-8/3}$ while for the other systems $T_{\text{merger}}\propto a_2^{-8/17}\propto {\chi }_{\text{eff}}^{-8/17}$ (Bavera et al. 2020). We therefore expect systems with χ_eff > 0 to have, on average, shorter merger timescales compared to systems with χ_eff = 0. This anti-correlation is the key to understanding the translation from the underlying BBH population (what one would observe with an infinitely sensitive detector) to the observed population. Current GW detectors are probing small redshifts (z ≲ 1) and cannot explore the peak of the Universe’s SFR at z ≃ 2 where most of the highly spinning BBHs form and merge (as they are preferentially formed in low metallicity environments, e.g. Stevenson et al. 2017).

2.2. Rate estimates

To compute the expected rate of detectable GW events, we need to convolve the redshift- and metallicity-dependent star-formation rate (SFR) with the selection effects of the detector array. To do this, we follow the approach shown in Appendix B of Bavera et al. (2020). We assume a flat ΛCDM cosmology with H₀ = 67.7 km s⁻¹ Mpc⁻¹ and Ω_m = 0.307 (Planck Collaboration XIII 2016), a cosmic SFR history as in Madau & Fragos (2017) and metallicities following a truncated log-normal distribution with standard deviation 0.5 dex around the empirical mean metallicity function derived by Madau & Fragos (2017). The log-normal distribution is truncated at the highest metallicity bin edge and the distribution is accordingly renormalised to ensure that ${\int }_{-\infty }^{Z_{\text{max}}}N({\log }_{10}(Z)|\mu (z),\sigma )\text{d}{\log }_{10}Z=1$, where Z_max is our highest metallicity edge bin. Portions of the distribution extending beyond the lower limit edge are included in the edge bin when integrating over metallicity. The population synthesis predictions are performed in finite time bins of Δt_i = 100 Myr and log-metallicity bins ΔZ_j. The detection rate of BBH mergers for a given detector network is calculated from the Monte Carlo simulations, cf. Eq. (13) in Bavera et al. (2020),

$$\begin{array}{c}\hfill R_{\det }={\displaystyle \sum _{\mathrm{\Delta }{t}_i}\sum _{\mathrm{\Delta }{Z}_j}\sum _k{f}_{\mathrm{corr}}\frac{\mathrm{fSFR}(z_{\mathrm{f},i})}{M_{\mathrm{sim},\mathrm{\Delta }{Z}_j}}4\pi c{D}_{\mathrm{c}}^2(z_{\mathrm{m},i,k})p_{\det ,i,k}\mathrm{\Delta }{t}_i{\mathrm{yr}}^{-1},}\end{array}$$(4)

where the argument of the summation is the cosmological weight contribution of the k-th binary born at redshift z_f, i with BH mass m_1, k and m_2, k, spin a_1, k and a_2, k and merging at redshift z_m, i, k. Furthermore, M_sim, ΔZj is the simulated mass per log-metallicity bin ΔZ_j and f_corr the normalisation constant which converts the simulated mass to the total stellar population (see Appendix A in Bavera et al. 2020). Here, D_c(z) is the comoving distance to the source, fSFR(z) is the SFR per log-metallicity range ΔZ_j and p_{det, i, k} ≡ p_det(z_m, i, k, m_1, k, m_2, k, a_1, k, a_2, k) accounts for the selection effects of the detector array.

In contrast to Bavera et al. (2020), we calculate the sensitivity of a GW detector to a source accounting for its network configuration as well as include the selection effects on the BH spins. We assume a 3-detector network configuration composed of LIGO–Hanford, LIGO–Livingston, and Virgo with simulated O3 sensitivity (mid high/late low from Abbott et al. 2018). Detector response functions are calculated using the PyCBC package (Nitz et al. 2020a). For each compact binary merger, we calculate the signal-to-noise ratio (S/R) as

$$\begin{array}{c}\hfill {\rho }^2=4\mathfrak{R}{\displaystyle {\int }_0^{\infty }\frac{{\stackrel{\sim }{h}}^{\star }(f)\stackrel{\sim }{h}(f)}{S_n(f)}\mathrm{d}f}\end{array}$$(5)

for each detector in the network, where S_n(f) is the one-sided power spectral density of the noise, and $\stackrel{\sim }{h}(f)$ is the GW strain, determined using the IMRPhenomPv2 waveform approximant (Hannam et al. 2014; Khan et al. 2016). The network S/R is the quadrature sum of the S/Rs in all the three detectors. Assuming a network detection threshold of ρ_det = 12, we Monte Carlo sample the sky location, inclination, and phase N times for each system and calculate ρ_net. The detection probably p_{det, i, k} is thus determined as

$$\begin{array}{c}\hfill p_{\det ,i,k}=\frac{1}{N}{\displaystyle \sum _{l=1}^N\mathcal{H}({\rho }_{\mathrm{net},l,i,k}-{\rho }_{\det }),}\end{array}$$(6)

where each l represents a random draw of extrinsic parameters and ℋ is the Heaviside step function; we perform N = 1000 realisations of extrinsic parameters for each system.

The total BBH merger rate density R_BBH(z) is the number of BBHs per comoving volume per year as a function of redshift. This quantity can be calculated knowing the contribution of each binary k placed at the centre of each formation time bin Δt_i in its corresponding metallicity bin ΔZ_j assuming p_{det, i, k} = 1,

$$\begin{array}{c}\hfill R_{\mathrm{BBHs}}(z_i)={\displaystyle \sum _{\mathrm{\Delta }{Z}_j}\sum _k{f}_{\mathrm{corr}}\frac{\mathrm{fSFR}(z_{\mathrm{f},i})}{M_{\mathrm{sim},\mathrm{\Delta }{Z}_j}}\frac{4\pi c{D}_{\mathrm{c}}^2(z_{\mathrm{m},i,k})}{\mathrm{\Delta }{V}_{\mathrm{c}}(z_i)}\mathrm{\Delta }{t}_i{\mathrm{Gpc}}^{-3}\hspace{0.17em}{\mathrm{yr}}^{-1},}\end{array}$$(7)

where ΔV_c(z) is the comoving volume shell corresponding to the cosmic time bin Δt_i,

$$\begin{array}{c}\hfill \mathrm{\Delta }{V}_{\mathrm{c}}(z_i)\equiv {\displaystyle {\int }_{\mathrm{\Delta }{z}_i}\frac{1}{1+z}\frac{\mathrm{d}{V}_{\mathrm{c}}}{\mathrm{d}z}\mathrm{d}z=\frac{4\pi c}{H_0}{\int }_{\mathrm{\Delta }{z}_i}\frac{D_{\mathrm{c}}^2(z)}{E(z)(1+z)}\mathrm{d}z,}\end{array}$$(8)

where $E(z)=\sqrt{{\mathrm{\Omega }}_m{(1+z)}^3+{\mathrm{\Omega }}_{\mathrm{\Lambda }}}$. Here Δz_i is the redshift interval corresponding to the cosmic time bin Δt_i centered at z_i ≡ z_f, i.

3. Results

We use our models to predict the distributions of some of the main GW observables: the effective inspiral spin parameter χ_eff, the chirp mass M_chirp and the binary mass ratio q. We investigate how these distributions change for different CE and accretion efficiencies. The distributions for SMT and CE channels are obtained by distributing the synthetic BBH population across the cosmic history of the Universe as described in Sect. 2.2.

Our models use detailed binary evolution simulations to determine the spin of the second-born BH, assuming that the first-born BH is formed with a negligible spin a₁ ≃ 0 because of the assumed efficient AM transport (Qin et al. 2018; Fuller & Ma 2019). If the second MT is stable the first-born BH can accrete material and spin up (Thorne 1974). Nevertheless, because in our reference models we assume Eddington-limited accretion efficiency onto compact objects, the accreted mass is small; this leads to small a₁ ≃ 0 also for the SMT channel. The Eddington limited accretion onto the BH is a crucial assumption for the existence of this channel. In Sect. 3.1.2 we show that if highly super-Eddington accretion onto the BH is allowed, the SMT channel contributes to a negligible part to the BBH rate density compared to the CE channel.

3.1. Underlying BBH population

3.1.1. Common envelope channel

In our CE channel, only the spin of the second-born BH is non-zero and, hence, contributes to the χ_eff parameter. The BH progenitor can be tidally spun up during binary evolution after the CE event. The efficiency of tides depends strongly on the orbital separation as the synchronisation timescale T_syn ∝ A^17/2, see Eq. (C.1). On the other hand, stellar winds can cause the binary to lose mass, widening the orbit and reducing or neutralising the effects of tides. He-stars have wind mass loss rates that strongly depend on metallicity. Hence, a second-born, tidally spun-up BH can only occur at low metallicity as shown in our detailed binary simulation (see Fig. D.1). A key point is that our detailed simulations do not show a dichotomy between tidally locked and non-spinning second-born BHs, but smoothly cover the whole range of a₂ ∈ [0, 1], (e.g. top panel of Fig. 1). This was pointed out by Qin et al. (2018) and Bavera et al. (2020) and is in contrast with results from semi-analytical models (Hotokezaka & Piran 2017; Zaldarriaga et al. 2018; Gerosa et al. 2018).

In the top panels of Fig. 2 we show the joint distribution of χ_eff and M_chirp for the underlying BBH population of the CE channel for the reference model with α_CE = 1 alongside the SMT channel and their combination. For the CE channel, we can see that the underlying BBH population has a non-negligible amount of positive χ_eff mergers due to the tidal spin-up of the second born BH’s progenitor. This is in agreement with our previous models (Bavera et al. 2020) obtained using the COMPAS code (Stevenson et al. 2017; Vigna-Gómez et al. 2018) which is based on the same stellar model fits but implements binary interactions differently (a comparison between the two codes is beyond the scope of this work). Because of the anti-correlation between the merger timescale T_merger and χ_eff (Sect. 2.1), these highly spinning systems merge soon after their formation. Current GW detectors probe only small redshifts (z ≲ 1), well below the peak of the cosmic SFR (z ∼ 2) where most of these systems are created and merge, as low metallicity environments are required for efficient tidal spin up.

Fig. 2. Model predictions for the underlying (intrinsic) BBH population (grey) and the O3 detected BBH population (orange) for our reference model with αCE = 1.0 and ηacc = 1. We show the joint distributions of chirp mass Mchirp and effective inspiral spin parameter χeff for the combined CE and SMT channels (left), CE formed BBHs (centre) and SMT formed BBHs (right). Lighter colours represent larger contour levels of 68%, 95% and 99%, respectively, constructed with pygtc module (Bocquet & Carter 2016). All histograms are plotted with 30 bins in the same range without any bin smoothing. We overlaid in grey the O1, O2 and O3a LVC GWTC-2 data with their 90% credible intervals. The 9 events of GWTC-2 in tension with our models are indicated in red (see Sect. 3.3), GW190521 is outside the plotted window.

CE efficiency has an important role in the determination of the post-CE orbital distribution. This is because α_CE correlates approximately linearly with the post-CE orbital separation A_postCE, see Eq. (B.2). Models assuming inefficient CE ejection (α_CE < 1.0) result in tighter post-CE orbits and more systems merging during CE compared to larger α_CE values (see Fig. 3). This occurs because the binaries need to deposit more orbital energy into the envelopes to successfully eject them. The opposite is true when assuming an efficient CE ejection (α_CE > 1.0). We expect, on average, larger χ_eff values for models with lower CE efficiency parameters, as more systems will undergo tidal spin up and small χ_eff for models assuming ultra-efficient CE ejection, as tides are weaker at larger orbital separations. Indeed, this trend is what we find. In Table 1, we report the median χ_eff, M_chirp and q with their 90% confidence interval (CI) for our different CE efficiencies models. For the underlying (intrinsic) BBH population, we observe a monotonic decrease of all these quantities for increasing α_CE (from 0.2 to 5.0). This trend is also found if we look at the relative fractions of massive and highly spinning BBHs, namely with χ_eff > 0.1 and M_chirp > 15 M_⊙, in Table 2. In both tables we see that on average models with small α_CE have larger M_chirp and q. This is because for the same orbital separation, massive binaries have a larger orbital energy reservoir compared to lighter systems and, hence, can deposit more energy into the envelope without shrinking to the point where they merge in the CE phase.

Fig. 3. Orbital separation of BH–He-star binaries post CE for αCE ∈ [0.2, 0.5, 1.0] represented with solid lines according to the legend. The histogram has units of M⊙ Gpc−3 and accounts for the total stellar mass formed per comoving volume integrated over the Universe cosmic history per log-orbital period bin. The grey shaded area represent the 90% CI of the systems forming merging BBHs with Tmerger < 14 Gyr for the αCE = 0.2 model (the other models have similar CIs). As the CE efficiency is lowered, the orbital separations become smaller and the distributions move to the left: for αCE = 0.2 the orbital separation peak enters the grey area boosting the merger rate.

In Table 1 we also report the local rate density for the cosmic time bin centreed at z = 0.01. The reference model with α_CE = 1.0 has a local rate density of 42.6 Gpc⁻³ yr⁻¹. If we increase α_CE, the post-CE orbital separations are larger, hence, the rate density decreases because fewer systems merge within the Hubble time. On the other hand, if we decrease α_CE, more systems merge during the CE event and the rate density decreases as well. This trend is not followed by the model with α_CE = 0.2 where the rate density suddenly jumps up to 113.0 Gpc⁻³ yr⁻¹. To understand the sudden increase in the rate density of this model, we need to carefully look at the post-CE binary orbital separations. In Fig. 3 we show a histogram of all BH–He-star orbital separations surviving CE for α_CE ∈ [0.2, 0.5, 1.0] (solid lines). The synthetic BBH population is weighted according to Eq. (B.10) of Bavera et al. (2020) which integrates the redshift- and metallicity-dependent SFR across the cosmic history of the Universe. In grey, we show the 90% CI of systems forming merging BBHs with T_merger < 14 Gyr for the model with α_CE = 0.2 (the other models have similar CIs). Systems with orbital separations smaller than the left boundary of the CI either form BH-NS binaries, merge during the BH–He-star evolution or widen the orbits (because of wind driven mass loss rate) past the point where they will merge within the Hubble time. Systems on the right of the CI form double compact objects with merging timescales larger than the Hubble time. In this figure we see that as the CE efficiency decreases, the orbital separations decrease. The total orbital separation distributions present a large peak of orbital separation preceded by a smaller flatter distribution of orbital separation. For α_CE = 0.2 we see that this large peak of orbital separation enters the merging BBH population (grey area of Fig. 3). This is the source of the sudden increase of the rate density.

The peak of BH–He-star orbital separations post CE is a metallicity product. All binaries going through CE are evolving during the He-burning phase. In our models, systems that are in the HG at onset of CE are considered to merge during the CE (because we assume the pessimistic CE scenario; see Belczynski et al. 2007). The maximum stellar radius of a star in the HG is metallicity dependent. Even though, on average stars with high metallicity have larger radii during this phase compared to lower metallicity stars, they have similar super-giant phase radii, see for example, Fig. 7 of Linden et al. (2010). This implies that binaries with high metallicities sample a smaller range of orbital separations at onset of CE, with the donor star having passed the HG phase, compared to binaries with lower metallicities. Therefore, low metallicity BH–He-star binaries sample a wide (approximately flat) A_preCE distribution which result to a wide (also approximately flat) range of A_postCE, while high metallicity systems sample a narrow A_preCE distribution which result in a narrow A_postCE range. Figure 3 shows the combination of these distributions for all metallicities. Since the average metallicity in the Universe is a monotonically increasing function, the yield of binaries at low metallicities is smaller than the yield at larger metallicities, hence the larger A_postCE peak.

3.1.2. Stable mass transfer channel

In our SMT channel, the spin of both BHs can be non-zero and hence affect the χ_eff parameter. Since in our reference model we assume Eddington-limited MT, the amount of accreted material is negligible compared to the BH mass and leads to small spins (a₁ ≃ 0.002 is the largest value in our population). In Table 1 we show that models with super-Eddington accretion efficiency limits result in larger first-born BH spins, as the BH is allowed to accrete at highly super-Eddington rates, and, hence, result in larger median χ_eff. After the second MT phase, tides can further spin up the second-born BH progenitor if the orbits are tight enough (this requires p < 1 day). Since the orbits cannot shrink as efficiently as in the CE channel, most of the systems formed through this evolutionary path will not undergo tidal spin up. Since the CE efficiency does not affect this evolutionary path, here we report only values from the model with α_CE = 1.0.

In Fig. 2 we show the underlying joint distribution of χ_eff and M_chirp for the SMT channel alongside the CE channel and their combination. We see that the underlying SMT BBH population presents a non-zero χ_eff contour at the 95% level. These are the systems accreting during the second MT event with tidal spin up during the subsequent phase. Even though the non-zero χ_eff distribution is a small part of the overall population (cf. median χ_eff in Table 1) this subpopulation has an astrophysical consequence. During the second MT, the accreting BHs are thought to form an accretion disk with strong X-ray emission. This partly explains the high-end of the luminosity function of stellar X-ray sources in galaxies (e.g. ultraluminous X-ray sources; Begelman 2002; Swartz et al. 2004; Kovlakas et al. 2020).

In Table 1 we report the rate density contribution of SMT channel to be 24.6 Gpc⁻³ yr⁻¹. This value is comparable to the contribution of the CE channel for α_CE ∈ [0.35, 5.0]. This result is consistent with other studies (e.g. Neijssel et al. 2019) and is strongly dependent on the assumed accretion efficiency limit onto the BH. If we allow for super-Eddington accretion the BBH rate density will decrease. The drop in the rate density when allowing for super-Eddington accretion occurs because conservative MT does not shrink the orbit as efficiently as non-conservative MT (in the case of Eddington-limited accretion) and thus the BBHs formed post-MT are more commonly too wide to merge within the Hubble time. This can be seen in Fig. 4 where we show a histogram of all BH–He-star orbital separation post SMT for different accretion efficiency limits η_acc ∈ [1, 10⁵, 10⁹] (in units of Eddington-limit), solid lines, and the 90% CI of systems forming merging BBHs with T_merger < 14 Gyr, in grey. For not Eddington-limited accretion (here arbitrarily limited to up to 10⁹ times the Eddington-limit) the rate density contribution of this channels drops by two orders of magnitudes down to 0.2 Gpc⁻³ yr⁻¹ while just decreasing the CE channel rate by ∼10%. This is a negligible contribution to the BBH merger rate compared to the yield of the CE channel.

Fig. 4. Orbital separation of BH–He-star binaries post SMT for accretion efficiency limit ηacc ∈ [1, 105, 109] (in units of Eddington-limit) represented with solid lines according to the legend. The histogram has units of M⊙ Gpc−3 and accounts for the total stellar mass formed per comoving volume integrated over the Universe cosmic history per log-orbital period bin. The grey shaded area represent the 90% CI of the systems forming merging BBHs with Tmerger < 14 Gyr for the model with ηacc = 1 (the other models have similar CIs). As the accretion efficiency is increased, the orbital separations become larger and the distributions move to the right: for highly super-Eddington accretion efficiency ηacc = 109 the tail almost exits the grey area, hence, decreasing the merger rate.

When allowing for super-Eddington accretion, the first-born BH accretes a non-negligible amount of matter leading to a different mass radio distribution compared to Eddington-limited SMT models. In Table 1 we can see that the median mass ratio decreases from 0.72 to 0.35 and 0.14 for increasing η_acc. Small mass ratios are also found by the BPASS team (Eldridge et al. 2017) who, in their models, allow for super-Eddington accretion onto BHs. As mentioned above, the first-born BH is spun-up when mass is accreted onto it. Indeed, for increasing η_acc, our SMT models predict larger median χ_eff.

3.2. Detected BBH population

O3a lasted 6 months and saw the detection of 37 GW signals from BBH mergers resulting in a total of 177.3 days of data suitable for coincident analysis. These detections translate to a detection rate of 76 yr⁻¹. In our model comparison with the data we include also the BBH detections of the first two observing runs O1 and O2 (Abbott et al. 2019) for a total of 47 events. Although evidence for additional signals in O1 and O2 data has been presented by other groups (e.g. Zackay et al. 2019a,b; Venumadhav et al. 2020; Nitz et al. 2020b), we do not include them in order to simplify our analysis and only consider a single treatment of selection effects. At the same time, adding a few additional, low significance events in the observed population is not expected to add significant discriminating power. For all events, we assume simulated O3 sensitivity (mid high and late low from Abbott et al. 2018) as the observable distributions of χ_eff and M_chirp are weakly dependent on the detector sensitivity for the channels considered (Bavera et al. 2020).

The detected joint distributions of χ_eff and M_chirp of our reference model with α_CE = 1.0 are presented in the bottom panels of Fig. 2. The figure shows CE and SMT channels alongside their combination. For a visual comparison with the observations, we add the LVC GW detections with their 90% credible intervals in grey. We can clearly see that the SMT channel only contributes with zero χ_eff and large M_chirp systems to the detected population. Hence, in our model, the only source of non-zero χ_eff in the detected BBH population comes from the CE channel.

In Fig. 5 we show the predicted two-dimensional distributions of χ_eff, M_chirp and q for the combined CE and SMT detected population at different CE efficiencies. We can clearly see how the models with inefficient CE (α_CE < 1) lead to overall larger χ_eff values compared to ultra efficient CE in the detected population (α_CE > 1). Next to each panel we also include the normalised one-dimensional histogram of each quantity where we also show the underlying BBH population for comparison. The one-dimensional M_chirp histograms show how the selection effect favouring higher BH masses changes the distribution, while in the one-dimensional χ_eff histograms we can see how the detectable population mostly probes low χ_eff. This happens because the GW detectors probes small redshifts (z ≲ 1) while highly spinning systems tend to form at high redshifts (z ∼ 2) and low metallicity environments, and merge quickly at a redshift close to the one of their formation (see the discussion about the anti-correlation between χ_eff and T_merger in Sect. 2.1). These high redshifts are outside current GW detection horizons.

Fig. 5. Model predictions for the O3 detected BBH population of CE and SMT channels combined for αCE ∈ [0.2, 0.35, 0.5, 1.0, 2.0, 5.0] and ηacc = 1, in orange. We show the joint distributions of chirp mass Mchirp, effective inspiral spin parameter χeff and binary mass ratio q. Lighter colours represent larger contour levels of 68%, 95% and 99%, respectively, constructed with pygtc (Bocquet & Carter 2016). All histograms are plotted with 30 bins without any bin smoothing. We overlaid in grey the LVC GWTC-2 data with their 90% credible intervals. The 9 events of GWTC-2 in tension with our models are indicated in red (see Sect. 3.3). We also show the one-dimensional projection of each quantity overplotted on the underlying (intrinsic) BBH population, in grey. For visualisation purposes, the χeff probability density function (pdf) is plotted in log-scale. For a visualisation of each model separately see Appendix G.

In Table 1 we report the detection rate of each model for O3 sensitivity. For the SMT channel our model predicts a detection rate of 86 yr⁻¹ while the detection rate for the CE model with α_CE = 1.0 is 108 yr⁻¹. As we increase or decrease CE efficiency we lower the detection rate to 56 yr⁻¹ for α_CE = 0.35 and 15 yr⁻¹ for α_CE = 5.0. The model with α_CE = 0.2 overpredicts the detections with 412 yr⁻¹ (see Sect. 3.1.1 for an explanation). On the other hand, the SMT model with highly super-Eddington accretion efficiency limit predicts a detection rate of 0.2 yr⁻¹ which is negligible compared to the CE channel contribution. Within the probed mass-transfer physics uncertainties, the combination of the two channels is roughly consistent with the observed rate. While our models are consistent with observations, the event rate does depend on many other uncertain evolutionary parameters (e.g. Dominik et al. 2015; Giacobbo et al. 2018; Barrett et al. 2018) and metallicity-specific star formation history (e.g. Chruslinska et al. 2019; Neijssel et al. 2019) which we have not explored. Therefore, while our results can illustrate the expected trend when varying CE efficiency, we cannot make definitive statements on the true value of α_CE without also considering the other evolutionary parameters.

3.3. Evidence for additional formation channels

In GWTC-2 there are BBH events that appear marginally consistent or inconsistent explained by our models of isolated binary evolution through CE or SMT. Using individual events to discriminate between models should be done with caution, as the information that individual events carry can be strongly affected by the choice of priors used in the analysis (e.g. Mandel & Fragos 2020; Zevin et al. 2020b; Fishbach & Holz 2020). Instead, one should attempt to derive conclusions based on the combined detected population. In this section we discuss such events which may originate from other active formation scenarios (see the discussion in Sect. 1), while in the following section we perform a model comparison based on the combined sample of events.

In the catalogue we have two high-significance events with asymmetric masses: GW190412 and GW190814. GW190412 has a binary mass ratio of $q=0.{28}_{-0.06}^{+0.13}$ (The LIGO Scientific Collaboration & the Virgo Collaboration 2020b) while GW190814 has $q=0.{112}_{-0.009}^{+0.008}$ (Abbott et al. 2020b). We find that these small mass ratios are consistent at the 90% level of BBHs formed though SMT with highly super-Eddington accretion. In these models the first-born BH accretes material during the second MT phase leading to unequal mass ratios. However, these models predict large χ_eff values as the first-born BH is spun up during accretion (Thorne 1974). The 90% CI of χ_eff in this model is not consistent with the ${\chi }_{\mathrm{eff}}=0.{25}_{-0.11}^{+0.08}$ and ${\chi }_{\mathrm{eff}}=-0.{002}_{-0.061}^{+0.060}$ of GW190412 and GW190814, respectively. If one assumes a different, astrophysically-motivated prior, such as a prior that assumes non-spinning primary BHs (different formation channels have different priors), rather than the one used in the LVC analysis, GW190412’s inferred mass ratio increases to 0.34 ≤ q ≤ 0.47 at the 90% level (Zevin et al. 2020b). The latter is marginally consistent with our models. The case of GW190814, with its lower-mass component being a 2.6 M_⊙ compact object, remains a challenge to explain with isolated BBH formation (Zevin et al. 2020a).

GW190521 is a GW signal with a BBH source with high component masses, ${85}_{-14}^{+21}{M}_{\odot }$ and ${66}_{-18}^{+17}{M}_{\odot }$ (Abbott et al. 2020c). Accounting for pulsational pair instability (PPI; Yoshida et al. 2016; Woosley 2017; Marchant et al. 2019) and a pair instability supernova (PISN; Fowler & Hoyle 1964; Rakavy & Shaviv 1967; Barkat et al. 1967) uncertainties, the primary mass falls in the mass gap predicted by PISN at around [40−65, 120] M_⊙ (Heger et al. 2003; Spera & Mapelli 2017; Giacobbo et al. 2018; Takahashi 2018; Woosley 2019; Marchant et al. 2019; Farmer et al. 2019; Marchant & Moriya 2020). Our work adopt fits to PPI and PISN models of Marchant et al. (2019) which for metallicity 0.1 Z_⊙ find that the maximum BH mass is ∼44 M_⊙. Hence, this system is a poor fit to our CE and SMT models. This conclusion is consistent with other studies; for example, van Son et al. (2020) investigated isolated binary evolution with super-Eddington accretion without finding any merging BBH systems with a total mass exceeding 100 M_⊙. Alternatively, Fishbach & Holz (2020) showed that if the event is analysed with a prior on the less massive BH of m_BH2 < 48 M_⊙ at 90% credibility, then, the primary BH has a 39% probability of being above the PISN gap. In our models we did not explore stellar and binary evolution above the PISN gap and hence, cannot rule out the formation through CE or SMT.

Amongst the O2 detections there is one event marginally consistent with our models: GW170729 (Abbott et al. 2019). This event has high chirp mass of $M_{\mathrm{chirp}}=35.4_{-4.8}^{+6.5}{M}_{\odot }$ and a high effective inspiral spin of ${\chi }_{\mathrm{eff}}=0.{37}_{-0.25}^{+0.21}$ (Abbott et al. 2019). Figure 5 shows that only models with low α_CE are consistent with this event at the 99% level. When we calculate the likelihood of this system originating from our CE and SMT models using Eq. (F.2), we find it to be extremely small compared to the other events.

Similarly, we can identify a groups of events in GWTC-2 with primary BH masses which support masses larger than our PPI maximum BH mass ∼44 M_⊙ (Marchant et al. 2019). These are GW190413_134308, GW190519_153544, GW190521, GW190602_175927, GW190620_030421, GW190701_203306, GW190706_222641 and GW190929_012149. These events, given our CE and SMT models, have extremely small likelihoods with respect to others, and hence cannot be readily explained by our models.

All the events discussed in this section, perhaps with the exception of GW190412 and GW190814, are likely not the outcome of isolated binary evolution though CE or SMT, given our models. Without considering explicitly models for alternative channels, which is outside the scope of our study, we cannot make a conclusive statement on their origins. However, if we assume that all these systems originated from a different formation channel we can put a lower bound on the contamination fraction in the detected BBH population from other channels to be 9/47 ≃ 0.2.

3.4. Models comparison

The real statistical power of model comparison lies in the combined information from all detected events. In Appendix F we explain how to compute the likelihood of observing N independent GW events, described by the physical parameters θ = (χ_eff, M_chirp, q), given an astrophysical model described by the set of parameters λ. How well the data are described by one model compared to another is described by the Bayes factor (BF), see Eq. (F.4). To compute the BFs, we use as our reference the model of CE and SMT with α_CE = 1.0 and Eddington-limited accretion onto the BH. For our model comparisons we remove the GW events discussed in the previous section (GW170729, GW190413_134308, GW190519_153544, GW190521, GW190602_175927, GW190620_030421, GW190701_203306, GW190706_222641, GW190929_012149) and consider them to not have been formed through the CE and SMT channels in this work. Of course, in the remaining population of events we cannot exclude contamination from other formation channels. A proper analysis would require a model comparison that includes all promising formation channels for BBHs and their branching fractions as model hyperparameters (Zevin et al. 2020c), this is beyond the scope of this paper as here we are only considering two formation channels.

To estimate which model describes best the events, and how sensitive this result is to the selection of events, we iterate the computation of the BF and remove at each iteration the event with lowest likelihood (with respect to the reference model) until the BF converges to a given value. In Fig. 6 we show the BFs of our reference model to the others as a function of sample size; this converge to a constant value after 5–6 events are removed. The BFs indicate moderate to strong evidence in favour of models with inefficient CE, namely α_CE < 1.0, while excluding the model with lowest α_CE = 0.2. If another model is chosen as the reference, then the order of events removed changes, but the end results is the same.

Fig. 6. BFs of CE and SMT models with respect to the reference model with αCE = 1.0 and Eddington-limited accretion efficiency as a function of sample size. The initial sample contains 38 GW events of GWTC-2 and exclude GW170729, GW190413_134308, GW190519_153544, GW190521, GW190602_175927, GW190620_030421, GW190701_203306, GW190706_222641, GW190929_012149. At each iteration the event with lowest likelihood with respect to the reference model is removed and indicated on the horizontal axis until 20 events are removed. By definition the BF of the model with αCE = 1.0 is equal to 1. The data show moderate to strong support for the models with low CE efficiency, αCE < 1.0. This result is robust because the BFs show a constant behaviour as a function of sample size.

Another question we could ask ourselves is which parameter out of the three (χ_eff, M_chirp, q) has the most discriminatory power in the BF analysis and, hence, carries most of the information about CE efficiency. To answer this question we repeat the analysis considering each parameter separately in θ. We find that the parameter carrying the least information is the mass ratio q which strongly disfavour only the model with α_CE = 0.2. Most of the information is contained in M_chirp; considering there is the greatest variation in BFs, for instance, for the model with α_CE = 5.0 the BF is initially disfavoured at ∼10⁻¹² and BF at ∼10⁻⁴ for the one with α_CE = 2.0. The only other model disfavoured by the M_chirp dimension is α_CE = 0.2 at BF ∼ 10⁻²⁰. The M_chirp parameter has the largest discriminating power because is the one affected the most by the variation of α_CE in the underlying distributions (cf. Tables 1 and 2). Larger α_CE leads to a smaller fraction of systems with M_chirp > 15 M_⊙ which is what most of the LVC data support. In contrast, the χ_eff dimension favours models with efficient CE ejection the most: models with α_CE = 2.0 at BF ∼ 10⁷–10⁴ and α_CE = 5.0 at BF ∼ 10^2.6–10⁰ while it strongly disfavours α_CE ∈ [0.2, 0.35]. This result is a direct consequence of the majority of LVC data supporting small χ_eff which in our models are achieved for ultra efficient α_CE. Finally, when we calculate the BFs with θ = (χ_eff, M_chirp) we find similar results to the original analysis with θ = (χ_eff, M_chirp, q). Moreover, we find that, in contrast to the three-dimensional BF analysis, the two-dimensional one shows rather constant BFs starting from the initial sample. The variation in the first 6 iterations of the three-dimensional analysis is caused by the q dimension. We should stress, however, that the model with α_CE = 0.2 significantly overpredicts the rate of events. Hence, overall an α_CE in the range 0.2 < α_CE < 1.0 is favoured.

The result of our model section analysis needs to be interpret with caution. Here we only explored one parameter of the models which is degenerate with others, for example MT stability and efficiency, BHs birth spin, etc. Moreover we found a non-negligible fraction (≥0.2) of BBHs originating from other formation channels and, hence, cannot rule out an even greater contamination in the studied population. Finally, other formation channels have been shown to predict BBH observable distributions degenerate with the CE and SMT channels (see references in Sect. 1). Hence, a proper analysis would require a model comparison of all promising formation scenarios (cf. Zevin et al. 2020c).

4. Model uncertainties

4.1. Initial binary properties

Lack of knowledge of the initial binary properties are a source of uncertainty in population synthesis studies. We assume that the primary mass follows the Kroupa (2001) initial-mass function (IMF). This is a power law which comes with uncertainties that can affect the rate estimate and to a lesser degree the observed distribution of BBHs (e.g. de Mink & Belczynski 2015). Moreover, the IMF may not be universal (e.g. Schneider et al. 2018; Farr & Mandel 2018). An additional source of uncertainty can be the initial binary mass fraction distribution and birth eccentricities, which we did not investigate here. We expect the impact of these uncertainties to be smaller compared to that of the IMF (de Mink & Belczynski 2015).

Another important assumption is the distribution of the birth orbital separation of the binaries. In our model, we assumed it follows an extended Sana et al. (2012) power law in log-orbital period p ∈ [0.4, 10^5.5] days, cf. Eq. (A.1). Here, we investigate the sensitivity of our results to this assumption. In Fig. 7 we show the histogram of 10⁵ initial orbital separations drawn from the assumed distribution, in blue, compared to a log-uniform distribution, in orange. We find that the rate density of the reference CE model with α_CE = 1.0 raises from 42.6 Gpc⁻³ yr⁻¹ to 73.3 Gpc⁻³ yr⁻¹ when assuming log-uniform orbital separation at birth. On the other hand, the rate density of SMT channel remains almost the same. This happens because the log-uniform distribution increases the yield of merging BBHs at large initial orbital periods (p ≳ 10² days), which are the binaries going through the CE evolutionary path, while it does not affect the yield of SMT BBHs as both initial orbital period distributions create roughly the same number of systems at p ∼ 10 days.

Fig. 7. Initial orbital separation of 105 binaries sampled from the extended Sana et al. (2012) distribution, cf. Eq. (A.1), in blue, and from a log-uniform distribution, in orange. Both distributions sample the range p ∈ [0.4, 105.5] days and are independent of metallicity. With a solid line, we show the same distributions after removing systems with Roche-lobe overflow at ZAMS for Z = 0.1 Z⊙.

Moreover, the sampled orbital period range can affect our results. In our model we lowered the smallest orbital period, compared to that measured by Sana et al. (2012), to include the portion of the parameter space leading to the chemical homogeneous formation of BBH (Mandel & de Mink 2016; Marchant et al. 2016; du Buisson et al. 2020). By default, we included systems overfilling their Roche lobe at ZAMS. Marchant et al. (2016) found that binaries that are already in contact at ZAMS can potentially survive and lead to the formation of BBHs through chemically homogeneous evolution. However, we only have observations of massive binaries when they are well-past their ZAMS, as in prior evolutionary phases they are still embedded in the formation clouds, undergoing accretion. Accretion onto a pre-MS binary significantly complicates its evolution (e.g. Sørensen et al. 2018), and thus including those binaries in our population models may be problematic. We now investigate how our results change if we exclude these systems. To exclude systems that are initially Roche lobe filling, we adopt the stellar radii fits of Tout et al. (1996). These fits are specific for ZAMS and offer more accuracy than the one of Hurley et al. (2000) which are meant to cover the entire stellar evolution and thus sacrifice some accuracy at ZAMS. By removing these binaries, we decrease the number of systems in small orbital periods as shown by the solid lines for both distributions in Fig. 7. The new way of drawing initial orbital periods is metallicity dependent because, in general, stars have larger radii at larger metallicities (Tout et al. 1996). In Table 3 we summarise the rates of all these models. As expected, we find for both distributions that excluding RLOF ZAMS increases the rates for both channels. We conclude that this uncertainty can affect our results by a factor of ≲2.

4.2. Mass-transfer stability

The critical mass ratio q_crit determines weather the MT is dynamically stable or unstable (cf. Appendix B.2). We chose our q_crit values to match the assumptions of our previous work (Bavera et al. 2020), which is based on the population synthesis model of Neijssel et al. (2019) obtained using the COMPAS code (Stevenson et al. 2017; Vigna-Gómez et al. 2018). In contrast to Neijssel et al. (2019), we implement the same q_crit fits to the GB and asymptotic AGB (Hjellming & Webbink 1987) but do not adopt Soberman et al. (1997) radial response to adiabatic mass loss for evolved stars beyond the HG (this option is not present in the current version of COSMIC). Despite this and other differences in the model assumptions (such as different CE λ fits) which might affect the results, we reached similar conclusions for the detected population of the CE channel with α_CE ∼ 1.0. Even though both models converge on similar detected population distributions, the two models have different mass-ratio distributions for the underlying BBH population. We suspect that this discrepancy is caused by the different q_crit assumption for GB and AGB stars as all our merging BBH systems are evolving through central He burning at onset of CE. In order to verify the source of the difference, a thorough code comparison is needed, which is outside the scope of this project.

To test the sensitivity of our results to this assumption, we run two additional models with different q_crit choices: model (C) with q_crit values corresponding to Claeys et al. (2014) and model (B) with q_crit values similar to Belczynski et al. (2008). We find model (C) to have similar rate density to our fiducial CE channel but almost five times larger rate density of the SMT channel. On the other hand, the model (B) CE rate density is half of our fiducial model and slightly larger rate densities for the SMT channel. These results are summarised in Table 4 where we also report for comparison the median mass ratio of the intrinsic and detected BBH populations. Both q_crit choices do not have a significant impact on the detected observable distributions even though they have different impact on the underlying populations. To determine q_crit, COSMIC uses the evolutionary type of the donor (as defined in Hurley et al. 2000) rather than the actual structure of its envelope. Recently, Klencki et al. (2021) showed that this tends to overpredict the number of systems that evolve through and survive a CE phase. This shows the limitation of parametric population synthesis codes. In fact, q_crit can be numerically determined by detailed binary evolution simulations, given the profile of the donor star. In the future, next-generation population synthesis tools based on detailed binary and stellar evolution simulations will remove this degree of freedom.

Finally, in our models we only explore the effects of MT accretion efficiency onto BHs, and did not investigate the effects of MT accretion efficiency between two non-degenerate objects. A recent study (Bouffanais et al. 2020) showed that, assuming Eddington-limited accretion onto BHs, stars need to accrete more than 30% of the mass lost by the donor stars during the first MT episode in order to explain O1 and O2 BBHs. MT efficiency between two non-degenerate stars depends strongly on the assumed specific AM of the material that reaches the surface of the mass-gaining star. In turns, this depends on the accretion disk physics and the coupling of the accretion disk to the star’s surface. The assumption that the accreted material carries the Keplerian specific AM of the accretor’s surface leads to a very efficient spin up of the mass-gaining star. The accretor then quickly reaches critical rotation and halts further accretion, leading to a highly inefficient mass transfer (a few to a few tens percent; e.g. de Mink et al. 2013; Langer et al. 2020). On the other hand, if one assumes that the AM is dissipated efficiently before it reached the accretor’s surface, and that that the material that is added onto the star has similar specific AM to that of its own surface layers, then MT can be significantly more efficient.

4.3. Effect of angular momentum transport and accretion efficiency on black hole spin

Our models assume efficient AM transport (Spruit 1999, 2002) which leads to the formation of non-spinning first-born BHs (Fragos & McClintock 2015; Qin et al. 2018). Although the Tayler–Spruit dynamo model helps to reproduce the flat rotation profile of our Sun (Fuller et al. 2014; Cantiello et al. 2014), as well as, neutron star and white dwarf spins (Heger et al. 2005; Suijs et al. 2008), it cannot reproduce the asteroseismic constrains for subgiants and red giants (Gehan et al. 2018). This would require an even higher efficiency in AM transport. Alternatively, a model with inefficient AM transport predicts highly spinning BHs (e.g. Arca Sedda & Benacquista 2019; Belczynski et al. 2020), which do not match current GW observations. If the second MT episode is stable, then the first-born BH can be spun up by accretion (Thorne 1974). If the MT accretion onto BHs is Eddignton-limited, the BH accretes a negligible amount of matter leading to a₁ ≃ 0. On the other hand, a super Eddington-limited MT accretion will result in larger spins. The extreme case of this would be highly super-Eddington accretion efficiency where the spin of the first-born BH can even approach to a₁ ≃ 1, but in this case the contribution to the merging BBH population of the SMT channel almost vanishes (Table 1).

The spin of the second-born BH is determined by the combined effects of stellar winds and tidal interactions during the BH–He-star binary evolution which we model in detail. During this evolutionary phase, the AM transport does not play an important role as the He-star is compact and will not expand during its final evolution (Bavera et al. 2020). The strength of the tidal interaction is primarily determined by the orbital separation during the BH–He-star evolutionary phase, see Eq. (C.1). In our SMT models, the orbital separation is determined by the accretion efficiency. Models where super-Eddington accretion is allowed will result in larger orbits than our reference model, decreasing further the small effect of tides on this evolutionary path.

4.4. Common-envelope prescription

In our CE models, the post-CE orbital separation, A_postCE, is linearly dependent on the CE parameterisation uncertainties as, approximately, δA_postCE/A_preCE ∝ δα_CE λ. Even though we did not explore different λ fits in our models, our parameter study of α_CE ∈ [0.2, 5.0] covers an uncertainty on δA_postCE/A_preCE ∝ 5.0/0.2 ≃ 25. Recently, Klencki et al. (2021) showed that λ fits similar to the one we used could underestimate the envelope binding energies of massive radiative-envelope giants, leading to an overestimation of the systems surviving CE. An additional free parameter in the calculation of λ which complicates a detailed treatment of CE is the assumed boundary down to which the envelope will be ejected. Unfortunately detailed stellar models cannot robustly predict this (e.g. Tauris & Dewi 2001) and hydrodynamic simulations of the CE phase are necessary. For example, Fragos et al. (2019) showed that for progenitors of binary neutron stars, a non-negligible fraction of hydrogen rich material will remain bound to the core after the successful ejection of the CE, that would in turn imply a relatively efficient ejection of the CE.

4.5. Other uncertainties

Our model may be limited by other uncertainties we did not explore which can alter the merger rate and, to a lesser degree, the predicted BBH property distributions. Uncertainties in the (i) physics of the supernova explosions, such as the kicks strength, can influence rates and affect the parameter distribution of BBH mergers (e.g. Dominik et al. 2013). When connecting the population synthesis code to our detailed MESA simulations we (ii) assumed the BH–He-star systems post second MT to be at ZAHeMS. This is not always the case as the binaries are evolving through central He burning at onset of the MT. This leads us to overestimate the lifetime of these He-stars. This overestimation is negligible as the binary only spends a few hundred thousand years in this state compared to its overall life of a few million years and much longer BBH inspiral. This overestimation might lead to less massive second-born BHs and smaller spins as winds act for a larger time window. However, we expect that the fraction of stars entering the MT on advanced He-burning phase is higher at low metallicities, as low-metallicity stars tend to expand later in their lives (Linden et al. 2010). At the same time, stellar winds in these stars are weaker due to the low metallicity, so the overall effect on the population of BBHs is expected to be limited. The uncertainty of the (iii) metallicity dependence of stellar winds for massive stars is another important ingredient of population synthesis studies which affects M_chirp distributions and the rates (cf. Barrett et al. 2018). The detection rate and density calculation is also affected by uncertainties in the (iv) redshift and metallicity dependent SFR (Chruslinska et al. 2019; Neijssel et al. 2019; Tang et al. 2020). A SFR favouring higher formation metallicities than the one assumed here would skew our results in favour of smaller χ_eff and lower rates as low metallicity systems are responsible for high χ_eff and the short merger timescales.

5. Conclusions

Mass-transfer physics is one of the most uncertain physical processes determining the observable properties of field binary black holes. The critical mass ratio q_crit determines the fraction of the parameter space going trough SMT and CE phases. Mass-accretion efficiency onto compact objects determines how efficiently binaries going through SMT will shrink, while CE efficiency α_CE determines post-CE orbital separations. In this work we investigated how the detected joint distributions of M_chirp, χ_eff and q of BBH formed through the CEe and SMT formation channels change given the uncertainties on these input physics. We investigated this by combining rapid binary population synthesis studies with detailed stellar and binary simulations. Rapid population synthesis studies allow us to obtain different BH–He-star populations for different input physics, while detailed simulations which take into account differential stellar rotation, tidal interaction, stellar winds and the evolution of the He-star stellar structure allow us to accurately determine the distributions of BBH observables. We then convolved the synthetic BBH population with the redshift- and metallicity-dependent star-formation rate, as well as selection effects from a 3-detector network to build a model capable of describing LIGO–Virgo detections. Our main findings are:

• We calculated the O3 detected joint distributions of χ_eff, M_chirp and q for CE and SMT channels. Assuming efficient AM transfer inside stars and Eddington-limited accretion efficiency, both channels lead to similar rate densities in the local Universe. We find that the CE channel is the only evolutionary path leading to non-zero χ_eff in the detected population as SMT channel cannot shrink the orbits enough for efficient tidal spin-up to take place.

• Inefficient CE (small α_CE values) leads to smaller orbital separation post CE. On average, these models lead to more systems being tidally spun up. However, the majority of these systems are not detected by current GW detectors because most of these systems are formed in low metallicity environments (otherwise stellar winds widen the binaries) and merge quickly at a redshift close to their formation (z ∼ 2 where the SFR peaks), far outside current GW detector horizons (z ≲ 1).

• Highly super-Eddington accretion efficiency onto compact objects reduces the rate densities of CE by ∼10%, while it reduces the contribution of SMT channel by two orders of magnitude, making the contribution of this channel almost negligible compared to the CE channel.

• The GW events GW170729, GW190413_134308, GW190519_153544, GW190521, GW190602_175927, GW190620_030421, GW190701_203306, GW190706_222641 and GW190929_012149 of GWTC-2 are not well explained by our models of isolated binary evolution through either CE or SMT. If we assume that these systems originate from other formation channels then we can put a lower bound on the detected branching fraction from other formation channels: 9/47 ≃ 0.2.

• We conducted a model comparison given the events of GWTC-2 consistent with our CE and SMT models to determine which CE efficiency is best supported by the data. The GW events show moderate to strong evidence in favour of the models with inefficient CE, 0.2 < α_CE < 1.0. We find this result to be robust considering different selections of events in our calculations. This analysis did not include rate estimates which the α_CE = 0.2 model significantly overpredicts.

We conclude that future works aiming to properly infer model parameters through model comparison will need to consider correlation between parameters as well as contamination from other formation channels in order to properly determine model parameters.

Acknowledgments

We would like to thank Simon Stevenson and Katerina Chatziioannou for useful discussions. This work was supported by the Swiss National Science Foundation Professorship grant (project number PP00P2 176868). This project has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie RISE action, grant agreement No 691164 (ASTROSTAT). MZ is supported by NASA through the NASA Hubble Fellowship grant HST-HF2-51474.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. CPLB is supported by the CIERA Board of Visitors Research Professorship. PM is supported by the FWO junior postdoctoral fellowship No. 12ZY520N. JJA and SC are supported by CIERA and AD, JGSP, and KAR are supported by the Gordon and Betty Moore Foundation through grant GBMF8477. KK received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement n. 617001. YQ acknowledges funding from the Swiss National Science Foundation under grant P2GEP2_188242. The computations were performed in part at the University of Geneva on the Baobab and Lesta computer clusters and at Northwestern University on the Trident computer cluster (the latter funded by the GBMF8477 grant). All figures were made with the free Python modules Matplotlib (Hunter 2007) and pygtc (Bocquet & Carter 2016). This research made use of Astropy (www.astropy.org), a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018).

References

Appendix A: Initial binary distributions

The parameters describing the initial conditions of a binary system are primary and secondary masses, m₁ and m₂, orbital period p, eccentricity e and metallicity Z of the stars at ZAMS. We assumed that the primary masses follow the initial mass function (IMF) of Kroupa (2001) which spans the mass range 0.01 M_⊙ ≤ m₁ ≤ 150 M_⊙. The upper limit is an extrapolation of the Kroupa (2001) IMF which is measured only up to 50 M_⊙. In our model, the lower limit represents the smallest theoretical mass for a star to support nuclear fusion (cf. Kumar 1963) while the arbitrary maximum stellar mass exclude BH formation above the upper mass-gap (e.g. Heger et al. 2002). The mass distribution of the less massive secondary star is given by m₂ = m₁q where the initial mass ratio q is drawn from a flat distribution (Sana et al. 2012) in the range q ∈ [0, 1]. Furthermore, we adopt a binary fraction of f_bin = 0.7 (Sana et al. 2012) and assume that at birth the distribution of log-orbital periods follow Sana et al. (2012) power law with coefficient π = −0.55 in the range p ∈ [10^0.15, 10^5.5] days and extrapolate down to the range p ∈ [0.4, 10^0.15] days assuming a log-flat distribution (as the power low is not defined for p < 10^0.15 days), that is

$$\begin{array}{c}\hfill f(p)=\{\begin{array}{cc}C\times 0.{15}^{-0.55}\hfill & 0.4\le p/\mathrm{days}<{10}^{0.15}\hfill \\ C{log}_{10}{(p/\mathrm{days})}^{-0.55}\hfill & {10}^{0.15}\le p/\mathrm{days}\le {10}^{5.5}\hfill \end{array},\end{array}$$(A.1)

where the normalisation constant C is determined by the condition ${\int }_{p_{\text{min}}}^{p_{\text{max}}}f(p)\text{d}p=1$. The lower limit of 0.4 days will ensure that we probe all the available parameter space, see for example, the detailed simulations of du Buisson et al. (2020) where they find binaries with initial orbital periods as small as 0.4 days forming BBHs. Finally, we assume that all binaries are born with circular orbits, namely with zero e. We assume that all these distributions are both independent of each other, as well as, independent of metallicity which in reality might not be the case (Moe & Di Stefano 2017). This oversimplification means that we also include systems overfilling the L1 and L2 Roche-lobe surfaces at ZAMS. In Sect. 4.1 we discuss how our BBH rate estimates are affected by omitting these systems or by assuming a log-uniform distribution. The omission of this portion of the parameter space does not lead to a qualitative difference in the predicted observable distributions.

In order to optimise the absolute number of binaries becoming BBHs per number of simulated binaries, we restrict the primary mass to the range 5 M_⊙ ≤ m₁ ≤ 150 M_⊙. For the metallicity, we divide uniformly the log-metallicity range Z ∈ [0.0001, 0.0309] in 30 bins. The largest metallicity bin was chosen to have a centre at 1.5 Z_⊙ where we adopt the solar reference Z_⊙ = 0.017 (Grevesse et al. 1996). This metallicity range is where the stellar model fits of Hurley et al. (2000) we use are defined. We evolve 5 million binaries per metallicity bin ΔZ for a total simulated mass per ΔZ of M_{sim,  ΔZ} ≃ 1.05 × 10⁸ M_⊙. Since we restricted the primary mass, we only model a fraction of the underlying stellar population. Hence, we need to re-normalise the simulated stellar mass M_{sim,  ΔZ} to obtain the total stellar population; the normalisation constant for our choice of initial binary properties is $f_{\text{corr}}^{-1}=4.72$, see Appendix A in Bavera et al. (2020) for a derivation of this quantity.

Appendix B: COSMIC population synthesis model

B.1. Single stellar models

To generate our BBH population models, we used COSMIC version v3.3.0. Stellar evolution in COSMIC (Breivik et al. 2020) is based on the analytical fits of Hurley et al. (2000, 2002) to the single stellar models of Pols et al. (1998). For O and B stars we adopt mass loss through stellar winds according to the prescription of Vink et al. (2001), which covers separately the temperature ranges 12 500 K < T_eff < 22 500 K and 27 500 K < T_eff < 50 000 K. Around T ≃ 25 000 K there is a bi-stability jump thatf leads to a mass-loss increase of a factor of about five. This jump is due to the recombination of the Fe IV to the Fe III ion in the lower part of the wind (Vink et al. 1999). For Wolf–Rayet stars mass-loss rate we adopt stellar winds as in Brott et al. (2011) who assume Hamann et al. (1995) reduced by a factor of 10 to correct for clumping and use a metallicity scaling of (Z/Z_⊙)^0.85 (Vink et al. 2001). For all these wind prescriptions we adopt the solar reference Z_⊙ = 0.017 (Grevesse et al. 1996).

During the post-carbon burning phase of massive stars, photons produced in the core are energetic enough to produce electron–positron pairs which softens the equation of state, diminishing the pressure support of the core (Woosley et al. 2007). This causes the core to rapidly contract and the temperature to increase, allowing for explosive oxygen burning (e.g. Woosley & Heger 2015). For He-core masses in the range ∼[30, 64] M_⊙ the released energy is insufficient to completely disrupt the star. This create a series of energetic pulses which eject material from the star before it collapses into a BH. This is the PPI (Yoshida et al. 2016; Woosley 2017; Marchant et al. 2019). If the He-core mass is in the range ∼[64, 133] M_⊙ the released energy is enough to reverse the collapse, unbinding and destroying the star. This event is a PISN (Fowler & Hoyle 1964; Rakavy & Shaviv 1967; Barkat et al. 1967). Similar to Stevenson et al. (2019), we adopt the fit to the grid of simulations from Marchant et al. (2019, see Table 1), which demonstrate a turnover in the relation between pre-supernova He-core mass and final mass. We use the 9th-order polynomial fit of Breivik et al. (2020), cf. their Eq. (4), to map the pre SN stellar mass in the range 31.99 ≤ M_preSN/M_⊙ ≤ 61.10 to the baryonic mass collapsing to form the BH.

B.2. Mass-transfer stability and common envelope

The stability of Roche-lobe overflow mass transfer is determined by the rate at which the Roche-lobe radius is changing as a result of mass-transfer dlog(R_L)/dlog(m) to the response of the radius of a star as its mass is changing dlog(R_*)/dlog(m). We use the approximation of Eggleton (1983) for the Roche-lobe radius while we approximate the radial response of the star depending on its stellar type. We adopt the values assumed in Neijssel et al. (2019), Bavera et al. (2020). The stability of the mass transfer can then be determined by solving this equation with respect to the critical mass ratio, defined as q_crit = m_donor/m_accretor. For MS stars we use dlog(R_*)/dlog(m) = 2.0 which correspond to q_crit ≃ 1.72 while for HG stars dlog(R_*)/dlog(m) = 6.5 which correspond to q_crit ≃ 3.83 (Ge et al. 2015). For stars on the GB and AGB we use fits from Hjellming & Webbink (1987). For stripped stars we adopt q_crit as in Claeys et al. (2014). Different choices of q_crit, especially for GB and AGB stars have an impact on the parameter space that leads to the formation of BBHs, hence on the merger rate, see Sect. 4.3 for a discussion of this uncertainty.

If the mass transfer is stable the companion star accretes a fraction of mass lost by the donor star. Any mass that is not accreted leaves the system instantaneously, taking away the specific AM of the accretor (Hurley et al. 2002). For our fiducial models we assume that the accretion of degenerate objects is Eddington-limited, this results in a highly non-conservative mass-transfer phase where degenerate objects accrete negligible amount of mass, hence the first-born BH will not spin up because of accretion during the second mass-transfer phase (Thorne 1974). For the Eddington-limited accretion efficiency onto a compact object, COSMIC uses the definition of Hurley et al. (2002),

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{edd}}=2.08\times {10}^{-3}{(1+X)}^{-1}(R_{\mathrm{acc}}/R_{\odot })M_{\odot }{\mathrm{yr}}^{-1}\equiv 1={\eta }_{\mathrm{acc}},\end{array}$$(B.1)

where X is the hydrogen mass fraction of the donor and R_acc is the accretion radius of the compact object which for a BH is chosen to be the Schwarzschild radius. In this work we define η_acc to be the MT accretion efficiency limit in units of the Eddington-limit. To investigate how our results depend on this limit, we explore different MT efficiency limits: η_acc ∈ [1, 10³, 10⁵, 10⁹].

If the mass transfer is unstable, the donor star will expand to form a CE of gas around the binary which can be expelled by the injection of orbital energy from the binary (Paczynski 1976). This is a complex phase and we parameterise it with the α_CE–λ formalism (see, e.g. Ivanova et al. 2013, for a review). In this parameterisation α_CE determines the efficiency factor for injecting orbital energy into the envelope, while λ characterises the binding energy of the envelope to its stellar core which depends on the structure of the donor’s envelope. Assuming that the pre-CE orbital energy is much smaller than the post-CE orbital energy, the initial and final orbital period of the CE event, A_preCE and A_postCE, are related by the following expression (Ivanova et al. 2013)

$$\begin{array}{c}\hfill A_{\mathrm{postCE}}\simeq \frac{1}{2}{\alpha }_{\mathrm{CE}}\lambda \frac{m_1{m}_{2,\mathrm{postCE}}}{m_{2,\mathrm{preCE}}{m}_{2,\mathrm{env}}}{\widehat{r}}_L{A}_{\mathrm{preCE}},\end{array}$$(B.2)

where ${\widehat{r}}_L$ is the dimensionless Roche-lobe radius, m₁ is the mass of the accretor, m_2, preCE and m_2, postCE are the donor star masses before and after the CE event and m_2, env the envelope mass. Previous studies of post-CE binaries have shown that the efficiency could be as low as α_CE = 0.2 (Zorotovic et al. 2010; Toonen & Nelemans 2013; Camacho et al. 2014) while other using detailed modelling of the CE phase for binary neutron star progenitors (Fragos et al. 2019), suggests CE efficiencies as high as α_CE = 5.0 (see also Giacobbo & Mapelli 2019). Approximately, uncertainties on α_CE or on λ linearly scale to uncertainties on the orbital separation post CE as δA_postCE/A_preCE ∝ δα_CE × λ + α_CE × δλ. For λ we adopt the fits as in Claeys et al. (2014).

Within the CE channel we distinguish and adopt a pessimistic scenario, in which unstable mass transfer from a donor star without a well-developed core-envelope structure, namely when a star finds itself in the HG, is always assumed to lead to a merger (Belczynski et al. 2007). An optimistic scenario which include these systems would result in an overestimation of the observed local BBHs merger rate density (Neijssel et al. 2019; Belczynski et al. 2020).

Appendix C: MESA detailed BH–He-star models

In our previous works (Qin et al. 2018; Bavera et al. 2020) we explored the evolution of tight BH–He-star systems computing grids of detailed binary evolution models using MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019). This is the late-end phase of the binary evolution of BBHs formed through the CE and SMT channels. Here we iterate on this work by adapting those models to the newer MESA version r11701, computing an even larger grid and modifying our stellar models to match the physical stellar assumptions made in the work of du Buisson et al. (2020), which expand on the work of Marchant et al. (2016) on chemical homogeneous evolution.

C.1. Single stellar physics

All mixing processes are treated as diffusive processes. Convection is modelled using the Ledoux criterion, adopting a mixing-length parameter of α_MLT = 1.5 (Böhm-Vitense 1958). Semiconvection is modelled according to Langer et al. (1983) with an efficiency parameter α_SEM = 1.0 (Langer 1991). We also take into account exponential core-overshooting and thermohaline mixing with α_TH = 1.0 (cf. Cantiello et al. 2009). Opacities are computed using CO-enhanced opacity tables from the OPAL project (Iglesias & Rogers 1996), computed using solar-scaled abundances based on Grevesse et al. (1996). As we are not interested in following the nucleosynthesis in detail, we use the simple networks provided with MESA basic.net which MESA automatically extend to co_burn.net to resolve carbon burning. Our models are stopped at C depletion. Finally, PPI and PISN are treated as in the single stellar models of COSMIC. We assume that the mass loss through PPI of any stellar model is lost from the surface of the star and only consider the AM of the remaining layers in the collapse of the star.

Rotational mixing and AM transport are also treated as diffusion processes (Heger et al. 2000, 2005), which involve the effects of Eddington–Sweet circulation, secular and dynamical shear mixing and the Goldreich–Schubert–Fricke instability with an efficiency parameter of f_c = 1/30 (Chaboyer & Zahn 1992; Heger et al. 2000). We include the effect of magnetic fields on the transport of AM assuming an efficient AM transport mechanism: the Tayler–Spruit dynamo (Spruit 1999, 2002). An efficient AM transport allows us to assume that all He-stars emerging from the CE or SMT phases are initially not rotating. This is because any initial or acquired rotation during the evolution of the secondary is erased by mass transfer and wind mass loss by the time it becomes a He-star. Assuming instead that the He-star is initially synchronized with the orbit, right after the formation of the BH–He-star binary, has been shown to have negligible effects in the final properties of the resulting BBH system (Qin et al. 2018).

We implement the same WR stars wind prescription as in COSMIC where we also include the enhancement of winds through rotation as in Heger & Langer (2000). When the rotation rate exceeds a given threshold, Ω/Ω_crit > 0.98, we implicitly compute the mass-loss rate required for the rotation rate to remain just below this value (Paxton et al. 2015). For a star with mass M and radius R, Ω_crit ≡ [(1−Γ)GM/R]^1/2 where Γ = κL/(4πcGM) is the Eddington factor and κ is the true flux-mean opacity coefficient (Langer et al. 1997).

C.2. Binary stellar physics

We use MESA single star module to first create a He-star with the desired mass and with abundances Y = 0.98 and Z = 0.02. We then load the model in the MESA binary module, relax the metallicity to the desired value and allow the star to relax until it reaches ZAHeMS. We define ZAHeMS to be the moment when the central luminosity becomes larger than ∼99% of the surface luminosity. To facilitate the relaxation of the star to ZAHeMS we adopt MLT++ treatment of convection (Paxton et al. 2013) which reduces the superadiabaticity in some radiation-dominated convective regions. Once the star reaches ZAHeMS we turn off MLT++. At ZAHeMS we check if the He-star overfill its Roche lobe, if this is the case we stop the evolution and assume the system to be non-physical. Whenever one component in the system attempts to overflow its Roche lobe during the binary evolution, we implicitly compute the mass-transfer rate using the Kolb scheme (Kolb & Ritter 1990). If the mass transfer exceed 10 M_⊙ yr⁻¹ we stop the run and assume that the binary will merge. Furthermore, we use the prescription of Misra et al. (2020) to check that the He-star does not overflow the L2 Roche volume of the binary. If this is the case, either at ZAHeMS or during the evolution of the binary, we stop the run and assume that it will lead to a merger.

Tidal forces are responsible for synchronising the spin of the He-star with the orbit (Zahn 1977; Hut 1981). Tidal effects are implemented as in Qin et al. (2019) for the case of stars with radiative envelope. This is a variation of the standard tidal prescription of MESA (Paxton et al. 2015) which synchronize the whole star. Here tides only operate on the radiative regions. This slight variation has a negligible impact on our results because the Tayler–Spruit dynamo guarantees strong coupling between the star’s layers. We assume that the orbits after the second MT phase are circular and the system remains circular during binary evolution. The strength of the interaction depends on the ratio of the stellar radius R to the orbital separation A, the He-star mass M, the binary mass ratio q and the moment of inertia I. The timescale for synchronization is given by (Hurley et al. 2002)

$$\begin{array}{c}\hfill \frac{1}{T_{\mathrm{syn}}}=3E_2{(1+q)}^{5/6}\frac{q^2}{r_{\mathrm{g}}^2}{\left(\frac{\mathit{GM}}{R^3}\right)}^{1/2}{\left(\frac{R}{A}\right)}^{17/2},\end{array}$$(C.1)

where r_g = [I/(MR²)]^1/2 is the dimensionless gyration radius of the He-star and E₂ is the second order tidal coefficient. As in our previous work, we use the fitting formula of E₂ for He-star from Qin et al. (2018).

Appendix D: Core collapse

BHs are formed during the core collapse of massive stars and, in some cases, their formation is associated with a supernova explosion. As in Bavera et al. (2020), we use Fryer et al. (2012) delayed supernova prescription to model how much baryonic remnant mass is left behind after the core collapse. This model avoids an enforced mass gap between neutron stars and BHs which are not consistent with current microlensing observations (Wyrzykowski & Mandel 2020) or with GW190814 (Abbott et al. 2020b). During the collapse of the star, asymmetric ejection of matter (Janka & Mueller 1994; Burrows 2013) or asymmetric emission of neutrinos (Bisnovatyi-Kogan 1993; Socrates et al. 2005) can provide a momentum kick to the newly formed BHs. We assume that the birth kicks of BHs follow a Maxwellian distribution with σ = 265 km s⁻¹ (Hobbs et al. 2005) rescaled by 1 − f_fb where f_fb is the fallback mass fraction (Fryer et al. 2012). For massive stars with carbon–oxygen core masses grater than 11 M_⊙, Fryer et al. (2012) prescription assumes a direct collapse, that is f_fb = 1. This implies that in our model all heavy BHs receive no natal kicks. The kicks impacted on the lighter BHs can tilt the orbit, which may generate a negative χ_eff, add eccentricity to the orbit or disrupt the binary. In practise systems with negative χ_eff are not statistically relevant in our models as the kicks are not strong enough to flip the orbits by more than 90°. Recently, Callister et al. (2020) showed that a σ ≃ 1000 km s⁻¹ is required in order to explain negative χ_eff in GWTC-2 data with CE BBHs. We take into account all these orbital changes, as well as, orbital changes due to neutrinos symmetric mass loss, following the analytical calculations of Kalogera (1996).

We estimate the spin of the resulting second-born BH following the framework presented in Bavera et al. (2020). Here, we adopt a different treatment of neutrino mass loss motivated by Zevin et al. (2020a) which prescribes that a collapsing proto-neutron star can lose up to 0.5 M_⊙ through neutrino emission. In order to calculate the final mass and spin of the BH resulting from the collapse, we need to consistently follow its accretion history soon after it is formed. We assume that the innermost 3 M_⊙ form a proto-neutron star, which collapse to form a BH of 2.5 M_⊙, while 0.5 M_⊙ are converted into neutrinos and leave the systems carrying away a fraction (0.5/3) of the proto-neutron star AM. We consider a collapsing star to be a collection of shells with mass m_shell and angular frequency Ω_shell that falls one by one onto the central BH. Once a shell reaches the BH’s event horizon, it is accreted by it. The amount of specific AM of the in-falling material, j(r, θ) = Ω_shell(r)r²sin(θ) where θ is the polar angle, determines the properties of the accretion flow. Low AM material collapses directly onto the BH transferring its entire mass and AM to the hole, while material with enough AM can create a disk around it. The mass M_BH and spin a of the accreting BH determine the threshold for disk formation and is given by the specific AM at the innermost stable circular orbit (ISCO; Bardeen et al. 1972)

$$\begin{array}{c}\hfill j_{\mathrm{ISCO}}=\frac{G{M}_{\mathrm{BH}}}{c}\frac{2}{3^{3/2}}[1+2{(\frac{3c^2{r}_{\mathrm{ISCO}}}{G{M}_{\mathrm{BH}}}-2)}^{1/2}],\end{array}$$(D.1)

where r_ISCO is the radius of the ISCO for prograde equatorial orbits,

$$\begin{array}{c}\hfill r_{\mathrm{ISCO}}=\frac{G{M}_{\mathrm{BH}}}{c^2}\{3+z_2-{\left[(3-z_1)(3+z_1+2z_2)\right]}^{1/2}\},\end{array}$$(D.2)

with z₁ = 1 + (1 − a²)^1/3[(1 − a)^1/3 + (1 + a)^1/3] and $z_2={(3a^2+z_1^2)}^{1/2}$. From the disk formation condition j(r, θ) > j_ISCO we can define the polar angle at which disk formation occurs as

$$\begin{array}{c}\hfill {\theta }_{\mathrm{disk}}=arcsin{\left(\frac{j_{\mathrm{ISCO}}}{{\mathrm{\Omega }}_{\mathrm{shell}}(r)r^2}\right)}^{1/2}.\end{array}$$(D.3)

The portion of the shell with θ < θ_disk will collapse directly onto the BH on a dynamical timescale, t_dyn ≃ [r³/(GM(r))]^1/2, transferring j(r, θ) to the hole, while the portion of the shell with θ ≥ θ_disk will form a disk and transfer only j_ISCO to the BH. The disk will be accreted on a viscous timescale $t_{\nu }\simeq {\alpha }_{\nu }^{-1}{(R_{\text{circ}}/H)}^2{t}_{\text{circ}}$ assumed to be much smaller than t_dyn (Batta & Ramirez-Ruiz 2019). Here H is the disk’s scale height, α_ν is the viscosity parameter and t_circ is the Keplerian orbital period at the accretion radius also known as circularisation radius $R_{\text{circ}}=j_{\text{shell}}^2/(G{M}_{\text{BH}})$. The collapsing shell contributes therefore to the AM of the BH by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill J_{\mathrm{shell}}\equiv J_{\mathrm{direct}}+J_{\mathrm{disk}}=& {\displaystyle {\int }_0^{{\theta }_{\mathrm{disk}}}{m}_{\mathrm{shell}}{\mathrm{\Omega }}_{\mathrm{shell}}(r)r^2{sin}^3(\theta )\mathrm{d}\theta +}\hfill \\ \hfill & +{\displaystyle {\int }_{{\theta }_{\mathrm{disk}}}^{\pi /2}{m}_{\mathrm{shell}}{j}_{\mathrm{ISCO}}sin(\theta )\mathrm{d}\theta .}\hfill \end{array}\end{array}$$(D.4)

The accretion disk has mass m_disk = m_shellcos(θ_disk) and the mass-energy accreted onto the BH from the disk is ΔM_disk = m_disk[1−2GM_BH/(3c²r_isco)]^1/2 (Bardeen 1970; Thorne 1974).

Fig. D.1. Examples of two-dimensional slices of the four-dimensional grid for 4 different metallicities, Z ∈ [0.00012, 0.00120, 0.00457, 0.01187], and 3 different BH masses mBH ∈ [6.3, 15.87, 29.38] M⊙. The final second-born BH spin value a2 is coloured for each successful track according to the legend. Each successful run stopped because of carbon depletion (square markers) while other termination flags are shown in the bottom legend.

Appendix E: Grids of detailed BH–He-star models

We use our detailed binary stellar models to cover the 4D parameter space defined by initial metallicity Z, black hole mass m_BH1, He-star mass m_{He − star} and orbital period p. We run grids for 30 different metallicities ranging from log₁₀(Z) = − 4.0 to log₁₀(Z) = log₁₀(1.5 Z_⊙) ≃ − 1.593 in steps of log₁₀(Z) ≃ 0.083. For each metallicity, we run 11 BH masses in the log-range [2.5, 54.4] M_⊙ in steps of log₁₀(m_BH1/M_⊙) ≃ 0.134. The lower limit is the smallest theoretical BH mass (in our model) while the maximum BH mass is chosen to have the second to last BH mass at 40 M_⊙ near the PISN cut. For each metallicity and BH mass we run 17 He-star masses in the log-range [8, 80] M_⊙ in steps of log₁₀(m_{He − star}/M_⊙) ≃ 0.063 and 20 binary periods in the log-range [0.09, 8] days in steps of log₁₀(p/days) ≃ 0.103. We verified that smaller He-star masses do not lead to BH formation for any metallicity we consider. The maximum He-star mass and smallest orbital period where chosen to include any BH–He-star system produced by our COSMIC models. The maximal orbital period range ensure that we cover the parameter space well past the point where the BBH systems are merging within the Hubble time.

In total, we calculated roughly 110 000 new binary evolution sequences, as compared to about 16 000 used in Bavera et al. (2020). The fraction of failed MESA runs vary from 3% to 0.1% depending on metallicity. To minimise the loss of information created by the failed runs we rerun those models with an He-star mass increased by 5% which reduced the failed runs by a factor of 3. In Fig. D.1 we show the spin of the second-born BH for different two-dimensional slices of the four-dimensional parameter space for four different metallicities, Z ∈ [0.00012, 0.00120, 0.00457, 0.01187], and three different BH masses, m_BH ∈ [6.3, 15.87, 29.38] M_⊙. We can see how the tides are more efficient at lower metallicities. This is because the stellar winds of He-stars are metallicity dependent ∝(Z/Z_⊙)^0.85 (Vink et al. 2001) and widen more efficiently the binaries at larger metallicities and hence reduce the impact of tides.

These grids were used to determine the final outcomes and final parameters of the late-end evolution stage of the binary systems through linear interpolation. Each metallicity is interpolated separately. We want to interpolate six quantities B_i: the He-star mass and its carbon-oxygen core mass before the supernova, the resulting BH mass, the orbital period before the supernova, the lifetime of the BH–He-star binary and the spin of the second-born BH. Before interpolating each quantity, we log-transformed it and re-scale it to the interval [ − 1, 1] to reduce the interpolation error. The interpolation itself relies on building a Delaunay triangulation of the input data points followed by barycentric linear interpolation over the vertices of the (hyper)triangle containing the location of interest. To test the accuracy of the interpolation: we computed around 1000 new MESA tracks for each of the four metallicities shown in Fig. D.1 and calculate the relative error of each transformed and rescaled quantity X_i ≡ log₁₀(B_i)^{[ − 1, 1]} as Δ_i = |X_true, i − X_interp, i|/X_true, i. In Fig. E.1 we show the median relative errors of these quantities as a function of metallicity. In the median calculation we exclude all the systems not becoming BBHs. Half of the quantities have median relative errors independent of metallicity with the exception of the second-born BH spin, He-star mass and orbital period before supernova. The decrease of median error for a₂ and p_preSN is explained by the fact that at high metallicities the orbits widen more, neutralising tides and resulting in systems with zero spin. On the other hand, the median error of m_{He − star} increases because at higher metallicities this quantity does not depend linearly on the initial He star mass which is caused by stellar winds that cause the He-star to lose a non-negligible amount of mass. The largest relative median error is the lifetime of the BH–He-star system. This is not a problem by itself as the delay time between the binary formation and merger is dominated by the GW inspiral which is many order of magnitudes larger than this timescale. The increase of interpolation accuracy compared to Bavera et al. (2020) is due to the larger data set used here, the fact that we use a regular grid and interpolate each metallicity independently.

Fig. E.1. Median relative error as a function of metallicity of the log-transformed and re-scaled six quantities Bi: the He-star mass (blue) and its carbon-oxygen core mass (orange) before the supernova, the orbital period before the supernova (red), the second-born BH mass (violet), the spin of the second-born black hole (brown), and the lifetime of the BH–He-star binary (green).

Appendix F: Model comparison

We use Bayesian hierarchical modelling to determine the likelihood of observing N independent GW events ${\left\{x_i\right\}}_{i=1}^N$ given an astrophysical model described by a set of parameters λ (e.g. Mandel et al. 2019). We assume that each GW event is an independent observation and is characterised by a set of physical parameters θ,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill p({\{x_i\}}_{i=1}^N|\mathit{\lambda })& ={\displaystyle \prod _{i=1}^N\frac{p(x_i|\mathit{\lambda })}{\xi (\mathit{\lambda })}=\prod _{i=1}^N\frac{1}{\xi (\mathit{\lambda })}\int p(x_i|\mathit{\theta })p(\mathit{\theta }|\mathit{\lambda })\mathrm{d}\mathit{\theta }=}\hfill \\ \hfill & ={\displaystyle \prod _{i=1}^N\frac{p(x_i)}{\xi (\mathit{\lambda })}\int \frac{p(\mathit{\theta }|x_i)}{p(\mathit{\theta })}p(\mathit{\theta }|\mathit{\lambda })\mathrm{d}\mathit{\theta },}\hfill \end{array}\end{array}$$(F.1)

where we have marginalised over the physical parameters of the individual events and used Bayes’ theorem to obtain the final line. Here p(θ) is the prior on the physical parameters that are used to generate the posterior samples. The normalisation factor ξ(λ) = ∫p_det(θ)p(θ|λ)dθ accounts for the overall probability of making an observation given a particular choice λ. Since in GWTC-2 we have samples drawn from the posterior p(θ | x_i), we can approximate posterior-weighted integrals as a sum over samples as

$$\begin{array}{c}\hfill p({\{x_i\}}_{i=1}^N|\mathit{\lambda })\simeq {\displaystyle \prod _{i=1}^N\frac{p(x_i)}{\xi (\mathit{\lambda })}\frac{1}{S}\sum _{k=1}^S\frac{p({\mathit{\theta }}_k|\mathit{\lambda })}{p({\mathit{\theta }}_k)},}\end{array}$$(F.2)

where S the number of posterior samples. Again, θ_k are the astrophysical parameters drawn from the GW posterior distribution and p(θ_k | λ) are their likelihood given by our astrophysical model (underlying distribution). Similarly, the normalisation constant ξ(λ) can be approximated as a sum over the weighed underlying BBH distribution p(θ|λ) as

$$\begin{array}{c}\hfill \xi (\mathit{\lambda })\simeq \frac{1}{T{\displaystyle {\sum }_{j=1}^T}{w}_j}\sum _{j=1}^T{p}_{\det }({\mathit{\theta }}_j)w_j,\end{array}$$(F.3)

where θ_j is the set of parameters describing the BBH j with a cosmological weight w_j given by the argument of the summation in Eq. (4) and T the total number of samples. Finally, to approximate the probability density functions of each event p(θ) and each model p(θ|λ) we use a three-dimensional kernel density estimator (KDE) where we graphically verify the accuracy by comparing random draws from the KDE with the real marginalised one-dimensional and two-dimensional sample distributions.

We can now compare two models, M₁ and M₂, described by λ₁ and λ₂, respectively. The amount by which the data supports a specific model is described by the BF, defined as:

$$\begin{array}{c}\hfill {\mathrm{BF}}_{12}\equiv \frac{p({\{x_i\}}_{i=1}^N|{\mathit{\lambda }}_1)}{p({\{x_i\}}_{i=1}^N|{\mathit{\lambda }}_2)}.\end{array}$$(F.4)

In the ratio of the two likelihoods the multiplicative constant p(x_i) is canceled out, leaving us with all the information required to compute the BF. This factor indicates whether any model is favoured or disfavoured by the data compared to another. Values larger than 1 favour the model M₁ while values smaller than 1 favours the model M₂. In our analysis, we adopt the convention of having M₂ as our reference model.

Appendix G: Model results

Here we present some extra figures which were not included in the paper because of their size. Figures G.1 and G.2 show the combined distributions of the main GW observables χ_eff, M_chirp and q of the CE channel for the detected and underlying BBH population, respectively, for different α_CE values. Similarly, Figs. G.3 and G.4 show the combined distributions of these observables of the SMT channel for the detected and underlying BBH population, respectively, for different η_acc values.

Fig. G.1. Model predictions for the O3 detected BBH population of the CE channel for different αCE values according to the legend, in orange. We show the joint distributions of chirp mass Mchirp, effective inspiral spin parameter χeff and binary mass ratio q. Lighter colours represent larger contour levels of 68%, 95% and 99%, respectively, constructed with pygtc (Bocquet & Carter 2016). All histograms are plotted with 30 bins without any bin smoothing.

Fig. G.2. Model predictions for the underlying (intrinsic) BBH population of the CE channel for different αCE values according to the legend, in grey. We show the joint distributions of chirp mass Mchirp, effective inspiral spin parameter χeff and binary mass ratio q. Lighter colours represent larger contour levels of 68%, 95% and 99%, respectively, constructed with pygtc (Bocquet & Carter 2016). All histograms are plotted with 30 bins without any bin smoothing.

Fig. G.3. Model predictions for the O3 detected BBH population of the SMT channel for different ηacc values according to the legend, in orange. We show the joint distributions of chirp mass Mchirp, effective inspiral spin parameter χeff and binary mass ratio q. Lighter colours represent larger contour levels of 68%, 95% and 99%, respectively, constructed with pygtc (Bocquet & Carter 2016). All histograms are plotted with 30 bins without any bin smoothing.

Fig. G.4. Model predictions for the underlying (intrinsic) BBH population of the SMT channel for different ηacc values according to the legend, in grey. We show the joint distributions of chirp mass Mchirp, effective inspiral spin parameter χeff and binary mass ratio q. Lighter colours represent larger contour levels of 68%, 95% and 99%, respectively, constructed with pygtc (Bocquet & Carter 2016). All histograms are plotted with 30 bins without any bin smoothing.

All Tables

All Figures

Fig. 1. Example of a two-dimensional slice of the four-dimensional grid showing their initial BH–He-star orbital period pi (in days) and initial He-star mass mHe − star (in M⊙) for log10(Z) ≃  − 3 and mBH1 = 29.39 M⊙. The final mass mBH2 and spin a2 of the second-born BH, as well as merger timescale Tmerger, are coloured according to the legend of each panel. All successful MESA simulation stopped at carbon depletion (square markers) while other termination flags are shown in the bottom legend. The merger timescale colour bar is capped at 14 Gyr.

In the text

Fig. 2. Model predictions for the underlying (intrinsic) BBH population (grey) and the O3 detected BBH population (orange) for our reference model with αCE = 1.0 and ηacc = 1. We show the joint distributions of chirp mass Mchirp and effective inspiral spin parameter χeff for the combined CE and SMT channels (left), CE formed BBHs (centre) and SMT formed BBHs (right). Lighter colours represent larger contour levels of 68%, 95% and 99%, respectively, constructed with pygtc module (Bocquet & Carter 2016). All histograms are plotted with 30 bins in the same range without any bin smoothing. We overlaid in grey the O1, O2 and O3a LVC GWTC-2 data with their 90% credible intervals. The 9 events of GWTC-2 in tension with our models are indicated in red (see Sect. 3.3), GW190521 is outside the plotted window.

In the text

Fig. 3. Orbital separation of BH–He-star binaries post CE for αCE ∈ [0.2, 0.5, 1.0] represented with solid lines according to the legend. The histogram has units of M⊙ Gpc−3 and accounts for the total stellar mass formed per comoving volume integrated over the Universe cosmic history per log-orbital period bin. The grey shaded area represent the 90% CI of the systems forming merging BBHs with Tmerger < 14 Gyr for the αCE = 0.2 model (the other models have similar CIs). As the CE efficiency is lowered, the orbital separations become smaller and the distributions move to the left: for αCE = 0.2 the orbital separation peak enters the grey area boosting the merger rate.

In the text

Fig. 4. Orbital separation of BH–He-star binaries post SMT for accretion efficiency limit ηacc ∈ [1, 105, 109] (in units of Eddington-limit) represented with solid lines according to the legend. The histogram has units of M⊙ Gpc−3 and accounts for the total stellar mass formed per comoving volume integrated over the Universe cosmic history per log-orbital period bin. The grey shaded area represent the 90% CI of the systems forming merging BBHs with Tmerger < 14 Gyr for the model with ηacc = 1 (the other models have similar CIs). As the accretion efficiency is increased, the orbital separations become larger and the distributions move to the right: for highly super-Eddington accretion efficiency ηacc = 109 the tail almost exits the grey area, hence, decreasing the merger rate.

In the text

Fig. 5. Model predictions for the O3 detected BBH population of CE and SMT channels combined for αCE ∈ [0.2, 0.35, 0.5, 1.0, 2.0, 5.0] and ηacc = 1, in orange. We show the joint distributions of chirp mass Mchirp, effective inspiral spin parameter χeff and binary mass ratio q. Lighter colours represent larger contour levels of 68%, 95% and 99%, respectively, constructed with pygtc (Bocquet & Carter 2016). All histograms are plotted with 30 bins without any bin smoothing. We overlaid in grey the LVC GWTC-2 data with their 90% credible intervals. The 9 events of GWTC-2 in tension with our models are indicated in red (see Sect. 3.3). We also show the one-dimensional projection of each quantity overplotted on the underlying (intrinsic) BBH population, in grey. For visualisation purposes, the χeff probability density function (pdf) is plotted in log-scale. For a visualisation of each model separately see Appendix G.

In the text

Fig. 6. BFs of CE and SMT models with respect to the reference model with αCE = 1.0 and Eddington-limited accretion efficiency as a function of sample size. The initial sample contains 38 GW events of GWTC-2 and exclude GW170729, GW190413_134308, GW190519_153544, GW190521, GW190602_175927, GW190620_030421, GW190701_203306, GW190706_222641, GW190929_012149. At each iteration the event with lowest likelihood with respect to the reference model is removed and indicated on the horizontal axis until 20 events are removed. By definition the BF of the model with αCE = 1.0 is equal to 1. The data show moderate to strong support for the models with low CE efficiency, αCE < 1.0. This result is robust because the BFs show a constant behaviour as a function of sample size.

In the text

	Fig. 7. Initial orbital separation of 105 binaries sampled from the extended Sana et al. (2012) distribution, cf. Eq. (A.1), in blue, and from a log-uniform distribution, in orange. Both distributions sample the range p ∈ [0.4, 105.5] days and are independent of metallicity. With a solid line, we show the same distributions after removing systems with Roche-lobe overflow at ZAMS for Z = 0.1 Z⊙.
In the text

Fig. D.1. Examples of two-dimensional slices of the four-dimensional grid for 4 different metallicities, Z ∈ [0.00012, 0.00120, 0.00457, 0.01187], and 3 different BH masses mBH ∈ [6.3, 15.87, 29.38] M⊙. The final second-born BH spin value a2 is coloured for each successful track according to the legend. Each successful run stopped because of carbon depletion (square markers) while other termination flags are shown in the bottom legend.

In the text

	Fig. E.1. Median relative error as a function of metallicity of the log-transformed and re-scaled six quantities Bi: the He-star mass (blue) and its carbon-oxygen core mass (orange) before the supernova, the orbital period before the supernova (red), the second-born BH mass (violet), the spin of the second-born black hole (brown), and the lifetime of the BH–He-star binary (green).
In the text

Fig. G.1. Model predictions for the O3 detected BBH population of the CE channel for different αCE values according to the legend, in orange. We show the joint distributions of chirp mass Mchirp, effective inspiral spin parameter χeff and binary mass ratio q. Lighter colours represent larger contour levels of 68%, 95% and 99%, respectively, constructed with pygtc (Bocquet & Carter 2016). All histograms are plotted with 30 bins without any bin smoothing.

In the text

Fig. G.2. Model predictions for the underlying (intrinsic) BBH population of the CE channel for different αCE values according to the legend, in grey. We show the joint distributions of chirp mass Mchirp, effective inspiral spin parameter χeff and binary mass ratio q. Lighter colours represent larger contour levels of 68%, 95% and 99%, respectively, constructed with pygtc (Bocquet & Carter 2016). All histograms are plotted with 30 bins without any bin smoothing.

In the text

Fig. G.3. Model predictions for the O3 detected BBH population of the SMT channel for different ηacc values according to the legend, in orange. We show the joint distributions of chirp mass Mchirp, effective inspiral spin parameter χeff and binary mass ratio q. Lighter colours represent larger contour levels of 68%, 95% and 99%, respectively, constructed with pygtc (Bocquet & Carter 2016). All histograms are plotted with 30 bins without any bin smoothing.

In the text

Fig. G.4. Model predictions for the underlying (intrinsic) BBH population of the SMT channel for different ηacc values according to the legend, in grey. We show the joint distributions of chirp mass Mchirp, effective inspiral spin parameter χeff and binary mass ratio q. Lighter colours represent larger contour levels of 68%, 95% and 99%, respectively, constructed with pygtc (Bocquet & Carter 2016). All histograms are plotted with 30 bins without any bin smoothing.

In the text