Gravitational-wave detection rates for compact binaries formed in isolation: LIGO/Virgo O3 and beyond

TL;DR

This paper estimates the detection rates of compact binary mergers formed in isolation for current and future gravitational-wave detectors, analyzing how stellar physics uncertainties impact these rates and discussing prospects with third-generation observatories.

Contribution

It provides new estimates of merger and detection rates incorporating uncertainties in stellar physics and explores future detector capabilities.

Findings

Detection rates increase with detector upgrades.

Stellar physics uncertainties significantly affect rate estimates.

Third-generation detectors will reveal detailed merger rate features.

Abstract

Using simulations performed with the population synthesis code MOBSE, we compute the merger rate densities and detection rates of compact binary mergers formed in isolation for second- and third-generation gravitational-wave detectors. We estimate how rates are affected by uncertainties on key stellar-physics parameters, namely common envelope evolution and natal kicks. We estimate how future upgrades will increase the size of the available catalog of merger events, and we discuss features of the merger rate density that will become accessible with third-generation detectors.

Tables

Table 1. Table 1: Catalog of MOBSE models considered in this study.

Model	${\sigma }_{\mathrm{CCSN}}$	$\alpha$
$\alpha 1$	$265\mathrm{km}/s$	$1$
$\alpha 3$	$265\mathrm{km}/s$	$3$
$\alpha 5$	$265\mathrm{km}/s$	$5$
CC15$\alpha 1$	$15\mathrm{km}/s$	$1$
CC15$\alpha 3$	$15\mathrm{km}/s$	$3$
CC15$\alpha 5$	$15\mathrm{km}/s$	$5$

Table 2. Table 2: Minimum and maximum detection rates ( yr − 1 superscript yr 1 {\rm yr}^{-1} ) across all models.

Detector	BNS	NSBH	BBH
O2	$0.028$ - $0.91$	$0.12$ - $1.1$	$27$ - $40$
O3	$0.11$ - $3.4$	$0.46$ - $3.9$	$94$ - $1.5\text{×}{10}^2$
AdLIGO	$0.27$ - $8.6$	$1.2$ - $9.3$	$2.2\text{×}{10}^2$ - $3.6\text{×}{10}^2$
A+	$0.88$ - $28$	$3.2$ - $26$	$5.6\text{×}{10}^2$ - $9.7\text{×}{10}^2$
A++	$2.3$ - $71$	$8.1$ - $63$	$1.3\text{×}{10}^3$ - $2.4\text{×}{10}^3$
Voyager	$32$ - $9.4\text{×}{10}^2$	$1.0\text{×}{10}^2$ - $7.8\text{×}{10}^2$	$9.7\text{×}{10}^3$ - $2.7\text{×}{10}^4$
ET-B	$1.1\text{×}{10}^3$ - $2.7\text{×}{10}^4$	$2.4\text{×}{10}^3$ - $2.2\text{×}{10}^4$	$4.9\text{×}{10}^4$ - $2.7\text{×}{10}^5$
CE	$1.6\text{×}{10}^4$ - $2.7\text{×}{10}^5$	$1.6\text{×}{10}^4$ - $1.4\text{×}{10}^5$	$8.6\text{×}{10}^4$ - $5.4\text{×}{10}^5$
Noiseless	$2.8\text{×}{10}^4$ - $4.5\text{×}{10}^5$	$2.0\text{×}{10}^4$ - $1.8\text{×}{10}^5$	$9.2\text{×}{10}^4$ - $5.7\text{×}{10}^5$

Equations

$$\mathcal{F}(m_1)\propto m_1^{-2.3}\mathrm{with}{m}_1\in [5-150]M_{\odot },$$

$$\mathcal{F}(q)\propto q^{-0.1}\mathrm{with}q\in [0.1-1].$$

$$\mathcal{F}(\mathcal{P})\propto {\mathcal{P}}^{-0.55}\mathrm{with}\mathcal{P}=\mathrm{lo}{\mathrm{g}}_{10}\left(\frac{P}{\mathrm{day}}\right)\in [0.15-5.5]$$

$$\mathcal{F}(e)\propto e^{-0.42}\mathrm{with}0\le e<1.$$

$$R_{\det }={\int }_0^{t_0}{p}_{\det }\mathcal{R}(z_m)\frac{d{V}_c}{d{t}_m}\frac{d{t}_m}{d{t}_{\det }}d{t}_m,$$

