# Metallicity-constrained merger rates of binary black holes and the stochastic gravitational wave background

Irina Dvorkin, Elisabeth Vangioni, Joseph Silk, Jean-Philippe Uzan, Keith A. Olive

► To cite this version:

Irina Dvorkin, Elisabeth Vangioni, Joseph Silk, Jean-Philippe Uzan, Keith A. Olive. Metallicity-constrained merger rates of binary black holes and the stochastic gravitational wave background. Monthly Notices of the Royal Astronomical Society, 2016, 461 (4), pp.3877 - 3885. 10.1093/mn-ras/stw1477. hal-01382790

HAL Id: hal-01382790

https://hal.sorbonne-universite.fr/hal-01382790v1

Submitted on 6 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Metallicity-constrained merger rates of binary black holes and the stochastic gravitational wave background

Irina Dvorkin,<fnref n="1" />★ Elisabeth Vangioni,<fnref n="1" /> Joseph Silk,<fnref n="1" />,<fnref n="2" />,<fnref n="3" />,<fnref n="4" /> Jean-Philippe Uzan<fnref n="1" /> and Keith A. Olive<fnref n="5" />

Accepted 2016 June 16. Received 2016 June 16; in original form 2016 April 22

## ABSTRACT

The recent detection of the binary black hole merger GW150914 demonstrates the existence of black holes more massive than previously observed in X-ray binaries in our Galaxy. This article explores different scenarios of black hole formation in the context of self-consistent cosmic chemical evolution models that simultaneously match observations of the cosmic star formation rate, optical depth to reionization and metallicity of the interstellar medium. This framework is used to calculate the mass distribution of merging black hole binaries and its evolution with redshift. We also study the implications of the black hole mass distribution for the stochastic gravitational wave background from mergers and from core-collapse events.

**Key words:** black hole physics – gravitational waves – binaries: general – stars: Population III.

## 1 INTRODUCTION

The recent detection by the Advanced Laser Interferometric Gravitational-wave Observatory (LIGO) of the gravitational wave source GW150914 (Abbott et al. 2016b) constitutes the first observational evidence for a merger of a binary black hole (BBH) system. The signal matches the waveform expected from a merger of two black holes (BHs) of masses $36_{-5}^{+5}M_{\odot}$ and $29_{-4}^{+4}M_{\odot}$ at a luminosity distance of $410_{-180}^{+160}Mpc$, corresponding to a redshift of $0.09_{-0.04}^{+0.03}$ (Abbott et al. 2016a,b).

One of the most interesting astrophysical questions raised by this discovery is how such 'heavy' BHs can form and what are the physical conditions required. The most massive X-ray binaries with reliably measured masses reach only $20 M_{\odot}$, making GW150914 the most massive stellar BBH ever observed (Abbott et al. 2016a). The masses of remnants that can form in a supernova (SN) have been studied extensively using various techniques (see Janka 2012, for a comprehensive review on SN explosion mechanisms). One of the first quantitative approaches (e.g. Woosley & Weaver 1995; Zhang, Woosley & Heger 2008) was to initiate the explosion artificially in a 1D (spherically symmetric) stellar model. These models typically find that the final mass of the collapsed remnant is sensitive to the explosion energy, the presupernova structure, the stellar mass and the metallicity. A reverse shock is generated when the SN

front shock travels through the hydrogen envelope of the star and can decelerate a significant amount of matter, further increasing the final remnant mass. BH masses obtained in these models are generally below $40 M_{\odot}$ with the maximum fraction of the mass of the progenitor star winding up in the BH of up to ~25 per cent for a solar-metallicity star and up to ~70 per cent for zero-metallicity star. Note that in order to calculate the mass distribution of BHs one needs to account for the stellar initial mass function (IMF) of the progenitor stars which peaks at low stellar masses (e.g. Kroupa 2001). As a result, piston-driven models tend to predict negligibly small number densities of heavy BHs.

