EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 635 (March 2020) A&A, 635 (2020) A97 Full HTML
Free Access
Issue
A&A
Volume 635, March 2020
Article Number A97
Number of page(s) 18
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/201936204
Published online 16 March 2020

© ESO 2020

1. Introduction

During the first and second observing runs O1/O2 of the advanced gravitational-wave (GW) detector network, Advanced LIGO (aLIGO; Aasi et al. 2015) and Advanced Virgo (Acernese et al. 2015) detected 10 GWs from binary black holes (BBHs). With the third observing run O3 that just started, this number is expected to increase. In the near future, sometime around 2020, the detectors will be upgraded to reach design sensitivity and we expect the detection of hundreds of such mergers (Abbott et al. 2019a).

To date, the origin of these BBHs is still an open scientific question. Various explanations of different formation channels for merging BBHs have entered into the scientific literature (see, e.g., Abbott et al. 2016a; Miller 2016; Mandel & Farmer 2018, for reviews). The most popular ones are as follows: isolated binary evolution where (i) the stars go through a common envelope (CE) phase due to an unstable mass transfer after the formation of the first-born black hole (BH) (e.g., Smarr & Blandford 1976; van den Heuvel 1976; Tutukov & Yungelson 1993; Kalogera et al. 2007; Postnov & Yungelson 2014; Belczynski et al. 2016a), (ii) massive stars with a nonextreme mass ratio after the formation of the first-born BH goes through stable mass transfer avoiding the CE phase (e.g., van den Heuvel et al. 2017; Pavlovskii et al. 2017; Inayoshi et al. 2017), (iii) massive stars in close orbits experiencing strong internal mixing go through chemically homogeneous evolution and produce massive BBHs (e.g., de Mink et al. 2009; Mandel & de Mink 2016; Marchant et al. 2016); dynamical formation (iv) in globular clusters and (v) galactic nuclear clusters where the BBHs are formed from stars not born in the same binary (e.g., Sigurdsson & Hernquist 1993; Portegies Zwart & McMillan 2000; Miller & Lauburg 2009; Rodriguez et al. 2015; Antonini et al. 2016); or (vi) Lidov–Kozai resonance bringing the inner binary to merge in hierarchical triple systems (e.g., Silsbee & Tremaine 2017). All of these scenarios possess some significant uncertainties in the prediction of merger rates due to the poorly constrained underlying physics or unconstrained distributions of initial conditions. The merger rate predictions for the isolated binary evolution via the CE phase are consistent (Abbott et al. 2017a) with the observed rate of BBH mergers of around ∼24−112 Gpc−3 yr−1 (Abbott et al. 2019b). The same holds for the stable mass transfer channel (Neijssel et al. 2019), while formation via chemically homogeneous evolution could yield tens of mergers per Gpc−3 yr−1 (Mandel & de Mink 2016; Marchant et al. 2016). Finally, predicted rates via the dynamical formation channel are closer to the lower end of the observed range (Fragione et al. 2018; Park et al. 2017); for example Rodriguez & Loeb (2018) find 4−18 Gpc−3 yr−1 from globular cluster in the local Universe.

Any astrophysical BH can be fully described by only two quantities: its mass M and its dimensionless spin parameter, a = cJ/(GM2), where J is the angular momentum of the BH. Using matched-filtering analysis, GW observations provide estimates for each of the above-mentioned quantities for both parent BHs. Although individual BH spin magnitudes and orientations are poorly constrained with present GW measurements, the effective inspiral spin parameter

χ eff = M 1 a 1 + M 2 a 2 M 1 + M 2 L ^ , $$ \begin{aligned} \chi _\mathrm{eff} = \frac{M_1 \mathbf a _1 + M_2 \mathbf a _2}{M_1+M_2} {\hat{\boldsymbol{L}}}, \end{aligned} $$(1)

the mass-weighted spin of the system projected onto the orbital angular momentum L, is reasonably well constrained (Abbott et al. 2019a). This is explained by the fact that the leading spin-orbit-coupling term in the post-Newtonian waveforms is dominated by this parameter (Santamaría et al. 2010). From the ten observed χeff, 8 are consistent with 0 within the 90% credible interval while the remaining two are determined with a positive value of χeff. Another important quantity characterizing the waveforms is the chirp mass

M chirp = ( M 1 M 2 ) 3 / 5 ( M 1 + M 2 ) 1 / 5 , $$ \begin{aligned} M_\mathrm{chirp} = \frac{(M_1 M_2)^{3/5}}{(M_1+M_2)^{1/5}}, \end{aligned} $$(2)

which, to first-order approximation, determines the frequency evolution of the GW signal emitted during the BBH’s inspiral phase (Cutler & Flanagan 1994). The ten observed Mchirp span the range of 7.9−35.7 M⊙ with a pile-up around 26 M⊙. In addition, the luminosity distance can be measured using the GW amplitude and, assuming a cosmological model, the cosmological redshift of the merger can be inferred. The distributions of these parameters for a population of merging BBHs can be used to distinguish between different formation channels. As pointed out in the literature, the effective inspiral spin parameter is sensitive to the evolutionary path of BBHs (see e.g., Rodriguez et al. 2016). For isolated field binary channels, the spins of the two BHs are expected to be preferentially aligned with the orbital angular momentum, whereas, assuming effective exchange interaction, the spin directions of BBHs formed in dynamical environments are expected to be randomly, isotropically distributed (Abbott et al. 2016a; Farr et al. 2017).

In this study, we focus on BBHs formed through classical isolated binary evolution that go through the CE phase. The main evolutionary phases of this pathway are now summarized. At the beginning, the stars are born in a relatively wide binary where the initially more massive star, called the “primary”, reaches the end of its main sequence first. At this stage the primary star expands its hydrogen-rich envelope past the Roche-lobe and begins transferring mass to the secondary until it loses its entire envelope, leaving a naked helium-burning star. Following wind-driven mass loss the primary collapses to form a BH. When the secondary reaches the end of its main sequence, the process repeats itself in reverse. This time, the mass transfer onto the black hole is unstable and this leads to the formation of a CE of gas around the binary (Paczynski 1976). The physical details of this phase are still not fully understood (Ivanova et al. 2013). The drag force on the BH from the envelope leads to a rapid inspiral and the dissipated orbital energy leads to the expulsion of the envelope and a decay of the orbital separation by more than two orders of magnitude. At this stage we are left with the immediate progenitor of the BBH system, namely a BH–He-star binary. Finally, the secondary eventually collapses into a BH and potential asymmetries in the core collapse may impart a kick on the newly formed BH and alter the orbit further. Eventually, due to energy and angular momentum loss from GW emission, the BBH system can coalesce into a single, more massive BH.

Previous theoretical works focused on the first few observed GW events suggest that these BBHs are consistent with having been formed through the CE formation channel (Stevenson et al. 2017; Giacobbo et al. 2018; Kruckow et al. 2018). These authors show how, at the respective appropriate metallicity regime, the observed BH masses are produced by their binary evolution models. Furthermore, their inferred merger rates are consistent with the one obtained from GW observations. In another study in favor of the CE formation channel, Belczynski et al. (2016b) carried out a detailed analysis of merger rates and found that BBHs formed though this channel should dominate the event rates in Advanced LIGO and Virgo.

In the CE formation channel, the physical process determining the spin of the first-born BH is the efficiency of angular momentum (AM) transport through the evolution of the progenitor star during the red supergiant phase. From observations of astro-seismology (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 is known that this mechanism must be efficient (Spruit 1999, 2002; Fuller et al. 2019). Thus, upon expansion, the initial angular momentum is mostly transported to the outer layers of the star which are subsequently lost due to Roche-lobe overflow mass transfer and wind mass loss. This leads to very slowly spinning BHs (a1 ≃ 0) as was shown by Qin et al. (2018; see also Fuller & Ma 2019). Assuming efficient AM transport, the angular momentum of the second-born BH is mainly determined by the net effect of the stellar winds and the tidal interaction of the BH–He-star binary system. This is because any initial or acquired rotation during the evolution of the secondary is erased through mass transfer and wind mass loss by the time it becomes a He-star. Several studies attempted to model the last evolutionary phase of this channel and derived constraints on the spin using analytical arguments and semi-analytical calculations (Kushnir et al. 2016; Hotokezaka & Piran 2017; Zaldarriaga et al. 2018). These studies found out that around half of the secondary BHs have zero spin and the other half are maximally spinning. When using detailed models to simulate the binary evolution and the stellar structure of the two components, Qin et al. (2018) did not reproduce this prediction of a bi-modal distribution of spins. Both Hotokezaka & Piran (2017) and Zaldarriaga et al. (2018) results are based on the approach outlined in Kushnir et al. (2016). Compared to the detailed binary simulations of Qin et al. (2018), these authors did not model self-consistently the orbit evolution of the binary due to the combined effects of tides and stellar winds, which in most cases leads to the widening of the orbit. Even when tides are initially efficient at synchronizing the spin of the helium star to the orbit, such wind-driven orbital widening can lead to tidal decoupling. Ignoring this effects underestimates the impact of stellar winds on the final spin of the second-born BH. Moreover they used approximate timescales for the process of tidal synchronization that do not take into account changes in the structure of the star during its lifetime and assumed that tides allow the He-star to remain tidally locked indefinitely. These approximations lead to results that disagree with what is found in our detailed binary simulations. Belczynski et al. (2017) used parametric binary population synthesis models, which share all the same approximations as the studies discussed above, to compare three different prescriptions for the efficiency of AM transport. They found that efficient AM transport is favored, as it results to distributions of χeff and BH masses qualitatively consistent with observations, while inefficient AM transport would lead to rapidly spinning BHs which are currently not observed by aLIGO.

In this paper, we present a model capable of predicting simultaneously the spin, mass and redshift distributions of coalescing BBHs formed from isolated field binaries that go through the CE phase. This aim is achieved by combining the parametric binary population syntheses code COMPAS with the detailed MESA stellar structure and binary evolution simulations. The study is structured as follows. In Sect. 2 we explain the methods used to generate a Monte-Carlo population of isolated field BBHs and how we take into account the redshift and metallicity dependence of the star formation rate and the observational selection effects. In Sect. 3 we present our main results and make detailed predictions for observing runs of the future GW detector network, distinguishing three regions of the BBH parameter space. Our model is able to successfully reproduce the observed χeff distribution providing strong support for the CE channel being the dominant formation channel. We then discuss our results in Sect. 4 where we compare our work to the current literature and demonstrate the importance of detailed binary evolution calculations. Finally, the conclusions of our work are given in Sect. 5.

2. Method

We use the parametric binary population synthesis code COMPAS (e.g., Stevenson et al. 2017, 2019; Barrett et al. 2017, 2018; Vigna-Gómez et al. 2018; Neijssel et al. 2019) to evolve isolated stellar binaries until the formation of BH–He-star systems, namely the immediate progenitors of BBHs. In Sect. 2.1 we briefly describe the COMPAS-model assumptions used in the simulation. Since we are not interested in a parameter study, we specifically picked a model capable of reproducing a BBHs merger rate which is in agreement with the observed one from Abbott et al. (2019a). For the last step of the evolution, which we consider to be the one that determines the spin distribution of the secondary BH, we use the stellar structure and evolution code MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019) to simulate the evolution of the binary systems. Assuming that the first-born black hole can be treated as a point-like particle, this approach allows us to track the angular momentum profile evolution of the He-star until the formation of the secondary BH (see Sect. 2.2). In Sect. 2.3, we explain in detail how we treat the collapse of the He-stars into BHs. Finally, taking into account the redshift dependence of metallicity, star-formation rate (SFR) and aLIGO sensitivity, we can distribute the population across cosmic time and compute the expected detection rate (see Sect. 2.5).

2.1. COMPAS model assumptions

We use the results from the simulations of Neijssel et al. (2019) and here we only highlight the main physical assumptions. We assume that the underlying stellar population spans the mass range 0.01 M⊙ <  m1 <  200 M⊙ following the initial mass function of Kroupa (2001). The mass distribution of the less massive secondary star is given by m2 = m1q0 where q0 is the initial mass ratio (0 <  q0 <  1) drawn from a flat distribution (Sana et al. 2012). We are interested in binaries with a primary star that ends up forming a BH, thus we restrict the initial mass distribution of primary masses between 5 M⊙ <  m1 <  150 M⊙. This means we only model a fraction fcorr of the underlying stellar population mass. We calculate this by assuming a binary fraction of 70% (Sana et al. 2012), see Appendix A. We assume that, at formation, binaries are distributed uniformly in log-orbital separation restricted to 0.1 <  A/AU <  1000 (Abt 1983) and have zero eccentricity. We assume that all these distributions are independent from each other as well as independent of metallicity. For the metallicity distribution of binaries we divide uniformly in 30 bins the log-metallicity range Z ∈ [0.0001, 0.0349]. We then evolve three million binaries per metallicity bin ΔZj with a total of star forming mass on the order of Msim, ΔZj = 6.5 × 107 M⊙.

COMPAS evolves stars according to the stellar models of Pols et al. (1998) and uses analytical fits of these models to rapidly evolve binaries (Hurley et al. 2000, 2002). We adopt wind mass loss rates as prescribed by Hurley et al. (2000) for stars with effective temperatures smaller than 12 500 K, and for hotter stars the winds of Vink et al. (2001) as implemented by Belczynski et al. (2008). If stars during their evolution cross the Humphrey-Davidson limit (Humphreys & Davidson 1994) and enter a region of the Hertzsprung–Russel diagram in which no stars are observed, we apply an additional wind mass loss rate of 1.5 × 10−4 M⊙ yr−1 (Belczynski et al. 2010).

When the primary star reaches the end of its main sequence, the star expands and loses its entire envelope through Roche-lobe overflow. In binaries where the first mass transfer episode is stable, the companion can accrete some mass with an efficiency that we assume depends on the ratio of the thermal timescales of the two stars (Hurley et al. 2002; Schneider et al. 2015; Neijssel et al. 2019), while the mass not accreted by the other star leaves the system carrying away the specific angular momentum of the accretor. Eventually, the envelope of the primary is stripped, leaving a naked helium burning star which, following wind-driven mass loss, collapses into a BH. The star collapses into a point-like BH following the “delayed” model of Fryer et al. (2012).

