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/00046361/201936204  
Published online  16 March 2020 
The origin of spin in binary black holes
Predicting the distributions of the main observables of Advanced LIGO
^{1}
Geneva Observatory, University of Geneva, Chemin des Maillettes 51, 1290 Versoix, Switzerland
email: simone.bavera@unige.ch
^{2}
Birmingham Institute for Gravitational Wave Astronomy and School of Physics and Astronomy, University of Birmingham, Birmingham B15 2TT, UK
^{3}
Monash Centre for Astrophysics, School of Physics and Astronomy, Monash University, Clayton, Victoria 3800, Australia
^{4}
ARC Centre of Excellence for Gravitational Wave Discovery – OzGrav, Australia
^{5}
Instituto Nacional de Astrofísica, Óptica y Electrónica, Tonantzintla, Puebla 72840, Mexico
^{6}
Consejo Nacional de Ciencia y Tecnología, Av. Insurgentes Sur 1582, Col. Crédito Constructor, Ciudad de México 03940, Mexico
^{7}
Niels Bohr Institute, University of Copenhagen, Blegdamsvej 17, 2100 Copenhagen, Denmark
^{8}
Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA), Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA
^{9}
Centre for Astrophysics and Supercomputing, Swinburne University of Technology, Hawthorn, VIC 3122, Australia
Received:
28
June
2019
Accepted:
21
January
2020
Context. After years of scientific progress, the origin of stellar binary black holes is still a great mystery. Several formation channels for merging black holes have been proposed in the literature. As more merger detections are expected with future gravitationalwave observations, population synthesis studies can help to distinguish between them.
Aims. We study the formation of coalescing binary black holes via the evolution of isolated field binaries that go through the common envelope phase in order to obtain the combined distributions of observables such as blackhole spins, masses and cosmological redshifts of mergers.
Methods. To achieve this aim, we used a hybrid technique that combines the parametric binary population synthesis code COMPAS with detailed binary evolution simulations performed with the MESA code. We then convolved our binary evolution calculations with the redshift and metallicitydependent starformation rate and the selection effects of gravitationalwave detectors to obtain predictions of observable properties.
Results. By assuming efficient angular momentum transport, we are able to present a model that is capable of simultaneously predicting the following three main gravitationalwave observables: the effective inspiral spin parameter χ_{eff}, the chirp mass M_{chirp} and the cosmological redshift of merger z_{merger}. We find an excellent agreement between our model and the ten events from the first two advanced detector observing runs. We make predictions for the third observing run O3 and for Advanced LIGO design sensitivity. We expect approximately 80% of events with χ_{eff} < 0.1, while the remaining 20% of events with χ_{eff} ≥ 0.1 are split into ∼10% with M_{chirp} < 15 M_{⊙} and ∼10% with M_{chirp} ≥ 15 M_{⊙}. Moreover, we find that M_{chirp} and χ_{eff} distributions are very weakly dependent on the detector sensitivity.
Conclusions. The favorable comparison of the existing LIGO/Virgo observations with our model predictions gives support to the idea that the majority, if not all of the observed mergers, originate from the evolution of isolated binaries. The firstborn black hole has negligible spin because it lost its envelope after it expanded to become a giant star, while the spin of the secondborn black hole is determined by the tidal spin up of its naked helium star progenitor by the firstborn black hole companion after the binary finished the commonenvelope phase.
Key words: stars: black holes / gravitational waves / binaries close / black hole physics
© ESO 2020
1. Introduction
During the first and second observing runs O1/O2 of the advanced gravitationalwave (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 firstborn 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 firstborn 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/(GM^{2}), where J is the angular momentum of the BH. Using matchedfiltering analysis, GW observations provide estimates for each of the abovementioned quantities for both parent BHs. Although individual BH spin magnitudes and orientations are poorly constrained with present GW measurements, the effective inspiral spin parameter
the massweighted 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 spinorbitcoupling term in the postNewtonian 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
which, to firstorder approximation, determines the frequency evolution of the GW signal emitted during the BBH’s inspiral phase (Cutler & Flanagan 1994). The ten observed M_{chirp} span the range of 7.9−35.7 M_{⊙} with a pileup 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 hydrogenrich envelope past the Rochelobe and begins transferring mass to the secondary until it loses its entire envelope, leaving a naked heliumburning star. Following winddriven 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–Hestar 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 firstborn BH is the efficiency of angular momentum (AM) transport through the evolution of the progenitor star during the red supergiant phase. From observations of astroseismology (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 Rochelobe overflow mass transfer and wind mass loss. This leads to very slowly spinning BHs (a_{1} ≃ 0) as was shown by Qin et al. (2018; see also Fuller & Ma 2019). Assuming efficient AM transport, the angular momentum of the secondborn BH is mainly determined by the net effect of the stellar winds and the tidal interaction of the BH–Hestar 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 Hestar. Several studies attempted to model the last evolutionary phase of this channel and derived constraints on the spin using analytical arguments and semianalytical 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 bimodal 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 selfconsistently 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 winddriven orbital widening can lead to tidal decoupling. Ignoring this effects underestimates the impact of stellar winds on the final spin of the secondborn 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 Hestar 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 MonteCarlo 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; VignaGómez et al. 2018; Neijssel et al. 2019) to evolve isolated stellar binaries until the formation of BH–Hestar systems, namely the immediate progenitors of BBHs. In Sect. 2.1 we briefly describe the COMPASmodel 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 firstborn black hole can be treated as a pointlike particle, this approach allows us to track the angular momentum profile evolution of the Hestar 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 Hestars into BHs. Finally, taking into account the redshift dependence of metallicity, starformation 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_{⊙} < m_{1} < 200 M_{⊙} following the initial mass function of Kroupa (2001). The mass distribution of the less massive secondary star is given by m_{2} = m_{1}q_{0} where q_{0} is the initial mass ratio (0 < q_{0} < 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_{⊙} < m_{1} < 150 M_{⊙}. This means we only model a fraction f_{corr} 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 logorbital 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 logmetallicity range Z ∈ [0.0001, 0.0349]. We then evolve three million binaries per metallicity bin ΔZ_{j} with a total of star forming mass on the order of M_{sim, ΔZj} = 6.5 × 10^{7} 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 HumphreyDavidson 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 Rochelobe 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 winddriven mass loss, collapses into a BH. The star collapses into a pointlike 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 Hestar 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 corotating 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 nonCE channel, a selfconsistent 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 Hertzsprunggap (HG). In the optimistic scenario we apply the usual α − λ prescription to evolve these systems. In the pessimistic scenario we assume that Hertzsprunggap stars have not yet developed a sufficiently sharp density gradient at the coreenvelope 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 lowmass BHs compared to the optimistic. This is because a greater fraction of the total postmainsequence expansion occurs during the HG for highmetallicity stars (Linden et al. 2010). Therefore, forbidding the channel with HG CE donors has a particularly strong effect at high metallicity, which yields lowermass BHs due to metallicityenhanced stellar winds.
In Fig. 1 we show the distributions of orbital separation, firstborn BH mass, Hestar mass and metallicity of our BH–Hestar binaries after the CE phase. These distributions are weighted by the redshift and metallicitydependent starformation 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 logdistributed and most of the firstborn BHs have masses smaller than 30 M_{⊙}. Moreover, we see that formation metallicities of progenitors of merging compactobject binaries follow a skewed lognormal distribution. This is because the mean metallicity of the Universe decreases as a function of the lookback time and the starformation rate peaks at a redshift ∼2, that is a lookback time of ∼10.5 Gyr, where most of the binaries are formed. These distributions are used as an initial condition for our detailed modeling.
Fig. 1. Parameter distributions of the binary population after the CE phase. These BH–Hestar 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, firstborn BH mass M_{BH}, Hestar mass M_{Hestar} and metallicity Z weighted by the integrated redshift and metallicitydependent 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 Hestar and tidal interaction between the BH and the Hestar. These simulations^{1} are based on the work of Qin et al. (2018) and are adapted for MESA r10398.
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 massloss rate of
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 Hestar 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 Hestar with the orbit. We assume that the CE ejection leaves a circular binary, and the system remains circular during Hestar 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
where the Hestar has mass M, radius R and moment of inertia I, r_{g} given by is the dimensionless gyration radius of the Hestar, q is the mass ratio of the BH mass to the Hestar mass, and E_{2} is the second order tidal coefficient. We take the new fitting formula of E_{2} as suggested by Qin et al. (2018) for Hestars
where R_{conv} is the radius of the convective core (see Appendix A in Qin et al. 2018, for an indepth discussion of E_{2}). 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 EddingtonSweet circulation, the GoldreichSchubertFricke instability, and secular as well as dynamical shear mixing. In addition, diffusion element mixing is included with an efficiency parameter of f_{c} = 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 Hestars 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 Hestar.
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 Hestar mass before the supernova, the carbonoxygen (CO) core mass presupernova, the resultant secondborn BH mass, the orbital period presupernova, the lifetime of the BH–Hestar binary system (from the expulsion of the CE to the collapse of the Hestar) and the spin of the secondborn BH, as a function of the initial parameters of the BH–Hestar binary: initial mass of the firstborn BH, m_{BH}, initial mass of the Hestar, m_{Hestar}, initial orbital period, p, and Hestar 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, m_{BH} ∈ [2.5 M_{⊙}, 60 M_{⊙}], m_{Hestar} ∈ [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 underperforming. 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 Hestar mass presupernova and resultant BH mass respectively, 0.20% for the CO core mass of the Hestar, 0.01% for the orbital period, 0.04% for the lifetime of the BH–Hestar binary and 0.41% for the logspin of the secondborn 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 nonfittable 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 Hestar mass presupernova and resultant BH mass, 8% of the CO core mass of the Hestar, 6% of the orbital period, 6% of the lifetime of the BH–Hestar binary and 17% of the logspin of the secondborn BH.
2.3. Corecollapse 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 shockwave 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, f_{fb}, 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 corecollapse 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 & RamirezRuiz (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 m_{shell} 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),
where r_{isco} is the radius at ISCO for prograde equatorial orbits,
with z_{1} = 1 + (1 − a^{2})^{1/3}((1 + a)^{1/3} + (1 − a)^{1/3}) and 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
where the disk formation angle is given by
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 m_{disk} = m_{shell} cos θ_{disk}. The massenergy accreted from the disk onto the BH is ΔM_{disk} = m_{disk}(1−2GM_{BH}/(3c^{2}r_{isco}))^{1/2} and the accreted angular momentum is ΔJ_{disk} = m_{disk}j_{isco} (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 pairinstability supernovae (PISNe) and pulsational pairinstability supernovae (PPISNe). Both events are caused by the production of electronpositron 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 Hecore 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 presupernovae stars with Hecore 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 M_{BH} ≳ 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 (BisnovatyiKogan 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 fallback mass fraction f_{fb} (Fryer et al. 2012). In the Fryer et al. (2012) that we adopt, this quantity depends on the carbonoxygen core mass m_{core} of the star before the collapse. For core masses grater than 11 M_{⊙}, f_{fb} = 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 secondborn 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
where m_{1} and m_{2} 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:
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 anticorrelated with the spin of the secondborn BH, a_{2}, 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 Hestar 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 . In the second relation we used Kepler’s third law with ω being the orbital frequency matching the Hestar’s angular frequency Ω and in the last one a_{1} = 0 as assumed in our model. Meanwhile, for the wider binaries partially synchronized by dynamical tides on a synchronization timescale , we recover small spins a_{2} ∼ Ω ∼ exp(1/T_{sync})∼1/T_{sync} ∼ A^{−17/2} (cf. Eq. (4)) and therefore . In Fig. 2 we show the combined distribution of the merger timescale and the effective spin parameter for a specific metallicity bin (centered at log_{10}(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).
Fig. 2. Combined distribution of the BBH merger timescale T_{merger} vs. the effective inspiral spin parameter χ_{eff} for our synthetic BBH population at metallicity log_{10}(Z) = −2.5, in gray. The lighter shades represent larger contour levels, 68, 95 and 99%, respectively. Tidally locked systems follow the relation , orange dashed line, while the others systems follow 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 starformation 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 and Ω_{m} = 0.307 (Planck Collaboration XIII 2016). Every simulated BBH k with BH masses m_{1, k} and m_{2, k}, born at redshift z_{f, i} and merging at redshift z_{m, i, k} set by the delay time of this binary contributes to the detection rate by the following weight
where subscripts f and m refer to formation and merger, respectively, and fSFR is the fractional starformation rate, that is the total mass of stars formed per comoving volume per year per metallicity interval ΔZ_{j}. 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 lognormally 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, M_{sim, ΔZj}/f_{corr} is the matter simulated in the metallicity bin ΔZ_{j} rescaled by the normalization factor f_{corr} (see Appendix A), is the comoving distance to the source and p_{det} 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
where we add the contribution of every binary placed at the center of each formation time bin Δt_{i} in its corresponding metallicity bin ΔZ_{j}. The population synthesis predictions are performed in finite time bins of Δt = 100 Myr and the logmetallicity 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 signaltonoise ratio (S/N). In our model we assume that signals are detected if their singledetector 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 p_{det} for a given set of parameters (m_{1}, m_{2}, z). The optimal S/N (for a faceon, 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 p_{det}.
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 secondborn 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 firstborn BH. As we already discussed earlier, here we assume that the spin of the firstborn BH is very low, a_{1} ≃ 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 Hestar. 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 stellarevolution 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 LIGOVirgo 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 metallicityspecific star formation history (e.g., Chruslinska et al. 2019; Neijssel et al. 2019).
In Fig. 3 we show the predicted twodimensional 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, M_{chirp} 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, z_{merger} vs. χ_{eff} and z_{merger} vs. M_{chirp}, 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.
Fig. 3. Model predictions for binary black hole observables: chirp mass M_{chirp}, effective inspiral spin parameter χ_{eff} and cosmological redshift of merger z_{merger} 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 LIGOVirgo Collaboration (LVC) data with their 90% credible intervals. 
To provide a qualitative comparison between our theoretical predictions of χ_{eff} and LIGOVirgo 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 a_{1, 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 zerocentered 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.
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 zerocentered 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. 
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.
Fig. 6. Model predictions for the chirp mass M_{chirp}, effective inspiral spin parameter χ_{eff}, cosmological redshift of merger z_{merger}, BBH merger timescale T_{merger}, 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 M_{chirp} vs. χ_{eff} with red lines into three regions at χ_{eff} = 0.1 and M_{chirp} = 15 M_{⊙}. RegionA contains 80% of the events, 10% are in RegionB and 10% are in RegionC. For illustrative purposes, all histograms are plotted with a smoothing scale of 0.8 bins with the exception of z_{formation} vs. χ_{eff}, z_{formation} vs. z_{merger}, q vs. χ_{eff} that have no smoothing. 
As pointed out previously, the twodimensional 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: RegionA with χ_{eff} < 0.1, RegionB with χ_{eff} ≥ 0.1 and M_{chirp} < 15 M_{⊙} and RegionC with χ_{eff} ≥ 0.1 and M_{chirp} ≥ 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 anticorrelation between the merging timescale of the BBHs T_{merger} and the spin of the secondborn BH a_{2} or equivalently the observed quantity χ_{eff} (see Sect. 2.4). RegionA contains systems with low spins which translate into long merger timescales. For reference, at formation of the secondborn 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 RegionB have short merger timescales since they have a high χ_{eff}: at the formation of the secondborn 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 RegionC have high spins and high chirp masses. These are binaries that formed with low metallicity and merge quickly, that is at z_{merger} ≃ z_{formation}. This part of the parameter space really probes the lowend 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 infinitelysensitive detector, in gray. The entire population of merging BBHs has a peak in the merger redshift at around z_{merger} ≃ 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 M_{chirp}. 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.
Fig. 7. Model predictions for the chirp mass M_{chirp}, effective inspiral spin parameter χ_{eff}, cosmological redshift of merger z_{merger}, BBH merger timescale T_{merger}, 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 firstborn BH. Meanwhile, the spin of the secondborn 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 subgiants 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 nonspinning helium star, we ran a grid of 5000 MESA simulations of BH–Hestar 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 secondborn 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 secondborn 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 Hestar 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 secondborn BH a_{2} from detailed BH–Hestar binaries simulations with masses M_{BH} = 30 M_{⊙}, M_{Hestar} = 35 M_{⊙} and metallicities Z_{⊙}, 10^{−1} Z_{⊙}, 10^{−2} Z_{⊙} assuming a uniform distribution in logorbital 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 winddriven mass loss causes the binaries to widen and the tidal interaction which spins up the Hestar gradually becomes less effective. In the same figure we also show our predicted a_{2} 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 a_{2} = 0 while the latter has a peak at around a_{2} = 0.2. Although it has a peak at a_{2} = 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 lognormal 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 metallicitydependent star formation rate and apply aLIGO selection effects. The anticorrelation between T_{merger} and a_{2} 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 noninclusion 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).
Fig. 8. Normalized distribution of the spin of the secondborn 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–Hestar binary systems with M_{BH} = 30 M_{⊙}, M_{Hestar} = 35 M_{⊙} and metallicities Z_{⊙} (orange), 10^{−1} Z_{⊙} (green) and 10^{−2} Z_{⊙} (red) under the assumption that the initial separation is loguniform 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 (a_{1, 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, , and high effective spin parameter, , which is marginally consistent with our model. If similar events with bettermeasured 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 dynamicallyformed 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, antialigned 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 M_{chirp} ≃ 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 Hestars after the CE phase are not necessarily zeroage 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 zeroage 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 lowmetallicity 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) redshiftdependent 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 RegionC in the χ_{eff} − M_{chirp} histogram of Fig. 6, namely BBHs with high M_{chirp} and high χ_{eff}, at the expense of systems in RegionA, 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–Hestar 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 Hestar, 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 starformation 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 populationsynthesis model assumptions and a specific choice of a metallicityspecific 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 M_{chirp}, χ_{eff} and z_{merger} 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 secondborn 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 M_{chirp} < 15 M_{⊙} are BBHs born in close orbits (with orbital periods of less than 1 day at the formation of the secondborn 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 M_{chirp} ≥ 15 M_{⊙} are BBHs formed at low metallicity at a range of redshifts; these systems again merge promptly. Efficient spinup 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 secondborn BH. Moreover, we find that the assumption of efficient angular momentum transport has a negligible impact on the spin of the secondborn BH.
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 SklodowskaCurie 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 PHY1607611. 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 communitydeveloped core Python package for Astronomy (Astropy Collaboration 2013, 2018).
References
 Aasi, J., Abbott, B. P., Abbott, R., et al. 2015, Class. Quant. Grav., 32, 074001 [CrossRef] [Google Scholar]
 Aasi, J., Abbott, B. P., Abbott, R., et al. 2016, Phys. Rev. D, 93, 042006 [NASA ADS] [CrossRef] [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016a, ApJ, 818, L22 [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016b, Phys. Rev. Lett., 116, 241102 [NASA ADS] [CrossRef] [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016c, Phys. Rev. X, 6, 041015 [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, Phys. Rev. Lett., 118, 221101 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017b, ApJ, 850, L40 [NASA ADS] [CrossRef] [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2018, Liv. Rev. Relativ., 21, 3 [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019a, Phys. Rev. X, 9, 031040 [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019b, ApJ, 882, L24 [NASA ADS] [CrossRef] [Google Scholar]
 Abt, H. A. 1983, ARA&A, 21, 343 [Google Scholar]
 Acernese, F., Agathos, M., Agatsuma, K., et al. 2015, Class. Quant. Grav., 32, 024001 [Google Scholar]
 Antonini, F., Chatterjee, S., Rodriguez, C. L., et al. 2016, ApJ, 816, 65 [NASA ADS] [CrossRef] [Google Scholar]
 Arca Sedda, M., & Benacquista, M. 2019, MNRAS, 482, 2991 [NASA ADS] [Google Scholar]
 Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Astropy Collaboration (PriceWhelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
 Baldry, I. K., Driver, S. P., Loveday, J., et al. 2012, MNRAS, 421, 621 [NASA ADS] [Google Scholar]
 Bardeen, J. M. 1970, Nature, 226, 64 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Bardeen, J. M., Press, W. H., & Teukolsky, S. A. 1972, ApJ, 178, 347 [NASA ADS] [CrossRef] [Google Scholar]
 Barkat, Z., Rakavy, G., & Sack, N. 1967, Phys. Rev. Lett., 18, 379 [NASA ADS] [CrossRef] [Google Scholar]
 Barrett, J. W., Mandel, I., Neijssel, C. J., Stevenson, S., & VignaGó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]
 Barrett, J. W., Gaebel, S. M., Neijssel, C. J., et al. 2018, MNRAS, 477, 4685 [NASA ADS] [CrossRef] [Google Scholar]
 Batta, A., & RamirezRuiz, E. 2019, ArXiv eprints [arXiv:1904.04835] [Google Scholar]
 Belczynski, K., Taam, R. E., Kalogera, V., Rasio, F. A., & Bulik, T. 2007, ApJ, 662, 504 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Kalogera, V., Rasio, F. A., et al. 2008, ApJS, 174, 223 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Bulik, T., Fryer, C. L., et al. 2010, ApJ, 714, 1217 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Holz, D. E., Bulik, T., & O’Shaughnessy, R. 2016a, Nature, 534, 512 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Repetto, S., Holz, D. E., et al. 2016b, ApJ, 819, 108 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Klencki, J., Meynet, G., et al. 2017, ArXiv eprints [arXiv:1706.07053] [Google Scholar]
 BisnovatyiKogan, G. S. 1993, Astron. Astrophys. Trans., 3, 287 [NASA ADS] [CrossRef] [Google Scholar]
 Bocquet, S., & Carter, F. W. 2016, J. Open Source Softw., 1 [Google Scholar]
 Bowler, R. A. A., Dunlop, J. S., McLure, R. J., et al. 2015, MNRAS, 452, 1817 [NASA ADS] [CrossRef] [Google Scholar]
 Broekgaarden, F. S., Justham, S., de Mink, S. E., et al. 2019, MNRAS, 490, 5228 [NASA ADS] [CrossRef] [Google Scholar]
 Burrows, A., & Hayes, J. 1996, Phys. Rev. Lett., 76, 352 [NASA ADS] [CrossRef] [Google Scholar]
 Cantiello, M., Mankovich, C., Bildsten, L., ChristensenDalsgaard, J., & Paxton, B. 2014, ApJ, 788, 93 [Google Scholar]
 Chaboyer, B., & Zahn, J. P. 1992, A&A, 253, 173 [NASA ADS] [Google Scholar]
 Chruslinska, M., Nelemans, G., & Belczynski, K. 2019, MNRAS, 482, 5012 [NASA ADS] [CrossRef] [Google Scholar]
 Cutler, C., & Flanagan, É. E. 1994, Phys. Rev. D, 49, 2658 [NASA ADS] [CrossRef] [Google Scholar]
 de Kool, M. 1990, ApJ, 358, 189 [NASA ADS] [CrossRef] [Google Scholar]
 de Mink, S. E., & Belczynski, K. 2015, ApJ, 814, 58 [NASA ADS] [CrossRef] [Google Scholar]
 de Mink, S. E., Cantiello, M., Langer, N., et al. 2009, A&A, 497, 243 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dewi, J. D. M., & Tauris, T. M. 2000, A&A, 360, 1043 [NASA ADS] [Google Scholar]
 Dominik, M., Belczynski, K., Fryer, C., et al. 2012, ApJ, 759, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Dominik, M., Belczynski, K., Fryer, C., et al. 2013, ApJ, 779, 72 [NASA ADS] [CrossRef] [Google Scholar]
 Dominik, M., Berti, E., O’Shaughnessy, R., et al. 2015, ApJ, 806, 263 [NASA ADS] [CrossRef] [Google Scholar]
 Farr, W. M., & Mandel, I. 2018, Science, aat361, 6506 [NASA ADS] [CrossRef] [Google Scholar]
 Farr, W. M., Stevenson, S., Miller, M. C., et al. 2017, Nature, 548, 426 [NASA ADS] [CrossRef] [Google Scholar]
 Figer, D. F. 2005, Nature, 434, 192 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Finkelstein, S. L., Ryan, Jr., R. E., Papovich, C., et al. 2015, ApJ, 810, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Finn, L. S., & Chernoff, D. F. 1993, Phys. Rev. D, 47, 2198 [NASA ADS] [CrossRef] [Google Scholar]
 Fowler, W. A., & Hoyle, F. 1964, ApJS, 9, 201 [Google Scholar]
 Fragione, G., Ginsburg, I., & Kocsis, B. 2018, ApJ, 856, 92 [NASA ADS] [CrossRef] [Google Scholar]
 Fragos, T., Andrews, J. J., RamirezRuiz, E., et al. 2019, ApJ, 883, L45 [NASA ADS] [CrossRef] [Google Scholar]
 Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, ApJ, 749, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Fuller, J., & Ma, L. 2019, ApJ, 881, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Fuller, J., Lecoanet, D., Cantiello, M., & Brown, B. 2014, ApJ, 796, 17 [Google Scholar]
 Fuller, J., Piro, A. L., & Jermyn, A. S. 2019, MNRAS, 485, 3661 [NASA ADS] [Google Scholar]
 Gehan, C., Mosser, B., Michel, E., Samadi, R., & Kallinger, T. 2018, A&A, 616, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gerosa, D., Berti, E., O’Shaughnessy, R., et al. 2018, Phys. Rev. D, 98, 084036 [NASA ADS] [CrossRef] [Google Scholar]
 Giacobbo, N., & Mapelli, M. 2018, MNRAS, 480, 2011 [NASA ADS] [CrossRef] [Google Scholar]
 Giacobbo, N., Mapelli, M., & Spera, M. 2018, MNRAS, 474, 2959 [NASA ADS] [CrossRef] [Google Scholar]
 Grazian, A., Fontana, A., Santini, P., et al. 2015, A&A, 575, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Grevesse, N., & Sauval, A. J. 1998, Space Sci. Rev., 85, 161 [NASA ADS] [CrossRef] [Google Scholar]
 Hamann, W. R., & Koesterke, L. 1998, A&A, 335, 1003 [NASA ADS] [Google Scholar]
 Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368 [NASA ADS] [CrossRef] [Google Scholar]
 Heger, A., Woosley, S. E., & Spruit, H. C. 2005, ApJ, 626, 350 [NASA ADS] [CrossRef] [Google Scholar]
 Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974 [NASA ADS] [CrossRef] [Google Scholar]
 Hotokezaka, K., & Piran, T. 2017, ApJ, 842, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Humphreys, R. M., & Davidson, K. 1994, PASP, 106, 1025 [NASA ADS] [CrossRef] [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543 [NASA ADS] [CrossRef] [Google Scholar]
 Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897 [NASA ADS] [CrossRef] [Google Scholar]
 Hut, P. 1981, A&A, 99, 126 [NASA ADS] [Google Scholar]
 Ilbert, O., McCracken, H. J., Le Fèvre, O., et al. 2013, A&A, 556, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Inayoshi, K., Hirai, R., Kinugawa, T., & Hotokezaka, K. 2017, MNRAS, 468, 5020 [NASA ADS] [CrossRef] [Google Scholar]
 Ishigaki, M., Kawamata, R., Ouchi, M., et al. 2015, ApJ, 799, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Ivanova, N., Justham, S., Chen, X., et al. 2013, A&ARv, 21, 59 [Google Scholar]
 Janka, H.T. 2013, MNRAS, 434, 1355 [Google Scholar]
 Janka, H. T., & Mueller, E. 1994, A&A, 290, 496 [NASA ADS] [Google Scholar]
 Janka, H. T., Langanke, K., Marek, A., MartínezPinedo, G., & Müller, B. 2007, Phys. Rep., 442, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Kajisawa, M., Ichikawa, T., Tanaka, I., et al. 2009, ApJ, 702, 1393 [NASA ADS] [CrossRef] [Google Scholar]
 Kalogera, V. 1996, ApJ, 471, 352 [NASA ADS] [CrossRef] [Google Scholar]
 Kalogera, V., Belczynski, K., Kim, C., O’Shaughnessy, R., & Willems, B. 2007, Phys. Rep., 442, 75 [NASA ADS] [CrossRef] [Google Scholar]
 Kalogera, V., Berry, C. P., Colpi, M., et al. 2019, Bull. Am. Astron. Soc., 51, 242 [Google Scholar]
 Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
 Kruckow, M. U., Tauris, T. M., Langer, N., Kramer, M., & Izzard, R. G. 2018, MNRAS, 481, 1908 [NASA ADS] [CrossRef] [Google Scholar]
 Kushnir, D., Zaldarriaga, M., Kollmeier, J. A., & Waldman, R. 2016, MNRAS, 462, 844 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, K.S., Ferguson, H. C., Wiklind, T., et al. 2012, ApJ, 752, 66 [NASA ADS] [CrossRef] [Google Scholar]
 LIGO Scientific Collaboration 2018, LIGO Algorithm Library – LALSuite, Free Software (GPL) [Google Scholar]
 Linden, T., Kalogera, V., Sepinsky, J. F., et al. 2010, ApJ, 725, 1984 [NASA ADS] [CrossRef] [Google Scholar]
 Livio, M., & Soker, N. 1988, ApJ, 329, 764 [NASA ADS] [CrossRef] [Google Scholar]
 Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415 [NASA ADS] [CrossRef] [Google Scholar]
 Madau, P., & Fragos, T. 2017, ApJ, 840, 39 [NASA ADS] [CrossRef] [Google Scholar]
 Mandel, I., & de Mink, S. E. 2016, MNRAS, 458, 2634 [NASA ADS] [CrossRef] [Google Scholar]
 Mandel, I., & Farmer, A. 2018, ArXiv eprints [arXiv:1806.05820] [Google Scholar]
 Mapelli, M., & Giacobbo, N. 2018, MNRAS, 479, 4391 [NASA ADS] [CrossRef] [Google Scholar]
 Marchant, P., Langer, N., Podsiadlowski, P., Tauris, T. M., & Moriya, T. J. 2016, A&A, 588, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marchant, P., Renzo, M., Farmer, R., et al. 2019, ApJ, 882, 36 [Google Scholar]
 McLeod, D. J., McLure, R. J., Dunlop, J. S., et al. 2015, MNRAS, 450, 3032 [NASA ADS] [CrossRef] [Google Scholar]
 McLeod, D. J., McLure, R. J., & Dunlop, J. S. 2016, MNRAS, 459, 3812 [NASA ADS] [CrossRef] [Google Scholar]
 Miller, M. C. 2016, Gen. Relativ. Grav., 48, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Miller, M. C., & Lauburg, V. M. 2009, ApJ, 692, 917 [NASA ADS] [CrossRef] [Google Scholar]
 Neijssel, C. J., VignaGómez, A., Stevenson, S., et al. 2019, MNRAS, 490, 3740 [NASA ADS] [CrossRef] [Google Scholar]
 Oesch, P. A., Bouwens, R. J., Illingworth, G. D., et al. 2015, ApJ, 808, 104 [NASA ADS] [CrossRef] [Google Scholar]
 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]
 Park, D., Kim, C., Lee, H. M., Bae, Y.B., & Belczynski, K. 2017, MNRAS, 469, 4665 [NASA ADS] [CrossRef] [Google Scholar]
 Pavlovskii, K., Ivanova, N., Belczynski, K., & Van, K. X. 2017, MNRAS, 465, 2092 [NASA ADS] [CrossRef] [Google Scholar]
 Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
 Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15 [Google Scholar]
 Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Peters, P. C. 1964, Phys. Rev., 136, B1224 [NASA ADS] [CrossRef] [Google Scholar]
 Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 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]
 Portegies Zwart, S. F., & McMillan, S. L. W. 2000, ApJ, 528, L17 [NASA ADS] [CrossRef] [Google Scholar]
 Postnov, K. A., & Kuranov, A. G. 2019, MNRAS, 483, 3288 [NASA ADS] [CrossRef] [Google Scholar]
 Postnov, K. A., & Yungelson, L. R. 2014, Liv. Rev. Relativ., 17, 3 [Google Scholar]
 Punturo, M., Abernathy, M., Acernese, F., et al. 2010, Class. Quant. Grav., 27, 194002 [Google Scholar]
 Qin, Y., Fragos, T., Meynet, G., et al. 2018, A&A, 616, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Qin, Y., Marchant, P., Fragos, T., Meynet, G., & Kalogera, V. 2019, ApJ, 870, L18 [NASA ADS] [CrossRef] [Google Scholar]
 Rakavy, G., & Shaviv, G. 1967, ApJ, 148, 803 [NASA ADS] [CrossRef] [Google Scholar]
 Rodriguez, C. L., & Loeb, A. 2018, ApJ, 866, L5 [NASA ADS] [CrossRef] [Google Scholar]
 Rodriguez, C. L., Morscher, M., Pattabiraman, B., et al. 2015, Phys. Rev. Lett., 115, 051101 [NASA ADS] [CrossRef] [Google Scholar]
 Rodriguez, C. L., Zevin, M., Pankow, C., Kalogera, V., & Rasio, F. A. 2016, ApJ, 832, L2 [NASA ADS] [CrossRef] [Google Scholar]
 Rodriguez, C. L., AmaroSeoane, P., Chatterjee, S., & Rasio, F. A. 2018, Phys. Rev. Lett., 120, 151101 [NASA ADS] [CrossRef] [Google Scholar]
 Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444 [Google Scholar]
 Santamaría, L., Ohme, F., Ajith, P., et al. 2010, Phys. Rev. D, 82, 064016 [NASA ADS] [CrossRef] [Google Scholar]
 Schneider, F. R. N., Izzard, R. G., Langer, N., & de Mink, S. E. 2015, ApJ, 805, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Schneider, F. R. N., Sana, H., Evans, C. J., et al. 2018, Science, 359, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Sigurdsson, S., & Hernquist, L. 1993, Nature, 364, 423 [NASA ADS] [CrossRef] [Google Scholar]
 Silsbee, K., & Tremaine, S. 2017, ApJ, 836, 39 [NASA ADS] [CrossRef] [Google Scholar]
 Smarr, L. L., & Blandford, R. 1976, ApJ, 207, 574 [NASA ADS] [CrossRef] [Google Scholar]
 Socrates, A., Blaes, O., Hungerford, A., & Fryer, C. L. 2005, ApJ, 632, 531 [NASA ADS] [CrossRef] [Google Scholar]
 Spruit, H. C. 1999, A&A, 349, 189 [NASA ADS] [Google Scholar]
 Spruit, H. C. 2002, A&A, 381, 923 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stevenson, S., VignaGómez, A., Mandel, I., et al. 2017, Nat. Commun., 8, 14906 [NASA ADS] [CrossRef] [Google Scholar]
 Stevenson, S., Sampson, M., Powell, J., et al. 2019, ApJ, 882, 121 [Google Scholar]
 Suijs, M. P. L., Langer, N., Poelarends, A. J., et al. 2008, A&A, 481, L87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Thorne, K. S. 1974, ApJ, 191, 507 [Google Scholar]
 Tutukov, A. V., & Yungelson, L. R. 1993, MNRAS, 260, 675 [NASA ADS] [CrossRef] [Google Scholar]
 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]
 van den Heuvel, E. P. J., Portegies Zwart, S. F., & de Mink, S. E. 2017, MNRAS, 471, 4256 [NASA ADS] [CrossRef] [Google Scholar]
 Venumadhav, T., Zackay, B., Roulet, J., Dai, L., & Zaldarriaga, M. 2019, ArXiv eprints [arXiv:1904.07214] [Google Scholar]
 VignaGómez, A., Neijssel, C. J., Stevenson, S., et al. 2018, MNRAS, 481, 4009 [NASA ADS] [CrossRef] [Google Scholar]
 Vink, J. S., & de Koter, A. 2005, A&A, 442, 587 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Webbink, R. F. 1984, ApJ, 277, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Woosley, S. E. 2017, ApJ, 836, 244 [NASA ADS] [CrossRef] [Google Scholar]
 Xu, X.J., & Li, X.D. 2010, ApJ, 716, 114 [NASA ADS] [CrossRef] [Google Scholar]
 Yoshida, T., Umeda, H., Maeda, K., & Ishii, T. 2016, MNRAS, 457, 351 [NASA ADS] [CrossRef] [Google Scholar]
 Zackay, B., Dai, L., Venumadhav, T., Roulet, J., & Zaldarriaga, M. 2019a, ArXiv eprints [arXiv:1910.09528] [Google Scholar]
 Zackay, B., Venumadhav, T., Dai, L., Roulet, J., & Zaldarriaga, M. 2019b, Phys. Rev. D, 100, 023007 [Google Scholar]
 Zahn, J. P. 1977, A&A, 500, 121 [NASA ADS] [Google Scholar]
 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 MonteCarlo simulation and generate a sample of half a million BH–Hestar 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 m_{A} < m_{primary} < m_{B} 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:
where the constant f_{0} is defined such that . Let f_{bin} 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
In the case when m_{A}, m_{B} > m_{2}, we only model the following fraction of all systems:
The mean mass of a binary system in our simulated population is;
Thus, the total modeled mass M_{*, model} represents only a fraction of the total stellar population mass M_{*}:
and we must renormalize by the inverse of f_{corr} 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, m_{1} = 0.08 M_{⊙}, m_{2} = 0.5 M_{⊙}, using the observed f_{bin} = 0.7 (Sana et al. 2012), arbitrarily choosing m_{min} = 0.01 M_{⊙} and m_{max} = 200 M_{⊙} as the minimum and the maximum stellar mass, and carrying out the simulation for primary masses in the range between m_{A} = 5 M_{⊙} and m_{B} = 150 M_{⊙} (Figer 2005), we obtain .
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 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),
We assume that the metallicities of the binaries follow a truncated lognormal distributed,
with standard deviation σ = 0.5 and mean where the mean metallicity is (Madau & Fragos 2017)
The lognormal distribution is truncated at the highest metallicity bin edge, Z_{max} = 0.034923, and the distribution is accordingly renormalized to ensure that . Portions of the distribution extending beyond the lower limit edge Z_{min} = 0.000091 are included in the edge bin when integrating over metallicity.
We compute the detection rate by integrating the cosmological merger rate 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
where the factor account for the difference in clock rates at merger and at the detector and p_{det} 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 , where
and D_{c}(z) is the comoving distance which is related to the luminosity distance as D_{L}(z) = D_{c}(z)(1 + z) and is computed as follows
where . The merger rate R(z_{m}) can be rewritten as the convolution of the starformation rate, SFR, that is the total mass of stars formed per comoving volume per year, and the number density of binaries per unit star forming mass M_{f} per unit masses m_{1}, m_{2} per unit logmetallicity interval Z per unit time delay τ = t_{m} − t_{f}:
where we used the compact notation z_{f} ≡ z(t_{f}). The time delay τ is mostly set by T_{merger}, since the GWdriven merger takes much longer than stellar evolution for BH progenitors. Performing the change of variable , the integral of Eq. (B.4) translates into the following MonteCarlo sum over the formation time intervals arbitrarily chosen as Δt_{i} = 100 Myr and 30 uniformlydistributed logmetallicity intervals for Z ∈ [Z_{min}, Z_{max}]
where M_{sim, ΔZj} is the total mass simulated per logmetallicity interval ΔZ_{j} and fSFR is the total mass of stars formed per comoving volume per year per logmetallicity interval ΔZ_{j},
where Z_{j} is the center of the logmetallicity bin ΔZ_{j} corresponding to the metallicity Z_{k} 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 ,
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 m_{BH} ∈ [2.5 M_{⊙}, 60 M_{⊙}], m_{Hestar} ∈ [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 Hestar mass and its CO core mass before the supernova, the resulting BH mass, the orbital period before the supernova, the lifetime of the BH–Hestar binary and the spin of the secondborn black hole. All physical quantities are logtransformed 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 X_{i} is computed as Δ_{i} = (X_{true, i} − X_{interp.,i})/X_{true, 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
where we arbitrarily limit the maximal combined relative error to , 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 nonfittable points. These stabilize at 0.01%, 0.20% and 0.04% for the Hestar 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–Hestar binary and 0.4% for the logspin of the secondborn 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 nonfittable 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 Hestar 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–Hestar binary and 17% for the logspin of the secondborn 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.
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 Rochelobe overflow (indicating a difference between COMPAS and MESA models). The right plot includes NaNs (obtained from nonfittable 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
Fig. 1. Parameter distributions of the binary population after the CE phase. These BH–Hestar 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, firstborn BH mass M_{BH}, Hestar mass M_{Hestar} and metallicity Z weighted by the integrated redshift and metallicitydependent 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 
Fig. 2. Combined distribution of the BBH merger timescale T_{merger} vs. the effective inspiral spin parameter χ_{eff} for our synthetic BBH population at metallicity log_{10}(Z) = −2.5, in gray. The lighter shades represent larger contour levels, 68, 95 and 99%, respectively. Tidally locked systems follow the relation , orange dashed line, while the others systems follow dictated by the tidal synchronization timescale, blue dashed line. Both lines are drawn at an arbitrary ordinate. 

In the text 
Fig. 3. Model predictions for binary black hole observables: chirp mass M_{chirp}, effective inspiral spin parameter χ_{eff} and cosmological redshift of merger z_{merger} 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 LIGOVirgo Collaboration (LVC) data with their 90% credible intervals. 

In the text 
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 zerocentered 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 
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 
Fig. 6. Model predictions for the chirp mass M_{chirp}, effective inspiral spin parameter χ_{eff}, cosmological redshift of merger z_{merger}, BBH merger timescale T_{merger}, 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 M_{chirp} vs. χ_{eff} with red lines into three regions at χ_{eff} = 0.1 and M_{chirp} = 15 M_{⊙}. RegionA contains 80% of the events, 10% are in RegionB and 10% are in RegionC. For illustrative purposes, all histograms are plotted with a smoothing scale of 0.8 bins with the exception of z_{formation} vs. χ_{eff}, z_{formation} vs. z_{merger}, q vs. χ_{eff} that have no smoothing. 

In the text 
Fig. 7. Model predictions for the chirp mass M_{chirp}, effective inspiral spin parameter χ_{eff}, cosmological redshift of merger z_{merger}, BBH merger timescale T_{merger}, 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 
Fig. 8. Normalized distribution of the spin of the secondborn 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–Hestar binary systems with M_{BH} = 30 M_{⊙}, M_{Hestar} = 35 M_{⊙} and metallicities Z_{⊙} (orange), 10^{−1} Z_{⊙} (green) and 10^{−2} Z_{⊙} (red) under the assumption that the initial separation is loguniform in 5 < A/R_{⊙} < 63. 

In the text 
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 Rochelobe overflow (indicating a difference between COMPAS and MESA models). The right plot includes NaNs (obtained from nonfittable 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 (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.