In recent years, other models of SN collapse were built using higher-dimensional numerical simulations and/or more accurate description of the neutrino physics (e.g. Fryer & Young 2007; O'Connor & Ott 2011; Lentz et al. 2012; Müller, Janka & Marek 2012; Ugliano et al. 2012; Kochanek 2014; Clausen, Piro & Ott 2015; Pejcha & Thompson 2015; Ertl et al. 2016). These studies confirmed that BHs form for progenitor masses above ~25 $M_{\odot}$ either via direct collapse or via fallback and elucidated the connections between the progenitor mass and metallicity and the remnant mass. However, it is challenging to explore a wide range of progenitor masses and metallicites using higher-dimensional simulations due to the high computational costs.

An alternative method is to use analytical tools to estimate the explosion energy, as was done in Fryer (2006) under the assumption that the energy reservoir is limited to the convective region bounded by the edge of the proto-neutron star and the SN shock.

<footnote n="1">¹Sorbonne Universités, UPMC Univ Paris 6 et CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="2">²AIM-Paris-Saclay, CEA/DSM/IRFU, CNRS, Univ Paris 7, F-91191, Gif-sur-Yvette, France</footnote>

<footnote n="3">³Department of Physics and Astronomy, The Johns Hopkins University, Baltimore, MD 21218, USA</footnote>

<footnote n="4">⁴BIPAC, University of Oxford, 1 Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="5">⁵William I. Fine Theoretical Physics Institute, School of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA</footnote>

<footnote>* E-mail: dvorkin@iap.fr</footnote>
This recipe was subsequently used in Fryer et al. (2012) to study the dependence of the compact remnant mass function on the delay between core bounce and explosion and to demonstrate that the difference between possible explosion mechanisms will be detectable by gravitational wave observatories.

If the progenitors of the BBH evolve in an isolated environment without dynamical interactions (see e.g. Ziosi et al. 2014, for an alternative scenario where BBH form through dynamical interactions in dense stellar clusters), then the masses of the remnants as well as the merger delay time will depend also on their binary interactions. Evolutionary models (e.g. Belczynski et al. 2010a; Spera, Mapelli & Bressan 2015) have shown that single stars can form BHs as massive as ~100 M⊙ and follow-up population synthesis codes of isolated binary evolution (e.g. Belczynski et al. 2010b) predict the existence of BBH systems that merge within the age of the Universe. In particular, it was recently shown (Belczynski et al. 2016; Eldridge & Stanway 2016) that the stellar progenitors of GW150914 had to form in a low-metallicity environment. Other models were developed that predict the formation of massive BHs under specific conditions, such as chemically homogeneous evolution of the binary (e.g. Marchant et al. 2016; Mandel & de Mink 2016) or the direct collapse of a single, fast-rotating star to a binary which subsequently merges (e.g. Loeb 2016; Woosley 2016).

It has also been proposed that the first stellar generation (Population III; PopIII) which was created from zero metallicity gas is responsible for a significant fraction or even the majority of BBH mergers observable with gravitational wave observatories (e.g. Bond & Carr 1984; Belczynski, Bulik & Rudak 2004; Kowalska, Bulik & Belczynski 2012). Kinugawa et al. (2014) found that the typical mass of PopIII BBH is ~30 M⊙ and the detection rate is expected to be as high as ~140 events per year (although note that the evolution of massive stars at low metallicity could differ substantially from that in metal-rich environments; see Szécsi et al. 2015). However the importance of PopIII stars as merging BBH progenitors in the context of realistic galaxy evolution models is still debated. Recently, Hartwig et al. (2016) concluded, based on their self-consistent cosmological semi-analytic model, that there is only an ~1 per cent probability that GW150914 is of PopIII origin and that the GW background of BBH mergers produced by PopIII stars is small compared with other contributions at $f \simeq 25$ Hz, where the LIGO network is most sensitive. On the other hand, Inayoshi et al. (2016) found that the GW background produced by PopIII remnants can dominate other populations and would enable tight constraints on the PopIII properties.

The purpose of this paper is to explore the merger rate of BBHs in the context of self-consistent cosmic metallicity evolution models. In particular, we focus on the differences in the expected mass distribution of merging BBHs between different models of BH formation. With more detections of coalescing BBHs by the Advanced LIGO and VIRGO observatories expected in the near future, it might be possible to constrain the SN explosion mechanism from the observed BBH mass spectrum. We also explore the stochastic gravitational wave background from merger and SN collapse events and show that it might be possible to disentangle these two contributions with observations in the range $f \gtrsim 400$ Hz.

In this work we compare two models of BH formation, the 1D model of Woosley & Weaver (1995) and the analytic description of Fryer et al. (2012) which give the remnant mass as a function of progenitor mass and metallicity for a wide range of masses and metallicities. Since we also use Woosley & Weaver (1995) for the stellar yields, this first choice constitutes a fully self-consistent description, while the model by Fryer et al. (2012) is based on

a physically motivated explosion mechanism and its description of BH formation is more realistic. We also explore the model developed by Kinugawa et al. (2014) for PopIII stars. In addition, we consider various prescriptions for the star formation rate (SFR) and IMF. While this choice of models is far from being exhaustive, our goal here is to show how different prescriptions can be discriminated using upcoming observations of gravitational waves combined with constraints from reionization and metallicity measurements.

The structure of this paper is as follows. In Section 2 we describe the BH formation models we employ and the chemical evolution model in which they are embedded. In particular, we discuss constraints from observations of metal absorption lines and optical depth to reionization on the various prescriptions for the SFR that we use. Possible contributions from PopIII stars are also discussed. In Section 3 we present the BBH merger rate as a function of the BH mass expected in different models. These results are used in Section 4 to calculate the stochastic gravitational wave background from mergers as well as from SN collapse. We conclude in Section 5.

## 2 BINARY BLACK HOLE FORMATION SCENARIOS

The formation of BHs occurs at the end of the nuclear burning phase in massive stars and can proceed via two routes. For the lower mass end of BH formation, a meta-stable proto-neutron star (NS) is produced, followed by a formation of a BH through accretion of the part of the stellar envelope that could not be expelled in the SN. Direct collapse (sometimes called failed SN) into a BH occurs in the case of the most massive stars. Thus, the mass of the remnant is determined mainly by the mass of the star at the moment of collapse, as well as the explosion energy.

Massive stars generally experience strong winds which cause them to shed a significant fraction of their envelopes during their lifetime (e.g. Meynet & Maeder 2007; Vink 2008; Georgy, Meynet & Maeder 2011). The strength of these winds depends on metallicity: stars at lower metallicities exhibit weaker winds due to reduced opacity and easier radiation transport. Other factors that influence the relationship between the zero-age main sequence (ZAMS) mass and the stellar core mass at the time of collapse are stellar rotation (e.g. de Mink et al. 2009) and the microphysics of stellar evolution (e.g. Jones et al. 2015; Meynet et al. 2015). The mass of the BH formed after the collapse of the core depends on the strength of the SN explosion which determines how much material is ejected. In addition, after the shocked material slows down some of it may be decelerated and fall back on to the proto-NS, adding to the remnant mass.

In this paper, we consider two models of BH formation. Model WWp follows Woosley & Weaver (1995) up to progenitor masses of $m = 40 M_{\odot}$. For stars with initial masses above 40 $M_{\odot}$ the remnant mass is extrapolated as follows:

$$ \frac{m_{\text{rem}}}{m} = A \left( \frac{m}{40 M_{\odot}} \right)^{\beta} \frac{1}{\left( \frac{Z}{0.01 Z_{\odot}} \right)^{\gamma} + 1} \quad (1) $$

where Z is the metallicity. This functional form was chosen so as to match the results of Woosley & Weaver (1995) at $m = 40 M_{\odot}$ for the range of metallicities they explored and the dependence of the remnant mass on metallicity given in Crocker et al. (2015) (see their fig. 4, which is based on Belczynski et al. 2010a). The fiducial values of our extrapolation are $A = 0.3$, $\beta = 0.8$ and $\gamma = 0.2$. The yields in this case are scaled from the tabulated values of Woosley & Weaver (1995) so as to ensure mass conservation.
The second model we consider, called *Fryer* in what follows, is based on the calculations of remnant masses taken from Fryer et al. (2012) (in particular their delayed model) with the yields taken from Woosley & Weaver (1995). This model is based on the assumption that the explosion is powered by a convection-enhanced, neutrino-driven engine and the explosion energy, mass loss and fallback are calculated analytically.

For each of these models we also consider the possibility of different BH formation scenarios from metal-poor stars. Compared with the present-day stellar population these stars are expected to be more massive, have smaller radii for the same mass, and less mass loss by stellar winds during their lifetime. Well before the actual detection of GW150914, Kinugawa et al. (2014) found that the typical mass of BHs formed from PopIII stars that would be merging today is ~30 M☉. Inspired by this prediction, we assume that stars below some metallicity limit $Z_{\text{limit}} = 10^{-3}Z_{\odot}$ produce BHs according to the relation found in Kinugawa et al. (2014). This prescription is based on the evolutionary models of Marigo et al. (2001) for zero-metallicity stars, which provide the stellar radius and core mass, and the fitting formulae from Belczynski, Kalogera & Bulik (2002) for the remnant mass. These models are named $WWp+K$ and Fryer+$K$, in which case this special prescription is applied to all stars below the chosen metallicity limit. In the following we explore different values of $Z_{\text{limit}}$.

The models described above provide the remnant mass as a function of ZAMS mass and metallicity. We then assume that the remnant is a BH for remnant masses above $2.5M_{\odot}$ and a neutron star for lower masses. While mergers of double neutron stars and BH-neutron star binaries are expected to be detectable by gravitational wave observatories, we do not include the neutron star population here and leave it to future work.

These models of BH formation are embedded in the cosmic chemical evolution model based on Daigne et al. (2004, 2006) and Rollinde et al. (2009). We assume a Salpeter stellar IMF with slope $x = 2.35$ in the mass range 0.1–100 M☉ (Salpeter 1955). We also consider the case of $x = 2.7$ (which can be the case for the high-mass tail of the IMF in dense and turbulent environments; see Chabrier, Hennebelle & Charlot 2014) which we denote as steep IMF (we note that this is an extreme scenario). For the cosmic SFR we use the functional form of Springel & Hernquist (2003):

$$ \psi(z) = \frac{a \exp[b(z - z_m)]}{a - b + b \exp[a(z - z_m)]} \quad (2) $$

where $z$ is the redshift. Our fiducial model is a fit to the observations of luminous galaxies compiled by Behroozi, Wechsler & Conroy (2013) and complemented by high-redshift observations from Bouwens et al. (2015) and Oesch et al. (2015). We use the fit parameters given in Vangioni et al. (2015), namely $\nu = 0.178 M_{\odot}yr^{-1}Mpc^{-3}$, $z_m = 2$, $a = 2.37$ and $b = 1.8$.

An alternative way to calibrate the SFR at high redshifts is by considering the rate of gamma-ray bursts (GRBs; Robertson & Ellis 2012; Wang 2013; Kistler, Yuksel & Hopkins 2013). We use the results of Trenti, Perna & Tacchella (2013) and Behroozi & Silk (2015) to obtain a fit combining low-redshift galaxy luminosity data and high-redshift GRB data. This choice corresponds to Model 2 in Vangioni et al. (2015) with the parameters $\nu = 0.146 M_{\odot}yr^{-1}Mpc^{-3}$, $z_m = 1.72$, $a = 2.8$ and $b = 2.46$. This set of models is dubbed GRB-based.

Finally, we also explore the possibility of an early PopIII component with the Salpeter IMF in the mass range 36–100 M☉ and SFR parameters $\nu = 0.002 M_{\odot}yr^{-1}Mpc^{-3}$, $z_m = 11.87$, $a = 13.8$ and $b = 13.36$. This set of models is named PopIII. The SFR as a

<figure><img src="image_2.png" /><figcaption>**Figure 1.** SFR as a function of redshift for the models explored here. The fiducial model (red line) is based on the data compilation by Behroozi et al. (2013) and high-redshift observations from Bouwens et al. (2015) and Oesch et al. (2015). Adding a PopIII component produces a peak at $z \sim z_m$ (the combined Fiducial+PopIII is shown in yellow). A model that is based on high-redshift GRB data and the normalization from Trenti et al. (2013) results in a higher SFR at high redshifts (purple line).</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Optical depth as a function of redshift for the set of SFR models shown in Fig. 1, compared with the constraints from Planck Collaboration XIII et al. (2015) using their Planck TT + low P combination. All the models considered here are consistent with the CMB measurement.</figcaption></figure>

function of redshift for a set of representative models is shown in Fig. 1. Note that the SFR calibrated to GRB data results in much higher SFR at high redshifts.

The viability of different SFR models can be studied using the constraints from the optical depth to reionization, constrained from the analysis of the cosmic microwave background (CMB). We calculate the evolution of the volume-filling fraction of ionized regions and the Thomson optical depth as in Greif & Bromm (2006), assuming an escape fraction of $f_{esc} = 0.2$. The number of ionizing photons for massive stars is calculated using the tables in Schaerer (2002). In Fig. 2 we compare the optical depth in a set of our representative models with the value obtained by Planck Collaboration XIII et al. (2015), using their Planck TT + low P combination: $\tau = 0.078 \pm 0.019$. It can be seen that all the models considered here are consistent with the CMB measurement. We note that the latest Planck analysis (Planck Collaboration XLVI, XLVII et al. 2016a,b), which appeared after the present article was submitted, results in a
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Metallicity evolution as a function of redshift for a few representative models, compared with DLA data from Rafelski et al. (2012, black crosses). Blue and red lines show the *Fryer* and *WWp* models, respectively, using our fiducial SFR. We also show the results for the *WWp* model with alternative SFR and IMF prescriptions. At lower redshifts, the different models produce very similar results; however at higher redshift the *PopIII* and *GRB-based* models predict much higher metallicities, which affects the masses of BHs formed in these scenarios. In the *Steep IMF* model, fewer massive stars are formed which leads to significantly lower metallicity (green line).</figcaption></figure>

significantly lower value of the optical depth ($\tau = 0.058 \pm 0.012$). This new value can be accommodated if the escape fraction is lowered to $f_{esc} = 0.1$, consistent with some recent measurements (see e.g. Matthee et al. 2016). We intend to explore the implications of these results in future work (see also the discussion in Vangioni et al. 2015, on the effect of the escape fraction).

Various studies (e.g. Belczynski et al. 2010b; Fryer et al. 2012) have shown that the remnant mass of a given progenitor star is very sensitive to the metallicity, mainly because stars at lower metallicity produce weaker winds. The cosmic evolution of metallicity is therefore an important part of any model that attempts to calculate BBH merger rates. Our model follows the chemical enrichment of the interstellar medium (ISM) in a self-consistent manner and reproduces the metallicities measured in high-redshift damped Lyα absorbers (DLAs), as shown in Fig. 3. At lower redshifts the different models are practically indistinguishable (except for the steep IMF model which produces much fewer massive stars); however at higher redshift the *PopIII* and *GRB-based* models predict much higher ISM metallicities, a direct result of the increased SFR in these models. Below we will explore the effect this evolution has on the efficiency of Kinugawa-like models. Note that the metallicity is slightly reduced in the fiducial *Fryer* model relative to the *WWp* model because in this case more mass remains locked in heavy BHs. Indeed, as the efficiency of producing heavy BHs increases the metal yield of the star decreases as is required by simple mass conservation.

The various models of BH formation, SFR and IMF are summarized in Table 1. As we have shown, our fiducial models are consistent with the observed SFR, optical depth to reionization and cosmic metallicity evolution. We now proceed to calculating the birth and merger rates of BBH.

### 3 BINARY BLACK HOLE MERGER RATES

The direct outcome of our calculation is the birthrate of BHs as a function of mass and redshift (or, equivalently, time) $R_{\text{birth}}(t, m_{\text{bh}})$

(in units of events per unit time per unit comoving volume per unit BH mass):

$$R_{\text{birth}}(t, m_{\text{bh}}) = \int \psi[t - \tau(m)]\phi(m)\delta(m - g_{\text{bh}}^{-1}(m_{\text{bh}}))dm \quad (3)$$

where $\tau(m)$ is the lifetime of a star of mass $m$ (taken from Schaerer 2002), $\phi(m)$ is the IMF, $\psi(t)$ is the SFR, $\delta(m)$ is the Dirac delta function and the ZAMS stellar mass and BH mass are related by some function $m_{\text{bh}} = g_{\text{bh}}(m)$ which is implicit in the equation above. $g_{\text{bh}}$ also depends on time through its metallicity dependence and is calculated according to the prescription of each of our models. However, the relevant quantity for gravitational waves is the merger rate which depends also on the binary fraction and the delay time between the binary formation and merger. This time delay is in most cases larger than the age of the Universe, and in general depends on the orbital parameters of the binary (semi-major axis and eccentricity) which are not resolved in our model. To circumvent this difficulty we assume a delay time distribution from the models of Belczynski et al. (2016) and convolve it with the birth rate $R_{\text{birth}}(t)$ given by our models as follows:

$$R_m(t, m) = N \int_{t_{\text{min}}}^{t_{\text{max}}} R_{\text{birth}}(t - t_d, m) P(t_d) dt_d, \quad (4)$$

where $t_d$ is the delay time whose distribution is $P(t_d) \propto 1/t_d$ for $t_{\text{min}} < t_d < t_{\text{max}}$ with $t_{\text{min}} = 50$ Myr and $t_{\text{max}}$ equal to the Hubble time. We note that, in principle, the time delay depends on the properties of the binary (e.g. masses and initial orbital parameters; see Peters 1964) as well as its environment if the binary is non-isolated (e.g. Ziosi et al. 2014; Rodriguez, Chatterjee & Rasio 2016). A complete treatment of the time delay distribution is beyond the scope of the present paper and will be discussed in our future work. To account for binaries that did not merge at all we normalize the total birth rate (sum over all masses) using the observed rate of $10^{-7}$ Mpc$^{-3}$ yr$^{-1}$ at $z=0$, where we used the estimate from Abbott et al. (2016c) that assumes a power-law distribution of BH masses. This normalization is expressed by the constant $N$. Note that, in principle, the delay time itself is a function of the masses of the BHs in the binary, in particular BHs with roughly equal mass are expected to merge faster, as well as those with larger total mass. The treatment of this effect is beyond the scope of the present paper and we plan to address it in future work.

The total merger rate as a function of redshift for a representative subset of our models is shown in Fig. 4. Note that all the models produce very similar total merger rates, which is a direct consequence of the similarity in the SFR. The GRB-based model predicts the highest BH birth and merger rates at high redshifts due to the enhanced SFR.

While it is impossible to distinguish between the different models using the total BBH merger rate shown in Fig. 4, the predictions for merger rates per unit BH mass vary among the different prescriptions. Fig. 5 shows the merger rate per unit BH mass at $z=0$ for a few of the *WWp* models, as well as two of the *Fryer* models. Here the differences between the different models become apparent above $m_{\text{bh}} \simeq 20M_{\odot}$. Clearly, the *Fryer* model is much more efficient in producing massive BHs. In contrast, the fiducial *WWp* model struggles to produce massive BHs in significant amounts: these models differ by more than two orders of magnitude for $M = 30 M_{\odot}$. The rate of production of massive BHs is increased if we use steep IMF or a Kinugawa-like prescription for metal-poor stars. The reason for the shift towards higher BH masses in the steep IMF model is that in this case the metallicity is significantly lower, as can be seen in Fig. 1, so that even though there are fewer massive progenitor stars,
<table><caption>**Table 1.** A summary of the models studied in this work. BH masses are taken from either *WWp* or *Fryer* and can be supplemented by *Kinugawa*-like prescriptions, named *WWp+K* and *Fryer+K*, respectively. The *PopIII* model for SFR is added to the *Fiducial* model to produce a bimodal SFR as shown in Fig. 1.</caption><thead><tr><th></th><th>Model name</th><th>Ref.</th><th>Parameters</th><th>Parameter values</th></tr></thead><tbody><tr><td rowspan="4">BH masses</td><td>WWp</td><td>Woodsley & Weaver (1995)</td><td>A, β, γ</td><td>0.3, 0.8, 0.2</td></tr><tr><td>Fryer</td><td>Fryer et al. (2012)</td><td>—</td><td>—</td></tr><tr><td>WWp+K</td><td rowspan="2">Kinugawa et al. (2014)</td><td rowspan="2">Z<sub>limit</sub>/Z<sub>☉</sub></td><td rowspan="2">0.001 or 0.01</td></tr><tr><td>Fryer+K</td></tr><tr><td rowspan="3">SFR</td><td>Fiducial</td><td rowspan="3">Vangioni et al. (2015)</td><td rowspan="3">v, z<sub>m</sub>, a, b</td><td>0.178, 2.00, 2.37, 1.8</td></tr><tr><td>PopIII</td><td>0.002, 11.87, 13.8, 13.36</td></tr><tr><td>GRB-based</td><td>0.146, 1.72, 2.8, 2.46</td></tr><tr><td rowspan="2">IMF</td><td>Fiducial</td><td>Salpeter (1955)</td><td rowspan="2">x</td><td>2.35</td></tr><tr><td>Steep IMF</td><td>Chabrier et al. (2014)</td><td>2.7</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Total merger rate of BBH as a function of redshift, normalized to $10^{-7}$ Mpc$^{-3}$ yr$^{-1}$ at $z = 0$. With the given normalization, the curves are in practice indistinguishable below $z = 3$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Merger rate per unit BH mass at $z = 0$ (the x-axis shows the mass of each component of the binary, assumed equal in our model). Even with the normalization of the total rate for all the models, the mass dependence is different.</figcaption></figure>

the masses of the BHs that form are, on average, higher. Note that if we include a bimodal IMF to account for an early burst of PopIII stars (yellow dotted line) the *Kinugawa* prescription is not effective at producing BBH that merge at $z = 0$. The reason is that in this case the metallicity rises quickly at high redshifts, as shown in Fig. 3 and the *Kinugawa* prescription is in fact not employed. While it pro-

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Merger rate per unit mass as a function of redshift. Different curves correspond to 1 $M_{☉}$-wide mass bins. Solid curves are for the *WWp* model, while the dashed curve represents the merger rate of 30 $M_{☉}$ BBH in the *Fryer* model. The merger rate is dominated by low-mass binaries, with most of the contribution coming from 3 $M_{☉}$ BHs (just above our chosen limit between neutrons stars and BHs). The same holds for the *Fryer* model (not shown in this plot).</figcaption></figure>

duces some BBH at higher redshift, their number density at $z = 0$ is negligible, and in fact slightly reduced relative to the fiducial *WWp* model. The same effect happens when we try to employ a *GRB-based* SFR (purple dotted line overlapping with the yellow dotted line) which increases the metallicity at high redshift. However, in the fiducial SFR model the *Kinugawa* prescription significantly increases the number of BBH mergers with masses above ~30 $M_{☉}$ (solid green and dashed blue lines). As expected, the outcome of our *Kinugawa*-like models depends strongly on the assumed metallicity floor $Z_{limit}$, below which we employ the *Kinugawa*-like prescription. When it is used for metallicities as high as 0.01 $Z_{⊙}$ the difference with the *Fryer* model is less than one order of magnitude for $M \gtrsim 35 M_{☉}$.

In contrast to the *WWp* family of models, all of the *Fryer* models are very similar, except for the steep IMF model, which shifts the local peak from ~20 $M_{☉}$ to ~25 $M_{☉}$. We note that this result might be sensitive to the mass ranges chosen for the fits in *Fryer* et al. (2012).

The evolution with redshift of the merger rate per unit BH mass is shown in Fig. 6 for the *WWp* model (solid lines) and the *Fryer* model (dashed line). We note that Fig. 5 corresponds to the $z = 0$ axis in Fig. 6, and we show only four mass bins to simplify the
plot. First, it can be seen that the merger rate of lower-mass BHs attains its maximal value at around the peak of the SFR, which is not surprising given that the chosen delay time distribution prefers very short delay times. The dominant contribution to the overall merger rate is from small BH masses, which is the case in all the models we considered. Events like GW150914 are therefore not expected to constitute the majority of the merger events. Note, however, that in view of the sensitivity of the LIGO detector, the majority of the observed events might well be similar to GW150914 (see e.g. fig. 3 in Belczyznski et al. 2016). In all the models we considered the dominant contribution was from BBHs just above our chosen limit between neutron stars and BBHs, as a consequence of the IMF which peaks at low masses. The WWp model discourages the formation of BHs with masses above ~20 M⊙, which only occurs at low metallicity. This explains the peak of the mergers of 30 M⊙ BHs which occurs at z ~ 5 in this model. For comparison we also show the merger rate of 30 M⊙ BHs in the Fryer model, in which case they continue to be formed and merge up to low redshift, with the peak occurring at z ~ 2. The difference between these two models in their predictions for the merger rate of 30 M⊙ BHs, shown in Fig. 5, is visible here for the whole redshift range.

We now turn to consider the implications of our models and the differences in the obtained BBH mass distributions for the stochastic gravitational wave background.

## 4 STOCHASTIC GRAVITATIONAL WAVE BACKGROUND

The stochastic background of gravitational waves is usually expressed in terms of the dimensionless density parameter:

$$ \Omega_{\text{gw}}(f) = \frac{1}{\rho_c} \frac{d\rho_{\text{gw}}}{d\ln f} \qquad (5) $$

where $\rho_{gw}$ is the gravitational energy density and $\rho_c$ is the critical energy density of the Universe. The density parameter is given by (Regimbau 2011):

$$ \Omega_{\text{gw}}(f_0) = \frac{8\pi G}{3c^2 H_0^3} f_0 \int d\theta p(\theta) \int dz \frac{R_{\text{source}}(z, \theta)}{(1+z)E_V(z)} \frac{dE_{\text{gw}}(\theta)}{df} \qquad (6) $$

where $f_0$ is the observed frequency, $f = (1+z)f_0$ is the frequency at emission, $p(\theta)$ is the distribution of source parameters $\theta$ (such as the source type, binary orbital parameters, etc.), $R_{\text{source}}(z, \theta)$ is the source rate density, $E_V(z) = \sqrt{\Omega_m(1+z)^3 + \Omega_\Lambda}$ is the cosmological volume parameter and $dE_{\text{gw}}/df$ is the gravitational spectral energy emitted.

The first contribution we consider is from the coalescence of two BHs and reflects the merger rate of BBH computed above. The merger event can be decomposed into three phases: inspiral, merger and ringdown, for which the spectrum is given approximately by (Zhu et al. 2011):

$$ \frac{dE_{\text{GW}}}{df_e} = \frac{(G\pi)^{2/3} M_c^{5/3}}{3} \begin{cases} f_e^{-1/3}, & f_e \le f_1 \\ \omega_1 f_e^{2/3}, & f_1 < f_e < f_2 \\ \omega_2 \left( \frac{f_e}{1 + \left( \frac{f_e}{\sigma/2} \right)^2} \right)^2, & f_2 \le f_e \le f_3 \end{cases} \qquad (7) $$

where $M_c = (m_1 m_2)^{3/5}/(m_1 + m_2)^{1/5}$ is the chirp mass. The set of parameters ($f_1, f_2, f_3, \sigma$), where $f_1, f_2$ correspond to the end of the inspiral and merger phases, respectively, is taken from Ajith et al.

(2008) for the case of non-spinning BHs for each set of masses (which we assume to be always equal). The constants $\omega_1 = f_1^{-1}$ and $\omega_2 = f_1^{-1} f_2^{-4/3}$ are chosen to make $dE_{\text{gw}}/df$ continuous.

The second contribution is from a collapse of a single star which reflects the BH birth rate. We assume, following Crocker et al. (2015), that most of the energy is dissipated via the ringdown of the $l = 2$ dominant quasi-normal mode whose frequency is given by (Echeverria 1989):

$$ f_* = \frac{\Delta(a)}{m_{\text{bh}}} \qquad (8) $$

where $m_{bh}$ is the mass of the BH and

$$ \Delta(a) = \frac{c^3}{2\pi G} [1 - 0.63(1 - a)^{0.3}] \qquad (9) $$

is a function of the dimensionless spin factor $a$. The energy spectrum of a single source is then given by:

$$ \frac{dE_{\text{GW}}}{df} = \epsilon m_{\text{bh}} c^2 \delta(f - f_*) \qquad (10) $$

where $\epsilon$ is the efficiency of GW production. We note that there could be other contributions which depend on the details of the post-core-bounce evolution of the SN (see e.g. Müller, Janka & Marek 2013; Ott et al. 2013).

Then the density parameter of gravitational waves can be calculated by plugging equations (8)–(10) into equation (6) and is given by:

$$ \Omega_{\text{gw}}(f_0) = \epsilon \frac{8\pi G}{3H_0^3} \int \frac{R_{\text{birth}}(z', m_{\text{bh}})}{(1+z')E_V(z')} m_{\text{bh}} dm_{\text{bh}} \qquad (11) $$

where $z'$ is a function of $m_{bh}$. The rate of formation of BHs per unit time, per unit volume and per unit mass $R_{\text{birth}}(z', m_{bh})$ is taken from our model described above where

$$ z' = \frac{\Delta(a)}{m_{\text{bh}} f_0} - 1 \qquad (12) $$

and we take the efficiency parameter $\epsilon = 10^{-5}$ from Crocker et al. (2015). We assume a constant spin parameter $a = 1$ for all the BHs. The resulting stochastic gravitational wave background is depicted in Fig. 7 (blue solids curve) calculated with the Fryer model, in particular using the BBH mass spectrum shown partly in Figs 5 and 6. The light blue shaded region corresponds to the upper and lower limits of the local merger rate, to which we normalize the total merger rate [the factor $N$ in equation (4)]: $1.02^{+1.98}_{-0.79} 10^{-7}$ Mpc$^{-3}$ yr$^{-1}$ (Abbott et al. 2016c). The solid red curve is the Fiducial model of The LIGO Scientific Collaboration & the Virgo Collaboration (2016) (taken from their fig. 1) calculated using $dE_{\text{GW}}/df_e$ from Ajith (2011) and assuming the merger rate is proportional to the cosmic SFR below metallicity $0.5Z_\odot$ and the same time delay distribution as used here. Note that in this case the normalization was taken to be $1.6^{+3.8}_{-1.3} 10^{-8}$ Mpc$^{-3}$ (Abbott et al. 2016c), since they assumed that all BBHs have masses identical to GW150914. The blue dashed line shows our calculation where we also assumed the same masses as in GW150914 for all the merger events, in which case our model coincides with the LIGO Fiducial model. The main difference between these models is the SFR: while we use our Fiducial SFR model, the red solid curve is calculated using the GRB-based model (see Fig. 1). The difference between these models has a negligible effect on the spectrum of the stochastic gravitational wave background. Another difference is that in our Fryer model with fixed mass the BBH birthrate is in practice proportional to the SFR (since all the BHs are assumed to be born with the same
<figure><img src="image_8.png" /></figure>

**Figure 7.** Gravitational wave stochastic background from merger events. The brown curve and the light orange area are taken from The LIGO Scientific Collaboration & the Virgo Collaboration (2016), where it was assumed that all the BBH have masses identical to GW150914 and the waveforms from Ajith (2011) were used. Dashed blue line is our calculation with the *Fryer* model under the same assumption and using the waveforms from Ajith et al. (2008) which produces a nearly identical result. The blue solid line corresponds to the *Fryer* model with the full mass spectrum taken into account. Note that the amplitude in this case is increased and there is additional power at high frequencies which is due to low-mass BHs. The light blue shaded area corresponds to the *Fryer* model taking into account the uncertainty in the measured local merging rate to which our calculation is normalized. The set of black curves shows our calculation using the waveforms from Ajith et al. (2011) which include PN corrections: dashed for fixed mass, solid for a mass distribution and dashed with squares for a mass distribution and a uniform spin distribution.

mass, the mass distributions shown in Figs 5 and 6 are irrelevant). The LIGO *Fiducial* model, on the other hand, assumes that the birthrate is proportional to the SFR at metallicity below $Z_{\odot}/2$, which, however, is proportional to the total SFR for most of the cosmic history (see Fig. 3). Since the same overall normalization to the local observed rate is used for both the red solid and black dashed curves, this difference is also unimportant.

On the other hand, it can be seen that if the whole mass distribution is taken into account the spectrum of the stochastic background shifts to higher values, since the BBH population is dominated by low-mass binaries. It is therefore important to correctly account for the mass distribution of BBHs discussed above. We note that there is also significant uncertainty due to the poorly constrained time delay distribution (see the discussion in The LIGO Scientific Collaboration & the Virgo Collaboration 2016) whose detailed treatment we leave to future work.

We also explore the effect of BH spins and post-Newtonian (PN) corrections by taking the spectrum from Ajith et al. (2011) (black lines in Fig. 7). Ajith et al. (2011) matched a PN description of the inspiral phase to a set of numerical relativity simulations to obtain an analytical inspiral-merger-ringdown waveform family. While PN corrections to the waveform are crucial in parameter estimation of merger events (e.g. Blanchet 2006), their contribution to the stochastic background is small in the case of both fixed masses (dashed black line) and a distribution of BH masses (solid black line). Moreover, adding a uniform distribution of BH spins in the range [-0.85, 0.85] (the range of validity of the models of Ajith et al. (2011); black dashed line with squares) does not significantly affect the result.

<figure><img src="image_9.png" /></figure>

**Figure 8.** Gravitational wave stochastic background from single-star collapse (dashed) and from merger events (solid). Blue curves correspond to the *Fryer* model (the blue solid line is the same as in the previous figure), whereas the pink lines correspond to the fiducial *WWp* model. As expected, the *Fryer* models predict higher amplitude and lower frequency of the peak owing to the larger BH masses. Note that for the *Fryer* models the contributions from the single and binary populations are separated in frequency and can be measured with experiments that are sensitive at $f \gtrsim 400$ Hz (although note the caveats discussed in the text). *WWp* models with a *Kinugawa*-like prescription for BH formation from metal-poor stars are shown in yellow and green for a metallicity limit of $Z_{\text{limit}} = 0.001Z_{\odot}$ and $Z_{\text{limit}} = 0.01Z_{\odot}$, respectively. For comparison, we also show one of the models from Crocker et al. (2015) (black dashed line, see text). All of the *single* models assume a dimensionless spin of $a = 1$.

Fig. 8 includes the contribution from SN collapse events (dashed lines) as well as from mergers (solid lines) and compare several of the models discussed above. Interestingly, in our fiducial *Fryer* model (blue) the two contributions are somewhat separated in frequency. Even though the amplitude of the contribution from SN collapse is at present highly uncertain, this separation might be observable with experiments sensitive in the $f \gtrsim 400$ Hz frequency domain, and when the rate of observed events becomes sufficiently large so as to reduce the uncertainty bands. Such an observation will provide a clear handle on the two different populations – single and binary – of BHs. We note, however, that there may be additional contributions in this frequency range, such as from binary neutron star mergers not discussed here. Furthermore, the mechanism of generation of gravitational waves during SN collapse is highly uncertain (see e.g. Müller et al. 2013; Ott et al. 2013).

The stochastic background from mergers and SN collapse in the fiducial *WWp* is shown by the pink solid and dashed lines in Fig. 8, respectively. The peak of the merger contribution is shifted towards higher frequency with respect to the *Fryer* model due to the lower BH masses in *WWp*. We also plot the background from SN collapse in *WWp* models that use a *Kinugawa*-like prescription for metal-poor stars (yellow and green dashed lines). This prescription produces a peak at around 150 Hz, but its amplitude is very small owing to the small number density of such SNe. For comparison we also show (black dashed curve) the calculation of the background from SN collapse from Crocker et al. (2015), where they assumed that 50 per cent of the progenitor mass goes into the BH (their model 1). We note that all the *Fryer* models (i.e. with steep IMF, PopIII stars, GRB-based SFR) produce very similar GW background and we do not expect the differences to be observable.
## 5 DISCUSSION

This article explores different scenarios of BH formation in the context of cosmic chemical evolution models which are consistent with measured SFRs, optical depth to reionization and metallicity evolution of the ISM. We have shown that the analytic model of Fryer et al. (2012) is much more efficient in producing heavy BBH, such as GW150914 recently observed by Advanced LIGO, than the piston-driven model of Woosley & Weaver (1995). While the sensitivity of the BBH mass distribution to the SN collapse model is not surprising (e.g. Fryer et al. 2012), our approach provides a convenient cosmological framework for the analysis of future observations.

We investigated various SFR prescriptions and found that models that produce large amounts of stars at high redshifts, such as GRB-based SFR and bimodal SFR which includes a contribution from PopIII stars, are not favourable to heavy BH production because of the accompanying rise in metallicity. On the other hand, models with very steep IMF which produce few massive stars also result in low cosmic metallicity which leads to higher BH masses. These results demonstrate the importance of self-consistent modelling. After this paper was submitted the Planck collaboration released new constraints on the optical depth to reionization. We intend to study the consequences of these new results in future work.

The role of PopIII stars in producing heavy BBH which merge within the age of the Universe is currently debated. Our results support the conclusion of Hartwig et al. (2016) in that the contribution of PopIII remnants is sub-dominant in realistic BBH formation models. In particular, we find that the inclusion of special prescriptions of BH formation from metal-poor stars, such as the model proposed by Kinugawa et al. (2014), does not affect the mass distribution of merging BBHs in realistic scenarios, such as our set of Fryer models. Moreover, we calculated the stochastic gravitational wave background and found that the contribution from PopIII stars is negligible in the entire frequency range we explored ($f \sim 10-1000$ Hz).

The analysis presented in this paper is far from being exhaustive and we plan to explore other BBH formation models in future work. We expect that future detections of merging BBH will enable us to discriminate between the different models of SN collapse and BBH evolution.

In this paper we normalized the total merger rate to the observed single event and used a universal distribution of delay times from stellar birth to merger. In reality, the delay time depends on the orbital parameters of each binary, and the overall normalization also depends on the binary fraction. We plan to include a realistic model of binary orbital parameters in a forthcoming paper.

Finally, we plan to apply the approach presented in this paper to a full galaxy evolution model in order to calculate the anisotropy of the stochastic gravitational wave background and the cross-correlation between the optical signal from galaxies and gravitational waves from BBH mergers.

## ACKNOWLEDGEMENTS

We thank Vuk Mandic and Chris Belczynski for their constructive comments on the manuscript and the referee Tania Regimbau for very helpful suggestions. The work of ID and JS was supported by the ERC Project No. 267117 (DARK) hosted by Université Pierre et Marie Curie (UPMC) - Paris 6, PI J. Silk. JS acknowledges the support of the JHU by NSF grant OIA-1124403. The work of KAO was supported in part by DOE grant DE-SC0011842 at the

University of Minnesota. This work has been carried out at the ILP LABEX (under reference ANR-10-LABX-63) supported by French state funds managed by the ANR within the Investissements d'Avenir programme under reference ANR-11-IDEX-0004-02.

## REFERENCES

* Abbott B. P. et al., 2016a, ApJ, 818, L22

* Abbott B. P. et al., 2016b, Phys. Rev. Lett., 116, 061102

* Abbott B. P. et al., 2016c, preprint (arXiv:1602.03842)

* Ajith P., 2011, Phys. Rev. D, 84, 084037

* Ajith P. et al., 2008, Phys. Rev. D, 77, 104017

* Ajith P. et al., 2011, Phys. Rev. Lett., 106, 241101

* Behroozi P. S., Silk J., 2015, ApJ, 799, 32

* Behroozi P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, 57

* Belczynski K., Kalogera V., Bulik T., 2002, ApJ, 572, 407

* Belczynski K., Bulik T., Rudak B., 2004, ApJ, 608, L45

* Belczynski K., Bulik T., Fryer C. L., Ruiter A., Valsecchi F., Vink J. S., Hurley J. R., 2010a, ApJ, 714, 1217

* Belczynski K., Dominik M., Bulik T., O'Shaughnessy R., Fryer C., Holz D. E., 2010b, ApJ, 715, L138

* Belczynski K., Holz D. E., Bulik T., O'Shaughnessy R., 2016, Nature, 534, 512

* Blanchet L., 2006, Living Rev. Relativ., 9

* Bond J. R., Carr B. J., 1984, MNRAS, 207, 585

* Bouwens R. J. et al., 2015, ApJ, 803, 34

* Chabrier G., Hennebelle P., Charlot S., 2014, ApJ, 796, 75

* Clausen D., Piro A. L., Ott C. D., 2015, ApJ, 799, 190

* Crocker K., Mandic V., Regimbau T., Belczynski K., Gladysz W., Olive K., Prestegard T., Vangioni E., 2015, Phys. Rev. D, 92, 063005

* Daigne F., Olive K. A., Vangioni-Flam E., Silk J., Audouze J., 2004, ApJ, 617, 693

* Daigne F., Olive K. A., Silk J., Stoehr F., Vangioni E., 2006, ApJ, 647, 773

* de Mink S. E., Cantiello M., Langer N., Pols O. R., Brott I., Yoon S.-C., 2009, A&A, 497, 243

* Echeverria F., 1989, Phys. Rev. D, 40, 3194

* Eldridge J. J., Stanway E. R., 2016, preprint (arXiv:1602.03790)

* Ertl T., Janka H.-T., Woosley S. E., Sukhbold T., Ugliano M., 2016, ApJ, 818, 124

* Fryer C. L., 2006, New Astron. Rev., 50, 492

* Fryer C. L., Young P. A., 2007, ApJ, 659, 1438

* Fryer C. L., Belczynski K., Wiktorowicz G., Dominik M., Kalogera V., Holz D. E., 2012, ApJ, 749, 91

* Georgy C., Meynet G., Maeder A., 2011, A&A, 527, A52

* Greif T. H., Bromm V., 2006, MNRAS, 373, 128

* Hartwig T., Volonteri M., Bromm V., Klessen R. S., Barausse E., Magg M., Stacy A., 2016, preprint (arXiv:1603.05655)

* Inayoshi K., Kashiyama K., Visbal E., Haiman Z., 2016, preprint (arXiv:1603.06921)

* Janka H.-T., 2012, Annu. Rev. Nucl. Part. Sci., 62, 407

* Jones S., Hirschi R., Pignatari M., Heger A., Georgy C., Nishimura N., Fryer C., Herwig F., 2015, MNRAS, 447, 3115

* Kinugawa T., Inayoshi K., Hotokezaka K., Nakauchi D., Nakamura T., 2014, MNRAS, 442, 2963

* Kistler M. D., Yuksel H., Hopkins A. M., 2013, preprint (arXiv:1305.1630)

* Kochanek C. S., 2014, ApJ, 785, 28

* Kowalska I., Bulik T., Belczynski K., 2012, A&A, 541, A120

* Kroupa P., 2001, MNRAS, 322, 231

* Lentz E. J., Mezzacappa A., Messer O. E. B., Liebendörfer M., Hix W. R., Bruenn S. W., 2012, ApJ, 747, 73

* Loeb A., 2016, ApJ, 819, L21

* Mandel I., de Mink S. E., 2016, MNRAS, 458, 2634

* Marchant P., Langer N., Podsiadlowski P., Tauris T. M., Moriya T. J., 2016, A&A, 588, A50

* Marigo P., Girardi L., Chiosi C., Wood P. R., 2001, A&A, 371, 152

* Matthee J., Sobral D., Best P., Khostovan A. A., Oteo I., Bouwens R., Röttgering H., 2016, preprint (arXiv:1605.08782)
Meynet G., Maeder A., 2007, A&A, 464, L11

Meynet G. et al., 2015, A&A, 575, A60

Müller B., Janka H.-T., Marek A., 2012, ApJ, 756, 84

Müller B., Janka H.-T., Marek A., 2013, ApJ, 766, 43

O'Connor E., Ott C. D., 2011, ApJ, 730, 70

Oesch P. A., Bouwens R. J., Illingworth G. D., Franx M., Ammons S. M., van Dokkum P. G., Trenti M., Labbé I., 2015, ApJ, 808, 104

Ott C. D. et al., 2013, ApJ, 768, 115

Pejcha O., Thompson T. A., 2015, ApJ, 801, 90

Peters P. C., 1964, Phys. Rev., 136, 1224

Planck Collaboration XIII et al., 2015, preprint (arXiv:1502.01589)

Planck Collaboration XLVII et al., 2016a, preprint (arXiv:1605.03507)

Planck Collaboration XLVI et al., 2016b, preprint (arXiv:1605.02985)

Rafelski M., Wolfe A. M., Prochaska J. X., Neeleman M., Mendez A. J., 2012, ApJ, 755, 89

Regimbau T., 2011, Res. Astron. Astrophys., 11, 369

Robertson B. E., Ellis R. S., 2012, ApJ, 744, 95

Rodriguez C. L., Chatterjee S., Rasio F. A., 2016, Phys. Rev. D, 93, 084029

Rollinde E., Vangioni E., Maurin D., Olive K. A., Daigne F., Silk J., Vincent F. H., 2009, MNRAS, 398, 1782

Salpeter E. E., 1955, ApJ, 121, 161

Schaerer D., 2002, A&A, 382, 28

Spera M., Mapelli M., Bressan A., 2015, MNRAS, 451, 4086

Springel V., Hernquist L., 2003, MNRAS, 339, 312

Szécsi D., Langer N., Yoon S.-C., Sanyal D., de Mink S., Evans C. J., Dermine T., 2015, A&A, 581, A15

The LIGO Scientific Collaboration the Virgo Collaboration 2016, Phys. Rev. Lett., 116, 131102

Trenti M., Perna R., Tacchella S., 2013, ApJ, 773, L22

Ugliano M., Janka H.-T., Marek A., Arcones A., 2012, ApJ, 757, 69

Vangioni E., Olive K. A., Prestegard T., Silk J., Petitjean P., Mandic V., 2015, MNRAS, 447, 2575

Vink J. S., 2008, New Astron. Rev., 52, 419

Wang F. Y., 2013, A&A, 556, A90

Woosley S. E., 2016, ApJL, 824, 10

Woosley S. E., Weaver T. A., 1995, ApJS, 101, 181

Zhang W., Woosley S. E., Heger A., 2008, ApJ, 679, 639

Zhu X.-J., Howell E., Regimbau T., Blair D., Zhu Z.-H., 2011, ApJ, 739, 86

Ziosi B. M., Mapelli M., Branchesi M., Tormen G., 2014, MNRAS, 441, 3703

This paper has been typeset from a *TEX/LATEX* file prepared by the author.