When the secondary reaches the end of its main sequence the process repeats in reverse and the mass transfer between the BH and the He-star can be either dynamically stable or unstable. Since we focus only on the subchannel that goes through the CE phase, we consider exclusively systems with dynamically unstable mass transfer which produces a non co-rotating CE of gas engulfing the binary. This represents only a subset of all merging binary black holes in the models of Neijssel et al. (2019). While uncertainties in the stability of mass transfer could reduce the importance of the non-CE channel, a self-consistent variation of the assumptions regarding mass transfer stability would also change the population of systems that evolve through a CE phase. This analysis is beyond the scope of this work, but could impact, in particular, our overall rate predictions for BBH mergers, which should be compared directly to Neijssel et al. (2019).

COMPAS uses the classical energy αCE − λ formalism (Webbink 1984; de Kool 1990; Dewi & Tauris 2000; Xu & Li 2010) to parameterize the uncertainties in the physics of the CE phase. During this phase the two stars spiral in due to friction with the envelope. The loss of orbital energy can heat up and expand the envelope. The efficiency of this energy transfer is parameterized by the αCE parameter which can vary (Livio & Soker 1988). We assume that all of the dissipated orbital energy goes into expelling the envelope, αCE = 1. The λ parameter, which characterizes the binding energy of the CE, depends on the structure of the donor’s envelope (de Kool 1990). We chose our λ according to the fits of Xu & Li (2010) as implemented by Dominik et al. (2012).

Within the CE subchannel, we distinguish two different scenarios for donor stars which are on the Hertzsprung-gap (HG). In the optimistic scenario we apply the usual α − λ prescription to evolve these systems. In the pessimistic scenario we assume that Hertzsprung-gap stars have not yet developed a sufficiently sharp density gradient at the core-envelope boundary to allow the inspiral during the CE to stop. Thus any CE event from donors in this evolutionary phase results in a merger which reduces the BBH merger rate (Belczynski et al. 2007). In this paper we present the results for the latter scenario. Both scenarios yield a similar distribution of spins, but the pessimistic one predicts fewer low-mass BHs compared to the optimistic. This is because a greater fraction of the total post-main-sequence expansion occurs during the HG for high-metallicity stars (Linden et al. 2010). Therefore, forbidding the channel with HG CE donors has a particularly strong effect at high metallicity, which yields lower-mass BHs due to metallicity-enhanced stellar winds.

In Fig. 1 we show the distributions of orbital separation, first-born BH mass, He-star mass and metallicity of our BH–He-star binaries after the CE phase. These distributions are weighted by the redshift- and metallicity-dependent star-formation rate (SFR) integrated over the cosmic time, see Eq. (B.10). We note that these distributions include all binaries, including those systems that are going to become BH–NS binaries and BBHs with GW inspiral timescales longer than the Hubble time. After the CE phase the orbital separations are no longer uniformly log-distributed and most of the first-born BHs have masses smaller than 30 M⊙. Moreover, we see that formation metallicities of progenitors of merging compact-object binaries follow a skewed log-normal distribution. This is because the mean metallicity of the Universe decreases as a function of the look-back time and the star-formation rate peaks at a redshift ∼2, that is a look-back time of ∼10.5 Gyr, where most of the binaries are formed. These distributions are used as an initial condition for our detailed modeling.

thumbnail Fig. 1.

Parameter distributions of the binary population after the CE phase. These BH–He-star binaries include systems that are going to form BH–NS binaries and BBHs with GW inspiral timescales bigger than the Hubble time. We show the distributions of orbital separation A, first-born BH mass MBH, He-star mass MHe-star and metallicity Z weighted by the integrated redshift- and metallicity-dependent SFR over the cosmic time (see Eq. (B.10)). The lighter shades represent larger contour levels, 68, 95 and 99%, respectively, constructed with pygtc (Bocquet & Carter 2016).

2.2. MESA model assumptions

We perform detailed stellar structure and binary evolution calculations that take into account wind mass loss, internal differential rotation of the He-star and tidal interaction between the BH and the He-star. These simulations1 are based on the work of Qin et al. (2018) and are adapted for MESA r-10398.

Stellar winds play an important role in binary evolution. Here, we take a slightly different approach compared to Qin et al. (2018), and we follow the wind prescriptions outlined in Belczynski et al. (2010), which is the same as the ones used in COMPAS. Namely, for helium stars we adopt a wind mass-loss rate of

d M d t = 10 − 13 ( L L ⊙ ) 1.5 ( Z Z ⊙ ) 0.86 M ⊙ yr − 1 , $$ \begin{aligned} \frac{\mathrm{d}M}{\mathrm{d}t} = 10^{-13}\left(\frac{L}{L_\odot }\right)^{1.5}\left(\frac{Z}{Z_\odot }\right)^{0.86}\, {M}_\odot \, \mathrm{yr} ^{-1} , \end{aligned} $$(3)

where L and Z are the star’s luminosity and metallicity, respectively. This prescription is a combination of Hamann & Koesterke (1998) and Vink & de Koter (2005) and takes into account He-star winds clumping and a strong dependence on the metallicity. Furthermore, we adopt Z⊙ = 0.017 as solar metallicity (Grevesse & Sauval 1998).

Tidal forces are responsible for synchronising the spin of the He-star with the orbit. We assume that the CE ejection leaves a circular binary, and the system remains circular during He-star evolution. It has been suggested that dynamical tides are dominant for stars with a radiative envelope and a convective core (Zahn 1977; Hut 1981). The strength of the interaction depends on the ratio of the stellar radius R to the orbital separation A. The timescale for synchronization is defined as

1 T sync = 3 E 2 ( 1 + q ) 5 / 6 q 2 r g 2 ( GM R 3 ) 1 / 2 ( R A ) 17 / 2 , $$ \begin{aligned} \frac{1}{T_\mathrm{sync} } = 3 E_2 \left(1+q\right)^{5/6} \frac{q^2}{r_\mathrm{g} ^2} \left( \frac{GM}{R^3} \right)^{1/2} \left( \frac{R}{A} \right)^{17/2} , \end{aligned} $$(4)

where the He-star has mass M, radius R and moment of inertia I, rg given by r g 2 =I/(M R 2 ) $ r_\mathrm{g}^2=I/(MR^2) $ is the dimensionless gyration radius of the He-star, q is the mass ratio of the BH mass to the He-star mass, and E2 is the second order tidal coefficient. We take the new fitting formula of E2 as suggested by Qin et al. (2018) for He-stars

E 2 = 10 − 0.93 ( R conv R ) 6.7 , $$ \begin{aligned} E_2 = 10^{-0.93} \left( \frac{R_{\rm conv}}{R} \right)^{6.7} , \end{aligned} $$(5)

where Rconv is the radius of the convective core (see Appendix A in Qin et al. 2018, for an in-depth discussion of E2). We highlight here that a variation of the implementation of tides is used (Qin et al. 2019). Instead of the standard tides prescription in MESA (Paxton et al. 2015) that synchronize the whole star, the tides here only operate on the radiative regions. However it has been verified that this slight variation has a very small impact on our results.

Rotational mixing and angular momentum transport are treated as diffusion processes (Heger et al. 2000, 2005), which mainly involve the effects of Eddington-Sweet circulation, the Goldreich-Schubert-Fricke instability, and secular as well as dynamical shear mixing. In addition, diffusion element mixing is included with an efficiency parameter of fc = 1/30 (Chaboyer & Zahn 1992; Heger et al. 2000) for all processes above. Furthermore, an efficient angular momentum transport mechanism (i.e., Tayler–Spruit dynamo: Spruit 1999, 2002) is included. For comparison, we also ran a small grid without the Tayler–Spruit dynamo and found that there is a negligible impact on our results. More details on this can be found in the discussion. Furthermore, efficient AM transport allows us to assume that all He-stars emerging from the CE phase 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.

Running these simulations for all binary systems computed by COMPAS is computationally too expensive. Therefore we run a grid that allows us to infer through interpolation the six parameters we are interested in, namely: the He-star mass before the supernova, the carbon-oxygen (CO) core mass pre-supernova, the resultant second-born BH mass, the orbital period pre-supernova, the lifetime of the BH–He-star binary system (from the expulsion of the CE to the collapse of the He-star) and the spin of the second-born BH, as a function of the initial parameters of the BH–He-star binary: initial mass of the first-born BH, mBH, initial mass of the He-star, mHe-star, initial orbital period, p, and He-star metallicity, Z. In order to optimally construct our grid, we first randomly generate 3000 points to cover the parameter space spanned by the binaries after the CE phase, namely, mBH ∈ [2.5 M⊙, 60 M⊙], mHe-star ∈ [2.5 M⊙, 89 M⊙], p ∈ [0.05 days, 8.5 days] and Z ∈ [0.0001, 0.0349], to which we add 1500 points drawn from a kernel density estimator (KDE) trained with the post CE phase parameters of the synthetic population. In Appendix C we explain how we try to minimize linear interpolation errors by running more simulations where the interpolator is under-performing. The accuracy of the linear interpolator at each step is verified conducting 50 “leave 5% of the sample out” tests. We run a total of 18 000 simulations and show how the median relative errors stabilize at 0.01% and 0.04% for the He-star mass pre-supernova and resultant BH mass respectively, 0.20% for the CO core mass of the He-star, 0.01% for the orbital period, 0.04% for the lifetime of the BH–He-star binary and 0.41% for the log-spin of the second-born BH. In the appendix we also show the spread of the relative errors. If in the 50 leave 5% of the sample out tests we also count non-fittable points, such as remote points at the boundary of the parameter space, we find the following percentages of test systems that have relative errors above 10% in the estimated quantities: 5% and 8% of the He-star mass pre-supernova and resultant BH mass, 8% of the CO core mass of the He-star, 6% of the orbital period, 6% of the lifetime of the BH–He-star binary and 17% of the log-spin of the second-born BH.

2.3. Core-collapse physics

Black holes are formed during the core collapse of massive stars and, in some cases, their formation may be accompanied by supernova explosions. The collapse occurs when the stellar core begins to contract under its own weight without being able to trigger any more nuclear burning in its iron core. This leads eventually to electron capture and dissociation of the core elements into alpha particles. The first process removes the degeneracy pressure support of the core while the second removes the thermal support. These two mechanisms combined accelerate the collapse until the core reaches nuclear densities and neutron degeneracy pressure halts the collapse. This sudden halt produces a bounce shock moving out of the core. The shock-wave moves outwards until it deposits all its energy into the surrounding layers. A supernova explosion occurs if the deposited energy can overcome the ram pressure of the infalling stellar material. A fraction, ffb, of the material ejected by the supernova then falls back onto the stellar remnant. If the remnant is massive enough, neutron degeneracy pressure fails to halt the collapse and a black hole is formed. Moreover, the most massive stars can directly overcome the neutron degeneracy pressure when the collapse starts and implode to form a black hole. For a thorough review of our current understanding of the core-collapse process see for example Janka et al. (2007).

We use Fryer et al. (2012) delayed supernova prescription to model how much baryonic remnant mass is left behind after the collapse of the secondary star. This differs from their rapid prescription which produces a mass gap between BHs and neutron stars by assuming a strong convection which allows instabilities to grow quickly after the core bounce, producing more energetic SN explosion. The two prescriptions are not expected to lead to significant differences in the detected BBH merger rate (Belczynski et al. 2016b) as the population is dominated, due to aLIGO’s selection effects, by more massive BHs.

Using the delayed prescription, we calculated the fraction of the star that collapses to form the BH, and we assume that any remaining outer layers that do not collapse are instantaneously ejected. In order to estimate the spin of the resulting BH we follow the framework described in Batta & Ramirez-Ruiz (2019). We assume that there is no pressure stopping or slowing down the collapse. We can think of the star mass distribution M(r) as a collection of shells with mass mshell and angular frequency Ωshell, that falls one by one onto the center of the star. We assume that at the center, the shells up to 2.5 M⊙ collapse directly to form a black hole conserving their angular momentum and mass. Once a shell reaches the BH’s event horizon, it is accreted by it. The amount of angular momentum of the infalling material determines the properties of the accretion flow. Low angular momentum material collapses directly onto the BH transferring its entire mass and angular momentum to the BH, while material with enough angular momentum can create a disk around it. The maximum amount of angular momentum the disk material can give to the BH is determined by the specific angular momentum at the innermost circular orbit (ISCO) around the BH (Bardeen et al. 1972),

j isco = G M BH c 2 3 3 / 2 ( 1 + 2 ( 3 c 2 r isco G M BH − 2 ) 1 / 2 ) , $$ \begin{aligned} j_{\rm isco} = \frac{GM_{\rm BH}}{c} \frac{2}{3^{3/2}} \left(1+2 \left(\frac{3c^2 r_{\rm isco}}{G M_{\rm BH}}-2 \right)^{1/2} \right) , \end{aligned} $$(6)

where risco is the radius at ISCO for prograde equatorial orbits,

r isco = G M BH c 2 ( 3 + z 2 − ( ( 3 − z 1 ) ( 3 + z 1 + 2 z 2 ) ) 1 / 2 ) , $$ \begin{aligned} r_{\rm isco} = \frac{GM_{\rm BH}}{c^2}\left(3+z_2 - \left( \left(3 - z_1 \right)\left(3 + z_1+ 2z_2 \right)\right)^{1/2}\right) , \end{aligned} $$(7)

with z1 = 1 + (1 − a2)1/3((1 + a)1/3 + (1 − a)1/3) and z 2 = (3 a 2 + z 1 2 ) 1/2 $ z_2 =(3a^2+z_1^2)^{1/2} $ where a is the spin of the BH. Assuming that the disk formed from the collapse of a shell is accreted before the next shell collapses, as the viscous timescale of the disk is shorter than the dynamical timescale of the collapsing shells, we can evolve the BH’s mass and spin as it accretes material through the accretion disk. Each mass shell then contributes to the angular momentum of the BH by

J shell = ∫ 0 θ disk m shell Ω shell ( r ) r 2 sin 3 θ d θ + ∫ θ disk π / 2 m shell j isco sin θ d θ , $$ \begin{aligned} J_{\rm shell} = \int _0^{\theta _{\rm disk}} m_\mathrm{shell} \Omega _\mathrm{shell} (r) r^2 \sin ^3 \theta \, \mathrm{d}\theta + \int _{\theta _{\rm disk}}^{\pi /2} m_\mathrm{shell} j_{\rm isco} \sin \theta \, \mathrm{d}\theta , \end{aligned} $$(8)

where the disk formation angle is given by

θ disk = arcsin ( ( j isco Ω shell ( r ) r 2 ) 1 / 2 ) $$ \begin{aligned} \theta _{\rm disk} = \arcsin \left( \left(\frac{j_{\rm isco}}{\Omega _\mathrm{shell} (r)r^2} \right)^{1/2} \right) \end{aligned} $$(9)