$$\mathcal{R}(z_m)$$

$$R_{\det }$$

$$\int \frac{dN}{d{t}_m}f(t_m)d{t}_m\to \sum _{i}f(t_m^i).$$

$$\frac{d{V}_c}{dz}(z)=4\pi \frac{c}{H_0}\frac{D_{\mathrm{c}}^2}{E(z)},$$

$$\log \frac{Z(z_f)}{Z_{\odot }}=\{\begin{array}{ll}-0.19z_f,& z_f\le 1.5\\ -0.22z_f,& z_f>1.5,\end{array}$$

$$\mathrm{sfr}(z)=\frac{0.015(1+z{)}^{2.7}}{1+[(1+z)/2.9{]}^{5.6}}{M}_{\odot }{\mathrm{Mpc}}^{-3}.$$

$$R_{\det }=\sum _i(s_i(t_f)\Delta t_f)\frac{p_{\det }}{1+z_m}\frac{d{V}_{\mathrm{c}}}{d{z}_m}\frac{d{z}_m}{d{t}_m},$$

$$s_i(z_f)\Delta t_f=f_{\mathrm{bin}}{f}_{\mathrm{IMF}}\frac{\mathrm{sfr}(z_f)}{M_{\mathrm{sim}}(Z_f)}\Delta t_f.$$

$${\rho }_0^2=4{\int }_0^{\infty }\frac{{\tilde{h}}^{\ast }(f)\tilde{h}(f)}{S_h(f)}df,$$

$$p_{\det }={\int }_{8/{\rho }_0}^{1}p(\omega )d\omega$$

Full text

I Introduction

The detection of gravitational waves (GWs) from $10$ binary black holes (BBHs) and a binary neutron star (BNS) in the first two LIGO/Virgo observing runs Abbott et al. (2018a), and the subsequent detections of numerous compact binary candidates in the third observing run, naturally lead to the question: how do these binaries form, and what is the physics that drives their evolution?

Advanced LIGO (AdLIGO) is expected to reach design sensitivity in the near future, the so-called A$+$ upgrade to current detectors was already approved for funding, and further upgrades (A$++$ and Voyager) are expected in the near future Abbott et al. (2018b); Adhikari (2014); Miller et al. (2015); LIGO Scientific Collaboration (2015); Abbott et al. (2017a). The GW community is also planning future, “third-generation” (3G) facilities, such as the Einstein Telescope (ET) ETw ; Punturo et al. (2010) and Cosmic Explorer (CE) Abbott et al. (2017a), which will extend the observable horizon to the very early Universe.

As GW detectors improve and the number of detections grows, we will gather information about the environments in which compact binaries form, and constrain the physical parameters that drive their evolution. Future GW detectors will measure compact binary parameters (such as masses and spins) within few per cent accuracy Vitale and Evans (2017), reconstructing fine details of distribution of these observables. They will observe sources up to redshifts as large as $z\sim 10^{2}$ Hall and Evans (2019), allowing us to study how the merger rate density evolves with redshift, and ultimately to constrain astrophysical models Fishbach et al. (2018); Vitale and Farr (2018); Abbott et al. (2018c). The large number of detections that comes with increased sensitivity will also reduce statistical errors on the parameters that describe compact binary populations to few per cent with $\sim 10^{3}$ observations  Barrett et al. (2018).

Compact-object binaries could form either in the field Belczynski et al. (2001); Hurley et al. (2002) or through dynamical interactions in young Banerjee et al. (2010); Ziosi et al. (2014); Mapelli (2016), nuclear Antonini and Rasio (2016); Hoang et al. (2018) or globular clusters Rodriguez et al. (2016); Askar et al. (2017). In this paper we present updated detection rates, and a roadmap of our prospects for constraining the astrophysics of compact binaries in the near future. We study how detection rates for binaries formed in isolation (“field binaries”) will evolve with future improvements of GW detectors, with the goal to understand if and when characteristic features of the astrophysical populations will become visible.

The plan of the paper is as follows. In Sec. II we present our astrophysical populations based on the MOBSE population-synthesis code (Mapelli et al., 2017; Giacobbo et al., 2018). In Sec. III we investigate how uncertainties in binary evolution affect the evolution of the merger rate density, and what new generation of detectors can tell us about this evolution. In Sec. IV we compute detection rates for each of the six models we consider and for different detector sensitivities. In Sec. V we summarize our findings and out line directions for future work. Appendix A gives details on how detection rates are computed from the MOBSE simulations. Throughout the paper we use the standard cosmological parameters determined by the Planck Collaboration Ade et al. (2016). We assume that a source is detected if the single-detector signal-to-noise ratio (SNR) $\rho$ is such that $\rho>8$.

II Astrophysical populations

We use simulations performed with the population-synthesis code MOBSE Giacobbo et al. (2018). MOBSE is an upgrade of the BSE code Hurley et al. (2002) which includes up-to-date prescriptions for the evolution of massive stars. The treatment of stellar winds accounts for the stellar metallicity and luminosity dependence of the mass loss. Compact objects are produced via different channels, including core-collapse, electron-capture and (pulsational) pair instability supernovae (SNe).

In our simulations, the primary star’s mass $m_{\mathrm{1}}$ is distributed according to the Kroupa mass function Kroupa (2001)

[TABLE]

while the mass ratio $q=m_{2}/m_{1}$ scales like Sana et al. (2012)

[TABLE]

The orbital period $P$ is drawn from

[TABLE]

and the eccentricity $e$ follows the distribution Sana et al. (2012)

[TABLE]

Among the many physical processes involved in the formation of compact binaries that can merge within a Hubble time, the so called common-envelope phase is believed to be critical Mapelli (2018); Mandel and Farmer (2018). When a star in a binary system overfills its Roche lobe, it starts transferring mass, and eventually forms a common envelope that engulfs the companion. The common envelope does not corotate with the stars or their cores, and this leads to a drag force. As a result, the stars spiral in and transfer their orbital energy to the envelope. The system will survive only if the energy transferred is sufficient to eject the envelope Ivanova et al. (2013a, b); Dominik et al. (2012). The efficiency of this mechanism constitutes a main uncertainty in compact-binary formation modelling.

Another important source of uncertainty are natal kicks. If a compact object forms from a supernova explosion, it is expected to receive a birth kick because of asymmetric mass ejection. A non-zero kick (the so-called Blaauw kick Blaauw (1961)) is expected even in the unlikely case where mass loss is symmetric, but the compact object is part of a binary system. This natal kick can disrupt the binary or substantially modify its orbit. Kicks set the fraction of stellar binaries which are unbound by the SN explosion and, consequently, play a major role in determining GW detection rates Belczynski et al. (2001); Mapelli et al. (2017); Gerosa et al. (2018).

As described by Giacobbo and Mapelli (2018) and summarized in Table 1, we consider six representative populations of merging binaries, aiming at bracketing the uncertainties in the physics of both common envelope and natal kicks. These two parameters might be the first to be constrained with GW data (see e.g. O’Shaughnessy et al. (2017); Barrett et al. (2018)).

The common envelope phase is treated using the so-called $\alpha\lambda$ formalism Webbink (1984); Ivanova et al. (2013b), where $\alpha$ quantifies the efficiency of energy transfer to the envelope and $\lambda$ represents the binding energy of the envelope. In this work we consider $\alpha$ as a free parameter, while $\lambda$ depends on the stellar type Dewi and Tauris (2000) and it is computed by using the prescriptions derived in Ref. Claeys et al. (2014). Kicks are extracted from a Maxwellian distribution with root-mean-square speed (rms) $\sigma_{\rm CCSN}$ for core-collapse SNe that produce neutron stars.111Neutron stars can also form through electron-capture SNe, which are less energetic, faster and do not develop large asymmetries. This is generally expected to lead to small kicks, and therefore we assume $\sigma_{\rm ECSN}=15\,{\rm km/s}$ (Giacobbo and Mapelli, 2019). For black holes, we reduce the kick velocity $v_{\rm BH}$ by taking into account fallback: $v_{\rm BH}=(1-f_{\rm fb})v_{\rm NS}$, where $v_{\rm NS}$ is the natal kick for neutron stars and $f_{\rm fb}$ parametrizes the amount of fallback on the proto-compact object Fryer et al. (2012).

Models CC15 produce natal kicks $\leq{}100$ km s*-1*, and therefore they are in tension with the proper motions of the fastest single Galactic neutron stars Hobbs et al. (2005). These models were chosen because they give a local merger rate density of binary neutron stars consistent with the one inferred from GW170817 Abbott et al. (2017b), without requiring exotic assumptions about common envelope.

MOBSE predicts the NS masses from $1.1$ to $2M_{\odot}$ where light (heavy) NSs are preferred during BNS (NSBH) mergers. On the other hand, NSBH mergers favor low BH masses ($<15M_{\odot}$) while BBH mergers could have BHs as heavy as $45M_{\odot}$ with most binaries having mass ratios close to unity Giacobbo and Mapelli (2018).

III Merger rate densities

The merger rate density ${\mathcal{R}}(z_{m})$ as a function of merger redshift $z_{m}$ tracks the distribution of merging binaries across cosmic time, and it depends on two factors:

• (i)

the rate of binary formation at a given redshift $z_{f}$, and

• (ii)

the distribution of time delays $t_{\rm delay}$ between the formation of the parent stars in the binary and the merger of their compact object remnants.

In turn, binary formation at $z_{f}$ depends on the star formation rate and the metallicity, both of which evolve over time. The time delay distribution is sensitive to the physics that drives binary evolution (see e.g. Mandel and de Mink (2016); Rodriguez et al. (2016); Di Carlo et al. (2019)).

In Fig. 1 we plot the evolution of the merger rate density for the six MOBSE models considered in this study. The low-redshift behavior is often parametrized as a power law: ${\mathcal{R}}(z)\approx{\mathcal{R}}_{0}(1+z)^{\lambda_{0}}$ Abbott et al. (2018c); Fishbach et al. (2018), where ${\mathcal{R}}_{0}$ is the local merger rate density and $\lambda_{0}$ is a model-dependent parameter that describes its evolution. The parameter $\lambda_{0}$ can be used to infer astrophysical information. The star formation rate is well approximated by $\lambda_{0}\simeq 2.4$ for $0.1<z<1$ Fishbach et al. (2018). Therefore, an observed $\lambda_{0}<2.4$ would imply that mergers peaked before the peak of star formation, which is only possible if compact-object binary formation is high at low metallicities and if the time delays are short enough Fishbach et al. (2018). Current detectors can only investigate the evolution of the merger rate at low redshift, but in the near future we will be able to trace the redshift evolution of the merger rate density.

Figure 1 shows that the BNS rate density follows quite closely the star formation rate, with a peak at slightly lower redshift (because of the short but finite time delays). Current observations favor models with low kicks and large $\alpha$: as shown by the red shaded region in the top panel of Fig. 1, only low-kick models with $\alpha=3$ and $\alpha=5$ can explain the high local merger rates resulting from the detection of GW170817 (Mapelli and Giacobbo, 2018; Giacobbo and Mapelli, 2018). Most BNS formation models have weak dependence on metallicity. Quite interestingly, some of them show a bimodal distribution, with a dip at $z_{m}\approx 5.6$ and a second peak at $z_{m}\approx 9$. Indeed, the efficiency in forming merging BNS has a minimum at intermediate metallicity $Z\sim 0.1Z_{\odot}$ (see e.g. Fig. 14 of Giacobbo and Mapelli (2018)). Stars at intermediate metallicities tend to develop larger radii, and this leads to the formation of wide BNS systems that either do not merge in a Hubble time, or are easily disrupted by a SN explosion (because of their large orbital separation). However, not all models that show a dip in the merger efficiency lead to a bimodal merger rate density. Since most detectors are not sensitive to binaries from such large redshifts, 3G detectors are needed to observe this behavior in the early Universe.

By contrast, BBH production is very efficient at low metallicities because of the impact of metallicity on stellar radii and evolutionary stages. At solar metallicity massive stars become Wolf-Rayet stars quite rapidly, after leaving the giant branch, because of stellar wind efficiency. Wolf-Rayet stars have small radii ($1-2~{}R_{\odot}$); thus, it is highly unlikely that such stars enter common envelope. Without common envelope, the binary star evolves into a BBH with a large orbital separation, which will not be able to merge within a Hubble time. In contrast, metal-poor massive stars can retain a large fraction of their hydrogen envelope and avoid the Wolf-Rayet stage, increasing the probability of undergoing mass transfer and entering common envelope. The rate density peaks at $z\gtrsim 2$, earlier than the peak of star formation, and the merger rate density at small redshifts is not as steep as the star formation rate (i.e., it has $\lambda_{0}<2.4$). We should soon be able to verify this trend with current detectors.

IV Detection rates

To study how detection rates will benefit from detector improvements, here we will consider noise power spectral densities for the AdLIGO design sensitivity noise Abbott et al. (2018b); planned upgrades to existing LIGO detectors (A+, A++ and Voyager Miller et al. (2015); Adhikari (2014); LIGO Scientific Collaboration (2015)); and 3G detectors, including CE Abbott et al. (2017a) and the Einstein Telescope (more specifically, ET-B ETw ). We approximate the detector noise for the O2 and O3 runs by rescaling the AdLIGO noise curve in such a way that the resulting BNS range is $90$ Mpc Abbott et al. (2018a) and $140$ Mpc Det , respectively. In Fig. 2 we plot the distribution of signal-to-noise ratios (SNRs) for these detectors using the low-kick model with $\alpha=5$. Most of the binaries with very large SNRs come from local Universe, so their distribution scales like $1/\rho{}^{4}$ Schutz (2011).222In the local Universe, the total number of binaries within luminosity distance $D_{*}$ is $N(D<D_{*})\propto D_{*}^{3}$, or equivalently $N(\rho>\rho_{*})\ \propto\rho_{*}^{-3}$, so the SNR probability distribution scales like $N(\rho_{*})=\frac{dN(\rho>\rho_{*})}{d\rho_{*}}\propto\rho_{*}^{-4}$. Since CE (and, for BBHs, also ET) will see past the peak of the merger rate density (cf. Fig. 1), the maximum detection redshift is not controlled by the detector capabilities, but by the physics that governs the merger rate density ${\mathcal{R}}(z_{m})$. Figure 3 shows the detection rates, $R_{\rm det}$ for different astrophysical models and different detectors, comparing them with the intrinsic merger rate in the Universe that would correspond to an ideal, noiseless detector (see Appendix A for details of the detection-rate calculations). According to our models, AdLIGO at design sensitivity could see $220-360$ BBH, up to $9$ NSBH and $9$ BNS mergers per year. Upgrading AdLIGO detectors to a configuration like A+ would increase the detection rates by a factor of $3$. With 3G detectors, BBH rates would increase by up to $2$–$3$ orders of magnitude, while NSBH and BNS detection rates would increase by up to $3$–$4$ orders of magnitude. CE would see at least $92\%$ of all BBH mergers in the Universe, compared to the $0.06$–$0.24\%$ seen by AdLIGO at design sensitivity. Current-generation detectors like AdLIGO have low BNS and NSBH detection rates, detecting only $10^{-5}$ ($\sim 10^{-4}$) of all BNS (NSBH) mergers in the Universe. By contrast, CE will see more than $50\%$ ($\sim 75\%$) of all BNS (NSBH) mergers.

It is also clear from Fig. 3 that $\alpha$ and $\sigma_{\rm CCSN}$ can affect detection rates of all compact binary systems by up to an order of magnitude. In particular, BBH and BNS rates are affected in different ways by the common-envelope efficiency parameter $\alpha$: lower values of $\alpha$ yield smaller rates for BNSs and larger rates for BBHs. This can be understood as follows. BBHs form from massive stars that can develop very large radii during their evolution, and therefore enter the common envelope phase with a wide orbital separation. If $\alpha>1$, the envelope will be ejected easily while the binary is still widely separated, and the outcome will be a wide binary that is unlikely to merge in a Hubble time Giacobbo and Mapelli (2018). In contrast, BNSs form from smaller stars, and the orbital separation at the beginning of the common envelope phase is smaller. Therefore high values of $\alpha$ lead to the formation of a close binary that can merge in a Hubble time, while small values of $\alpha$ cause a premature merger of the system.

Low kicks (CC$15\alpha 1$, CC$15\alpha 3$, CC$15\alpha 5$) lead to higher detections rates for BNS and NSBH mergers, because strong kicks are efficient at disrupting these binaries. On the other hand, most BBH progenitors undergo direct collapse in the models presented here: nearly all of the star’s mass falls back onto the compact object, and kicks are suppressed. For this reason, BBH detection rates are nearly insensitive to natal kicks.333BBH merger rates are found to strongly depend on SN kicks if fallback is suppressed (Mapelli et al., 2017; Wysocki et al., 2018; Gerosa et al., 2018).

Local NSBH merger rates for low-kick models are larger than high-kick models by a factor of 3–10. If we assume low (high) SN kicks, the NSBH merger rate increases (decreases) with $\alpha$. This is because large SN kicks tend to unbind the binary. If the natal kick is high, a small value of $\alpha{}$ increases the probability that the system merges, because if $\alpha$ is small the system’s semi-major axis shrinks considerably during CE, after the first supernova. Thus, if the kick is high a small value of $\alpha{}$ increases the NSBH merger rate. In contrast, if the kick is low, a small value of $\alpha{}$ might trigger the premature merger of the binary, before the second compact object has formed. Thus, if the kick is low, the highest NSBH merger rate is achieved for a rather large value of $\alpha{}$, as already explained in Mapelli and Giacobbo (2018).

We list minimum and maximum rates across all models in Table 2.

V Conclusions

We studied the detection rates and redshift evolution of BNS, NSBH and BBH merger rate densities. The redshift distribution of the merger rates contains important clues about the physics that drives the evolution of these compact objects (see also the companion papers Mapelli et al. (2017); Mapelli and Giacobbo (2018); Mapelli et al. (2019)). The merger rate history of compact-object binaries is obtained by convolving their formation history with the time-delay distribution. The formation rate depends on both star formation rate and metallicity. The formation of BNSs depends only mildly on metallicity, and therefore their formation across cosmic time follows quite closely the star formation rate (but it is shifted to slightly lower redshifts, because of finite delay times). Therefore for BNSs we expect $\lambda_{0}\gtrsim 2.4$, i.e. the merger rate peak occurs after, but very close to the peak of star formation. Current detectors have small BNS horizons, so they will mainly see binaries that formed in the local Universe, where metallicity is high, but 3G detectors should allow us to observe large-redshift BNSs and to verify this prediction. In contrast, BBH production (and, marginally, NSBH production) is very efficient at low metallicities. Most BBHs form at $z\gtrsim 2$, before the peak of star formation, and their merger rate density evolves slowly compared to BNSs: most BBHs and NSBHs formed before the peak of star formation, yielding $\lambda_{0}<2.4$. Only CE (and, in the case of BBHs, ET) will allow us to see beyond the merger rate peak of compact object binaries.

We also investigated how these rates are affected by common-envelope efficiency and natal kicks, considering both second- and third-generation detectors. We found that a lower common envelope efficiency leads to smaller BNS detection rates, and larger BBH detection rates. This is because lower efficiency causes a longer inspiral of the stellar cores, leading to BNS progenitors that merge prematurely, before they can collapse into a neutron star. By contrast, BBH progenitors are much larger, and their orbits are wider compared to BNS progenitors. Natal kick assumptions affects only BNS and NSBH mergers in our models: high kicks can more easily disrupt binaries and usually lead to lower detection rates. On the other hand, BBH kicks are suppressed because of the large amount of material that falls back onto the compact object after the supernova explosion.

In Fig. 4 we plot the growth of the GW catalog size as detectors improve, based on the rate calculations of Fig. 3. We assume 1 year of observations for O3, which started in 2019. The observing run O4 for AdLIGO at design sensitivity is expected to start in 2021, and it should last for $\sim 2$ years, followed by 1 year of commissioning period for upgrades to A+ (which is currently targeted to be operational by 2024 LIGO Scientific Collaboration (2018)). We assume the operational time for A+ to be 6 years LIGO Scientific Collaboration (2017), with further upgrades to “A++” in 2027. By the beginning of the 2030s, when new detectors – Voyager in the existing LIGO facilities, and CE/ET in separate facilities – may start operations, we could have a GW catalog of up to $10^{4}$ events. In Fig. 4 we assume a 5-year observation period before Voyager is superseded by CE.

As the detectors improve, the rapid growth of the GW catalog should allow us to place stringent constraints on the population parameters that influence the final stages of the lives of massive stars.

Acknowledgments

MM and YB acknowledge financial support by the European Research Council for the ERC Consolidator grant DEMOBLACK, under contract no. 770017. EB and VB are supported by NSF Grant No. PHY-1841464, NSF Grant No. AST-1841358, NSF-XSEDE Grant No. PHY-090003, and NASA ATP Grant No. 17-ATP17-0225. This work has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 690904. The authors would like to acknowledge networking support by the COST Action GWverse CA16104. Computational work was performed on the University of Birmingham’s BlueBEAR cluster and at the Maryland Advanced Research Computing Center (MARCC).

Appendix A Detection rate calculations

The detection rate is given by Dominik et al. (2015); Belczynski et al. (2016)

[TABLE]

where $t_{0}$ is the age of universe and $p_{\rm det}$ is the probability of detecting a given binary, defined in Eq. (15) below. The factor ${dt_{m}}/{dt_{\rm det}}={1}/{(1+z_{m})}$ accounts for the different clock rates at merger and at the detector. The source-frame merger rate density at redshift $z_{m}$ is

[TABLE]

where the star-formation rate is ${\rm sfr}(z_{f})\equiv\frac{dM_{f}}{dV_{c}dt_{f}}$. The second term in the integrand accounts for the number of binaries per unit star-forming mass that form at $t_{f}$ and merge at $t_{m}$. Here, we have marginalized over the distribution of component masses and time delays. We can rewrite Eq. (5) (after switching the order of the integrals over $t_{f}$ and $t_{m}$) as

[TABLE]

In the second line above, we converted the integral over a distribution to a Monte-Carlo sum,

[TABLE]

In practice, the term in parentheses is evaluated by Monte Carlo integrations, where the samples $t^{i}_{m}$ are generated from the distribution ${dN}/{dt_{m}}$. The comoving volume element $dV_{\rm c}/dz$ is given by

[TABLE]

where $E(z)$ is the function that describes the evolution of Hubble parameter, i.e. $H(z)=H_{0}E(z)$, and $D_{\rm c}$ is comoving distance Hogg (1999). The factor of $4\pi$ takes into account the angular integration over the sky.

In practice, at a given metallicity $Z_{f}$, MOBSE starts with a given total mass $M_{\rm sim}$ and outputs a distribution of binaries. For each set of free parameters in Table 1, we have $12$ simulations of $10^{7}$ binaries each, with metallicities $Z=0.01$–$1\,Z_{\odot}$. We simulate a set of compact-object binaries formed at different times $t_{f}$ inside bins of $\Delta t_{f}=10{\rm~{}Myr}$. At the time of formation $t_{f}$, we assume that the metallicity is given by

[TABLE]

i.e. we follow the metallicity evolution of Ref. Rafelski et al. (2012), but we rescale it so that $Z(0)=Z_{\odot}$. Each formation time bin is assigned one the $12$ metallicities according to Eq. (10). However, since the MOBSE simulation started with total binary mass, $M_{\rm sim}$, we need to rescale this mass according to the star formation in that particular time bin. We have adopted the following fit for star formation rate Madau and Dickinson (2014):

[TABLE]

These binaries are then evolved in time until they merge at $t_{m}$. This produces a catalog of binaries that form at $t_{f}$ and merge at $z_{m}$. The integral in Eq. (7) can be now be written as

[TABLE]

where all terms except the first are evaluated at the merger redshift $z_{m}$. The first term is the number density of binaries formed at redshift $z_{f}$,

[TABLE]

The factors $f_{\rm bin}=0.5$ and $f_{\rm IMF}=0.285$ take into account the fact that MOBSE only simulates binaries with primary mass larger than $5M_{\odot}$.

Finally, a binary is assumed to be detected if it has the signal-to-noise ratio (SNR) $\rho=\rho_{0}w>8$, where $\rho_{0}$ is the SNR assuming that the binary is optimally oriented and located in the sky, while $0\leq w\leq 1$ is the projection factor that depends on the binary’s sky position and orientation. The optimal SNR is calculated as

[TABLE]

where $h(f)$ is the frequency-domain GW signal and $S_{h}(f)$ is the detector noise power spectral density Moore et al. (2015); Sathyaprakash and Schutz (2009). The horizon $z_{h}$ is the farthest redshift for which a binary with component masses $m_{1}$ and $m_{2}$ can be detected, i.e. $\rho_{0}(m_{1},m_{2},z_{h})=8$. The quantity $\rho_{0}$ determines the probability of detecting a binary that lies within the detector’s horizon (i.e. $\rho_{0}>8$, or equivalently $z<z_{h}$):

[TABLE]

where $p(w)$ is the probability distribution function of $\omega$ Finn and Chernoff (1993). Detection rates only depend on $p_{\rm det}$, hence $\rho_{0}$. We calculate the signal-to-noise ratio of BBH mergers using the waveform approximant IMRPhenomD, while for NSBH and BNS mergers we use TaylorF2. Since MOBSE does not have any prescriptions to evolve the spins, we assume black holes and neutron stars to be non-spinning. Spins are expected to impact detection rates within a factor 1.5 Gerosa et al. (2018), which should be added to the error budget of our estimates.

Note that in Fig 2, where we looked at the distribution of $\rho=\rho_{0}w$, we sample $p(\omega)$ for each binary in the catalogs mentioned above and assign the SNR accordingly.