and if the argument of arcsin exceeds 1, there is no disk formation. The first term in Eq. (8) represents material with low angular momentum that collapses directly onto the BH, while the second term corresponds to the material that forms the accretion disk with mass mdisk = mshell cos θdisk. The mass-energy accreted from the disk onto the BH is ΔMdisk = mdisk(1−2GMBH/(3c2risco))1/2 and the accreted angular momentum is ΔJdisk = mdiskjisco (Bardeen 1970; Thorne 1974). When treating the accretion of the portions of the shell that collapse directly onto the BH, we take into account 10% of baryonic mass loss through neutrinos (Fryer et al. 2012).

In our population synthesis study we neglect the effects of pair-instability supernovae (PISNe) and pulsational pair-instability supernovae (PPISNe). Both events are caused by the production of electron-positron pairs in the cores of very massive stars. In a PISN, pair production leads to a drop in the radiation pressure support in the core, causing the core to contract and the core temperature to increase. This results in explosive oxygen burning which reverses the collapse, unbinding the star. A PPISN is similar but the release of 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. PISNe cause massive stars with He-core masses between approximately 60 and 150 M⊙ to be completely disrupted. Therefore, PISNe put a second theoretical mass gap into the distribution of BH masses (Fowler & Hoyle 1964; Rakavy & Shaviv 1967; Barkat et al. 1967). On the other hand, PPISNe affect pre-supernovae stars with He-core masses between around 35 and 60 M⊙ enhancing the loss of mass before the supernova event and resulting in less massive BHs (Yoshida et al. 2016; Woosley 2017; Marchant et al. 2019). Neglecting these two phenomena leads us to overestimate the mass of the most massive BHs with MBH ≳ 35 M⊙. For a recent population synthesis study of the effect of PISNe and PPISNe on the population of coalescing BBHs using the same code as in this work see Stevenson et al. (2019).

During a supernova, the asymmetric ejection of matter (Janka & Mueller 1994; Burrows & Hayes 1996; Janka 2013) or asymmetric emission of neutrinos (Bisnovatyi-Kogan 1993; Socrates et al. 2005) can provide a momentum kick to the newly formed compact object. Here we assume that the birth kicks of BHs follow a Maxwellian distribution with σ = 265 km s−1 (Hobbs et al. 2005), which is then rescaled by one minus the fall-back mass fraction ffb (Fryer et al. 2012). In the Fryer et al. (2012) that we adopt, this quantity depends on the carbon-oxygen core mass mcore of the star before the collapse. For core masses grater than 11 M⊙, ffb = 1, which means that in our model all heavy black holes receive no natal kicks. These kicks can tilt the orbit of the BBH, which may generate a negative χeff, add eccentricity to the orbit or disrupt the binary. We take into account all these orbital changes, as well as orbital changes due to symmetric mass loss, following the analytical calculations of Kalogera (1996) and Abbott et al. (2017b).

2.4. Inspiral due to gravitational waves

After the birth of the second-born BH, GW emission removes energy and angular momentum from the orbit, shrinking it, and eventually leading to the merger of the two compact objects. The merger timescale for eccentric BBHs is computed as

T merger = 15 304 c 5 G 3 1 m 1 m 2 ( m 1 + m 2 ) A 4 f ( e ) , $$ \begin{aligned} {T}_{\rm merger} = \frac{15}{304} \frac{c^5}{G^3} \frac{1}{m_1m_2(m_1+m_2)} A^4 f(e) , \end{aligned} $$(10)

where m1 and m2 are the masses of the BHs, A is the orbital separation (Peters 1964) and f(e) is a numerical factor that account for the orbital eccentricity:

f ( e ) = ( 1 − e 2 ) 4 ∫ 0 e e ′ 29 / 19 ( 1 + 121 304 e ′ 2 ) 1181 / 2299 ( 1 − e ′ 2 ) 3 / 2 d e ′ e 48 / 19 ( 1 + 121 304 e 2 ) 3480 / 2299 · $$ \begin{aligned} f(e) = \frac{ \left(1-e^2\right)^4 \int _0^e \frac{e^{\prime 29/19}(1+\frac{121}{304}e{^\prime 2})^{1181/2299}}{(1-e^{\prime 2})^{3/2}} \mathrm{d}e^{\prime } }{e^{48/19} (1+\frac{121}{304}e^2)^{3480/2299}} \cdot \end{aligned} $$(11)

There is an important point to make here. As was already explained by other authors (e.g., Kushnir et al. 2016; Zaldarriaga et al. 2018; Qin et al. 2018), the merger timescale is anti-correlated with the spin of the second-born BH, a2, or the observed quantity χeff. This is because in order to form a fast rotating BH, tides should be strong and therefore the orbital separation between the parent He-star and the BH companion should be small. Since the merger timescale scales as the fourth power of the orbital separation, for tidally locked systems we can recover the following proportionality T merger ∼ A 4 ∼ ω −8/3 ∼ a 2 −8/3 ∼ χ eff −8/3 $ T_\mathrm{merger} \sim A^4 \sim \omega^{-8/3} \sim a_2^{-8/3} \sim \chi_\mathrm{eff}^{-8/3} $. In the second relation we used Kepler’s third law with ω being the orbital frequency matching the He-star’s angular frequency Ω and in the last one a1 = 0 as assumed in our model. Meanwhile, for the wider binaries partially synchronized by dynamical tides on a synchronization timescale T sync − 1 = | Ω ˙ | / Ω $ T_{\mathrm{sync}}^{-1} = |\dot{\Omega}|/\Omega $, we recover small spins a2 ∼ Ω ∼ exp(1/Tsync)∼1/Tsync ∼ A−17/2 (cf. Eq. (4)) and therefore T merger ∼ A 4 ∼ a 2 −8/17 ∼ χ eff −8/17 $ T_\mathrm{merger} \sim A^4 \sim a_2^{-8/17} \sim \chi_\mathrm{eff}^{-8/17} $. In Fig. 2 we show the combined distribution of the merger timescale and the effective spin parameter for a specific metallicity bin (centered at log10(Z) = −2.5) of our synthetic BBH population, namely not integrating over redshift and not accounting for any selection effect, in gray. Indeed, systems with high χeff follow the scaling relation for tidally locked systems (orange dashed line), while those with low χeff follow the scaling dictated by the tidal synchronization timescale (blue dashed line).

thumbnail Fig. 2.

Combined distribution of the BBH merger timescale Tmerger vs. the effective inspiral spin parameter χeff for our synthetic BBH population at metallicity log10(Z) = −2.5, in gray. The lighter shades represent larger contour levels, 68, 95 and 99%, respectively. Tidally locked systems follow the relation log 10 ( T merger ) ∼ − 8 3 log 10 ( χ eff ) $ \log_{10} ({T}_{\mathrm{merger}}) \sim - \frac{8}{3} \log_{10}(\chi_{\mathrm{eff}}) $, orange dashed line, while the others systems follow log 10 ( T merger ) ∼ − 8 17 log 10 ( χ eff ) $ \log_{10} ({T}_{\mathrm{merger}}) \sim - \frac{8}{17} \log_{10}(\chi_{\mathrm{eff}}) $ dictated by the tidal synchronization timescale, blue dashed line. Both lines are drawn at an arbitrary ordinate.

2.5. Detection rate

To compute the expected rate of detectable GW events, we need to convolve the star-formation rate (SFR) and metallicity distribution at different redshift epochs with the selection effects of the detectors. To do this we follow a similar approach to the one used in Belczynski et al. (2016b). Here we briefly summarize our approach, which is described in detail in Appendix B.

In our cosmological calculation we adopt the flat ΛCDM model with H 0 = 67.7 km s − 1 Mpc $ H_0=67.7 \, \frac{\text{km}\,\text{s} ^{-1}}{\text{ Mpc}} $ and Ωm = 0.307 (Planck Collaboration XIII 2016). Every simulated BBH k with BH masses m1, k and m2, k, born at redshift zf, i and merging at redshift zm, i, k set by the delay time of this binary contributes to the detection rate by the following weight

w i , j , k = fSFR ( z f , i ) M sim , Δ Z j f corr 4 π c D c 2 ( z m , i , k ) p det ( z m , i , k , m 1 , k , m 2 , k ) Δ t i $$ \begin{aligned} { w}_{i,j,k} = \frac{\text{ fSFR}(z_{\mathrm{f},i})}{M_{\mathrm{sim} ,\Delta Z_j}} \, f_{\rm corr} \, 4 \pi c\, D_{\rm c}^2(z_{\mathrm{m},i,k})\, p_{\rm det}(z_{\mathrm{m},i,k}, m_{1,k}, m_{2,k}) \, \Delta t_i \end{aligned} $$(12)

where subscripts f and m refer to formation and merger, respectively, and fSFR is the fractional star-formation rate, that is the total mass of stars formed per comoving volume per year per metallicity interval ΔZj. We adopt the SFR and metallicity distribution of Madau & Fragos (2017). The SFR formula we adopt is computed from UV and infrared surveys and is an updated version of Madau & Dickinson (2014) that better reproduces recent results at high redshifts 4 ≤ z ≤ 10 (Bowler et al. 2015; Finkelstein et al. 2015; Ishigaki et al. 2015; McLeod et al. 2015, 2016; Oesch et al. 2015). The metallicities are log-normally distributed with standard deviation 0.5 dex around the mean metallicity function of Madau & Fragos (2017). The mean metallicity function is obtained fitting observations assuming that the galaxy mass–metallicity relation holds at any redshift (Kajisawa et al. 2009; Baldry et al. 2012; Lee et al. 2012; Ilbert et al. 2013; Grazian et al. 2015). Furthermore, Msim, ΔZj/fcorr is the matter simulated in the metallicity bin ΔZj rescaled by the normalization factor fcorr (see Appendix A), D c (z)=c/ H 0 ∫ 0 z ( Ω m (1+ z ′ ) 3 + Ω Λ ) −1/2 d z ′ $ D_{\rm c}(z) = c/H_0 \smallint_0^z (\Omega_{\rm m} (1+z^{\prime})^3+\Omega_\Lambda)^{-1/2} {\rm d}z^{\prime} $ is the comoving distance to the source and pdet accounts for the selection effects of the detector. The total rate of detectable BBH mergers for a given detector network is calculated from the Monte Carlo simulations as a sum

R det = ∑ Δ t i ∑ Δ Z j ∑ k w i , j , k , $$ \begin{aligned} R_{\rm det} = \sum _{\Delta t_i} \sum _{\Delta Z_j} \sum _k { w}_{i,j,k}, \end{aligned} $$(13)

where we add the contribution of every binary placed at the center of each formation time bin Δti in its corresponding metallicity bin ΔZj. The population synthesis predictions are performed in finite time bins of Δt = 100 Myr and the log-metallicity range Z ∈ [0.0001, 0.0349] is divided in 30 bins.

The sensitivity of a GW detector to a source depends on the distance to the source, its orientation and position relative to the detector and on its physical characteristics. The detectability of a signal depends on its signal-to-noise ratio (S/N). In our model we assume that signals are detected if their single-detector S/N exceeds a threshold value of 8 (Aasi et al. 2016). For the two observing runs O1/O2 of aLIGO, we assumed a detector sensitivity equal to the target “early high sensitivity” (Abbott et al. 2018). This simplification is motivated by the fact that the sensitivity of O2 was close to that of O1 (Abbott et al. 2018). We assume a target “late low sensitivity” for the third observing run O3 and for design sensitivity the corresponding one (Abbott et al. 2018). We follow the methodology and implementation of Barrett et al. (2018) (see their Sect. 3.2) to compute the detection probability pdet for a given set of parameters (m1, m2, z). The optimal S/N (for a face-on, i.e. zero inclination, directly overhead source) is computed for a single detector using the sensitivity above with GW waveforms from lalsuite (LIGO Scientific Collaboration 2018). This S/N is then convolved with the antenna pattern function distribution (Finn & Chernoff 1993) in order to efficiently sample over the four angles involved, two for the sky location and two for the source orientation, which allows us to estimate the probability of detection. In our simplification we ignored the impact of BH spin on detectability, although high χeff may slightly enhance pdet.

3. Results

We use our model to predict the distributions of the three main observables inferred from GW detections: the chirp mass, the effective inspiral spin parameter and the cosmological merger redshift (Abbott et al. 2016b). Every binary in our population synthesis model contributes to the total distributions of every observable quantity with the weight given in Eq. (12).

Our detailed binary evolution models give predictions about the spin of the second-born BH and its misalignment with the orbit. However, in order to estimate the observable χeff, we need to also have information about the spin of the first-born BH. As we already discussed earlier, here we assume that the spin of the first-born BH is very low, a1 ≃ 0. This is due to two reasons. First, while the progenitor of the BH evolves through the red supergiant phase, most of the angular momentum is transported to the outer layers of the star upon expansion (because of the assumed efficient AM transport). This depletes the angular momentum of the core and eventually, due to mass transfer or stellar winds which remove the outer layers of the star, leads to a slowly rotating naked He-star. Second, the initial orbital separation of the two stars is quite large compared to the later stage of the evolution. Thus, even if tides can efficiently synchronize the rotation of the star to the orbit, the angular frequency of the envelope is too low to efficiently spin up the core. To quantitatively check these assumptions, Qin et al. (2018) used detailed stellar-evolution simulations to show that main sequence stars with initial angular rotations up to Ωinitial ≲ 0.5 Ωcritical evolve to yield BHs with negligible spins at all metallicities, even when assuming that the angular momentum of the core is conserved in the collapse. A small subset of the most rapidly spinning stars undergo efficient internal mixing and evolve chemically homogeneously. These stars never expand to become giant stars and hence do not evolve through the standard CE formation channel.

3.1. aLIGO O1, O2, and O3 observing runs

The first and second observing runs of aLIGO (and, for parts of it, Advanced Virgo) lasted for 4 and 9 months, respectively, resulting in a total of 166 days of data suitable for coincident analysis (Abbott et al. 2016c). Ten GW signals from BBH mergers were detected (Abbott et al. 2019a). These ten detections translate to a rate of 22 BBH mergers per year. In our model comparison to the data we only include these ten detections from the LIGO-Virgo Collaboration’s catalog, although evidence for additional signals in O1/O2 data has been presented by Zackay et al. (2019a,b) and Venumadhav et al. (2019). Our goal here is to model the combined distributions of observable quantities of the CE formation channel, with a special focus on the spins of the BHs which we investigate in detail. We intentionally picked a population synthesis model that approximately matches the observed rate for this study. Using Eq. (13), our model predicts for O1/O2 a detection rate of 27 yr−1, while for the ongoing observing run O3, we predict a detection rate of around 100 yr−1. However, the predicted event rate depends sensitively on a number of uncertain evolutionary model 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).

In Fig. 3 we show the predicted two-dimensional distributions of chirp mass, effective inspiral spin parameter and cosmological merger redshift for O1/O2 in orange and O3 in green. Lighter colors delineate larger contour levels of 68, 95 and 99%, respectively. For a comparison with the observations, we overlay the ten GW detections with their 90% credible intervals in black. These detections agree visually very well with our model prediction. In the first histogram, Mchirp vs. χeff, we see that the selection effects of the detectors at different sensitivities do not significantly affect the chirp mass and the effective inspiral spin parameter distributions. Meanwhile in the other two histograms, zmerger vs. χeff and zmerger vs. Mchirp, we see that at higher detector sensitivity we are able to detect events at higher cosmological redshift, namely at further distances, and that more massive sources can be observed out to a higher redshift, as one might trivially expect.

thumbnail Fig. 3.

Model predictions for binary black hole observables: chirp mass Mchirp, effective inspiral spin parameter χeff and cosmological redshift of merger zmerger distributions. We represent O1/O2 observing runs in orange and O3 in green, while lighter colors represent larger contour levels of 68, 95 and 99%, respectively. We overlaid in black the O1/O2 LIGO-Virgo Collaboration (LVC) data with their 90% credible intervals.

To provide a qualitative comparison between our theoretical predictions of χeff and LIGO-Virgo Collaboration (LVC) measurements, we conducted a visual cumulative distribution function (cdf) graphical check and a Bayesian model comparison between our model and the LVC prior (uniform spin magnitudes a1, 2 ∈ [0, 1] and isotropic spin directions). For the graphical check, shown in Fig. 4, we generate 5000 sets of 10 mock events from a KDE of our O1/O2 model predictions (solid orange line in Fig. 4) and compute the corresponding cdf. In the left panel of Fig. 4 we plot the median cdf (dashed orange line) and 90% credible interval (shaded orange region) of these sets of mock observations without any measurement uncertainty. Our model mostly predicts positive χeff and only a few slightly negative χeff but cannot reproduce χeff ≪ 0. In the right panel of Fig. 4 we plot the same quantities, but this time we add uncertainties to the 5000 sets generated from our model. These uncertainties are generated from the zero-centered LVC likelihoods. We show the cdfs from 5000 sets of 10 samples, one each from the LVC data likelihoods of the 10 observed events, in gray in both panels of Fig. 4. These likelihoods are obtained by weighting the posteriors of the ten O1/O2 GW observations by the inverse of the average projected LVC prior, which is found by combining the samples of all ten priors. We see that, when accounting for observational uncertainties, we can also reproduce negative values of χeff as in the observed cdf tail. We conclude that, visually, our model agrees well with the data. Of course, there is no single statistical check to unambiguously test the goodness of a model. In addition to the graphical check described above, we could ask, for example, whether the existing observations are statistically consistent with a model that predicts a negligible number of events with negative χeff. They are indeed consistent. All individual observations allow for zero or positive values of χeff well within their 90% credible intervals. However, if future GW observations have significant negative χeff inconsistent with zero, this will be an indication that those systems have been formed through alternative channels, such as dynamical formation. In Fig. 5 we show the probability density function of χeff as predicted by our O1/O2 model, in orange, and the average projection of the LVC prior onto χeff. For reference we added the ten LVC GW observations with their 90% credible intervals at an arbitrary vertical position. We carried out a Bayesian model comparison between our model and the LVC prior given the observational data. This test results in a Bayes factor of 15.7 which favors our model.

thumbnail Fig. 4.

Cumulative distribution functions (cdfs) of the effective inspiral spin parameter χeff as predicted for O1/O2 by our model (solid orange line). Left: from our model we generate 5000 sets of 10 measurements and plot the median cdf (dashed orange line) together with the 90% credible interval (shaded orange region). Right: we generate again 5000 sets of 10 measurements from our model to which we now add mock measurement uncertainties generated from the zero-centered LVC data likelihoods. We plot the median cdf (dashed orange line) and the 90% credible interval (shaded orange region). For both graphs, in gray, we plot the median cdf (dashed line) and the 90% credible interval (shaded region) of 5000 samples from the 10 actual observations generated from their respective LVC data likelihoods.

thumbnail Fig. 5.

Probability density function of χeff predicted by our O1/O2 model (orange) and the average LIGO prior (blue). At an arbitrary vertical position, we plot in black the LVC data with their respective 90% credible intervals. The Bayes factor between our model and the LIGO prior is 15.7.

3.2. aLIGO design sensitivity

We use the target design sensitivity curve of Abbott et al. (2018) and predict a BBH merger rate of around 400 yr−1 for advanced detectors operating at design sensitivity. In Fig. 6, we show the expected properties of the BBH population detectable at aLIGO design sensitivity: the effective inspiral spin parameter, chirp mass and cosmological merger redshift, as well as the BBH merger timescale, metallicity and binary mass ratio.

thumbnail Fig. 6.

Model predictions for the chirp mass Mchirp, effective inspiral spin parameter χeff, cosmological redshift of merger zmerger, BBH merger timescale Tmerger, metallicity Z and the binary mass ratio q distributions of binary black holes observables at design sensitivity. Lighter colors represent larger contour levels of 68, 95 and 99%. We arbitrarily divide the two dimensional histogram Mchirp vs. χeff with red lines into three regions at χeff = 0.1 and Mchirp = 15 M⊙. Region-A contains 80% of the events, 10% are in Region-B and 10% are in Region-C. For illustrative purposes, all histograms are plotted with a smoothing scale of 0.8 bins with the exception of zformation vs. χeff, zformation vs. zmerger, q vs. χeff that have no smoothing.

As pointed out previously, the two-dimensional distribution of the effective inspiral spin parameter vs. the chirp mass does not change with different detector sensitivities. We arbitrarily divide the parameter space and identify three regions: Region-A with χeff <  0.1, Region-B with χeff ≥ 0.1 and Mchirp <  15 M⊙ and Region-C with χeff ≥ 0.1 and Mchirp ≥ 15 M⊙. Our model predicts that 80, 10 and 10% of detectable BBH mergers fall into these three regions, respectively.

To understand these different regions of the parameter space, we recall that there is an anti-correlation between the merging timescale of the BBHs Tmerger and the spin of the second-born BH a2 or equivalently the observed quantity χeff (see Sect. 2.4). Region-A contains systems with low spins which translate into long merger timescales. For reference, at formation of the second-born BH most of the orbital periods are between 1 and 5 days. These systems might have formed at redshifts up to 10 and they probe a wide range of metallicities and chirp masses. BBHs in Region-B have short merger timescales since they have a high χeff: at the formation of the second-born BH they are in close orbits with periods smaller than 1 day. These BBHs are formed in the local Universe, at redshifts close to zero, where the average metallicity is high and the chirp mass is low because high metallicity massive stars tend to lose a lot of mass due to stellar winds and thus the resulting BH masses are lower. Finally, systems in Region-C have high spins and high chirp masses. These are binaries that formed with low metallicity and merge quickly, that is at zmerger ≃ zformation. This part of the parameter space really probes the low-end tail of the metallicity distribution out to the observational redshift horizon of aLIGO. These are intrinsically rare systems but are amplified by aLIGO’s higher sensitivity to high BH masses.

In Fig. 7 we further investigate GW selection effects that favor high BH masses. We show, in blue, the model prediction for aLIGO at design sensitivity against the overall BBH underlying population that one would observe with an infinitely-sensitive detector, in gray. The entire population of merging BBHs has a peak in the merger redshift at around zmerger ≃ 2. While aLIGO is not sensitive to mergers at such high redshifts, future third generation GW detectors, such as the Einstein Telescope (Punturo et al. 2010; Kalogera et al. 2019), are able to observe them. The selection effects in favor of higher BH masses are clearly visible in the distribution of Mchirp. Our observable predictions show a bimodal distribution of chirp masses with peaks at 11 M⊙ and 33 M⊙, while the underlying population has only one peak at around 10 M⊙. Selection effects allow us to observe massive BHs formed at high redshifts where the mean metallicity is lower than today. GW detectors preferentially observe BBHs that do not merge quickly, namely have wider orbits and slower spin (see the peak at χeff ≃ 0 in the χeff histogram). We note that in our treatment of the selection effects we did not take into account the potentially greater sensitivity to GWs from BBHs with high χeff; this may influence the tail of the effective inspiral spin parameter distribution, accentuating the second peak at χeff = 0.4.

thumbnail Fig. 7.

Model predictions for the chirp mass Mchirp, effective inspiral spin parameter χeff, cosmological redshift of merger zmerger, BBH merger timescale Tmerger, metallicity Z and the binary mass ratio q distributions of binary black holes observables at design sensitivity, in blue, vs. the underlying population of merging BBHs one would observe with a GW detector with infinite sensitivity, in gray. Lighter shades represent larger contour levels of 68, 95 and 99%, respectively.

4. Discussion

4.1. Angular momentum efficiency

Our results are obtained assuming efficient angular momentum transport (Spruit 1999, 2002; Fuller et al. 2019) which plays an important role in determining the spin of the first-born BH. Meanwhile, the spin of the second-born BH is mainly determined by the combined effects of stellar wind and tidal interaction during the binary evolution. 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 sub-giants and red giants (Gehan et al. 2018), which would require an even higher efficiency in angular momentum transport. Alternatively, a model with inefficient angular momentum transport predicts highly spinning BHs, χeff ≃ 1 (Belczynski et al. 2017; Arca Sedda & Benacquista 2019), which do not match current GW observations. To test that angular momentum transport efficiency affects only the spin of the first born BH, and perhaps the initial rotation of the helium star after the common envelope, but not how tides can spin up an initially non-spinning helium star, we ran a grid of 5000 MESA simulations of BH–He-star binaries with inefficient angular momentum transport, namely without the Tayler–Spruit dynamo. We found that there is a negligible impact on the spin of the second-born BH.

4.2. Comparison with other studies

When comparing our results with other studies of the CE channel, we find some discrepancies. For example Zaldarriaga et al. (2018) found a bimodal distribution of spins of the second-born BH, with around half of the BHs having spin zero and the other half maximally spinning. When accounting for stellar winds and tidal interaction between the BH and He-star in a detailed binary evolution calculation we find that this bimodal distribution is an oversimplification. In Fig. 8 we show the normalized distribution of the spin of the second-born BH a2 from detailed BH–He-star binaries simulations with masses MBH = 30 M⊙, MHe-star = 35 M⊙ and metallicities Z⊙, 10−1 Z⊙, 10−2 Z⊙ assuming a uniform distribution in log-orbital separation. Indeed, we find an approximately bimodal distribution of spins (similar to Fig. 4 from Zaldarriaga et al. 2018) at low metallicity where stellar winds are weak and do not affect the orbital evolution. However, at higher metallicity the wind-driven mass loss causes the binaries to widen and the tidal interaction which spins up the He-star gradually becomes less effective. In the same figure we also show our predicted a2 distribution from our model of BBH mergers observable at the aLIGO design sensitivity (blue shaded region) vs. the underlying population of merging BBHs (gray shaded region). The former one has a peak at around a2 = 0 while the latter has a peak at around a2 = 0.2. Although it has a peak at a2 = 0, it is much flatter and both do not show the second peak at high spins. A key reason is that we did not assume a log-normal distribution of orbital separation after the CE phase, as was assumed by Zaldarriaga et al. (2018), but used the predictions of our population synthesis model (cf. Fig 1). Moreover, we also take into account the redshift- and metallicity-dependent star formation rate and apply aLIGO selection effects. The anti-correlation between Tmerger and a2 means that binary systems with high spins merge quickly. Thus, their merger redshift distribution follows the SFR evolution with a peak around z ≃ 2, beyond aLIGO sensitivity. This further reduces the number of rapidly spinning BHs detectable by aLIGO. We believe that the non-inclusion of detailed binary evolution calculations that carefully track the angular momentum evolution due to tidal interaction and stellar winds is also the reason for overestimated χeff distributions derived in other studies (Gerosa et al. 2018; Postnov & Kuranov 2019).

thumbnail Fig. 8.

Normalized distribution of the spin of the second-born BH as predicted by our model for GWs observable at aLIGO design sensitivity, in light blue, vs. the underlying population of merging BBHs, in gray. For comparison, we show the same distribution obtained with MESA simulations of BH–He-star binary systems with MBH = 30 M⊙, MHe-star = 35 M⊙ and metallicities Z⊙ (orange), 10−1 Z⊙ (green) and 10−2 Z⊙ (red) under the assumption that the initial separation is log-uniform in 5 <  A/R⊙ <  63.

4.3. Comparison with other formation scenarios

We now compare our findings with theoretical results from other formation scenarios. All isolated binary evolution channels produce BH spins which are expected to be mostly aligned. It was shown by Marchant et al. (2016) that the chemically homogeneous evolution channel mostly produce highly spinning BHs (a1, 2 >  0.4), which are currently not observed. An indicator that would rule in favor of this scenario are the detections of high effective spin parameters, say χeff >  0.8, which are not predicted by the CE channel. Zackay et al. (2019b) recently reported the finding of a new BBH merger by reanalyzing the publicly available raw data from O1, using an independently developed pipeline. Assuming a flat χeff prior they found an event with high chirp mass, M chirp = 31 − 3 + 2 M ⊙ $ {M}_{\mathrm{chirp}}=31^{+2}_{-3}\,M_\odot $, and high effective spin parameter, χ eff = 0 . 81 − 0.21 + 0.15 $ \chi_{\mathrm{eff}}=0.81^{+0.15}_{-0.21} $, which is marginally consistent with our model. If similar events with better-measured parameters are found in the future, such BBHs would probably have been formed through the chemically homogeneous evolution channel, since the high χeff is outside our model prediction range for the CE channel. For dynamically-formed BBHs in dense star cluster, Rodriguez et al. (2018) found a symmetric distribution of χeff with a peak at χeff = 0, regardless of the BH birth spins. Therefore, anti-aligned systems (χeff <  0) would be a key indicator of the dynamical formation channel, as these are not predicted for either the CE channel or the chemically homogeneous evolution channel.

4.4. Uncertainties

Our model may be limited by some uncertainties which can alter the merger rate and, to a lesser degree, the predicted BBH property distributions. Uncertainties in (i) how the CE phase is accounted for, namely different choices of the αCE parameter which characterizes the efficiency of transferring orbital energy into unbinding the CE might lead to different rate predictions (see e.g., Giacobbo & Mapelli 2018). For example, Fragos et al. (2019) report a very efficient (high) αCE for a specific binary neutron stars (BNSs) system analyzed with 1D hydrodynamic simulations. However, in their estimate of αCE they do not include the envelope’s thermal energy in the calculation of the envelope’s binding energy, and, most importantly, their results may not carry over to BBH formation, which tends to happen at lower metallicities, with more similar donor and accretor masses at common envelope onset. Indeed, Mapelli & Giacobbo (2018) showed that uncertainties in αCE correspond to a variation of a factor of around 1.5 in the BBH merger rate estimates while a factor of 10 for BNS merger rates. Another example are uncertainties in the (ii) physics of the supernova explosions, such as the kicks strength, which can influence rates and affect the parameter distribution of BBH mergers (Dominik et al. 2013). To test how the delayed SN mechanism affects our results, we relaxed the model to account for direct collapse of the second born BH (still assuming a possible 10% of mass defect). We found similar distributions to the one of the delayed collapse with only a slight increase around Mchirp ≃ 11 M⊙ and χeff ≃ 0.2 for detectable binaries and an increase of the merger rate by ∼10% due to the survival of BBHs disrupted by natal kicks in the Fryer et al. (2012) model and the detectability of the slightly more massive BBHs are greater distances. Furthermore, in our binary population, (iii) the He-stars after the CE phase are not necessarily zero-age helium main sequence stars, as assumed in the second step of our detailed binary evolution calculations. This is because some of the progenitors of the helium stars overflowed their Roche lobes and entered the common envelope phase after helium burning was initiated in their core. The remaining lifetime of these stars is shorter than the duration of their zero-age helium main sequence. They lose less mass through stellar winds in their remaining lifetimes, which cause the orbits to not widen as much and result in more massive BHs with higher spins. However, we expect that the fraction of stars that enter the CE while burning helium in their core is higher at low metallicities, as low-metallicity stars tend to expand later in their lives. 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. Moreover, our detection rate calculation is affected by (iv) uncertainties in the redshift dependent SFR, (v) redshift-dependent metallicity distribution and (vi) the initial mass function (de Mink & Belczynski 2015; Chruslinska et al. 2019; Neijssel et al. 2019) which may not be universal (e.g., Kroupa 2001; Schneider et al. 2018), but see Farr & Mandel (2018). For example there is an uncertainty in the UV and IR data used to infer the SFR and mean metallicity at high redshift z >  4 (Madau & Dickinson 2014). A SFR favoring lower formation metallicities than the one assumed here would skew our result in favor of systems in Region-C in the χeff − Mchirp histogram of Fig. 6, namely BBHs with high Mchirp and high χeff, at the expense of systems in Region-A, namely BBHs with low χeff.

5. Conclusions

One of the biggest open questions in GW astrophysics today is how merging binary black holes are formed. Isolated binaries that go through the CE phase are one of the main proposed formation channels for BBHs. In this work we investigate the combined distributions of masses, spins and merger redshifts of a population of BBHs formed through this channel that would be detectable by advanced GW detectors. We combine binary population synthesis studies with detailed stellar structure and binary evolution simulations. Rapid population synthesis allows us to obtain a population of BBH progenitors: BH–He-star binaries. Meanwhile, the detailed simulations that take into account the effects of differential stellar rotation, tidal interactions, wind mass loss and the evolution of the structure of the He-star, allow us to accurately predict the distribution of the properties of BBH systems at their formation. We then take into account the redshift and metallicity dependence of the star-formation rate together with the selection effects of the detectors to build a model capable of reproducing all observable properties of the current sample of 10 BBH mergers. We also predict what future GW experiments are likely to observe. Our main findings can be summarized in the following points:

  • Our model is the first one to use detailed stellar structure and binary evolution simulations to successfully reproduce the observed χeff population: most with χeff ≃ 0 and a few with positive χeff. Hence, it provides strong support for the CE channel as the dominant formation channel for the observed BBH mergers.

  • We find that the ten O1/O2 GW detections are consistent with having formed through the CE channel. We predict a detection rate of 27 yr−1 for a particular set of population-synthesis model assumptions and a specific choice of a metallicity-specific star formation history, which is consistent with the 10 GW detections found in 167 days of total coincident observing time during the first two advanced detector observing runs.

  • We predict the combined distributions of Mchirp, χeff and zmerger for the current O3 observing run and for future data at design sensitivity.

  • We distinguish three different regions of observable BBH mergers. At design sensitivity, we expect around 80% of events with χeff <  0.1 and a wide range of chirp masses: these systems formed in relatively wide orbits (mostly with periods of 1–5 days at the formation of the second-born BH) and might have formed at redshifts up to 10, probing a wide range of metallicities. Around 10% of events with χeff ≥ 0.1 and Mchirp <  15 M⊙ are BBHs born in close orbits (with orbital periods of less than 1 day at the formation of the second-born BH) in the local Universe at redshift close to 0 where the metallicity is high. These systems merge promptly. Finally, the remaining 10% of events with χeff ≥ 0.1 and Mchirp ≥ 15 M⊙ are BBHs formed at low metallicity at a range of redshifts; these systems again merge promptly. Efficient spin-up of the secondary, which yields high χeff, requires the BBH to be born in a close orbit, which then allows for a prompt merger through GW emission.

  • We find that the total population of merging BBHs, namely the one that would be observed by a GW detector with infinite sensitivity, has a peak in the merger redshift at around 2, far beyond aLIGO sensitivity. This peak is set by a combination of the star formation rate, which peaks at a redshift of 2; the metallicity of star formation, which is lower in the early Universe and favors more efficient BBH formation; and the delay time distribution until merger.

  • We show that in order to understand the distribution of BBH spins, population synthesis studies of isolated field binary formation channels should include detailed binary evolution calculations that carefully track the angular momentum evolution due to the tidal interaction and stellar winds which are the origin of the spin of the second-born BH. Moreover, we find that the assumption of efficient angular momentum transport has a negligible impact on the spin of the second-born BH.

1

The detailed list of parameters used for the simulations can be found at http://mesastar.org/results

Acknowledgments

We would like to thank Pablo Marchant, Michael Zevin, Vicky Kalogera and Christopher Berry for helpful suggestions regarding Sect. 2.4 and Appendix B. 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). This work was performed in part at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611. The computations were performed in part at the University of Geneva on the Baobab and Lesta computer clusters and at the University of Birmingham. AB is funded by the program Cátedras CONACYT para Jóvenes Investigadores and by the Danish National Research Foundation (project number DNRF132). SS is supported by the Australian Research Council Centre of Excellence for Gravitational Wave Discovery (OzGrav), through project number CE170100004. All figures were made with the free Python modules Matplotlib (Hunter 2007) and pygtc (Bocquet & Carter 2016). This research made use of Astropy (http://www.astropy.org) a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018).

References

  1. Aasi, J., Abbott, B. P., Abbott, R., et al. 2015, Class. Quant. Grav., 32, 074001 [CrossRef] [Google Scholar]
  2. Aasi, J., Abbott, B. P., Abbott, R., et al. 2016, Phys. Rev. D, 93, 042006 [NASA ADS] [CrossRef] [Google Scholar]
  3. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016a, ApJ, 818, L22 [Google Scholar]
  4. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016b, Phys. Rev. Lett., 116, 241102 [NASA ADS] [CrossRef] [Google Scholar]
  5. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016c, Phys. Rev. X, 6, 041015 [Google Scholar]
  6. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, Phys. Rev. Lett., 118, 221101 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  7. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017b, ApJ, 850, L40 [NASA ADS] [CrossRef] [Google Scholar]
  8. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2018, Liv. Rev. Relativ., 21, 3 [Google Scholar]
  9. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019a, Phys. Rev. X, 9, 031040 [Google Scholar]
  10. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019b, ApJ, 882, L24 [NASA ADS] [CrossRef] [Google Scholar]
  11. Abt, H. A. 1983, ARA&A, 21, 343 [Google Scholar]
  12. Acernese, F., Agathos, M., Agatsuma, K., et al. 2015, Class. Quant. Grav., 32, 024001 [Google Scholar]
  13. Antonini, F., Chatterjee, S., Rodriguez, C. L., et al. 2016, ApJ, 816, 65 [NASA ADS] [CrossRef] [Google Scholar]
  14. Arca Sedda, M., & Benacquista, M. 2019, MNRAS, 482, 2991 [NASA ADS] [Google Scholar]
  15. Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
  17. Baldry, I. K., Driver, S. P., Loveday, J., et al. 2012, MNRAS, 421, 621 [NASA ADS] [Google Scholar]
  18. Bardeen, J. M. 1970, Nature, 226, 64 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  19. Bardeen, J. M., Press, W. H., & Teukolsky, S. A. 1972, ApJ, 178, 347 [NASA ADS] [CrossRef] [Google Scholar]
  20. Barkat, Z., Rakavy, G., & Sack, N. 1967, Phys. Rev. Lett., 18, 379 [NASA ADS] [CrossRef] [Google Scholar]
  21. Barrett, J. W., Mandel, I., Neijssel, C. J., Stevenson, S., & Vigna-Gómez, A. 2017, in Astroinformatics, eds. M. Brescia, S. G. Djorgovski, E. D. Feigelson, G. Longo, & S. Cavuoti, IAU Symp., 325, 46 [NASA ADS] [Google Scholar]
  22. Barrett, J. W., Gaebel, S. M., Neijssel, C. J., et al. 2018, MNRAS, 477, 4685 [NASA ADS] [CrossRef] [Google Scholar]
  23. Batta, A., & Ramirez-Ruiz, E. 2019, ArXiv e-prints [arXiv:1904.04835] [Google Scholar]
  24. Belczynski, K., Taam, R. E., Kalogera, V., Rasio, F. A., & Bulik, T. 2007, ApJ, 662, 504 [NASA ADS] [CrossRef] [Google Scholar]
  25. Belczynski, K., Kalogera, V., Rasio, F. A., et al. 2008, ApJS, 174, 223 [NASA ADS] [CrossRef] [Google Scholar]
  26. Belczynski, K., Bulik, T., Fryer, C. L., et al. 2010, ApJ, 714, 1217 [NASA ADS] [CrossRef] [Google Scholar]
  27. Belczynski, K., Holz, D. E., Bulik, T., & O’Shaughnessy, R. 2016a, Nature, 534, 512 [NASA ADS] [CrossRef] [Google Scholar]
  28. Belczynski, K., Repetto, S., Holz, D. E., et al. 2016b, ApJ, 819, 108 [NASA ADS] [CrossRef] [Google Scholar]
  29. Belczynski, K., Klencki, J., Meynet, G., et al. 2017, ArXiv e-prints [arXiv:1706.07053] [Google Scholar]
  30. Bisnovatyi-Kogan, G. S. 1993, Astron. Astrophys. Trans., 3, 287 [NASA ADS] [CrossRef] [Google Scholar]
  31. Bocquet, S., & Carter, F. W. 2016, J. Open Source Softw., 1 [Google Scholar]
  32. Bowler, R. A. A., Dunlop, J. S., McLure, R. J., et al. 2015, MNRAS, 452, 1817 [NASA ADS] [CrossRef] [Google Scholar]
  33. Broekgaarden, F. S., Justham, S., de Mink, S. E., et al. 2019, MNRAS, 490, 5228 [NASA ADS] [CrossRef] [Google Scholar]
  34. Burrows, A., & Hayes, J. 1996, Phys. Rev. Lett., 76, 352 [NASA ADS] [CrossRef] [Google Scholar]
  35. Cantiello, M., Mankovich, C., Bildsten, L., Christensen-Dalsgaard, J., & Paxton, B. 2014, ApJ, 788, 93 [Google Scholar]
  36. Chaboyer, B., & Zahn, J. P. 1992, A&A, 253, 173 [NASA ADS] [Google Scholar]
  37. Chruslinska, M., Nelemans, G., & Belczynski, K. 2019, MNRAS, 482, 5012 [NASA ADS] [CrossRef] [Google Scholar]
  38. Cutler, C., & Flanagan, É. E. 1994, Phys. Rev. D, 49, 2658 [NASA ADS] [CrossRef] [Google Scholar]
  39. de Kool, M. 1990, ApJ, 358, 189 [NASA ADS] [CrossRef] [Google Scholar]
  40. de Mink, S. E., & Belczynski, K. 2015, ApJ, 814, 58 [NASA ADS] [CrossRef] [Google Scholar]
  41. de Mink, S. E., Cantiello, M., Langer, N., et al. 2009, A&A, 497, 243 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Dewi, J. D. M., & Tauris, T. M. 2000, A&A, 360, 1043 [NASA ADS] [Google Scholar]
  43. Dominik, M., Belczynski, K., Fryer, C., et al. 2012, ApJ, 759, 52 [NASA ADS] [CrossRef] [Google Scholar]
  44. Dominik, M., Belczynski, K., Fryer, C., et al. 2013, ApJ, 779, 72 [NASA ADS] [CrossRef] [Google Scholar]
  45. Dominik, M., Berti, E., O’Shaughnessy, R., et al. 2015, ApJ, 806, 263 [NASA ADS] [CrossRef] [Google Scholar]
  46. Farr, W. M., & Mandel, I. 2018, Science, aat361, 6506 [NASA ADS] [CrossRef] [Google Scholar]
  47. Farr, W. M., Stevenson, S., Miller, M. C., et al. 2017, Nature, 548, 426 [NASA ADS] [CrossRef] [Google Scholar]
  48. Figer, D. F. 2005, Nature, 434, 192 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  49. Finkelstein, S. L., Ryan, Jr., R. E., Papovich, C., et al. 2015, ApJ, 810, 71 [NASA ADS] [CrossRef] [Google Scholar]
  50. Finn, L. S., & Chernoff, D. F. 1993, Phys. Rev. D, 47, 2198 [NASA ADS] [CrossRef] [Google Scholar]
  51. Fowler, W. A., & Hoyle, F. 1964, ApJS, 9, 201 [Google Scholar]
  52. Fragione, G., Ginsburg, I., & Kocsis, B. 2018, ApJ, 856, 92 [NASA ADS] [CrossRef] [Google Scholar]
  53. Fragos, T., Andrews, J. J., Ramirez-Ruiz, E., et al. 2019, ApJ, 883, L45 [NASA ADS] [CrossRef] [Google Scholar]
  54. Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, ApJ, 749, 91 [NASA ADS] [CrossRef] [Google Scholar]
  55. Fuller, J., & Ma, L. 2019, ApJ, 881, L1 [NASA ADS] [CrossRef] [Google Scholar]
  56. Fuller, J., Lecoanet, D., Cantiello, M., & Brown, B. 2014, ApJ, 796, 17 [Google Scholar]
  57. Fuller, J., Piro, A. L., & Jermyn, A. S. 2019, MNRAS, 485, 3661 [NASA ADS] [Google Scholar]
  58. Gehan, C., Mosser, B., Michel, E., Samadi, R., & Kallinger, T. 2018, A&A, 616, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Gerosa, D., Berti, E., O’Shaughnessy, R., et al. 2018, Phys. Rev. D, 98, 084036 [NASA ADS] [CrossRef] [Google Scholar]
  60. Giacobbo, N., & Mapelli, M. 2018, MNRAS, 480, 2011 [NASA ADS] [CrossRef] [Google Scholar]
  61. Giacobbo, N., Mapelli, M., & Spera, M. 2018, MNRAS, 474, 2959 [NASA ADS] [CrossRef] [Google Scholar]
  62. Grazian, A., Fontana, A., Santini, P., et al. 2015, A&A, 575, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Grevesse, N., & Sauval, A. J. 1998, Space Sci. Rev., 85, 161 [NASA ADS] [CrossRef] [Google Scholar]
  64. Hamann, W. R., & Koesterke, L. 1998, A&A, 335, 1003 [NASA ADS] [Google Scholar]
  65. Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368 [NASA ADS] [CrossRef] [Google Scholar]
  66. Heger, A., Woosley, S. E., & Spruit, H. C. 2005, ApJ, 626, 350 [NASA ADS] [CrossRef] [Google Scholar]
  67. Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974 [NASA ADS] [CrossRef] [Google Scholar]
  68. Hotokezaka, K., & Piran, T. 2017, ApJ, 842, 111 [NASA ADS] [CrossRef] [Google Scholar]
  69. Humphreys, R. M., & Davidson, K. 1994, PASP, 106, 1025 [NASA ADS] [CrossRef] [Google Scholar]
  70. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [NASA ADS] [CrossRef] [Google Scholar]
  71. Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543 [NASA ADS] [CrossRef] [Google Scholar]
  72. Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897 [NASA ADS] [CrossRef] [Google Scholar]
  73. Hut, P. 1981, A&A, 99, 126 [NASA ADS] [Google Scholar]
  74. Ilbert, O., McCracken, H. J., Le Fèvre, O., et al. 2013, A&A, 556, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  75. Inayoshi, K., Hirai, R., Kinugawa, T., & Hotokezaka, K. 2017, MNRAS, 468, 5020 [NASA ADS] [CrossRef] [Google Scholar]
  76. Ishigaki, M., Kawamata, R., Ouchi, M., et al. 2015, ApJ, 799, 12 [NASA ADS] [CrossRef] [Google Scholar]
  77. Ivanova, N., Justham, S., Chen, X., et al. 2013, A&ARv, 21, 59 [Google Scholar]
  78. Janka, H.-T. 2013, MNRAS, 434, 1355 [Google Scholar]
  79. Janka, H. T., & Mueller, E. 1994, A&A, 290, 496 [NASA ADS] [Google Scholar]
  80. Janka, H. T., Langanke, K., Marek, A., Martínez-Pinedo, G., & Müller, B. 2007, Phys. Rep., 442, 38 [NASA ADS] [CrossRef] [Google Scholar]
  81. Kajisawa, M., Ichikawa, T., Tanaka, I., et al. 2009, ApJ, 702, 1393 [NASA ADS] [CrossRef] [Google Scholar]
  82. Kalogera, V. 1996, ApJ, 471, 352 [NASA ADS] [CrossRef] [Google Scholar]
  83. Kalogera, V., Belczynski, K., Kim, C., O’Shaughnessy, R., & Willems, B. 2007, Phys. Rep., 442, 75 [NASA ADS] [CrossRef] [Google Scholar]
  84. Kalogera, V., Berry, C. P., Colpi, M., et al. 2019, Bull. Am. Astron. Soc., 51, 242 [Google Scholar]
  85. Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
  86. Kruckow, M. U., Tauris, T. M., Langer, N., Kramer, M., & Izzard, R. G. 2018, MNRAS, 481, 1908 [NASA ADS] [CrossRef] [Google Scholar]
  87. Kushnir, D., Zaldarriaga, M., Kollmeier, J. A., & Waldman, R. 2016, MNRAS, 462, 844 [NASA ADS] [CrossRef] [Google Scholar]
  88. Lee, K.-S., Ferguson, H. C., Wiklind, T., et al. 2012, ApJ, 752, 66 [NASA ADS] [CrossRef] [Google Scholar]
  89. LIGO Scientific Collaboration 2018, LIGO Algorithm Library – LALSuite, Free Software (GPL) [Google Scholar]
  90. Linden, T., Kalogera, V., Sepinsky, J. F., et al. 2010, ApJ, 725, 1984 [NASA ADS] [CrossRef] [Google Scholar]
  91. Livio, M., & Soker, N. 1988, ApJ, 329, 764 [NASA ADS] [CrossRef] [Google Scholar]
  92. Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415 [NASA ADS] [CrossRef] [Google Scholar]
  93. Madau, P., & Fragos, T. 2017, ApJ, 840, 39 [NASA ADS] [CrossRef] [Google Scholar]
  94. Mandel, I., & de Mink, S. E. 2016, MNRAS, 458, 2634 [NASA ADS] [CrossRef] [Google Scholar]
  95. Mandel, I., & Farmer, A. 2018, ArXiv e-prints [arXiv:1806.05820] [Google Scholar]
  96. Mapelli, M., & Giacobbo, N. 2018, MNRAS, 479, 4391 [NASA ADS] [CrossRef] [Google Scholar]
  97. Marchant, P., Langer, N., Podsiadlowski, P., Tauris, T. M., & Moriya, T. J. 2016, A&A, 588, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  98. Marchant, P., Renzo, M., Farmer, R., et al. 2019, ApJ, 882, 36 [Google Scholar]
  99. McLeod, D. J., McLure, R. J., Dunlop, J. S., et al. 2015, MNRAS, 450, 3032 [NASA ADS] [CrossRef] [Google Scholar]
  100. McLeod, D. J., McLure, R. J., & Dunlop, J. S. 2016, MNRAS, 459, 3812 [NASA ADS] [CrossRef] [Google Scholar]
  101. Miller, M. C. 2016, Gen. Relativ. Grav., 48, 95 [NASA ADS] [CrossRef] [Google Scholar]
  102. Miller, M. C., & Lauburg, V. M. 2009, ApJ, 692, 917 [NASA ADS] [CrossRef] [Google Scholar]
  103. Neijssel, C. J., Vigna-Gómez, A., Stevenson, S., et al. 2019, MNRAS, 490, 3740 [NASA ADS] [CrossRef] [Google Scholar]
  104. Oesch, P. A., Bouwens, R. J., Illingworth, G. D., et al. 2015, ApJ, 808, 104 [NASA ADS] [CrossRef] [Google Scholar]
  105. Paczynski, B. 1976, in Structure and Evolution of Close Binary Systems, eds. P. Eggleton, S. Mitton, & J. Whelan, IAU Symp., 73, 75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  106. Park, D., Kim, C., Lee, H. M., Bae, Y.-B., & Belczynski, K. 2017, MNRAS, 469, 4665 [NASA ADS] [CrossRef] [Google Scholar]
  107. Pavlovskii, K., Ivanova, N., Belczynski, K., & Van, K. X. 2017, MNRAS, 465, 2092 [NASA ADS] [CrossRef] [Google Scholar]
  108. Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
  109. Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [NASA ADS] [CrossRef] [Google Scholar]
  110. Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15 [Google Scholar]
  111. Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  112. Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  113. Peters, P. C. 1964, Phys. Rev., 136, B1224 [NASA ADS] [CrossRef] [Google Scholar]
  114. Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  115. Pols, O. R., Schröder, K.-P., Hurley, J. R., Tout, C. A., & Eggleton, P. P. 1998, MNRAS, 298, 525 [NASA ADS] [CrossRef] [Google Scholar]
  116. Portegies Zwart, S. F., & McMillan, S. L. W. 2000, ApJ, 528, L17 [NASA ADS] [CrossRef] [Google Scholar]
  117. Postnov, K. A., & Kuranov, A. G. 2019, MNRAS, 483, 3288 [NASA ADS] [CrossRef] [Google Scholar]
  118. Postnov, K. A., & Yungelson, L. R. 2014, Liv. Rev. Relativ., 17, 3 [Google Scholar]
  119. Punturo, M., Abernathy, M., Acernese, F., et al. 2010, Class. Quant. Grav., 27, 194002 [Google Scholar]
  120. Qin, Y., Fragos, T., Meynet, G., et al. 2018, A&A, 616, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  121. Qin, Y., Marchant, P., Fragos, T., Meynet, G., & Kalogera, V. 2019, ApJ, 870, L18 [NASA ADS] [CrossRef] [Google Scholar]
  122. Rakavy, G., & Shaviv, G. 1967, ApJ, 148, 803 [NASA ADS] [CrossRef] [Google Scholar]
  123. Rodriguez, C. L., & Loeb, A. 2018, ApJ, 866, L5 [NASA ADS] [CrossRef] [Google Scholar]
  124. Rodriguez, C. L., Morscher, M., Pattabiraman, B., et al. 2015, Phys. Rev. Lett., 115, 051101 [NASA ADS] [CrossRef] [Google Scholar]
  125. Rodriguez, C. L., Zevin, M., Pankow, C., Kalogera, V., & Rasio, F. A. 2016, ApJ, 832, L2 [NASA ADS] [CrossRef] [Google Scholar]
  126. Rodriguez, C. L., Amaro-Seoane, P., Chatterjee, S., & Rasio, F. A. 2018, Phys. Rev. Lett., 120, 151101 [NASA ADS] [CrossRef] [Google Scholar]
  127. Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444 [Google Scholar]
  128. Santamaría, L., Ohme, F., Ajith, P., et al. 2010, Phys. Rev. D, 82, 064016 [NASA ADS] [CrossRef] [Google Scholar]
  129. Schneider, F. R. N., Izzard, R. G., Langer, N., & de Mink, S. E. 2015, ApJ, 805, 20 [NASA ADS] [CrossRef] [Google Scholar]
  130. Schneider, F. R. N., Sana, H., Evans, C. J., et al. 2018, Science, 359, 69 [NASA ADS] [CrossRef] [Google Scholar]
  131. Sigurdsson, S., & Hernquist, L. 1993, Nature, 364, 423 [NASA ADS] [CrossRef] [Google Scholar]
  132. Silsbee, K., & Tremaine, S. 2017, ApJ, 836, 39 [NASA ADS] [CrossRef] [Google Scholar]
  133. Smarr, L. L., & Blandford, R. 1976, ApJ, 207, 574 [NASA ADS] [CrossRef] [Google Scholar]
  134. Socrates, A., Blaes, O., Hungerford, A., & Fryer, C. L. 2005, ApJ, 632, 531 [NASA ADS] [CrossRef] [Google Scholar]
  135. Spruit, H. C. 1999, A&A, 349, 189 [NASA ADS] [Google Scholar]
  136. Spruit, H. C. 2002, A&A, 381, 923 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  137. Stevenson, S., Vigna-Gómez, A., Mandel, I., et al. 2017, Nat. Commun., 8, 14906 [NASA ADS] [CrossRef] [Google Scholar]
  138. Stevenson, S., Sampson, M., Powell, J., et al. 2019, ApJ, 882, 121 [Google Scholar]
  139. Suijs, M. P. L., Langer, N., Poelarends, A. J., et al. 2008, A&A, 481, L87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  140. Thorne, K. S. 1974, ApJ, 191, 507 [Google Scholar]
  141. Tutukov, A. V., & Yungelson, L. R. 1993, MNRAS, 260, 675 [NASA ADS] [CrossRef] [Google Scholar]
  142. van den Heuvel, E. P. J. 1976, in Structure and Evolution of Close Binary Systems, eds. P. Eggleton, S. Mitton, & J. Whelan, IAU Symp., 73, 35 [NASA ADS] [CrossRef] [Google Scholar]
  143. van den Heuvel, E. P. J., Portegies Zwart, S. F., & de Mink, S. E. 2017, MNRAS, 471, 4256 [NASA ADS] [CrossRef] [Google Scholar]
  144. Venumadhav, T., Zackay, B., Roulet, J., Dai, L., & Zaldarriaga, M. 2019, ArXiv e-prints [arXiv:1904.07214] [Google Scholar]
  145. Vigna-Gómez, A., Neijssel, C. J., Stevenson, S., et al. 2018, MNRAS, 481, 4009 [NASA ADS] [CrossRef] [Google Scholar]
  146. Vink, J. S., & de Koter, A. 2005, A&A, 442, 587 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  147. Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  148. Webbink, R. F. 1984, ApJ, 277, 355 [NASA ADS] [CrossRef] [Google Scholar]
  149. Woosley, S. E. 2017, ApJ, 836, 244 [NASA ADS] [CrossRef] [Google Scholar]
  150. Xu, X.-J., & Li, X.-D. 2010, ApJ, 716, 114 [NASA ADS] [CrossRef] [Google Scholar]
  151. Yoshida, T., Umeda, H., Maeda, K., & Ishii, T. 2016, MNRAS, 457, 351 [NASA ADS] [CrossRef] [Google Scholar]
  152. Zackay, B., Dai, L., Venumadhav, T., Roulet, J., & Zaldarriaga, M. 2019a, ArXiv e-prints [arXiv:1910.09528] [Google Scholar]
  153. Zackay, B., Venumadhav, T., Dai, L., Roulet, J., & Zaldarriaga, M. 2019b, Phys. Rev. D, 100, 023007 [Google Scholar]
  154. Zahn, J. P. 1977, A&A, 500, 121 [NASA ADS] [Google Scholar]
  155. Zaldarriaga, M., Kushnir, D., & Kollmeier, J. A. 2018, MNRAS, 473, 4174 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Mass renormalization of the population synthesis simulation

We used COMPAS to run a Monte-Carlo simulation and generate a sample of half a million BH–He-star binaries. When we perform binary population synthesis simulations, we only model binary systems, neglecting the population of single stars. Furthermore, to save on computational costs, we restrict the mass of the primary star to a suitable range mA <  mprimary <  mB so that we only consider initial binaries that can be progenitors of the systems we want to study (this is a basic version of the importance sampling approach described by Broekgaarden et al. 2019). This means that we model only a fraction of the underlying stellar population. Here we show how to renormalize the population synthesis simulation to the total stellar mass of the underlying stellar population.

Let us consider a stellar population of total mass M* with an initial mass function (IMF) of single star masses:

f ( m ) = { f 0 ( m m min ) − α 1 m min ≤ m ≤ m 1 f 0 ( m 1 m min ) − α 1 ( m m 1 ) − α 2 m 1 ≤ m ≤ m 2 f 0 ( m 1 m min ) − α 1 ( m 2 m 1 ) − α 2 ( m m 2 ) − α 3 m 2 ≤ m ≤ m max $$ \begin{aligned} f(m) = \left\{ \begin{array}{ll} f_0 \left(\frac{m}{m_\mathrm{min} }\right)^{-\alpha _1}&m_\mathrm{min} \le m \le m_1 \\ f_0 \left(\frac{m_1}{m_\mathrm{min} }\right)^{-\alpha _1}\left(\frac{m}{m_{1}}\right)^{-\alpha _2}&m_1 \le m \le m_2 \\ f_0 \left(\frac{m_1}{m_\mathrm{min} }\right)^{-\alpha _1}\left(\frac{m_2}{m_{1}}\right)^{-\alpha _2}\left(\frac{m}{m_{2}}\right)^{-\alpha _3}&m_2 \le m \le m_\mathrm{max} \end{array} \right. \end{aligned} $$(A.1)

where the constant f0 is defined such that ∫ m min m max f (m) dm=1 $ \smallint_{m_\mathrm{min}}^{m_\mathrm{max}}f(m) \, {\rm d}m = 1 $. Let fbin be the fraction of stars in binaries and assume that the distribution of mass ratios in binaries is flat between 0 and 1, that is g(q) = 1. Then, the mean mass of a stellar system in the population is

m ¯ ⋆ = ( 1 − f bin ) ∫ m min m max m f ( m ) d m + f bin ∫ m min m max ∫ 0 1 [ f ( m ) g ( q ) ( m + q m ) ] d q d m = f 0 ( 1 + f bin 2 ) [ 1 2 − α 1 ( m 1 2 − α 1 − m min 2 − α 1 m min − α 1 ) + 1 2 − α 2 ( m 1 m min ) − α 1 ( m 2 2 − α 2 − m 1 2 − α 2 m 1 − α 2 ) + 1 2 − α 3 ( m 1 m min ) − α 1 ( m 2 m 1 ) − α 2 ( m max 2 − α 3 − m 2 2 − α 3 m 2 − α 3 ) ] . $$ \begin{aligned} \begin{split} \bar{m}_{\star }&= (1-f_\mathrm{bin} )\int _{m_\mathrm{min} }^{m_\mathrm{max} } m f(m)\,\mathrm{d}m + f_\mathrm{bin} \int _{m_\mathrm{min} }^{m_\mathrm{max} }\int _{0}^{1}\left[f(m){ g}(q)\left(m+qm\right)\right] \, \mathrm{d}q\mathrm{d}m \\&= f_0 \left(1+\frac{f_\mathrm{bin} }{2}\right) \left[ \frac{1}{2-\alpha _1} \left(\frac{m_{1}^{2-\alpha _1}-m_\mathrm{min} ^{2-\alpha _1}}{m_\mathrm{min} ^{-\alpha _1}} \right) + \frac{1}{2-\alpha _2}\left(\frac{m_1}{m_\mathrm{min} }\right)^{-\alpha _1} \left(\frac{m_{2}^{2-\alpha _2}-m_{1}^{2-\alpha _2}}{m_{1}^{-\alpha _2}} \right) +\frac{1}{2-\alpha _3}\left(\frac{m_1}{m_\mathrm{min} }\right)^{-\alpha _1} \left(\frac{m_2}{m_{1}}\right)^{-\alpha _2} \left(\frac{m_\mathrm{max} ^{2-\alpha _3}-m_{2}^{2-\alpha _3}}{m_{2}^{-\alpha _3}} \right)\right] . \end{split} \end{aligned} $$(A.2)

In the case when mA, mB >  m2, we only model the following fraction of all systems:

f model = f bin ∫ m A m B f ( m ) d m = f bin f 0 1 1 − α 3 ( m 1 m min ) − α 1 ( m 2 m 1 ) − α 2 ( m B 1 − α 3 − m A 1 − α 3 m 2 − α 3 ) . $$ \begin{aligned} f_{\mathrm{model} } = f_\mathrm{bin} \int _{m_A}^{m_B}f(m)\,\mathrm{d}m = f_\mathrm{bin} f_0\frac{1}{1-\alpha _3}\left(\frac{m_1}{m_\mathrm{min} }\right)^{-\alpha _1} \left(\frac{m_2}{m_{1}}\right)^{-\alpha _2} \left(\frac{m_{B}^{1-\alpha _3}-m_{A}^{1-\alpha _3}}{m_{2}^{-\alpha _3}} \right) . \end{aligned} $$(A.3)

The mean mass of a binary system in our simulated population is;

m ¯ ∗ , model = 1 ∫ m A m B f ( m ) d m ∫ m A m B ∫ 0 1 f ( m ) g ( q ) ( m + q m ) d q d m = 3 2 1 − α 3 2 − α 3 ( m B 2 − α 3 − m A 2 − α 3 m B 1 − α 3 − m A 1 − α 3 ) . $$ \begin{aligned} \bar{m}_{*, \mathrm{model} }&= \frac{1}{\int _{m_{A}}^{m_{B}}f(m)dm} \int _{m_A}^{m_B} \int _{0}^1 f(m) { g}(q) (m+qm) \, \mathrm{d}q \mathrm{d}m = \frac{3}{2} \frac{1-\alpha _3}{2-\alpha _3} \left(\frac{m_{B}^{2-\alpha _3}-m_{A}^{2-\alpha _3}}{m_{B}^{1-\alpha _3}-m_{A}^{1-\alpha _3}} \right). \end{aligned} $$(A.4)

Thus, the total modeled mass M*, model represents only a fraction of the total stellar population mass M*:

f corr = M ∗ , model M ∗ = f model m ¯ ∗ , model m ¯ ∗ $$ \begin{aligned} f_\mathrm{corr} = \frac{M_{*,\mathrm{model} }}{M_*} = f_{\mathrm{model} } \frac{\bar{m}_{*,\mathrm{model} }}{\bar{m}_{*}} \end{aligned} $$(A.5)

and we must renormalize by the inverse of fcorr in order to return to the population we intended to simulate.

Adopting the Kroupa (2001) IMF, namely α1 = 0.3, α2 = 1.3, α3 = 2.3, m1 = 0.08 M⊙, m2 = 0.5 M⊙, using the observed fbin = 0.7 (Sana et al. 2012), arbitrarily choosing mmin = 0.01 M⊙ and mmax = 200 M⊙ as the minimum and the maximum stellar mass, and carrying out the simulation for primary masses in the range between mA = 5 M⊙ and mB = 150 M⊙ (Figer 2005), we obtain f corr −1 =4.78 $ f_\mathrm{corr}^{-1}=4.78 $.

Appendix B: Detection rate

To compute the detected BH merger rate by GW detectors, we follow a similar procedure to the one of Belczynski et al. (2016b) which is a refined version of Dominik et al. (2015). In our cosmological calculation, we adopt the flat ΛCDM model with H 0 = 67.7 km s − 1 Mpc $ H_0=67.7 \, \frac{\text{km}\,\text{s}^{-1}}{\text{ Mpc}} $ and Ωm = 0.307 (Planck Collaboration XIII 2016). We follow the model of Madau & Fragos (2017) for the star formation rate (SFR) model as a function of redshift, which is an updated version of Madau & Dickinson (2014),

SFR ( z ) = 0.01 · ( 1 + z ) 2.6 1 + ( ( 1 + z ) / 3.2 ) 6.2 M ⊙ yr − 1 Mpc − 3 . $$ \begin{aligned} \text{ SFR}(z) = \frac{0.01 \cdot (1 + z)^{2.6}}{1 + ((1 + z) / 3.2)^{6.2}} \, {M}_{\odot } \, \text{ yr}^{-1} \, \text{ Mpc}^{-3}. \end{aligned} $$(B.1)

We assume that the metallicities of the binaries follow a truncated log-normal distributed,

N ( log 10 ( Z ) | μ ( z ) , σ ) ≡ d P d log 10 ( Z ) ( z ) = 1 σ 2 π exp ( − ( log 10 ( Z ) − μ ( z ) ) 2 2 σ 2 ) $$ \begin{aligned} \mathcal{N} (\log _{10}(Z) \, | \, \mu (z), \sigma ) \equiv \frac{\mathrm{d}P}{\mathrm{d} \log _{10}(Z)}(z) = \frac{1}{\sigma \sqrt{2\pi }} \exp \left(- \frac{\left( \log _{10}(Z)-\mu (z) \right)^2}{2\sigma ^2} \right) \end{aligned} $$(B.2)

with standard deviation σ = 0.5 and mean μ = log 10 ( Z ¯ ( z ) ) − ln ( 10 ) 2 σ 2 $ \mu = \log_{10}\left(\bar{Z}(z)\right) - \frac{\ln(10)}{2} \sigma^2 $ where the mean metallicity is (Madau & Fragos 2017)

Z ¯ ( z ) = Z ⊙ × 10 0.153 − 0.074 · z 1.34 . $$ \begin{aligned} \bar{Z}(z) = Z_\odot \times 10^{0.153 - 0.074 \cdot z ^{1.34}}. \end{aligned} $$(B.3)

The log-normal distribution is truncated at the highest metallicity bin edge, Zmax = 0.034923, and the distribution is accordingly renormalized to ensure that ∫ −∞ log 10 Z max N ( log 10 (Z)|μ(z),σ)d log 10 Z=1 $ \smallint_{-\infty}^{\log_{10}Z_\mathrm{max}} \mathcal{N}(\log_{10}(Z)\, |\mu(z), \sigma)\,d\log_{10}Z = 1 $. Portions of the distribution extending beyond the lower limit edge Zmin = 0.000091 are included in the edge bin when integrating over metallicity.

We compute the detection rate by integrating the cosmological merger rate R ( z m ) = d N d m 1 d m 2 d V c d t m $ R(z_{\mathrm{m}}) = \frac{\mathrm{d}N}{\mathrm{d}m_1 \mathrm{d}m_2 \mathrm{d}V_{\mathrm{c}} \mathrm{d}t_{\mathrm{m}}} $ per unit masses, per unit comoving volume, per unit time as measured in the source frame at the redshift of the merger as in Eq. (5) of Dominik et al. (2015), that is

R det = ∫ ∫ ∫ R ( z m ) d t m d t det p det ( z m , m 1 , m 2 ) d m 1 d m 2 d V c , $$ \begin{aligned} R_{\rm det} = \int \int \int R(z_{\rm m}) \frac{\mathrm{d}t_{\rm m}}{\mathrm{d}t_{\rm det}} p_\mathrm{det} (z_{\rm m},m_{1},m_2) \, \mathrm{d}m_1 \mathrm{d}m_2 \mathrm{d}V_{\rm c} , \end{aligned} $$(B.4)

where the factor d t m d t det = 1 1 + z m $ \frac{\mathrm{d}t_{\mathrm{m}}}{\mathrm{d}t_{\mathrm{det}}}=\frac{1}{1+z_{\mathrm{m}}} $ account for the difference in clock rates at merger and at the detector and pdet is the detection probability accounting for the detector’s selection effects. The integration over the comoving volume can be calculated with a change of variable over the redshift of merger, namely d V c = d V c d z m d z m $ \mathrm{d}V_{\mathrm{c}} = \frac{\mathrm{d}V_{\mathrm{c}}}{\mathrm{d}z_{\mathrm{m}}} \mathrm{d} z_{\mathrm{m}} $, where

d V c d z m = 4 π c H 0 D c 2 ( z m ) E ( z ) , $$ \begin{aligned} \frac{\mathrm{d}V_{\rm c}}{\mathrm{d}z_{\rm m}}=\frac{4\pi c}{H_0} \frac{D^2_{\rm c}(z_{\rm m})}{E(z)}, \end{aligned} $$(B.5)

and Dc(z) is the comoving distance which is related to the luminosity distance as DL(z) = Dc(z)(1 + z) and is computed as follows

D c ( z ) = c H 0 ∫ 0 z d z ′ E ( z ′ ) $$ \begin{aligned} D_{\rm c}(z)=\frac{c}{H_0} \int _0^{z}\frac{\mathrm{d}z^{\prime} }{E(z^{\prime} )} \end{aligned} $$(B.6)

where E ( z ) = Ω m ( 1 + z ) 3 + Ω Λ $ E(z) = \sqrt{\Omega_{\mathrm{m}} (1+z)^3+\Omega_\Lambda} $. The merger rate R(zm) can be rewritten as the convolution of the star-formation rate, SFR = d 2 M f d V c d t f $ =\frac{\mathrm{d}^2 M_{\mathrm{f}}}{\mathrm{d}V_{\mathrm{c}} \mathrm{d}t_{\mathrm{f}}} $, that is the total mass of stars formed per comoving volume per year, and the number density of binaries per unit star forming mass Mf per unit masses m1, m2 per unit log-metallicity interval Z per unit time delay τ = tm − tf:

R ( z m ) ≡ R ( z ( t m ) ) = ∫ 0 t m ∫ d 2 M f d V c d t f ( z f ) d 5 N d M f d m 1 d m 2 d τ d log 10 Z ( τ ) N ( log 10 ( Z ) | μ = log 10 ( Z ¯ ( z f ) ) − ln ( 10 ) 2 σ 2 , σ = 0.5 ) d log 10 Z d t f , $$ \begin{aligned} R(z_{\rm m}) \equiv R(z(t_{\rm m})) = \int ^{t_{\rm m}}_0 \int \frac{\mathrm{d}^2 M_{\rm f}}{\mathrm{d}V_{\rm c} \mathrm{d}t_{\rm f}} (z_{\rm f}) \, \frac{\mathrm{d}^5 N}{\mathrm{d}M_{\rm f} \mathrm{d}m_1 \mathrm{d}m_2 \mathrm{d}\tau \mathrm{d}\log _{10}Z} (\tau )\, \mathcal{N} \left( \log _{10}(Z) \, \bigg | \, \mu = \log _{10} \left( \bar{Z} (z_{\mathrm{f}})\right) - \frac{\ln (10)}{2} \sigma ^2 , \sigma = 0.5 \right) \, \mathrm{d}\log _{10}Z \, \mathrm{d}t_{\rm f}, \end{aligned} $$(B.7)

where we used the compact notation zf ≡ z(tf). The time delay τ is mostly set by Tmerger, since the GW-driven merger takes much longer than stellar evolution for BH progenitors. Performing the change of variable d z m = d z m d t m d t m = H 0 ( 1 + z m ) E ( z m ) d t m $ \mathrm{d}z_{\mathrm{m}} = \frac{\mathrm{d}z_{\mathrm{m}}}{\mathrm{d}t_{\mathrm{m}}} \mathrm{d}t_{\mathrm{m}} = H_0 (1+z_{\mathrm{m}}) E(z_{\mathrm{m}}) \mathrm{d}t_{\mathrm{m}} $, the integral of Eq. (B.4) translates into the following Monte-Carlo sum over the formation time intervals arbitrarily chosen as Δti = 100 Myr and 30 uniformly-distributed log-metallicity intervals for Z ∈ [Zmin, Zmax]

R det = ∑ Δ t i ∑ Δ Z j ∑ k fSFR ( z f , i ) M sim , Δ Z j f corr 4 π c D c 2 ( z m , i , k ) p det ( z m , i , k , m 1 , k , m 2 , k ) Δ t i , $$ \begin{aligned} R_{\rm det} = \sum _{\Delta t_i} \sum _{\Delta Z_j} \sum _k \frac{\text{ fSFR}(z_{\mathrm{f},i})}{M_{\mathrm{sim} ,\Delta Z_j}} \, f_\mathrm{corr} \, 4 \pi c\, D_{\rm c}^2(z_{m,i,k})\, p_\mathrm{det} (z_{m,i,k}, m_{1,k}, m_{2,k}) \, \Delta t_i , \end{aligned} $$(B.8)

where Msim, ΔZj is the total mass simulated per log-metallicity interval ΔZj and fSFR is the total mass of stars formed per comoving volume per year per log-metallicity interval ΔZj,

fSFR ( z f , i ) = ∫ Δ Z j SFR ( z f , i ) N ( log 10 ( Z ) | μ = log 10 ( Z ¯ ( z f ) ) − ln ( 10 ) 2 σ 2 , σ = 0.5 ) d log 10 Z = SFR ( z f , i ) [ CDF ( log 10 ( Z j ) + Δ Z j 2 ) − CDF ( log 10 ( Z j ) − Δ Z j 2 ) ] M ⊙ yr − 1 Mpc − 3 , $$ \begin{aligned} \begin{split} \text{ fSFR}(z_{\mathrm{f},i})&= \int _{\Delta Z_j} \text{ SFR} \left(z_{\mathrm{f},i} \right)~\mathcal{N} \left( \log _{10}(Z) \, \bigg | \, \mu = \log _{10} \left( \bar{Z} (z_{\mathrm{f}})\right) - \frac{\ln (10)}{2} \sigma ^2 , \sigma = 0.5 \right) d\log _{10}Z \\&= \text{ SFR}(z_{\mathrm{f},i}) \left[ \mathrm{CDF}\bigg (\log _{10}(Z_j)+ \frac{\Delta Z_j}{2} \bigg )\,{-}\,\mathrm{CDF} \bigg (\log _{10}(Z_j)-\frac{\Delta Z_j}{2} \bigg ) \right] \, {M}_{\odot } \, \text{ yr}^{-1} \, \text{ Mpc}^{-3} \, , \end{split} \end{aligned} $$(B.9)

where Zj is the center of the log-metallicity bin ΔZj corresponding to the metallicity Zk of the binary k. Meanwhile, the integrated SFR (iSFR) over the cosmic time used to obtain the weighted distributions of parameters after the CE phase is computed with the change of variable d t = d t d z d z = ( H 0 ( 1 + z ) E ( z ) ) − 1 d z $ \mathrm{d}t = \frac{\mathrm{d}t}{\mathrm{d}z} \mathrm{d}z = \left(H_0(1+z)E(z)\right)^{-1} \mathrm{d}z $,

iSFR ( Z ) = ∫ 0 ∞ SFR ( z ) N ( log 10 ( Z ) | μ = log 10 ( Z ¯ ( z f ) ) − ln ( 10 ) 2 σ 2 , σ = 0.5 ) d t d z d z M ⊙ Mpc − 3 , $$ \begin{aligned} \text{ iSFR}(Z)&= \int _0^{\infty } \text{ SFR}\left(z\right)\, \mathcal{N} \left( \log _{10}(Z) \, \bigg | \, \mu \,{=}\,\log _{10} \left( \bar{Z} (z_{\mathrm{f}})\right)\,{-}\,\frac{\ln (10)}{2} \sigma ^2 , \sigma \,{=}\,0.5 \right) \frac{\mathrm{d}t}{\mathrm{d}z} \mathrm{d}z \, {M}_{\odot } \, \text{ Mpc}^{-3}\!, \end{aligned} $$(B.10)

which gives the total mass of stars formed per comoving volume at a given metallicity Z.

Appendix C: Linear interpolation of the MESA simulations

Running MESA simulations on the entire simulated binary population is too computationally expensive. Instead, we use linear interpolation over a simulated grid to estimate the physical observables of the binaries that we are interested in. To generate the first simulations we sample stochastically in the logarithmic parameter space of initial masses, orbital periods and metallicities. We generate 3000 initial points with mBH ∈ [2.5 M⊙, 60 M⊙], mHe-star ∈ [2.5 M⊙, 89 M⊙], p ∈ [0.05 days, 8.5 days] and Z ∈ [0.0001, 0.0349]. We add to these 1500 points drawn from a kernel density estimator (KDE) of the parameter distribution of the synthetic binary population.

We cover these points by running MESA binary simulations as described in Sect. 2.2. We want to interpolate six quantities: the He-star mass and its CO 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 black hole. All physical quantities are log-transformed and rescaled to the interval [−1, 1] before going through the interpolation algorithm. The interpolation itself relies on building a Delaunay triangulation of the input data points followed by barycentric linear interpolation over the vortices of the (hyper)triangle containing the location of interest. The relative error on each quantity Xi is computed as Δi = (Xtrue, i − Xinterp.,i)/Xtrue, i (in the original units, except for spins, where the relative error is computed on the spin logarithm to avoid excessive sensitivity to true values close to zero). We then combine the relative errors of the quantities to obtain the combined relative error

Δ = min ( ∑ i = 1 6 Δ i 2 , 6 ) , $$ \begin{aligned} \Delta = \min \left( \sqrt{\sum _{i=1}^6 \Delta _i^2} \, , \sqrt{6} \right), \end{aligned} $$(C.1)

where we arbitrarily limit the maximal combined relative error to 6 $ \sqrt{6} $, that is a point with all relative errors equal to |Δi| = 1.

We check the accuracy of the linear interpolator by conducting 50 leave 5% of the sample out tests. We use the combined relative errors as weights to sample an additional 500 points where the interpolator is performing the worst. We iterate this procedure 21 times for a total of 10,500 simulations, stopping because almost all the new points generated through this procedure in the 22nd iteration would be on the boundaries. The triangulation scheme can still fail near the parameter space boundaries; in this case, we find that 2.5% of synthetic population systems cannot be interpolated, and we run 3000 simulations to bring the parameter space coverage to 100%.

Figure C.1 shows the relative errors in the interpolated quantities over the series of 50 leave 5% of the sample out tests. The left panel shows the median percentage relative errors excluding non-fittable points. These stabilize at 0.01%, 0.20% and 0.04% for the He-star mass before supernova, its CO core mass, and resultant BH mass; 0.01% for the orbital period; 0.04% for the lifetime of the BH–He-star binary and 0.4% for the log-spin of the second-born BH. The log spin is the parameter which shows the biggest relative errors because it can have very large negative values for spins close to zero. The right panel of Fig. C.1 shows the fraction of relative errors larger than 10% as a function of the number of simulations for the different interpolated quantities where we also count non-fittable points, such as points at the boundary of the parameter space or isolated regions of the parameter space. At the last iteration the mean of the 50 leave 5% of the sample out tests shows the following fraction of relative errors greater than 10%: 5% for the He-star mass before the supernova, 8% for its CO core mass and the resulting BH mass, 6% for the orbital period, 6% for the lifetime of the BH–He-star binary and 17% for the log-spin of the second-born BH. The apparent increase in the fraction of relative errors larger than 10% with the number of simulations happens because with the last iterations we are sampling mostly the boundaries of the parameter space and the test picks up more points that cannot be interpolated (we note that the median relative errors does not show this trend and stabilizes). The last simulations used to bring the coverage of the parameter space to 100% are run in disconnected and remote regions of the parameter space, and the “leave out” tests pick up the newly added samples, artificially increasing the apparent fraction of large relative errors.

thumbnail Fig. C.1.

Median relative error of the interpolation expressed as a percentage (left) and the fraction of relative errors above 10% (right) from all iterations of the 50 leave 5% of the sample out tests for six interpolated quantities. The points on each plot, moving from left to right, represent the different iterations; we exclude all simulations that stopped due to initial Roche-lobe overflow (indicating a difference between COMPAS and MESA models). The right plot includes NaNs (obtained from non-fittable points, e.g., points at the boundary of the parameter space) when counting relative errors larger than 10%, while the left plot excludes them.

All Figures

thumbnail Fig. 1.

Parameter distributions of the binary population after the CE phase. These BH–He-star binaries include systems that are going to form BH–NS binaries and BBHs with GW inspiral timescales bigger than the Hubble time. We show the distributions of orbital separation A, first-born BH mass MBH, He-star mass MHe-star and metallicity Z weighted by the integrated redshift- and metallicity-dependent SFR over the cosmic time (see Eq. (B.10)). The lighter shades represent larger contour levels, 68, 95 and 99%, respectively, constructed with pygtc (Bocquet & Carter 2016).
In the text
thumbnail Fig. 2.

Combined distribution of the BBH merger timescale Tmerger vs. the effective inspiral spin parameter χeff for our synthetic BBH population at metallicity log10(Z) = −2.5, in gray. The lighter shades represent larger contour levels, 68, 95 and 99%, respectively. Tidally locked systems follow the relation log 10 ( T merger ) ∼ − 8 3 log 10 ( χ eff ) $ \log_{10} ({T}_{\mathrm{merger}}) \sim - \frac{8}{3} \log_{10}(\chi_{\mathrm{eff}}) $, orange dashed line, while the others systems follow log 10 ( T merger ) ∼ − 8 17 log 10 ( χ eff ) $ \log_{10} ({T}_{\mathrm{merger}}) \sim - \frac{8}{17} \log_{10}(\chi_{\mathrm{eff}}) $ dictated by the tidal synchronization timescale, blue dashed line. Both lines are drawn at an arbitrary ordinate.
In the text
thumbnail Fig. 3.

Model predictions for binary black hole observables: chirp mass Mchirp, effective inspiral spin parameter χeff and cosmological redshift of merger zmerger distributions. We represent O1/O2 observing runs in orange and O3 in green, while lighter colors represent larger contour levels of 68, 95 and 99%, respectively. We overlaid in black the O1/O2 LIGO-Virgo Collaboration (LVC) data with their 90% credible intervals.
In the text
thumbnail Fig. 4.

Cumulative distribution functions (cdfs) of the effective inspiral spin parameter χeff as predicted for O1/O2 by our model (solid orange line). Left: from our model we generate 5000 sets of 10 measurements and plot the median cdf (dashed orange line) together with the 90% credible interval (shaded orange region). Right: we generate again 5000 sets of 10 measurements from our model to which we now add mock measurement uncertainties generated from the zero-centered LVC data likelihoods. We plot the median cdf (dashed orange line) and the 90% credible interval (shaded orange region). For both graphs, in gray, we plot the median cdf (dashed line) and the 90% credible interval (shaded region) of 5000 samples from the 10 actual observations generated from their respective LVC data likelihoods.
In the text
thumbnail Fig. 5.

Probability density function of χeff predicted by our O1/O2 model (orange) and the average LIGO prior (blue). At an arbitrary vertical position, we plot in black the LVC data with their respective 90% credible intervals. The Bayes factor between our model and the LIGO prior is 15.7.
In the text
thumbnail Fig. 6.

Model predictions for the chirp mass Mchirp, effective inspiral spin parameter χeff, cosmological redshift of merger zmerger, BBH merger timescale Tmerger, metallicity Z and the binary mass ratio q distributions of binary black holes observables at design sensitivity. Lighter colors represent larger contour levels of 68, 95 and 99%. We arbitrarily divide the two dimensional histogram Mchirp vs. χeff with red lines into three regions at χeff = 0.1 and Mchirp = 15 M⊙. Region-A contains 80% of the events, 10% are in Region-B and 10% are in Region-C. For illustrative purposes, all histograms are plotted with a smoothing scale of 0.8 bins with the exception of zformation vs. χeff, zformation vs. zmerger, q vs. χeff that have no smoothing.
In the text
thumbnail Fig. 7.

Model predictions for the chirp mass Mchirp, effective inspiral spin parameter χeff, cosmological redshift of merger zmerger, BBH merger timescale Tmerger, metallicity Z and the binary mass ratio q distributions of binary black holes observables at design sensitivity, in blue, vs. the underlying population of merging BBHs one would observe with a GW detector with infinite sensitivity, in gray. Lighter shades represent larger contour levels of 68, 95 and 99%, respectively.
In the text
thumbnail Fig. 8.

Normalized distribution of the spin of the second-born BH as predicted by our model for GWs observable at aLIGO design sensitivity, in light blue, vs. the underlying population of merging BBHs, in gray. For comparison, we show the same distribution obtained with MESA simulations of BH–He-star binary systems with MBH = 30 M⊙, MHe-star = 35 M⊙ and metallicities Z⊙ (orange), 10−1 Z⊙ (green) and 10−2 Z⊙ (red) under the assumption that the initial separation is log-uniform in 5 <  A/R⊙ <  63.
In the text
thumbnail Fig. C.1.

Median relative error of the interpolation expressed as a percentage (left) and the fraction of relative errors above 10% (right) from all iterations of the 50 leave 5% of the sample out tests for six interpolated quantities. The points on each plot, moving from left to right, represent the different iterations; we exclude all simulations that stopped due to initial Roche-lobe overflow (indicating a difference between COMPAS and MESA models). The right plot includes NaNs (obtained from non-fittable points, e.g., points at the boundary of the parameter space) when counting relative errors larger than 10%, while the left plot excludes them.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.