Issue 
A&A
Volume 619, November 2018



Article Number  A77  
Number of page(s)  19  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/201833025  
Published online  08 November 2018 
Impact of intercorrelated initial binary parameters on double black hole and neutron star mergers
^{1}
Department of Astrophysics/IMAPP, Radboud University, PO Box 9010, 6500 GL Nijmegen, The Netherlands
email: j.klencki@astro.ru.nl
^{2}
Steward Observatory, University of Arizona, 933 N. Cherry Ave., AZ 85721 Tucson, USA
^{3}
Astronomical Observatory, Warsaw University, Ujazdowskie 4, 00478 Warsaw, Poland
^{4}
Enrico Fermi Institute, Department of Physics, Department of Astronomy and Astrophysics, and Kavli Institute for Cosmological Physics, University of Chicago, IL 60637 Chicago, USA
^{5}
Kavli Institute for Particle Astrophysics & Cosmology and Physics Department, Stanford University, Ujazdowskie 4, CA 94305 Stanford, USA
^{6}
Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, ul. Bartycka 18, 00716 Warsaw, Poland
Received:
14
May
2018
Accepted:
7
September
2018
The distributions of the initial mainsequence binary parameters are one of the key ingredients in obtaining evolutionary predictions for compact binary (BH–BH/BH–NS/NS–NS) merger rates. Until now, such calculations were done under the assumption that initial binary parameter distributions were independent. For the first time, we implement empirically derived intercorrelated distributions of initial binary parameters primary mass (M_{1}), mass ratio (q), orbital period (P), and eccentricity (e). Unexpectedly, the introduction of intercorrelated initial binary parameters leads to only a small decrease in the predicted merger rates by a factor of ≲2–3 relative to the previously used noncorrelated initial distributions. The formation of compact object mergers in the isolated classical binary evolution favours initial binaries with stars of comparable masses (q ≈ 0.5–1) at intermediate orbital periods (log P (days) = 2–4). New distributions slightly shift the mass ratios towards lower values with respect to the previously used flat q distribution, which is the dominant effect decreasing the rates. New orbital periods (∼1.3 more initial systems within log P (days) = 2–4), together with new eccentricities (higher), only negligibly increase the number of progenitors of compact binary mergers. Additionally, we discuss the uncertainty of merger rate predictions associated with possible variations of the massivestar initial mass function (IMF). We argue that evolutionary calculations should be normalized to a star formation rate (SFR) that is obtained from the observed amount of UV light at wavelength 1500 Å (an SFR indicator). In this case, contrary to recent reports, the uncertainty of the IMF does not affect the rates by more than a factor of ∼2. Any change to the IMF slope for massive stars requires a change of SFR in a way that counteracts the impact of IMF variations on compact object merger rates. In contrast, we suggest that the uncertainty in cosmic SFR at low metallicity can be a significant factor at play.
Key words: binaries: general / stars: black holes / stars: neutron / stars: massive / gravitational waves
© ESO 2018
1. Introduction
The Laser Interferometer Gravitationalwave Observatory (LIGO) began its first upgraded observational run (O1) in September 2015; the first ever detection of a gravitational wave signal from a binary black hole (BH–BH) coalescence came shortly afterwards: i.e. GW150914 (Abbott et al. 2016b,c). Since then, four additional BH–BH mergers and, most recently, a double neutron star (NS–NS) merger, were detected and reported to the community: i.e. GW151226 (BH–BH; Abbott et al. 2016a), GW170104 (BH–BH; Abbott et al. 2017a), GW170608 (BH–BH; Abbott et al. 2017b), GW170814 (BH–BH; Abbott et al. 2017c), and GW170817 (NS–NS; Abbott et al. 2017d). The last three events were observed during the second observational run (O2); the Advanced Virgo detector (Acernese et al. 2015) joined the run on August 1, 2017 and contributed to the analysis of GW170814 and GW170817.
The LIGO discovery marks the beginning of the gravitationalwave era. Detections of the coalescence signals from compact binary mergers (CBM) are of utmost astrophysical significance as, among other applications, they will constrain potential formation scenarios, stellar evolution models, and other assumptions associated with theoretical predictions (e.g. Stevenson et al. 2015; O’Shaughnessy et al. 2017; Barrett et al. 2018). Various formation scenarios for CBM have been proposed. Those most widely discussed include the isolated binary evolution channel involving a common envelope (CE) phase (Eldridge & Stanway 2016; Belczynski et al. 2016b; Chruslinska et al. 2018; Mapelli et al. 2017; Giacobbo et al. 2018; Stevenson et al. 2017) or stable mass transfer (van den Heuvel et al. 2017), isolated evolution of field triples (e.g. Antonini & Rasio 2016), dynamical evolution in dense stellar environments such as globular clusters (GCs; Rodriguez et al. 2016a,b; Askar et al. 2017; Park et al. 2017), nuclear star clusters (Miller & Lauburg 2009; Antonini & Rasio 2016) or even discs of active galactic nuclei (Stone et al. 2017), and the formation of compact objects in very close and tidally locked binaries through chemically homogeneous evolution (de Mink & Mandel 2016; Mandel & de Mink 2016; Marchant et al. 2016). We note that while it is still possible to distinguish between different types of mergers (i.e. BH–BH, BH–NS, or NS–NS) in a modelindependent way (Mandel et al. 2015, 2017) using their gravitationalwave signatures, we still lack a reliable way to determine the formation channel of an observed merger. This is especially true in case of the LIGO/Virgo NS–NS merger (Belczynski et al. 2018). In the case of BH–BH mergers there is some hope connected to the measurement of the BH–BH spinorbit misalignments (Stevenson et al. 2017; Farr et al. 2017, 2018), although this may not work if the BH spins are intrinsically very small (Belczynski et al. 2017).
Regardless of the formation scenario, the theoretical predictions for the compact binary merger rates are burdened with significant uncertainties due to numerous assumptions and models with poorly constrained parameters, for example the infamous CE phase (Dominik et al. 2012) or the BH/NS natal kicks (Repetto & Nelemans 2015; Belczynski et al. 2016c). One of the key ingredients in the calculations are the initial conditions for the simulations: the birth properties of stellar clusters in the case of dynamical scenarios and the characteristics of primordial binaries in the case of isolated binary evolution channels.
Recently, de Mink & Belczynski (2015) incorporated the updated primordial binary parameter distributions obtained by Sana et al. (2012) from spectroscopic measurements of massive Otype stars in very young (∼2 Myr) open star clusters and associations. The updated distributions included a much stronger bias towards close binary orbits with respect to the previously adopted Öpik’s law, i.e. a flat in logarithm distribution (Öpik 1924; Abt 1983). Intuitively, this change should favour interacting binaries (including those undergoing the CE phase) and possibly cause a notable increase in merger rates. However, de Mink & Belczynski (2015) found only a very small increase of less than a factor of 2. This is because the distributions obtained by Sana et al. (2012) also show a heavy bias towards low eccentricities with respect to the thermalequilibrium distribution of Heggie (1975); this distribution results in a nearly unchanged distribution of periastron separations, which is the essential separation regulating the onset of mass transfer.
The sample of binaries observed by Sana et al. (2012) suffers from a significant limitation: it is restricted to systems with log P (days) < 3.5 (spectroscopically detectable binaries) and dominated by very shortperiod orbits with P < 20 days (hence the huge fraction of circularized systems). Since the BH–BH mergers can originate from primordial binaries of up to log P ∼ 5.5 (de Mink & Belczynski 2015) the binary parameter distributions obtained by Sana et al. (2012) need to be extrapolated to longer periods, which automatically assumes no intrinsic correlations between parameters. However, the joint probability density function cannot be decomposed into independent distribution functions of the individual parameters, i.e. f(M_{1}, q, P, e) ≠ f(M_{1})f(q)f(P)f(e). Observational studies have hinted at probable correlations (Abt et al. 1990; Duchêne & Kraus 2013), but hitherto the selection biases have been too large to accurately quantify the intrinsic interrelations. Recently, Moe & Di Stefano (2017; hereafter MD17) analysed more than 20 surveys of massive binary stars, corrected for their respective selection effects, combined the data in a homogeneous manner, and fit analytic functions to the corrected distributions. These authors confirm that many of the physical binary star parameters are indeed correlated at a statistically significant level.
de Mink & Belczynski (2015) de Mink & Belczynski (2015) concluded that the most significant variations of merger rates associated with the initial parameters are due to uncertainties of the initial mass function (IMF) powerlaw slope for massive stars (merger rates going up and down by a factor of six in the case of BH–BH). Cosmological calculations of the BH–BH merger rates (e.g. Dominik et al. 2013, 2015; Belczynski et al. 2016b; Kruckow et al. 2018; Mapelli et al. 2017) are performed based on the assumption of a universal IMF across the cosmic time. Often assumed is a socalled “canonical” IMF, which is a multipart powerlaw distribution dN/dM = ξ(M) ∝ M^{−αi}, where α_{1} = 1.3 for M/ M_{⊙} ∈ [0.08,0.5] and α_{2} = α_{3} = 2.3 for M/ M_{⊙} ∈ [0.5,1.0] and [1.0,150.0], respectively (Kroupa 2001).
Although clear evidence for strong IMF variations with environmental conditions is still lacking, there are a growing number of results hinting at departures from the IMF universality (see reviews by Bastian et al. 2010; Kroupa et al. 2013). Notably, a recent spectroscopic survey of massive stars in the 30 Doradus starforming region in the Large Magellanic Cloud has led to a discovery of an excess of stars with masses above 30 M_{⊙} with respect to the canonical IMF (Schneider et al. 2018); the bestfitted single powerlaw exponent for M > 1 M_{⊙} is (although see a technical comment on the data analysis from Farr & Mandel 2018, suggesting somewhat larger values for α_{3}).
The unknowns associated with the massivestar IMF are often considered to be one of the significant contributors to uncertainty of compact binary merger rates calculated based on population synthesis. In this work, we argue that by normalizing the simulated stellar population to the total amount of farUV light that it emits, rather than to its total mass, one can significantly reduce the uncertainty associated with possible variations of the IMF slope in different environments; see Sect. 6.2. As an example of such a variation and its impact on merger rate calculations, we study in detail the case of a possible correlation between the massivestar IMF slope and metallicity Z.
With the exception of abovementioned results of Schneider et al. (2018), numerous observations of OB associations and clusters in the Local Group did not reveal any significant deviations from the canonical IMF slope for massive stars α_{3} ≈ 2.3 (Massey 2003). These included surveys of the Milky Way (Z ≈ Z_{⊙} = 0.02; Daflon & Cunha 2004)^{1}, the Small Magellanic Cloud (Z ≈ 0.004; Korn et al. 2000) and the Large Magellanic Cloud (Z ≈ 0.008; Korn et al. 2000), indicating that for Z ≥ 0.004 ([Fe/H] ≳ −0.7) the highmass end of IMF does not significantly depend on metallicity. Moe & Di Stefano (2013) showed that the same is true for the parameters of close binaries with massive stars. Even at solar metallicity, the discs of massive protostars are already highly prone to gravitational instability and fragmentation, explaining why the close binary fraction of massive stars is so large (Kratter & Matzner 2006; Kratter et al. 2008, 2010; Moe & Di Stefano 2017). Further reducing the metallicity can therefore only marginally increase the close binary fraction of massive stars (Tanaka & Omukai 2014; Moe et al. 2018). However, as we show, the vast majority of BH–BH mergers evolving from the CE channel are expected to originate from Z < 0.004 environments, for which there is no direct observational evidence for the persistence of the canonical IMF.
From a theoretical point of view, both the Jeansmass formalism (Jeans 1902; Larson 1998; Bate & Bonnell 2005; Bonnell et al. 2006) and the model of stellar formation as a selfregulated balance between the rate of accretion from the protostellar envelope and the radiative feedback from the forming star (Adams & Fatuzzo 1996) predict that at a certain sufficiently low metallicity the IMF becomes top heavy^{2}. In the first case, the prediction comes from the fact that the Jeans mass has to be larger at lower Z owing to less effective cooling of the protostellar cloud (M_{J} ∝ ρ^{−1/2} T^{3/2}). The radiative feedback, on the other hand, is also metallicity dependent since photons couple less effectively to gas of lower metallicity. Hydrodynamical simulations of the formation of Population III stars (Z/ Z_{⊙} < 10^{−6}) demonstrate that the first generation of stars were almost exclusively massive OBtype mainsequence (MS) stars (Bromm et al. 1999; Yoshida et al. 2006; Clark et al. 2011). However, by redshift z ≈ 10, the mean metallicity of actively forming stars was Z/ Z_{⊙} ≈ 10^{−3} (Tornatore et al. 2007; Madau & Dickinson 2014). For such Population II stars with intermediate metallicities, the IMF is expected to be only moderately top heavy compared to the canonical IMF (Fang & Cen 2004; Daigne et al. 2006; Greif et al. 2008).
From the observational side, there is indirect evidence for topheavy IMF variations at cosmological times. The relative paucity of metalpoor Gdwarf stars in the Milky Way (the socalled Gdwarf problem; Pagel 1989) can be explained by applying an IMF, which is increasingly deficient in lowmass stars the earlier the star formation took place (Larson 1998). By modelling the abundances of Lymanbreak and submillimetre galaxies in the ΛCDM cosmology, Baugh et al. (2005) found that the observations can be reproduced only if some episodes of star formation with a topheavy IMF were also present in addition to the canonical IMF. Nagashima et al. (2005) used the same model with two modes of star formation to explain the elemental abundances in the intracluster medium of galaxy clusters. Wilkins et al. (2008) pointed out that the local stellar mass density is significantly lower than the value obtained from integrating the cosmic star formation history with a Salpeter IMF (single slope, α = 2.35; Salpeter 1955). At low redshifts (z < 0.5) they manage to match the observations using a slightly topheavy (α_{3} = 2.15) IMF. For higher redshifts, however, Wilkins et al. (2008) argued that no universal IMF can reproduce the measured stellar mass densities.
The only quantitative calibration of the relation between the highmass IMF and metallicity obtained thus far has been made possible thanks to the observations of some GCs in the Milky Way. Djorgovski et al. (1993) noticed that the higher metallicity GCs tend to have a bottomlight present day mass function (i.e. a deficit of lowmass stars). De Marchi et al. (2007) later found that these GCs are also the least concentrated sources in their observed sample. Such a trend contradicts basic cluster evolution theory – GCs are expected to be losing lowmass stars from their outer regions owing to their collapse into dense, highly concentrated clusters. An explanation was put forth by Marks et al. (2008), who proposed that lowmass stars could be unbound from masssegregated GCs during expulsion of their residual gas. This introduces a dependence on metallicity, since the process of gas expulsion is expected to be enhanced in metalrich environments (Marks & Kroupa 2010). Finally, Marks et al. (2012) concluded that in order to provide enough radiative feedback to expel the residual gas and match the characteristics of the observed GCs, their IMF had to be top heavy with the value of α_{3} decreasing with cluster metallicity.
It should be noted, however, that the model of residual gas expulsion applied by Marks et al. (2012) relies on simplified assumptions concerning the radiative feedback. First, the amount of energy deposited from stars into the ISM is fitted to stellar models of only three different masses (35, 60, and 85 M_{⊙}; Baumgardt et al. 2008), and second there is no metallicity dependence. Additionally, Marks et al. (2012) assumed that the energy needed for residual gas removal is provided by the stellar winds and radiation only (e.g. no feedback from supernova taken into account) and that all the energy radiated by the stars is deposited into the ISM. Thus, the exact relation between the highmass IMF slope and metallicity remains highly uncertain.
The purpose of this study is twofold. First, we aim to incorporate the interrelated initial binary parameter distributions and multiplicity statistics obtained by MD17 to determine their effects on the predicted rates and properties of CBM. Second, we investigate the importance of possible IMF variations for the merger rate calculations using the Marks et al. (2012) calibration of the IMF dependency on metallicity as an example of such a variation. To achieve this, we perform comparative population synthesis where we use the works of de Mink & Belczynski (2015) and Belczynski et al. (2016b) as references (hereafter dMB15 and B16, respectively). In Sect. 2 we describe the new initial conditions and compare these with the previously used distributions. In Sect. 3 we describe our computational method. In Sect. 4 we compare the distributions of the initial parameters of double compact merger progenitors in our simulations. In Sect. 5 we present the impact of the incorporated changes on the cosmological merger rates. In Sect. 6 we discuss the metallicity distribution of BH–BH mergers in our simulations and the significance of the topheavy IMF in lowZ environments for the LIGO predictions. We conclude in Sect. 7.
2. Initial distributions
2.1. Intercorrelated initial binary parameters
In the present study, we account for intrinsic correlations between the initial binary star physical parameters. We used the distribution functions presented in Sect. 9 of MD17. The correlations in the MS binary initial conditions are thoroughly discussed in that paper, but we summarize in this section the main results pertinent to the formation of compact objects mergers. dMB15 found that the majority of compact object mergers derive from MS binaries with initial orbital periods log P (days) = 2–4. Although for our simulations we generated companions across all orbital periods according to the distribution functions provided in MD17, in this section we focus the discussion on the companion star properties across intermediate periods log P (days) = 2–4.
Binaries versus triples. dMB15 and B16 assumed all companions to massive stars (M > 10 M_{⊙}) were in binaries. In contrast, MD17 have modelled a companion star fraction and period distribution that includes true binaries as well as inner binaries and outer tertiaries in hierarchical triples. With increasing orbital period, the likelihood that a companion is an outer tertiary must increase. In fact, essentially all wide companions (log P > 5) to massive stars are the outer tertiaries in hierarchical triples/quadruples (Sana et al. 2014; MD17). We model the probability that a companion is in the inner binary in a manner that reproduces the multiplicity statistics derived in MD17. Specifically, for M_{1} = 8 M_{⊙}, we find the probability decreases from for log P = 2 to for log P = 4. For M_{1} = 30 M_{⊙}, the triple/quadruple star fraction is larger, and so we adopt for log P = 2 and for log P = 4. We illustrate in Fig. 1.
Companion star fraction. dMB15 and B16 implemented a binary star fraction such that F_{logP = 2–4;bin} ≈ 28% of massive primaries had binary star companions with periods log P (days) = 2–4 and mass ratios q = 0.1–1.0. The companion star fraction increases dramatically with primary mass (Abt et al. 1990; Raghavan et al. 2010; Sana et al. 2012; Chini et al. 2012; Duchêne & Kraus 2013; Moe & Di Stefano 2013, 2017). According to MD17, the intermediateperiod companion star fraction increases from F_{logP} = 2–4;comp = 0.46 for M_{1} = 8 M_{⊙} to F_{logP = 2–4;comp} = 0.66 for M_{1} = 30 M_{⊙}. Approximately twothirds of Otype primaries have a companion with log P (days) = 2–4 and q > 0.1. A significant fraction of these systems, especially those with longer periods log P (days) = 3.5–4.0, are actually outer tertiaries in triples (see above and Fig. 1).
Period distribution. dMB15 and B16 adopted a powerlaw binary period distribution f_{logP;bin} ∝ (log P)^{−0.55} motivated by spectroscopic observations of massive binaries (Sana et al. 2012). These works normalize this period distribution such that integration reproduces the binary fraction across log P (days) = 2–4. In the current study, we adopt a companion star period distribution f_{logP;comp}(M_{1}, P) based on the analytic fits in MD17 that are shown in the lower panel of their Fig. 37. We then derive the inner binary period distribution . In Fig. 1, we show f_{logP;bin} from dMB15 and B16, f_{logP;comp} for M_{1} = 8 M_{⊙} and 30 M_{⊙} taken directly from MD17, and f_{logP;bin} for M_{1} = 8 M_{⊙} and 30 M_{⊙} implemented in the current study. While companions in general are weighted towards longer periods, the inner binary period distribution is skewed towards shorter periods as implemented in dMB15 and B16 and found in Sana et al. (2012).
Fig. 1. Frequency f_{logP} of companions with q > 0.1 per decade of orbital period across log P (days) = 2–4. We compare the distribution f_{logP;comp} of all companions (binaries and tertiaries; dashed) to the distribution of only the inner binaries (solid), where is the probability that the companion is a member of the inner binary. The overall binary fraction and binary period distribution across log P (days) = 2–4 are similar between those adopted by dMB15 and B16 (green) and the values for M_{1} = 8 M_{⊙} (red) and M_{1} = 30 M_{⊙} (blue) primaries based on the multiplicity statistics presented in MD17. 

Open with DEXTER 
Eccentricity distribution. dMB15 and B16 incorporated a powerlaw eccentricity distribution p_{e} ∝ e^{η} with exponent η = −0.4 across the domain e = 0.0–0.9, as motivated by Sana et al. (2012). However, the Sana et al. (2012) sample of spectroscopic binaries are dominated by systems with P < 20 days, therefore the powerlaw slope η = –0.4 is only appropriate for such shortperiod systems. MD17 have found that eccentricities become weighted towards higher values with increasing orbital period, which they model with two parameters. First, MD17 have defined the domain of the eccentricity distribution across the interval e = 0.0–e_{max}, where e_{max} (P) is the maximum eccentricity possible without substantially filling the Roche lobes of the binary (Eq. (3) in MD17). For P = 5 days, the binaries must have e < e_{max} ≈ 0.5 to be initially nonRochelobe filling. Meanwhile, binaries with log P (days) = 2–4 can extend up to e_{max} ≈ 0.95. Second, MD17 have found the powerlaw slope η also increases with increasing period. For massive binaries, they fit η ≈ 0.8 for log P (days) = 2–4. We adopt the powerlaw slopes η(M_{1}, P) presented in MD17 across the domain e = 0–0.8e_{max}. At high eccentricities e = (0.8–1.0)e_{max}, we assume the probability distribution function turns over according to a decreasing linear function such that p_{e}(e = e_{max}) = 0. In Fig. 2, we compare the cumulative distribution function of the eccentricity distribution adopted in dMB15 and B16 to the updated distributions at log P (days) = 2 and 4. As expected, the eccentricity distribution based on the Sana et al. (2012) sample is skewed towards lower values.
Fig. 2. Cumulative distribution function of eccentricities. dMB15 and B16 adopted an eccentricity distribution p_{e} ∝ e^{−0.4} (green) based on a sample of massive spectroscopic binaries (Sana et al. 2012) that are dominated by shortperiod systems P < 20 days and therefore weighted towards lower eccentricities 〈e〉 = 0.3. MD17 (and references therein) found that massive binaries with intermediate (blue) and long (red) orbital periods are weighted towards higher eccentricities 〈e〉 = 0.6 and are sufficiently modelled by a powerlaw distribution p_{e} ∝ e^{0.8} with a turnover above e > 0.8. 

Open with DEXTER 
Mass ratio distribution. dMB15 and B16 adopted a uniform massratio distribution f_{q} = q^{0}. Again, this was based on a sample of massive spectroscopic binaries (Sana et al. 2012) that is dominated by shortperiod systems P < 20 days. Based on a series of observational evidence (long baseline interferometry, companions to Cepheids, eclipsing binaries, adaptive optics, Hubble imaging, etc.), MD17 have demonstrated that slightly wider massive binaries are weighted considerably towards smaller mass ratios. Across intermediate periods log P (days) = 2–4, MD17 have found the massratio distribution is accurately described by a twocomponent powerlaw p_{q} ∝ q^{γ} with slopes γ_{smallq} across small mass ratios q = 0.1–0.3 and γ_{largeq} across large mass ratios q = 0.3–1.0. For massive binaries and log P (days) = 2, MD17 have fit γ_{smallq} = 0.0 and γ_{largeq} = –1.4, while for log P (days) = 4, they have measured γ_{smallq} = –0.7 and γ_{largeq} = –2.0 (see their Fig. 35). In Fig. 3, we compare the cumulative massratio distributions used by dMB15 and B16 to the updated distributions measured at log P (days) = 2–4.
Fig. 3. Cumulative distribution function of mass ratios. dMB15 and B16 adopted a uniform massratio distribution p_{q} ∝ q^{0} (green) based on a sample of massive spectroscopic binaries (Sana et al. 2012) that are dominated by shortperiod systems P < 20 days. MD17 (and references therein) found that massive binaries with intermediate periods (blue and red) are weighted significantly towards smaller mass ratios. 

Open with DEXTER 
Very massive binaries. At present, there are no reliable observational measurements of the properties of very massive binaries with intermediate periods log P (days) = 2–4 and primary masses M_{1} > 40 M_{⊙}. At these intermediate separations, the massratio and eccentricity distributions do not significantly change across primary masses 8 M_{⊙} < M_{1} < 40 M_{⊙}. In addition, the inner binary fraction and inner binary period distribution across log P (days) = 2–4 varies only slightly across the primary mass interval 8 M_{⊙} < M_{1} < 40 M_{⊙} (see Fig. 1). We therefore assume that systems with M_{1} > 40 M_{⊙} have the same period, eccentricity, massratio distributions, and binary fractions as those systems with M_{1} = 40 M_{⊙}.
Triplestar evolution. We model the evolution of all inner binaries with StarTrack. For triples, we ignore the outer tertiaries except for the following case. If the inner binary initially has P < 5 days and q < 0.4, it most likely will undergo unstable Case A Roche lobe overflow (RLOF), merge on the MS, and effectively evolve as a single star. If the tertiary has log P_{outer} (days) ≲ 4, then such a triple may lead to the formation of compact object merger. We model these triples with inner binaries P_{bin} < 5 days and q_{bin} < 0.4 as simple binaries, where the original inner binary merges to become the effective primary with combined mass M_{1} + M_{2} and the tertiary becomes the effective binary companion with mass M_{3}. For simplicity, we assume no rejuvenation of the merger product, i.e. evolve it as a star with the same age as its companion with mass M_{3}. About 10% of massive inner binaries with P < 5 days have outer tertiaries in tight, hierarchical configurations log P_{outer} (days) < 4 (MD17). This triplestar channel certainly does not dominate, but can contribute an additional ≈5–10% to the overall compact merger rate.
2.2. IMF and its dependence on metallicity
We adopt the same IMF that was used by B16, as guided by the recent observations (Bastian et al. 2010),
where for our standard model α_{3} = 2.3; we note that dMB15 adopt α_{3} = 2.7 but also investigate α_{3} = 2.2 and 3.2.
Additionally we calculate a model in which the IMF slope for massive stars depends on metallicity, as described below. In the following, we assume a conversion between [Fe/H] and Z given by Bertelli et al. (1994), appropriate for stars with solarlike fraction of iron among all the metal components.
Based on the observations of stellar clusters in the Milky Way and Nbody simulations (Sect. 1) Marks et al. (2012) calibrated that α_{3;Marks} = 2.63 + 0.66[Fe/H] for [Fe/H] < −0.5, which corresponds to Z ≲ 0.0063. We note that the above relation comes from a rather uncertain fit (see Fig. 4 of Marks et al. 2012) and is based on a model with several caveats (see discussion at the end of Sect. 1). Because the observations of OB associations and clusters in the Local Group do not support any significant deviations from α = 2.3 for metallicities Z ≥ 0.004 (Massey 2003), we modify the relation for α_{3} proposed by Marks et al. (2012) such that the IMF dependence towards lower metallicity does not occur until Z < 0.004, and α_{3} = 2.3 for Z ≥ 0.004. Marks et al. (2012) observational data extends down to about Z = 0.0001 ([Fe/H] ≈ –2.3, see their Fig. 4), and we decided to limit the further decrease of α_{3} for lower metallicities. Explicitly we adopt the following formula:
For some example values of Z, this formula yields α_{3}(Z = 0.1 Z_{⊙}) = 2.07 and α_{3}(Z = 0.01 Z_{⊙}) = 1.31.
We note that we treat the IMF dN/dM (all stars) and primary mass function dN/dM_{1} (single stars and primaries in binaries) interchangeably. Across large masses M > 1 M_{⊙}, both the IMF and primary mass function have the same slope α_{3} = 2.3 (Kroupa et al. 2013). The IMF and primary mass function only differ slightly across small masses M < 1 M_{⊙}. Since α_{3} (and its dependence on metallicity) is the most important parameter for determining merger rates, the assumption that dN/dM equals dN/dM_{1} is justified.
3. Computational method
3.1. Physical assumptions – the StarTrack code
We utilized the STARTRACK population synthesis code (Belczynski et al. 2002, 2008; Dominik et al. 2012). STARTRACK is developed based on analytical formulae for the evolution of a nonrotating star obtained by Hurley et al. (2000) from the fits to the grid of evolutionary tracks calculated by Pols et al. (1998). The original Hurley et al. (2000) models were updated with the prescriptions for the wind mass loss from OB type stars (Vink et al. 2001), WolfRayet stars (Hamann & Koesterke 1998; Vink & de Koter 2005), and enhanced mass loss rates for luminous blue variables (Belczynski et al. 2010a). For corecollapse supernovae we adopt the convection driven, neutrino enhanced “rapid” supernova engine from Fryer et al. (2012), which reproduces the observed mass gap between NSs and BHs (Belczynski et al. 2012). The key parameter of this model, dependent on the mass of the CO core at the time of explosion, is the fraction of material ultimately falling back onto the proto compact object (f_{b}) and contributing to the final remnant mass. The supernova kick velocity is drawn from a Maxwellian distribution with σ_{1D} = 265 km s^{−1} (Hobbs et al. 2005) and lowered proportionally to the amount of fallback, i.e. V_{NK} = V_{NK;Hobbs} (1 − f_{b}). This effectively means that the most massive BHs in our simulations (up to ∼15 M_{⊙} for Z_{⊙} and ∼40 M_{⊙} for 0.1 Z_{⊙}) are also typically those that receive a relatively small or zero natal kick (direct collapse).
We calculated the CE evolution in one step by applying the energy balance prescription of Webbink (1984). We adopt α_{CE} = 1 for the efficiency of energy transfer. The envelope binding parameter λ is taken from the fits of Xu & Li (2010) and depends on the radius and the evolutionary state of the donor. There are large uncertainties concerning the stability of mass transfer initiated by a Hertzsprung gap (HG) donor star. Pavlovskii et al. (2017) recently reported that for a large range of donor radii and masses such mass transfer may actually be stable, rather than leading to a CE phase. While revised criteria for a range of different metallicities and other parameters are still under preparation (Ivanova, priv. comm.), we opted to not follow the evolution of systems with HG donors in which our standard criteria for the mass transfer stability (Eq. (49) of Belczynski et al. 2008) indicate a dynamically unstable event.
Up until this point our physical assumptions are the same as in submodel B of dMB15 and standard model M1 of B16. The only recent update is the addition of mass loss due to pairinstability pulsations (model M10; Belczynski et al. 2016a). Such pulsations can affect stars with massive helium cores M_{He} between ∼40–45 M_{⊙} and ∼60–65 M_{⊙} and deplete them of a significant amount of mass (M_{ejecta} ∼5–20 M_{⊙}; Heger & Woosley 2002). We modelled this by assuming that stars with M_{He} = 45–65 M_{⊙} are subject to pairinstability pulsations and lose all their mass above 45 M_{⊙}. We note that the inclusion of pairinstability mass loss does not affect our predictions for detections of double compact mergers with the LIGO O2 observational run sensitivity (Belczynski et al. 2016a), and is consistent with existing data (Fishbach & Holz 2017).
3.2. Cosmological calculations of the merger rates
We placed our simulated systems into a cosmological background by populating the Universe up to z = 15. We modelled binaries across 32 different metallicities covering the range from 0.005 Z_{⊙} to 1.5 Z_{⊙}. In our standard model (I1) the IMF does not depend on metallicity, and so for each Z we use the same sample of 2 × 10^{6} systems (single stars, binaries, triples) generated according to the multiplicity statistics described in Sect. 2 and a primary mass function with α_{3} = 2.3 across 5 M_{⊙} < M_{1} < 150 M_{⊙}. We only evolve binaries (∼72% of systems with M_{1} > 5 M_{⊙}) and hierarchical triples (∼3.6%), in which the inner binary is likely to merge during an early case A mass transfer (Sect. 2.1). For submodel I2 we assume the metallicitydependent primary mass function according to Eq. (2), and for each Z we generate a corresponding sample of 10^{6} initial systems.
We calculate the merger rate density as a function of redshift by integrating through the cosmic star formation rate (SFR) history as described in see Sect. 4 of Dominik et al. (2015) and Sect. 2.2 of Belczynski et al. (2016c). For each redshift bin Δz = 0.1 up to z = 15 we use the cosmic SFR
from Madau & Dickinson (2014) as implemented in B16. The contribution from each of the metallicities in our simulations at a given redshift of star formation z_{SF} is then calculated as a fraction of SFR(z_{SF}) according to the mean metallicity evolution model from Madau & Dickinson (2014), increased by 0.5 dex to better fit observational data (Vangioni et al. 2015), i.e.
We adopt a return fraction R = 0.27 (fraction of stellar mass returned into the interstellar medium), a net metal yield of y = 0.019 (mass of metals ejected into the medium by stars per unit mass locked in stars), a baryon density with Ω_{b} = 0.045 and h_{0} = 0.7, a SFR from Eq. (3), and in a flat cosmology with Ω_{Λ} = 0.7, Ω_{M} = 0.3, Ω_{k} = 0, and H_{0} = 70.0 km s^{−1} Mpc^{−1}. We assume a log–normal distribution of metallicity with σ = 0.5 dex around the mean value at each redshift (Dvorkin et al. 2015).
Madau & Dickinson (2014) obtained their cosmic SFR based on the farUV (FUV; 1500 Å) and IR (8– 1000 μm) luminosity functions. Both are regarded as good tracers of star formation as they are dominated by the contribution from shortlived massive stars: FUV directly and IR through reradiation of dustabsorbed UV. Specifically, Madau & Dickinson (2014) converted strength of the 1500 Å line in the UV spectrum L_{ν}(FUV) into a cosmic SFR by applying an adequate conversion factor , such that (their Eq. (10)). To calculate they assume a Salpeter IMF with slope α = 2.35 across masses M = 0.1 – 100 M_{⊙} (Salpeter 1955). It is therefore inconsistent to implement the above Eq. (3) in combination with a more realistic Kroupalike IMF that turns over below M ≲ 0.5 M_{⊙}. We introduce a conversion factor for a given IMF such that
where SFR_{Salpeter}(z) is given by Eq. (3). In order to calculate for various IMFs we compute the corresponding UV spectra via STARBURST99 code, designed to model spectrophotometric properties of starforming galaxies (Leitherer et al. 1999, 2014; Vázquez & Leitherer 2005). In other words, when changing the IMF, we normalize the SFR from Madau & Dickinson (2014) in such a way that the amount of UV observed at 1500 Å stays the same. Details are given in Appendix A, together with a tabularized set of values at solar and subsolar metallicities for different Kroupalike IMFs.
In the model with metallicitydependent IMF at each redshift of star formation z_{SF} we use the mean metallicity of that redshift Z_{mean}(z_{SF}) to calculate a highmass IMF slope (Eq. (2)). We then use an appropriate value of from Table A.1 to correct the cosmic SFR (Eq. (5)). For example, the conversion to Kroupalike IMF (Kroupa et al. 1993) with α_{3} = 2.3 (Bastian et al. 2010), the IMF which was assumed in a the recent studies (Belczynski et al. 2016a,b), requires multiplying the cosmic SFR (Eq. (3)) by ^{3}.
The output of STARTRACK and the cosmological calculations described above is in the form of CBM happening at different redshifts. Each event is assigned with a normalization factor calculated as a convolution of the cosmic SFR and the metallicity distribution at the redshift of the binary formation (Eq. (7) of Belczynski et al. 2016c). Based on that we compute the source frame merger rate density (Gpc^{−3} yr^{−1}).
3.3. Calculations of predictions for the LIGO/Virgo detectors
Finally, we need to account for the detector sensitivity to produce predictions for the LIGO/Virgo observations. For each of our mergers we modelled the full inspiralmergerringdown waveform using the IMRPhenomD gravitational waveform templates (Khan et al. 2016; Husa et al. 2016). We adopt a fiducial O2 noise curve (“midhigh”) from Abbott et al. (2016d). Following Belczynski et al. (2016b) we consider a merger to be detectable if the signaltonoise ratio in a single detector is above a threshold value S/N = 8. We estimate detection rates as described in Belczynski et al. (2016c). This includes the calculation of detection probability p_{det}, which takes into account the detector antenna pattern (LIGO samples a “peanutshaped” volume rather than an entire spherical volume enclosed by the horizon redshift; Chen et al. 2017). For increased accuracy with respect to Belczynski et al. (2016c), where we used an analytic fit to obtain p_{det} (Eq. (A2) of Dominik et al. 2015), in this work we interpolate the numerical data for the cumulative distribution function available online^{4}. This improvement leads to a few per cent increase in detection rates.
3.4. Differences with respect to B16 and Belczynski et al. (2016a)
In terms of physical assumptions (Sect. 3.1), in the current study, we adopt exactly the model M10 from Belczynski et al. (2016a). Model M1 from B16 is also very similar, only differing by the lack of mass loss due to pairinstability pulsations.
There are, however, two differences between our current calculations and these two previous studies. We correct the assumed SFR (Eq. (3)) derived from a Salpeter IMF for a more realistic Kroupalike IMF (see Appendix A), which lowers the predicted merger rates by a factor of ∼0.51 in the case of α_{3} = 2.3. In both Belczynski et al. (2016a) and B16 there was a mistake in the way the simulated systems were normalized to match the entire stellar population (i.e. the calculation of the simulated mass M_{sim}). As a result, the calculated merger rates were about ∼1.82 smaller. The net effect of these two changes is a slight decrease of the merger rates (by a factor of 0.926). This is why the numbers we present for model M 10 in the following sections are 0.926 times smaller than the corresponding values presented in Belczynski et al. (2016a). We note that in the most recent work, Belczynski et al. (2017), the correction factor of 0.926 is already applied to all the models.
4. Results
4.1. Birth properties of compact binary mergers
Starting from the zeroage MS (ZAMS), we evolved ∼1.5 × 10^{6} binaries with primaries of at least 5 M_{⊙} drawn from the initial distributions of MD17. According to the assumed IMF and multiplicity statistics, our sample constitutes about 29% of the mass of all the stars forming at ZAMS, and so our simulations are normalized to a starforming mass of M_{sim} = 1.687 × 10^{8} M_{⊙}. The formation efficiencies of CBM, depend significantly on metallicity, but are generally between 10^{−5} and 10^{−7} (M_{⊙}^{−1}) (see Table 1); we define CMBs as the number of merging systems of a given type in our simulations per unit of starforming mass, i.e. .
We compare the birth properties of CBM progenitors between the two simulated samples: (1) one assuming the initial MS binary distributions of (Sana et al. 2012; old simulations from Belczynski et al. 2016a) and the other adopting the recent results of MD17. We discuss the NS–NS mergers together with BH–BHs because these two merger types originate from systems with similar mass ratios and their distributions of initial pericentre separations and primary masses are clearly distinguishable. The BH–NS progenitors are presented separately.
Fig. 4. Birth distributions of initial ZAMS binary parameters of the progenitors of BH–BH and NS–NS mergers. The histograms are normalized to unity for an easier comparison. NS–NS mergers dominate at Z = Z_{⊙} (upper panels) while BH–BH mergers dominate at Z = 0.1 Z_{⊙} and Z = 0.01 Z_{⊙} (lower panels). 

Open with DEXTER 
Comparison of the formation efficiencies per unit of starforming mass of CBM: BH–BH, BH–NS, and NS–NS between the two models with different initial distributions: either Sana et al. (2012, old) or Moe & Di Stefano (2017, new).
4.1.1. Progenitors of BH–BH and NS–NS mergers
In Fig. 4, we show the birth properties of systems that evolve to form BH–BH and NS–NS mergers. We only include double compact binaries that are expected to merge within the Hubble time. The histograms are normalized to unity for an easier comparison; see Table 1 for the relative formation efficiencies.
Primary masses. The distributions of primary masses are noticeably different in different metallicities. The peak around 10 M_{⊙} at Z = 0.02 corresponds to NS–NS progenitors, which are the dominant type of CBM at solar metallicity. At Z = 0.002 the formation of BH–BH begins to dominate and happens for primary masses M_{1} > 20 M_{⊙}, peaking around 25 and 45 M_{⊙}. This bimodal shape is a direct consequence of the bimodal distribution of CO core masses M_{CO} of supernova progenitors expected to undergo a direct BH formation. In the rapid supernova engine we adopted (Fryer et al. 2012), BHs receive no natal kick on formation if M_{CO} is either between 6 and 7 M_{⊙} or above 11 M_{⊙} (see their Eq. (16)). At Z = 0.01 Z_{⊙} the higher mass peak corresponding to M_{CO} > 11 M_{⊙} dominates.
The differences between the new and old initial distributions are very small. The NS–NS as well as lower mass BH–BH corresponding to M_{CO} = 6–7 M_{⊙} are most preferably formed from systems with very high initial mass ratios q > 0.8, which is why there is relatively less of such systems in the simulations with MD17 initial conditions.
Orbital periods, eccentricities, and periastron separations. The new initial binary period distribution of MD17 is in general very similar in shape to that proposed by Sana et al. (2012) in the crucial range of log P (days) between 2 and 4 (especially for systems with massive primaries ∼30 M_{⊙}; see Fig. 1). However, there is a noticeable shift towards larger periods (and separations) among compact binary merger progenitors in the new simulations. This can be explained in connection with the eccentricity distribution of these systems. In the case of MD17 the eccentricity distribution is slightly skewed towards e > 0.5, which is very different from that of Sana et al. (2012) showing a strong preference for nearly circular orbits. As a result the periastron separation distributions in the new and old simulations are very similar. This showcases the fact that the pericentre separation d_{per} determines the evolutionary path of a binary with given component masses and decides whether or not it will form a compact binary merger. In the evolutionary channel involving a CE phase, only a fixed range of d_{per}, that is independent of the choice of initial properties, results in the formation of CBM.
At Z = Z_{⊙}, where the formation of NS–NS dominates, the periastron separations are centred around ∼100R_{⊙}. At Z = 0.1 Z_{⊙} BH–BH become dominant and shift the d_{per} distribution towards higher values of over ∼100R_{⊙}. At this metallicity the distribution is clearly bimodal with peaks around ∼750R_{⊙} and ∼4000R_{⊙}. This comes from the fact that at Z ∼ 0.1 Z_{⊙} there are two relevant formation channels of merging BH–BH: the dominant channel (∼70% of the systems), in which there is only one CE phase taking place after the first BH has already formed (the peak at d_{per} ∼ 750R_{⊙}); and another channel (∼20% of the systems), in which there are two CE phases, one before each BH formation. The systems originating from the latter channel are those responsible for the second peak at d_{per} ∼ 4000R_{⊙}. These systems need to have larger initial d_{per} to prevent any RLOF during the HG evolutionary stage of the primary; we note that we do not allow for the CE phase from HG donors. Instead, the primary needs to initiate a RLOF and the CE phase later during the core helium burning phase. This channel becomes ineffective in our simulations for extremely low metallicities (Z ≲ 0.0005), where stars do not increase their radii sufficiently after the HG stage. For that reason at Z = 0.01 Z_{⊙}, the d_{per} distribution is no longer clearly bimodal. It is also slightly shifted towards smaller separations with respect to the distribution at Z = 0.1 Z_{⊙}, again because of smaller sizes of the stars. Overall, the BH–BH formation channel involving two different CE phases is only relevant (i.e. at least 5% of all the BH–BH formed this way) at metallicity range Z = 0.0005–0.004, contributing up to at most ∼22% of the BH–BH formation at metallicity Z ≈ 0.0025.
We note that the secondary peak in pericentre separations around ∼4000 R_{⊙} at Z = 0.1 Z_{⊙} is noticeably smaller in the simulations with the new initial distributions. This is because of the fact that, according to the multiplicity statistics of MD17, the number of inner binaries with massive primaries (≳ 30 M_{⊙}) having orbital periods of about log P (days) ≈ 3.5–4.0 is smaller with respect to the old distributions of Sana et al. (2012); see Fig. 1.
Mass ratios. Merging NS–NS and BH–BH in our simulations originate almost exclusively from binaries with high mass ratios: q_{NS–NS} ≳ 0.6 and q_{BH–BH} ≳ 0.5, respectively (at Z = 0.1 Z_{⊙} the systems with initial q < 0.5 are all BH–NS). The initial mass ratio distribution of Sana et al. (2012) is flat in q. Meanwhile, the new results of MD17 include a distribution decreasing quickly with q (see Fig. 3). This is reflected in the changes of the mass ratio distribution of CBM progenitors, which are slightly shifted towards q ∼ 0.5 at Z = 0.1 Z_{⊙} and Z = 0.01 Z_{⊙} in the new simulations.
4.1.2. Progenitors and formation of BH–NS mergers
In Fig. 5, we show the birth properties of systems that evolve to form BH–NS mergers. This time, for simplicity, we focus only on the three most important parameters: primary mass, mass ratio, and periastron separation.
The formation of BH–NS mergers is very different in the three different metallicites shown: Z_{⊙}, 0.1 Z_{⊙}, and 0.01 Z_{⊙}. In the solar metallicity there are hardly any BH–NS systems formed. The histograms for Z = Z_{⊙} are based on only around a dozen binaries, and, as a result, these histograms are not very meaningful. The reason why the BH–NS formation is so ineffective in our simulations at solar metallicity is that the stars with Z = Z_{⊙} already expand to their nearly maximum radius during the HG phase. This means that, in the majority of cases, the RLOF is initiated from the primary while it is a HG star. Because BH–NS systems typically originate from binaries with relatively small initial mass ratios of q ≲ 0.6, then, according to the standard STARTRACK criteria, such a mass transfer would most likely be flagged as unstable, leading to a CE phase from a HG donor. In this work, we decided to not follow the evolution of such systems (see Sect. 3.1 for more details).
At subsolar metallicity Z = 0.1 Z_{⊙} the formation of BH–NS mergers is much more effective, and the initial distributions are representative for the majority of BH–NS mergers formed across all of the 32 metallicites we compute. The primary expands much more after the HG phase than in the case of solar metallicity (Fig. 2 of Belczynski et al. 2010b) and, as a corehelium burning giant with a convective envelope, initiates an unstable RLOF and a CE phase. After the BH formation, the secondary eventually evolves off the MS, expands, and initiates a second CE episode. As a result, the binary is very close, the natal kick velocity gained by a newly formed NS does not disrupt the system and the binary merges within the Hubble time. We note that this is the exact same evolutionary route (i.e. two CE phases) that also operates in the case of BH–BH mergers in intermediate metallicites Z = 0.0005–0.004; see the paragraph on orbital periods in Sect. 4.1.1.
The distribution of primary masses at Z = 0.1 Z_{⊙} is slightly shifted towards lower values (peak at ∼ 40 M_{⊙}) with respect to the same distribution for BH–BH progenitors (peak at ∼50 M_{⊙}). The dominant range of mass ratios is between 0.3 and 0.6 (the contribution from q > 0.95 is described below). In the simulation with the updated initial distributions from MD17, the mass ratio distribution is slightly shifted towards smaller values. This, in turn, causes a small shift towards higher M_{1} values, so that the secondary mass distribution (not explicitly shown here) stays roughly the same. The optimal periastron separation is about 3000 R_{⊙} similar to the case of BH–BH mergers evolving through the same channel involving two CE phases; see the right peak in the periastron separation distribution for Z = 0.1 Z_{⊙} in Fig. 4. We note that in the case of new simulations with the initial conditions of MD17 there are fewer binaries with initial periods above log P (days) ≈ 3.5 with respect to the old distributions of Sana et al. (2012), which is why the peak in the periastron separations is slightly shifted towards smaller values (also similarly to the BH–BH case at Z = 0.1 Z_{⊙}).
In the case of very low metallicities (such as Z = 0.01 Z_{⊙}) the main formation channel of BH–NS mergers described above is no longer efficient. This is because the lower the metallicity the more massive the BHs formed. For comparison, the least massive BHs in our simulations are of about 5.5 M_{⊙} at Z = 0.1 Z_{⊙} and about 7.0 M_{⊙} at metallicity Z = 0.01 Z_{⊙}. As a result, at Z = 0.01 Z_{⊙}, the mass ratio of the components at the point when the secondary expands and initiates a RLOF is not sufficiently high to cause an unstable mass transfer and a CE phase.
Fig. 5. Birth distributions of initial ZAMS binary parameters of the progenitors of BH–NS mergers. The histograms are normalized to unity for an easier comparison. We note that the upper panels (Z = Z_{⊙}) are made based on a very small number of systems (<20) that form in solar metallicity. 

Open with DEXTER 
Fig. 6. Formation !efficiency (i.e. number of systems formed per unit of starforming mass, ) of different types of CBM at different metallicities. We compute a grid of 32 different metallicities (indicated with filled squares at the BH–BH curves), ranging from 0.005 Z_{⊙} to 1.5 Z_{⊙}. We indicate three different models: M10 in blue (initial distributions of Sana et al. 2012), I1 in green (initial distributions of MD17) and I2 in red (same as I1 but with a topheavy IMF at low metallicities; see Eq. (2)). We note that model I2 is only different from I1 at Z < 0.004 = 0.2 Z_{⊙}. 

Open with DEXTER 
At Z = 0.01 Z_{⊙} another formation channel dominates, which is also present, yet not so important at higher metallicities. This channel involves binaries with a mass ratio close to unity (q ≳ 0.9), comprised of stars with masses that are close to the threshold between the NS and BH formation (M ∼ 20 M_{⊙}). Because the binary components have similar masses, the secondary evolves transfer between the two corehelium burning stars, which reverses the mass ratio of the components. The primary eventually collapses to form a BH in a direct event with no mass eject and no natal kick. The “Rapid” BH formation model of Fryer et al. (2012) predicts a 100% mass fallback for the least massive BH progenitors with CO core masses between 6 and 7 M_{⊙}; see their Eq. (16). The rejuvenated secondary continues to expand and initiates a second RLOF. This time the mass ratio is lower and the mass transfer becomes unstable, leading to a CE phase. The binary becomes close enough to survive the supernova and NS formation, and later becomes a BH–NS merger. This scenario requires a specific range of initial binary parameters and is much less effective than the more standard evolutionary route dominating the BH–NS formation at intermediate metallicities.
4.2. Formation of compact binary mergers across metallicity
In Fig. 6 we show the relation between the formation efficiency X_{CBM} () and metellicity for different types of CBM. We note that similar trends were obtained recently by Giacobbo & Mapelli (2018). It is evident that the formation of BH–BH mergers becomes most effective at low metallicities Z ≲ 0.1 Z_{⊙}. This is a wellknown result (e.g. Belczynski et al. 2010b, 2016b; Dominik et al. 2012; Eldridge & Stanway 2016; Mapelli et al. 2017; Giacobbo & Mapelli 2018) for mainly two reasons: first, the fact that at lower metallicities the BH progenitors are relatively more massive and, therefore, more BHs are formed through direct collapse with zero natal kick, and second there is a higher chance for a CE event to be initiated by a corehelium burning donor, rather than by a HG donor, because of smaller radii of HG stars (Fig. 2 of Belczynski et al. 2010b).
According to stellar tracks from Hurley et al. (2000), at metallicities close to the solar Z ≃ 0.02, there is only a small range of separations at which a mass transfer from corehelium burning donors of > 20 M_{⊙} is initiated. Much more likely, the giant causes a RLOF earlier during its rapid HG evolution, in which case we choose to not allow for a CE evolution. Apart from strong stellar winds, this is what suppresses the formation of BH–BH and BH–NS compact binaries at metallicities close to the solar Z ≃ 0.02 in our simulations.
The formation efficiency of BH–NS mergers X_{BHNS}() does not increase in a similar fashion as X_{BHBH} () does at very low metallicities. There are two major factors at play. The main formation channel of BH–NS involving two CE phases becomes less effective in very low metallicities (Sect. 4.1.2). Additionally, natal kicks for NSs formed in core collapse are significant even at very low metallicities.
The formation of NS–NS mergers is much less metallicitydependent, as was similarly found by Giacobbo & Mapelli (2018). It is slightly more efficient at metallicities similar to solar (Z ≃ 0.02) than it is at lower metallicities Z ≲ 0.005. There is no one single reason for this behaviour, and it is most likely a combination of small effects such as different radii and different wind mass loss rates of stars with different metallicities. We thus consider this result to be much less robust than the increasing formation efficiency of BH–BH mergers towards lower metallicity.
The updated initial binary distributions of MD17 (green lines, I1) result in smaller formation efficiencies of CBM across all the metallicities with respect to the previously adopted distributions of Sana et al. (2012; blue lines, M10). The only exception is the case of BH–NS mergers at Z ≈ Z_{⊙}, where the formation efficiency X_{BHNS} () is at its lowest and there are few such systems in our simulations (fewer than 20). On the other hand, the inclusion of a topheavy IMF for Z < 0.004 = 0.2 Z_{⊙} (red lines, I2) significantly increases the formation efficiency of BH–BH and BH–NS mergers by up to nearly an order of magnitude for the lowest Z = 0.0001 = 0.005 Z_{⊙} in the case of BH–BH. The NS–NS mergers originate from binaries with less massive primaries, which is why there is little change between models I1 and I2 in this case.
Fig. 7. Source frame BH–BH merger rate density Gpc^{−3} yr^{−1} as a function of redshift for our three models: M10 (binary distributions of Sana et al. 2012), I1 (binary distributions of MD17), and I2 (same as I1 but with a metallicitydependent IMF). Upper panel: peak around z ≈ 2 in all the models corresponds to the maximum in the SFR (see text). Lower panel: range of redshifts accessible with the current and future sensitivity of LIGO detectors. We indicate the horizon redshifts for observation of an optimally inclined BH–BH merger with a total mass of ∼ 80 M_{⊙} (roughly the most massive systems in our simulations) for the sensitivities of O1 and O2 observing runs, and for the Advanced LIGO (AdLIGO). The detection distances for NS–NS mergers (Chen et al. 2017) were about d_{NSNS} = 70 Mpc for O1 and O2 sensitivity (Abbott et al. 2018). The local BH–BH merger rate density (within z < 0.02) in our models are about 203 (M10), 89 (I1), and 181 Gpc^{−3} yr^{−1} (I2). The current limits imposed by the existing LIGO detections are 12–219 Gpc^{−3} yr^{−1} (Abbott et al. 2017a). 

Open with DEXTER 
5. Merger rates and LIGO/Virgo predictions
5.1. Sourceframe merger rate density
In Fig. 7 we show the source frame BH–BH merger rate density (Gpc^{−3} yr^{−1}) as a function of redshift for our three models. We note that the results for M10 are also corrected with respect to Belczynski et al. (2016a) for a typo in the calculation of simulation normalization, yielding a small net decrease of merger rates by a factor of 0.926 (Sect. 3.4).
The behaviour of the merger rate density relation with redshift (upper plot) closely resembles the star formation history (Madau & Dickinson 2014) and peaks around z ≈ 2. This results from the fact that the distribution of merger delay times of close BH–BH binaries follows a powerlaw (Dominik et al. 2012; Belczynski et al. 2016b), and most of the systems merge relatively shortly after their formation (a few hundred Myr). We note that the maximum of BH–BH merger rate density in model I2 is very slightly (Δz ≈ 0.2) shifted towards smaller redshifts. This is because in this variation, at each redshift, we assume a different IMF based on the mean metallicity at that redshift, and consequently apply the adequate IMFSFR correction factor (see Appendix A for details). As a result, the location of the SFR maximum is slightly different.
The lower panel shows the range of redshifts up to z = 3.0. We also indicate the horizon redshifts for observation of an optimally inclined BH–BH merger with a total mass of ∼80 M_{⊙} (roughly the most massive systems in our simulations) for the sensitivities of O1 and O2 observing runs (z_{horizon} ≈ 0.7), and for the design sensitivity of Advanced LIGO (z_{horizon} ≈ 2.0; Abbott et al. 2018). For comparison, the detection distances for NS–NS mergers (d_{NSNS}; i.e. the radius of a sphere with an equal volume to the peanutshaped LIGO response function; see Chen et al. 2017) were about d_{NSNS} = 70 Mpc for O1 and O2 sensitivity. The local BH–BH merger rate density (within z < 0.02) in our models are about 203 (M10), 89 (I1), and 181 Gpc^{−3} yr^{−1} (I2), which are all within the most uptodate range determined by the existing LIGO detections (12–219 Gpc^{−3} yr^{−1}; Abbott et al. 2017a).
5.2. LIGO/Virgo detection rates
In Table 2, we summarize our predictions for the local (i.e. within z < 0.02) merger rate densities, detection rates (R_{det} yr^{−1}) with the sensitivity of O1/O2 LIGO/Virgo observing runs, and corresponding numbers of detections assuming 120 days of coincident data acquisition. The inclusion of the updated initial MS binary distributions of MD17 (model I1) has led to a decrease in both the local BH–BH merger rate density and the detection rate during O2 by a factor of ∼2.3 (with respect to model M10). The assumption of a metallicitydependent IMF (model I2, see Eq. (2)), on the other hand, causes an increase in the local BH–BH merger rate density by a factor of ∼2, and an increase in the O2 detection rate by a larger factor of ∼2.5 (with respect to model I1). The difference between these two factors arises from the fact that the topheavy IMF for low metallicities assumed in I2 favours more massive BHs, which are easier to detect. We discuss the comparison of our results with LIGO/Virgo detection rates in Sect. 6.4.
5.3. Metallicity distribution of BH–BH progenitors
In Fig. 8, we show the metallicity distributions of different types of CBM weighted by their detection rates with O1/O2 sensitivity. In the case of BH–BH mergers, the huge preference for very low metallicities (Z ≲ 0.001 = 0.05 Z_{⊙}) with respect to metallicities close to solar (Z ≳ 0.01 = 0.5 Z_{⊙}) is already present when looking at the formation efficiencies across Z alone (see Fig. 6). At that stage, which does not account for any cosmological calculations yet (i.e. SFR, cosmic metallicity distribution, BH–BH merger delay times), there are about 2 orders of magnitude more BH–BH mergers formed with low Z ≲ 0.001 compared to Z ∼ Z_{⊙} (solely due to binary evolution effects). Notably, this difference grows much higher in the detector frame, with the values of dR_{det}/dZ being about 4 order of magnitude larger for Z ≲ 0.001 than for Z ∼ Z_{⊙}. This is a combination of two effects. First, at low metallicities BHs are born more massive and their mergers are easier to detect with the LIGO/Virgo sensitivity. Second, in our approach to population synthesis on cosmological scales (Sect. 3.2) the metallicity distribution of stars forming at redshifts around the maximum of the cosmic SFR at z ≈ 2, from which most of the BH–BH merger progenitors originates (see Fig. 2 of B16), is centred at Z < 0.1 Z_{⊙} (see Fig. 6 of Extended Data of B16).
Fig. 8. Detectorframe distribution of metallicities of BH–BH, BH–NS, and NS–NS merger progenitors. More than 90% of BH–BH merger detections in all our models are systems formed with metallicites lower than 10% of the solar metallicity, i.e. Z < 0.1 Z_{⊙}. Such a strong dependence of compact binary merger formation with metallicity (see also Fig. 6), especially in the case of BH–BH systems, indicates the significance of SFRs at different metallicites across the Universe for the compact binary merger rates (see discussion in Sect. 6.3). 

Open with DEXTER 
Merger rate densities and detection rates for the LIGO/Virgo O2 run sensitivity for three models: M10 (initial binary distributions from Sana et al. 2012), I1 (updated distributions from MD17), and I2 (I1 plus topheavy IMF for low metallicities from Marks et al. 2012).
6. Discussion
6.1. Understanding the impact of the new initial distributions
Within the classical binary evolution scenario, our predictions for the merger rates of double compact objects based on the updated initial MS binary parameters from MD17 are ∼2–2.5 times smaller than the results previously obtained by dMB15 and B16 from the distributions of Sana et al. (2012). This can be explained by looking at the differences between the two distributions within the areas of the parameter space that are most relevant for the formation of double compact object mergers. Let us first take a look at BH–BH and NS–NS progenitors.
(i) In agreement with dMB15, we find that the majority of BH–BH and NS–NS mergers originate from binaries with initial periods within log P (days) = 2–4. The fraction of primaries residing in such systems is about F_{logP = 2–4;bin} = 28% in the case of Sana et al. (2012) distributions. Meanwhile, this value is higher by about ∼1.35 in the case of statistics obtained by MD17 (Fig. 1): i.e. F_{logP = 2–4;bin} = 37% for M_{1} = 8 M_{⊙} primaries, and F_{logP = 2–4;bin} = 40% for M_{1} = 30 M_{⊙} primaries (regimes for NS–NS and BH–BH progenitors, respectively).
(ii) The eccentricity distribution based on the Sana et al. (2012) sample is skewed towards lower values with 〈e〉 = 0.3 compared to the statistics reported by MD17 with 〈e〉 = 0.6 in the case of massive binaries (Fig. 2). The progenitors of CBM derive from binaries with periastron separations log d_{per} (R_{⊙}) = 2–3.5, depending slightly on metallicity and the exact formation channel. By increasing the average initial eccentricity from 〈e〉 = 0.3 to 0.6, then the average orbital separation of the progenitors increases by a factor of (1 –0.3)/(1 –0.6) = 1.75. This corresponds to a shift in logarithmic orbital period of Δlog P ≈ 0.2. Considering the progenitors derive from a broad parameter space of orbital periods log P (days) = 2–4, a slight shift of ΔlogP ≈ 0.2 due to the updated eccentricity distribution does not affect the predicted rates of compact object mergers.
(iii) Similar to dMB15, we find that nearly all BH–BH and NS–NS mergers derive from initial MS binaries with mass ratios q = 0.5–1.0. According to the old massratio distribution from Sana et al. (2012), ≈55% of earlytype binaries have mass ratios across this interval. According to the updated perioddependent massratio distribution of MD17, however, only (17–28)% of massive binaries with log P (days) = 2–4 have mass ratios q = 0.5–1.0, which results in about ∼2.0–3.2 times fewer potential BH–BH and NS–NS merger progenitors.
Altogether, the new mass ratio distribution that is shifted towards smaller values with respect to the old distribution (iii) is the most significant update in the context of BH–BH and NS–NS mergers, and the differences between period distributions are of secondary importance (i). Based on the estimates presented above, we could expect that there should be roughly between 1.5 and 2.4 fewer BH–BH and NS–NS mergers in the updated simulations. This explains the results of our computations (see Tables 1 and 2).
In the case of BH–NS mergers, the binary parameters of their progenitors depend strongly on metallicity (Fig. 5). Around solar metallicity, there are few to no BH–NS systems formed and our sample is not statistically significant. At subsolar metallicities of ∼0.1 Z_{⊙}, the main formation channel involves two CE phases, and this channel dominates when summed across all metallicities. It requires the initial mass ratios between 0.3 and 0.6 (peaking around 0.4) and, more importantly, pericentre separations of log d_{per} (R_{⊙}) = 3.2–4, preferably over 3.5. The updated initial distributions predict fewer massive binaries this wide. Additionally, in most cases the NS is formed from a star with mass at the higher end of the mass range of NS progenitors (15 ≲ M_{2}/M_{⊙} ≲ 22), which helps to produce a high enough mass ratio at a BH–MS binary stage, so that a dynamically unstable mass transfer occurs and a CE phase is triggered. The shift towards smaller mass ratios in the updated distributions causes a corresponding shift towards slightly more massive primaries. This additionally decreases the BH–NS formation efficiency due to a declining slope of the IMF.
At very low Z ∼ 0.01 Z_{⊙}, the dominant formation channel involves binaries with components of initially similar masses that evolve into doublegiant systems and experience a stable mass transfer instead of the first CE (Sect. 4.1.2). The requirement of q > 0.9 means that there are about three times fewer potential progenitors for this channel in the updated distributions. With all the metallicities combined, all the factors described above result in the merger rate of BH–NS systems being about two times smaller in the new simulations (Table 2).
6.2. Small impact of a topheavy IMF
The uncertainty in! the IMF is often considered to have a substantial impact on the compact binary merger rates, by a factor of a few dMB15 up to an order of magnitude (see Fig. 20 of Kruckow et al. 2018). One could thus expect that the topheavy IMF for low metallicities (Z < 0.2 Z_{⊙}; see Eq. (2)) we implemented in model I2 would result in a significant increase of BH–BH merger rates. According to Fig. 8, the vast majority of BH–BH mergers detectable with LIGO/Virgo were formed in Z < 0.1 Z_{⊙} metallicity environments. In model I2, the IMF powerlaw exponent for massive stars is α_{3} ≈ 2.07 for Z = 0.1 Z_{⊙} decreasing down to as low as α_{3} ≈ 1.31 for Z = 0.01 Z_{⊙}. For a fixed amount of stellar mass formed (i.e. fixed SFR) and with respect to the standard value α_{3} = 2.3 as in models M10 and I1, this would correspond to an increase in the number of primaries within M_{1} = 45–55 M_{⊙} (highest bin in the primary mass distribution of BH–BH progenitors, Fig. 4) by a factor of about ∼2.2 at Z = 0.1 Z_{⊙} and ∼5 at Z = 0.01 Z_{⊙}. Meanwhile, the merger rate density of BH–BH in model I2 (with a topheavy IMF) is higher than in model I1 by only a factor of ∼2.
This showcases the significance of keeping the consistency between an adopted cosmic SFR and an assumed IMF. Because the cosmic SFR is measured based on UV observations (Madau & Dickinson 2014), we normalize our simulations in such a way that the total amount of UV light at 1500 Å stays roughly constant, irrespective of the assumed IMF (see Appendix A). As most of the UV light is produced by massive stars, a topheavy IMF is characterized by a higher UV luminosity per unit of stellar mass formed and must result in a rescaling of the SFR to a smaller value. A smaller SFR means fewer stars formed, which counteracts the fact that with a topheavy IMF a higher fraction of these stars are massive. A similar but inverted argument could be made for a toplight IMF (i.e. α_{3}>2.3). Thus, the total number of massive stars formed stays roughly constant irrespective of the IMF variations, as it is constrained by the UV observations.
We show this quantitatively in Fig. 9, where, as a function of α_{3}, we plot the relative number ℛ(α_{3}) of stars formed in a given mass bin M_{low} – M_{high} (e.g. between 8 and 13 M_{⊙}, and so on) defined as
where ξ(α_{3}, M) is a Kroupalike IMF in the form given by Eq. (1) for which we vary the highmass slope α_{3}, and M_{total;1} and M_{total;2} are the total amounts of stellar mass formed in models with different IMFs. The IMFs are normalized to unity, i.e. for any slope α_{3}. The upper panel represents simulations normalized to a fixed amount of stellar mass formed, irrespective of the IMF (the usually made assumption), so M_{total;1} = M_{total;2}. The lower panel, in turn, assumes a fixed UV luminosity of a starforming population, so that . The amount of UV light per unit of starforming mass is inversely proportional to coefficients in Eq. (A.1) and Table A.1.
We note that although the relative number of primaries does not change by more than a factor of ∼2 in any of the mass ranges and IMFs shown, the ratio of progenitors of massive BH–BH to NS–NS mergers changes more significantly^{5}.
6.3. Population synthesis on cosmological scales
The discovery of GW150914 revealed the existence of stellar BHs with masses of 30–40 M_{⊙} (Abbott et al. 2016c). Black holes this massive are only expected to form from lowmetallicity stars Z/Z_{ ⊙ } ≲ 0.3 (Belczynski et al. 2010a), and this conclusion is strengthened by a recent update to WolfRayet star winds (Vink 2017). Additionally, Fig. 6 shows that the formation efficiency of BH–BH mergers is very strongly dependent on metallicity. For that reason, predictions for the rates of such events are sensitive to the applied distribution of metallicity of star formation throughout the Universe and they cannot be made by simply extrapolating results for a Milky Waylike galaxy, as still done by some authors.
We have adopted the mean metallicityredshift relation from the chemical evolution model of Madau & Dickinson (2014, see Eq. (4)) with the level increased by 0.5 dex and with a Gaussian spread of σ = 0.5 dex, and used this relation to weight the contributions from different metallicities to the SFR at different redshifts. This cosmological information, combined with the LIGO/Virgo sensitivity, allowed us to transform our simulated CBM into detectorframe predictions.
Fig. 9. Relative numbers of massive stars ℛ formed within multiple initial mass ranges as a function of powerlaw exponent α_{3} for the highmass IMF component (M > 1 M_{⊙}). See Eq. (6) for the exact definition of ℛ. Upper panel: same SFR, irrespective of the IMF changes. The number of stars is thus calculated with respect to a fixed stellar mass formed. Lower panel: SFR corrected for the IMF changes (see Eq. (5)). The number of stars is thus calculated with respect to a fixed amount of UV radiation at 1500 Å wavelength that is constrained by the observations (Madau & Dickinson 2014). 

Open with DEXTER 
However, the mean metallicity of the Universe (defined as the total mass of heavy elements ever produced per mass of baryons in the Universe) in general does not equal the mean metallicity of star formation at a given redshift. This simplification is a caveat of our model to keep in mind. The mean metallicity at which star formation takes place is likely higher, since the highest SFRs are revealed by massive galaxies (LaraLópez et al. 2013), which are also relatively metal rich (e.g. Tremonti et al. 2004). Lowmetallicity star formation in the local Universe, on the other hand, occurs mostly in lowmass dwarf galaxies (Andrews & Martini 2013), which have relatively little star formation (Boogaard et al. 2018).For that reason, we are likely somewhat overestimating the current SFR of metalpoor massive stars. For example, according to our model, about 18% of massive stars that recently formed in the local Universe (z = 0) have Z < 0.1 Z_{⊙}. Observations of nearby dwarf and spiral galaxies suggest only a few percent, perhaps at most 10%, of massive stars form locally at such low metallicity (see Appendix B for a detailed calculation). While this signifies the existence of a problem, the issue is more complicated and requires a coherent and observationally based model of cosmic star formation at different metallicities and across different redshifts (Chruslinska et al., in prep.).
A different approach to population synthesis on a cosmological scale is to distribute simulated systems over individual galaxies. This could be obtained in a semianalytical manner by applying a galaxy mass function (Elbert et al. 2018) or fits to galaxy trees from cosmological simulations (Lamberts et al. 2016), and combining these with observationally inferred galaxy scaling relations for SFRs and metallicity distributions (e.g. Behroozi et al. 2013). Alternatively, results from cosmological simulations could be applied directly to obtain similar information about star formation in individual galaxies (Schneider et al. 2017; Mapelli et al. 2017). Interestingly, Mapelli et al. (2017) predict that the distribution of metallicity of BH–BH merger progenitors peaks at around ∼ 0.1 Z_{⊙} and declines towards lower Z, which conflicts with our Fig. 8. The most likely explanation of this discrepancy is that we are overestimating the contribution from very low metallicities at close redshifts of star formation. The majority of the lowmetallicity content of the recent Universe is locked in dwarf galaxies with little star formation, for which our simplification that the mean metallicity of the Universe ≈ mean metallicity of star formation is likely less accurate. This interpretation also agrees with the fact that Mapelli et al. (2017) did not find the clearly bimodal distribution of massive BH–BH merger birth redshifts as shown in the case of our simulations in Fig. 2 of B16, but rather only the single maximum around the peak of cosmic star formation at z ≈ 2. However, even if we were overpredicting the star formation in low metallicities Z < 0.1 Z_{⊙} at redshifts z < 2 to the point at which we would have to discard all our mergers formed in these conditions, this would still not affect our predictions for the merger rates by more than a factor of 2.
Although making predictions for CBM based on cosmological simulations as carried out by Schneider et al. (2017) or Mapelli et al. (2017; or, in a simplified way, also by Belczynski et al. 2018) may be burdened with additional biases, this ultimately seems to be the superior approach that hopefully will allow us to confront models with various galactic properties in the case of mergers with identified hosts (such as GW170817, Abbott et al. 2017d).
Finally, we wish to highlight that the small impact of IMF variations on merger rates (Sect. 6.2) is a general result that should also appear in the predictions based on cosmological simulations. Even though such simulations often provide SFR values as their output, changing the assumptions on the massivestar IMF without renormalizing the SFR make their results inconsistent with star formation tracers such as UV observations.
6.4. Comparison with LIGO/Virgo detection rates and other studies
We show the merger (Gpc^{−3} yr^{−1}) and detection rates (yr^{−1}) from our simulations in Table A.1. While our models perform relatively well in retrieving the BH–BH merger rate as constrained by the LIGO/Virgo observations R_{BHBH} = 12–213 Gpc^{−3} yr^{−1} (Abbott et al. 2017a), they underpredict the rate of NS–NS mergers that was inferred after the recent detection of GW170817: (Abbott et al. 2017d)^{6}. Recently, Chruslinska et al. (2018) analysed NS–NS merger rates within a large suite of STARTRACK models, varying multiple binary evolution assumptions, and found that it is, in fact, very difficult to obtain R_{NSNS} consistent with the LIGO/Virgo constraints. The only three models marginally consistent with the reported NS–NS merger rates, which have R_{NSNS} values of about 380, 450, and 630 Gpc^{−3} yr^{−1}, overpredict the merger rate of double BHs, resulting in R_{BHBH} of about 1070, 310, and 700 Gpc^{−3} yr^{−1}, respectively. This discrepancy could perhaps be a hint that the stability criteria for mass transfer in Rochelobe overflowing highmass Xray binaries were different in the cases of BH and NS accretors, resulting in a stable mass transfer (i.e. avoiding CE) for a large space of binary parameters in the case of BH accretors (see also Pavlovskii et al. 2017).
Similarly, using code C_{OM}B_{IN}E, Kruckow et al. (2018) obtained NS–NS merger rates of around 10 Gpc^{−3} yr^{−1} in their default model; this is well below the LIGO/Virgo limits. Their most optimistic model yields a value of 159 Gpc^{−3} yr^{−1}, which could be increased further up to 400 Gpc^{−3} yr^{−1} to be marginally consistent with the observational constraints if natal kick velocities were reduced by half. Contrary to Chruslinska et al. (2018), the optimistic model of Kruckow et al. (2018) does not overpredict the BH–BH merger rates. However, it was only calculated for a single metallicity value Z = 0.0088 = 0.44 Z_{⊙} above the lowmetallicity regime at which the BH–BH formation is the highest.
Interestingly, using code MOBSA combined with the results of Illustris simulations (Nelson et al. 2015), Mapelli & Giacobbo (2018), were able to obtain both NS–NS and BH–BH simultaneously consistent with the LIGO/Virgo constraints. It should be noted that they did have to assume both very low natal kicks (σ ≃ 15 km s^{−1}) for all of their corecollapse supernovae and a rather high efficiency of the CE ejection (α_{CE} = 5), which is a value that usually calls for some additional process aiding the CE ejection. The former assumption can to some extent be justified by the fact that NS progenitors in binary systems tend to get their envelopes stripped in mass transfer episodes, which leads to lessenergetic supernovae and smaller mass ejections, and possibly natal kicks as small as ∼ 15 km s^{−1}, although likely only in the case of ultrastripped stars that also lose their helium envelope to a compact object accretor (Tauris et al. 2015). As not all NS progenitors in NS–NS systems are expected to be ultrastripped (especially those that form the first NS), a single σ ≃ 15 km s^{−1} value for all the core collapse supernovae may be an oversimplification (for comparison, see the model of Kruckow et al. 2018). Observationally, some Galactic closeorbit NS–NS systems show evidence of rather large natal kicks of over >100 km s^{−1}, while others require small kick velocities <50 km s^{−1} (see Sect. 6.4 of Tauris et al. 2017).
7. Summary
Observational studies have long hinted at probable correlations between the initial binary parameters primary mass (M_{1}), mass ratio (q), orbital period (P), and eccentricity (e) (Abt et al. 1990; Duchêne & Kraus 2013), but hitherto the selection biases have been too large to accurately quantify the intrinsic relations. A recently published paper, Moe & Di Stefano (2017), analysed results from more than 20 massive star surveys and, for the first time, obtained a joint probability density function f(M_{1}, q, P, e) for the initial ZAMS binary parameters (Sect. 2). We implement this result and analyse its impact on the predictions of merger rates of double compact object detectable with the LIGO/Virgo interferometers that originate from the isolated evolution scenario involving a CE phase. Using the STARTRACK rapid binary evolution code (Sect. 3.1; Belczynski et al. 2002, 2008), we evolve a large population of massive MS binaries across a range of 32 metallicities from 0.005 Z_{⊙} to 1.5 Z_{⊙} with their initial parameters drawn from the interrelated distribution of Moe & Di Stefano (2017). We distribute our systems at cosmological scales according to an observationally inferred cosmic SFR and a mean metallicityredshift relation obtained within a chemical evolution model by Madau & Dickinson (2014, see Sect. 3.2). We compare our results to those obtained by de Mink & Belczynski (2015) and Belczynski et al. (2016b) who performed similar simulations but with the noncorrelated initial binary distributions from Sana et al. (2012).
Additionally, we discuss the level of uncertainty of compact binary merger rate predictions associated with possible variations of the massive star IMF. We make an indepth study of one such variation, according to which the IMF becomes more top heavy (i.e. flatter powerlaw exponent for massive stars) the lower the metallicity of star formation (Marks et al. 2012, see Sects. 1 and 2.2). We make sure that along with the changes of the IMF, the amount of farUV light at 1500 Å wavelength (based on which the cosmic SFR is measured Madau & Dickinson 2014) stays the same.
We note that our merger rates are in good agreement with the LIGO/Virgo constraints in the case of BH–BH systems, but they are strikingly too small when it comes to NS–NS mergers (Sect. 5.2, see also Chruslinska et al. 2018). However, we predict that our results of the relative impact of initial distributions and IMF variations will hold in the general case of isolated evolution channels involving a CE phase, regardless of the absolute numbers of mergers obtained.
We arrive at the following conclusions:
(a) The introduction of the updated initial ZAMS binary distributions from MD17 (model I1) decreases the formation efficiency and, consequently, merger and detection rates of all types of CBM by a factor of about 2.1–2.6 with respect to the old simulations adopting Sana et al. (2012) distributions (model M10); see Sect. 5. This is a very small change in comparison to uncertainties associated with the binary evolution (e.g. Eldridge & Stanway 2016; Mapelli et al. 2017; Stevenson et al. 2017; Chruslinska et al. 2018; Kruckow et al. 2018).
In the case of BH–BH and NS–NS mergers, the major factor in play is the difference in the initial mass ratio distributions. Sana et al. (2012) obtained their binary statistics based on the spectroscopic measurement of massive Otype stars. For that reason their sample is limited to log P (days) < 3.5 (spectroscopically detectable binaries) and dominated by very shortperiod orbits with P < 20 days, to which they fit a flat mass ratio distribution. However, MD17 have shown that wider binaries are weighted towards considerably smaller mass ratios (Fig. 3). According to their updated distributions, across intermediate periods log P (days) = 2–4, there are ∼2.5 fewer systems with q > 0.5, which is the range of orbital periods and mass ratios at which BH–BH and NS–NS mergers are formed (see Fig. 4).
Of secondary importance is the difference between the orbital period distributions (up to 40% more system in the orbital period range log P (days) = 2–4, depending on metallicity), whereas the changes in the eccentricity statistics have negligible influence on merger rates even though an increase from 〈e〉 = 0.3 to 〈e〉 = 0.6 may seem significant at first.
(b) The introduction of a topheavy IMF at low metallicities (Sect. 2.2; Marks et al. 2012) results in a small increase of BH–BH merger rates by only a factor of ∼2, and even less significant differences in the case of BH–NS and NS–NS mergers (Table 2). This might seem surprising because for the values of powerlaw exponent α_{3} for the IMF of massive stars of α_{3} ≈ 1.8 at Z = 0.05 Z_{⊙} and even α_{3} ≈ 1.3 at Z = 0.01 Z_{⊙}, we could expect an increase in the number of massive BH–BH progenitors (for a fixed stellar mass formed) by a factor of ∼4–5 with respect to the standard value α_{3} = 2.3. Interestingly, recent spectroscopic observations of 247 massive stars in the 30 Doradus star formation region reveal an IMF slope of in the range 15–200 M_{⊙} (Schneider et al. 2018). However, because the formation of CBM likely occurs at cosmological scales (e.g. Dominik et al. 2013; Mapelli et al. 2017), we argue that the simulations should be normalized to a fixed amount of UV light at wavelength 1500 Å (used as an SFR indicator) produced by a starforming population. A more topheavy IMF produces more UV light per unit of stellar mass, which forces a rescale of the cosmic SFR to a lower value and prevents the BH–BH merger rates from increasing by more than a factor of ∼2. This is true even for very topheavy IMFs with α_{3} value as low as 1.0 (Fig. 9). We conclude that the previously reported uncertainty by a factor of 6 up and down of the merger rate predictions due to possible variations in the IMF (de Mink & Belczynski 2015) is, in fact, not more significant than a factor of ∼2.
Acknowledgments
We thank the anonymous referee for very useful remarks. JK gratefully acknowledges comments and pointers from JJ Eldridge, Michela Mapelli, and Gijs Nelemans. KB acknowledges support from the Polish National Science Center (NCN) grants: Sonata Bis 2 DEC2012/07/E/ST9/01360, Meastro 2015/18/A/ST9/00746, and Opus (2015/19/B/ST9/01099). MM acknowledges financial support from NASA’s Einstein Postdoctoral Fellowship programme PF5160139. DEH was partially supported by NSF grant PHY1708081. They were also supported by the Kavli Institute for Cosmological Physics at the University of Chicago through an endowment from the Kavli Foundation.
Appendix A: Conversion coefficients between SFR and different IMFs
As we describe in Sect. 3.2, any cosmic SFR is associated with an assumed IMF and it is inconsistent to use it in combination with a different IMF. We therefore introduce a conversion factor to correct the cosmic SFR of Madau & Dickinson (2014 Eq. (3)) obtained for the Salpeter IMF (Salpeter 1955) for the purpose of combining it with a more realistic Kroupalike IMF, i.e.
where SFR_{Salpeter}(z) is given by Eq. (3). We calculate as a relative strength of the 1500 Å line in the UV spectrum of a starforming region with SFR = 1.0 M_{⊙} yr^{−1} and the Salpeter IMF with respect to the same SFR but with another given IMF.
In order to compute the UV spectra for different IMFs, we utilize the STARBURST99 code, designed to model spectrophotometric properties of starforming galaxies (Leitherer et al. 1999, 2014; Vázquez & Leitherer 2005). We configurate STARBURST99 to use the newest evolutionary tracks published by the Geneva group: for solar (Z = 0.014; Ekström et al. 2012) and subsolar metallicity (Z = 0.002; Georgy et al. 2013), which cover the mass range M/M_{⊙} ∈ [0.8, 120]. In this paper our population synthesis simulations covered a much larger grid of 32 different metallicities from Z = 0.0001 up to Z = 0.03. However, the differences in UVtoSFR conversion factors for different metallicities are not very significant (see Fig. 3 of Madau & Dickinson 2014). Thus, we apply calculated for Z = 0.014 Geneva stellar tracks in the case of the STARTRACK models with Z ≥ 0.006, and apply calculated for Z = 0.002 to all the other STARTRACK models for Z < 0.006. The upper mass limit 120 M_{⊙} is a limitation of all the stellar tracks implemented in STARBURST99. We calculate the UV spectra up to 10^{8} yr after the beginning of a constant star formation. By that time the 1500 Å line is expected to have reached its asymptotic strength already (see Fig. 2 of Madau & Dickinson 2014).
In Table A.1 we list the values of computed at solar and subsolar metallicity for multiple powerlaw exponents of the highmass end of a Kroupalike IMF. Even if the highmass end of the IMF is significantly top heavy (e.g. α_{3} = 1.0), we find there are only ∼30% more stars with M_{1} = 45–55 M_{⊙} to match the same 1500 Å emission strength.
We wish to mention that the STARBURST99 code does not incorporate the evolution of binaries. Interacting binaries, for instance, are expected to affect the emission properties of a starforming population such as its ionizing flux or a UV luminosity (e.g. Stanway et al. 2016). Depending on metallicity and for a Kroupalike IMF truncated at 100 M_{⊙}, Eldridge et al. (2017) have calculated that binary interactions increase the amount of radiation ν L_{ν} at 1500 Å from a starforming population by a factor of ∼1.2–1.25 relative to a corresponding population of single stars (see their Table 4). We expect that the values of our conversion coefficients in Table A.1 could differ by a similar factor if the effect of binaries on spectral properties was taken into account.
Conversion coefficients for a cosmic SFR between the Salpeter and different Kroupalike IMFs.
Appendix B: Local SFR at metallicities Z < 0.1 Z_{⊙}
We use observational results to do a rough estimate of the local (i.e. z ∼ 0) SFR density at metallicities below 0.1 Z_{⊙}. At close redshifts, star formation at Z < 0.1 Z_{⊙} occurs in galaxies with mass < 10^{7} M_{⊙} (Andrews & Martini 2013). The local number density of starforming galaxies with masses between 10^{6} and 10^{7} M_{⊙} can be estimated as 10^{−0.8} Mpc^{−3} (from extrapolating the blue fit in Fig. 15 of Baldry et al. 2012), while their average SFR is measured to be around 10^{−2.5} M_{⊙} yr^{−1} Boogaard et al. 2018. This implies a local SFR density at Z < 0.1 Z_{⊙} of 10^{−3.3} M_{⊙} Mpc^{−3} yr^{−1}, which constitutes about ∼ 3% of the total SFR density at redshift z = 0 (from Eq. (3)). While the above estimate is highly uncertain, it seems unlikely that this fraction could be higher than 10%.
Appendix C: Completeness of simulations
Number of mergers of respective types, BH–BH, BH–NS, or NS–NS, in our samples simulated for models I1 (MD17 initial binary conditions) and I2 (additionally a metallicitydependent IMF; values given in parenthesis) for each of the 32 metallicities we computed.
Throughout this study, we adopt Z_{⊙} = 0.02 (Villante et al. 2014).
Such a correction is missing in all our previous studies, except for Chruslinska et al. (2018). However, its effect w ould be of secondary importance compared to uncertainties arising from evolutionary models (e.g. natal kicks), so thei r conclusions remain unchanged.
The simplified way we correct the SFR explains why the BH–BH merger rate increases by a factor 2 between model s I1 and I2 while none of the curves in lower panel of Fig. 9 actually reaches the value of 2. In practice, at each re dshift bin we assume that the amount of UV light per unit of starforming mass is such as for the median metallicity ( and IMF corresponding to it) at that redshift. However, given the uncertainty of the contributions from different meta llicities to the cosmic SFR, this simplification has no meaningful impact on the results.
References
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016a, Phys. Rev. Lett., 116, 241103 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016b, Phys. Rev. Lett., 116, 061102 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016c, Phys. Rev. Lett., 116, 241102 [NASA ADS] [CrossRef] [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016d, Liv. Rev. Rel., 19, 1 [NASA ADS] [CrossRef] [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, 851, L35 [NASA ADS] [CrossRef] [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017c, Phys. Rev. Lett., 119, 141101 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017d, Phys. Rev. Lett., 119, 161101 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2018, Liv. Rev. Rel., 21, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Abt, H. A. 1983, ARA&A, 21, 343 [NASA ADS] [CrossRef] [Google Scholar]
 Abt, H. A., Gomez, A. E., & Levy, S. G. 1990, ApJS, 74, 551 [NASA ADS] [CrossRef] [Google Scholar]
 Acernese, F., Agathos, M., Agatsuma, K., et al. 2015, Class. Quant. Grav., 32, 024001 [NASA ADS] [CrossRef] [Google Scholar]
 Adams, F. C., & Fatuzzo, M. 1996, ApJ, 464, 256 [NASA ADS] [CrossRef] [Google Scholar]
 Andrews, B. H., & Martini, P. 2013, ApJ, 765, 140 [NASA ADS] [CrossRef] [Google Scholar]
 Antonini, F., & Rasio, F. A. 2016, ApJ, 831, 187 [NASA ADS] [CrossRef] [Google Scholar]
 Askar, A., Szkudlarek, M., GondekRosi´nska, D., Giersz, M., & Bulik, T. 2017, MNRAS, 464, L36 [NASA ADS] [CrossRef] [Google Scholar]
 Baldry, I. K., Driver, S. P., Loveday, J., et al. 2012, MNRAS, 421, 621 [NASA ADS] [Google Scholar]
 Barrett, J. W., Gaebel, S. M., Neijssel, C. J., et al. 2018, MNRAS, 477, 4685 [NASA ADS] [CrossRef] [Google Scholar]
 Bastian, N., Covey, K. R., & Meyer, M. R. 2010, ARA&A, 48, 339 [NASA ADS] [CrossRef] [Google Scholar]
 Bate, M. R., & Bonnell, I. A. 2005, MNRAS, 356, 1201 [NASA ADS] [CrossRef] [Google Scholar]
 Baugh, C. M., Lacey, C. G., Frenk, C. S., et al. 2005, MNRAS, 356, 1191 [NASA ADS] [CrossRef] [Google Scholar]
 Baumgardt, H., Kroupa, P., & Parmentier, G. 2008, MNRAS, 384, 1231 [NASA ADS] [CrossRef] [Google Scholar]
 Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Kalogera, V., & Bulik, T. 2002, ApJ, 572, 407 [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. 2010a, ApJ, 714, 1217 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Dominik, M., Bulik, T., et al. 2010b, ApJ, 715, L138 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K.,Wiktorowicz, G., Fryer, C. L., Holz, D. E., & Kalogera, V. 2012, ApJ, 757, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Heger, A., Gladysz, W., et al. 2016a, A&A, 594, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belczynski, K., Holz, D. E., Bulik, T., & O’Shaughnessy, R. 2016b, Nature, 534, 512 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Repetto, S., Holz, D. E., et al. 2016c, ApJ, 819, 108 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Klencki, J., Meynet, G., et al. 2017, ArXiv eprints [arXiv:1706.07053] [Google Scholar]
 Belczynski, K., Askar, A., ArcaSedda, M. et al. 2018, A&A, 615, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bertelli, G., Bressan, A., Chiosi, C., Fagotto, F., & Nasi, E. 1994, A&AS, 106, B94 [Google Scholar]
 Bonnell, I. A., Clarke, C. J., & Bate, M. R. 2006, MNRAS, 368, 1296 [NASA ADS] [CrossRef] [Google Scholar]
 Boogaard, L. A., Brinchmann, J., Bouché, N., et al. 2018, A&A, 619, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bromm, V., Coppi, P. S., & Larson, R. B. 1999, ApJ, 527, L5 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Chen, H. Y., Holz, D. E., Miller, J., et al. 2017, ArXiv eprints [arXiv:1709.08079] [Google Scholar]
 Chini, R., Hoffmeister, V. H., Nasseri, A., Stahl, O., & Zinnecker, H. 2012, MNRAS, 424, 1925 [NASA ADS] [CrossRef] [Google Scholar]
 Chruslinska, M., Belczynski, K., Klencki, J., & Benacquista, M. 2018, MNRAS, 474, 2937 [NASA ADS] [CrossRef] [Google Scholar]
 Clark, P. C., Glover, S. C. O., Klessen, R. S., & Bromm, V. 2011, ApJ, 727, 110 [NASA ADS] [CrossRef] [Google Scholar]
 Daflon, S., & Cunha, K. 2004, ApJ, 617, 1115 [NASA ADS] [CrossRef] [Google Scholar]
 Daigne, F., Olive, K. A., Silk, J., Stoehr, F., & Vangioni, E. 2006, ApJ, 647, 773 [NASA ADS] [CrossRef] [Google Scholar]
 De Marchi, G., Paresce, F., & Pulone, L. 2007, ApJ, 656, L65 [NASA ADS] [CrossRef] [Google Scholar]
 de Mink, S. E., & Belczynski, K. 2015, ApJ, 814, 58 [NASA ADS] [CrossRef] [Google Scholar]
 de Mink, S. E., & Mandel, I. 2016, MNRAS, 460, 3545 [NASA ADS] [CrossRef] [Google Scholar]
 Djorgovski, S., Piotto, G., & Capaccioli, M. 1993, AJ, 105, 2148 [NASA ADS] [CrossRef] [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]
 Duchêne, G., & Kraus, A. 2013, ARA&A, 51, 269 [NASA ADS] [CrossRef] [Google Scholar]
 Dvorkin, I., Silk, J., Vangioni, E., Petitjean, P., & Olive, K. A. 2015, MNRAS, 452, L36 [NASA ADS] [CrossRef] [Google Scholar]
 Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Elbert, O. D., Bullock, J. S., & Kaplinghat, M. 2018, MNRAS, 473, 1186 [NASA ADS] [CrossRef] [Google Scholar]
 Eldridge, J. J., & Stanway, E. R. 2016, MNRAS, 462, 3302 [NASA ADS] [CrossRef] [Google Scholar]
 Eldridge, J. J., Stanway, E. R., Xiao, L., et al. 2017, PASA, 34, e058 [NASA ADS] [CrossRef] [Google Scholar]
 Fang, T., & Cen, R. 2004, ApJ, 616, L87 [NASA ADS] [CrossRef] [Google Scholar]
 Farr, W. M., & Mandel, I. 2018, Science, 361, aat6506 [NASA ADS] [CrossRef] [Google Scholar]
 Farr, W. M., Stevenson, S., Miller, M. C., et al. 2017, Nature, 548, 426 [NASA ADS] [CrossRef] [Google Scholar]
 Farr, B., Holz, D. E., & Farr, W. M. 2018, ApJ, 854, L9 [NASA ADS] [CrossRef] [Google Scholar]
 Fishbach, M., & Holz, D. E. 2017, ApJ, 851, L25 [NASA ADS] [CrossRef] [Google Scholar]
 Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, ApJ, 749, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Georgy, C., Ekström, S., Eggenberger, P., et al. 2013, A&A, 558, A103 [NASA ADS] [CrossRef] [EDP Sciences] [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]
 Greif, T. H., Schleicher, D. R. G., Johnson, J. L., et al. 2008, in LowMetallicity Star Formation: From the First Stars to Dwarf Galaxies, eds. Hunt L. K., Madden S. C., & Schneider R., IAU Symp., 255, 33 [NASA ADS] [CrossRef] [Google Scholar]
 Hamann, W.R., & Koesterke, L. 1998, A&A, 335, 1003 [NASA ADS] [Google Scholar]
 Heger, A., & Woosley, S. E. 2002, ApJ, 567, 532 [NASA ADS] [CrossRef] [Google Scholar]
 Heggie, D. C. 1975, MNRAS, 173, 729 [NASA ADS] [CrossRef] [Google Scholar]
 Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974 [NASA ADS] [CrossRef] [Google Scholar]
 Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543 [NASA ADS] [CrossRef] [Google Scholar]
 Husa, S., Khan, S., Hannam, M., et al. 2016, Phys. Rev. D, 93, 044006 [NASA ADS] [CrossRef] [Google Scholar]
 Jeans, J. H. 1902, Phil. Trans. R. Soc. London, Ser. A, 199, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Khan, S., Husa, S., Hannam, M., et al. 2016, Phys. Rev. D, 93, 044007 [NASA ADS] [CrossRef] [Google Scholar]
 Korn, A. J., Becker, S. R., Gummersbach, C. A., & Wolf, B. 2000, A&A, 353, 655 [NASA ADS] [Google Scholar]
 Kratter, K. M., & Matzner, C. D. 2006, MNRAS, 373, 1563 [NASA ADS] [CrossRef] [Google Scholar]
 Kratter, K. M., Matzner, C. D., & Krumholz, M. R. 2008, ApJ, 681, 375 [NASA ADS] [CrossRef] [Google Scholar]
 Kratter, K. M., Matzner, C. D., Krumholz, M. R., & Klein, R. I. 2010, ApJ, 708, 1585 [NASA ADS] [CrossRef] [Google Scholar]
 Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
 Kroupa, P., Tout, C. A., & Gilmore, G. 1993, MNRAS, 262, 545 [NASA ADS] [CrossRef] [Google Scholar]
 Kroupa, P., Weidner, C., PflammAltenburg, J., et al. 2013, in The Stellar and SubStellar Initial Mass Function of Simple and Composite Populations, eds. Oswalt T. D., & Gilmore G., 115 [Google Scholar]
 Kruckow, M. U., Tauris, T. M., Langer, N., Kramer, M., & Izzard, R. G. 2018, MNRAS, 481, 1908 [NASA ADS] [CrossRef] [Google Scholar]
 Lamberts, A., GarrisonKimmel, S., Clausen, D. R., & Hopkins, P. F. 2016, MNRAS, 463, L31 [NASA ADS] [CrossRef] [Google Scholar]
 LaraLópez, M. A., Hopkins, A. M., LópezSánchez, A. R., et al. 2013, MNRAS, 434, 451 [NASA ADS] [CrossRef] [Google Scholar]
 Larson, R. B. 1998, MNRAS, 301, 569 [NASA ADS] [CrossRef] [Google Scholar]
 Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Leitherer, C., Ekström, S., Meynet, G., et al. 2014, ApJS, 212, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415 [NASA ADS] [CrossRef] [Google Scholar]
 Mandel, I., & de Mink, S. E. 2016, MNRAS, 458, 2634 [NASA ADS] [CrossRef] [Google Scholar]
 Mandel, I., Haster, C.J., Dominik, M., & Belczynski, K. 2015, MNRAS, 450, L85 [NASA ADS] [CrossRef] [Google Scholar]
 Mandel, I., Farr, W. M., Colonna, A., et al. 2017, MNRAS, 465, 3254 [NASA ADS] [CrossRef] [Google Scholar]
 Mapelli, M., & Giacobbo, N. 2018, MNRAS, 479, 4391 [NASA ADS] [CrossRef] [Google Scholar]
 Mapelli, M., Giacobbo, N., Ripamonti, E., & Spera, M. 2017, MNRAS, 472, 2422 [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]
 Marks, M., & Kroupa, P. 2010, MNRAS, 406, 2000 [NASA ADS] [Google Scholar]
 Marks, M., Kroupa, P., & Baumgardt, H. 2008, MNRAS, 386, 2047 [NASA ADS] [CrossRef] [Google Scholar]
 Marks, M., Kroupa, P., Dabringhausen, J., & Pawlowski, M. S. 2012, MNRAS, 422, 2246 [NASA ADS] [CrossRef] [Google Scholar]
 Massey, P. 2003, ARA&A, 41, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Miller, M. C., & Lauburg, V. M. 2009, ApJ, 692, 917 [NASA ADS] [CrossRef] [Google Scholar]
 Moe, M., & Di Stefano, R. 2013, ApJ, 778, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Moe, M., & Di Stefano, R. 2017, ApJS, 230, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Moe, M., Kratter, K. M., & Badenes, C. 2018, ApJ, submitted, [arXiv:1808.02116] [Google Scholar]
 Nagashima, M., Lacey, C. G., Baugh, C. M., Frenk, C. S., & Cole, S. 2005, MNRAS, 358, 1247 [NASA ADS] [CrossRef] [Google Scholar]
 Nelson, D., Pillepich, A., Genel, S., et al. 2015, Astron. Comput., 13, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Öpik, E. 1924, Publ. Tartu Astrofiz. Observatory, 25, 1 [Google Scholar]
 O’Shaughnessy, R., Gerosa, D., & Wysocki, D. 2017, Phys. Rev. Lett., 119, 011101 [NASA ADS] [CrossRef] [Google Scholar]
 Pagel, B. E. J. 1989, in Evolutionary Phenomena in Galaxies, eds. Beckman J. E., & Pagel B. E. J., 201 [Google Scholar]
 Park, D., Kim, C., Lee, H. M., Bae, Y.B., & Belczynski, C. 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]
 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]
 Raghavan, D., McAlister, H. A., Henry, T. J., et al. 2010, ApJS, 190, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Repetto, S., & Nelemans, G. 2015, MNRAS, 453, 3341 [NASA ADS] [CrossRef] [Google Scholar]
 Rodriguez, C. L., Chatterjee, S., & Rasio, F. A. 2016a, Phys. Rev. D, 93, 084029 [NASA ADS] [CrossRef] [Google Scholar]
 Rodriguez, C. L., Haster, C.J., Chatterjee, S., Kalogera, V., & Rasio, F. A. 2016b, ApJ, 824, L8 [NASA ADS] [CrossRef] [Google Scholar]
 Salpeter, E. E. 1955, ApJ, 121, 161 [NASA ADS] [CrossRef] [Google Scholar]
 Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Sana, H., Le Bouquin, J.B., Lacour, S., et al. 2014, ApJS, 215, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Schneider, R., Graziani, L., Marassi, S., et al. 2017, MNRAS, 471, L105 [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]
 Stanway, E. R., Eldridge, J. J., & Becker, G. D. 2016, MNRAS, 456, 485 [NASA ADS] [CrossRef] [Google Scholar]
 Stevenson, S., Ohme, F., & Fairhurst, S. 2015, ApJ, 810, 58 [NASA ADS] [CrossRef] [Google Scholar]
 Stevenson, S., VignaGómez, A., Mandel, I., et al. 2017, Nat. Commun., 8, 14906 [NASA ADS] [CrossRef] [Google Scholar]
 Stone, N. C., Metzger, B. D., & Haiman, Z. 2017, MNRAS, 464, 946 [NASA ADS] [CrossRef] [Google Scholar]
 Tanaka, K. E. I., & Omukai, K. 2014, MNRAS, 439, 1884 [NASA ADS] [CrossRef] [Google Scholar]
 Tauris, T. M., Langer, N., & Podsiadlowski, P. 2015, MNRAS, 451, 2123 [NASA ADS] [CrossRef] [Google Scholar]
 Tauris, T. M., Kramer, M., Freire, P. C. C., et al. 2017, ApJ, 846, 170 [NASA ADS] [CrossRef] [Google Scholar]
 Tornatore, L., Borgani, S., Dolag, K., & Matteucci, F. 2007, MNRAS, 382, 1050 [NASA ADS] [CrossRef] [Google Scholar]
 Tremonti, C. A., Heckman, T. M., Kauffmann, G., et al. 2004, ApJ, 613, 898 [NASA ADS] [CrossRef] [Google Scholar]
 van den Heuvel, E. P. J., Portegies Zwart, S. F., de Mink, S. E., et al. 2017, MNRAS, 471, 4256 [NASA ADS] [CrossRef] [Google Scholar]
 Vangioni, E., Olive, K. A., Prestegard, T., et al. 2015, MNRAS, 447, 2575 [NASA ADS] [CrossRef] [Google Scholar]
 Vázquez, G. A., & Leitherer, C. 2005, ApJ, 621, 695 [NASA ADS] [CrossRef] [Google Scholar]
 Villante, F. L., Serenelli, A. M., Delahaye, F., & Pinsonneault, M. H. 2014, ApJ, 787, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Vink, J. S. 2017, A&A, 607, L8 [NASA ADS] [CrossRef] [EDP Sciences] [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]
 Wilkins, S. M., Hopkins, A. M., Trentham, N., & Tojeiro, R. 2008, MNRAS, 391, 363 [NASA ADS] [CrossRef] [Google Scholar]
 Xu, X.J., & Li, X.D. 2010, ApJ, 716, 114 [NASA ADS] [CrossRef] [Google Scholar]
 Yoshida, N., Omukai, K., Hernquist, L., & Abel, T. 2006, ApJ, 652, 6 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Comparison of the formation efficiencies per unit of starforming mass of CBM: BH–BH, BH–NS, and NS–NS between the two models with different initial distributions: either Sana et al. (2012, old) or Moe & Di Stefano (2017, new).
Merger rate densities and detection rates for the LIGO/Virgo O2 run sensitivity for three models: M10 (initial binary distributions from Sana et al. 2012), I1 (updated distributions from MD17), and I2 (I1 plus topheavy IMF for low metallicities from Marks et al. 2012).
Conversion coefficients for a cosmic SFR between the Salpeter and different Kroupalike IMFs.
Number of mergers of respective types, BH–BH, BH–NS, or NS–NS, in our samples simulated for models I1 (MD17 initial binary conditions) and I2 (additionally a metallicitydependent IMF; values given in parenthesis) for each of the 32 metallicities we computed.
All Figures
Fig. 1. Frequency f_{logP} of companions with q > 0.1 per decade of orbital period across log P (days) = 2–4. We compare the distribution f_{logP;comp} of all companions (binaries and tertiaries; dashed) to the distribution of only the inner binaries (solid), where is the probability that the companion is a member of the inner binary. The overall binary fraction and binary period distribution across log P (days) = 2–4 are similar between those adopted by dMB15 and B16 (green) and the values for M_{1} = 8 M_{⊙} (red) and M_{1} = 30 M_{⊙} (blue) primaries based on the multiplicity statistics presented in MD17. 

Open with DEXTER  
In the text 
Fig. 2. Cumulative distribution function of eccentricities. dMB15 and B16 adopted an eccentricity distribution p_{e} ∝ e^{−0.4} (green) based on a sample of massive spectroscopic binaries (Sana et al. 2012) that are dominated by shortperiod systems P < 20 days and therefore weighted towards lower eccentricities 〈e〉 = 0.3. MD17 (and references therein) found that massive binaries with intermediate (blue) and long (red) orbital periods are weighted towards higher eccentricities 〈e〉 = 0.6 and are sufficiently modelled by a powerlaw distribution p_{e} ∝ e^{0.8} with a turnover above e > 0.8. 

Open with DEXTER  
In the text 
Fig. 3. Cumulative distribution function of mass ratios. dMB15 and B16 adopted a uniform massratio distribution p_{q} ∝ q^{0} (green) based on a sample of massive spectroscopic binaries (Sana et al. 2012) that are dominated by shortperiod systems P < 20 days. MD17 (and references therein) found that massive binaries with intermediate periods (blue and red) are weighted significantly towards smaller mass ratios. 

Open with DEXTER  
In the text 
Fig. 4. Birth distributions of initial ZAMS binary parameters of the progenitors of BH–BH and NS–NS mergers. The histograms are normalized to unity for an easier comparison. NS–NS mergers dominate at Z = Z_{⊙} (upper panels) while BH–BH mergers dominate at Z = 0.1 Z_{⊙} and Z = 0.01 Z_{⊙} (lower panels). 

Open with DEXTER  
In the text 
Fig. 5. Birth distributions of initial ZAMS binary parameters of the progenitors of BH–NS mergers. The histograms are normalized to unity for an easier comparison. We note that the upper panels (Z = Z_{⊙}) are made based on a very small number of systems (<20) that form in solar metallicity. 

Open with DEXTER  
In the text 
Fig. 6. Formation !efficiency (i.e. number of systems formed per unit of starforming mass, ) of different types of CBM at different metallicities. We compute a grid of 32 different metallicities (indicated with filled squares at the BH–BH curves), ranging from 0.005 Z_{⊙} to 1.5 Z_{⊙}. We indicate three different models: M10 in blue (initial distributions of Sana et al. 2012), I1 in green (initial distributions of MD17) and I2 in red (same as I1 but with a topheavy IMF at low metallicities; see Eq. (2)). We note that model I2 is only different from I1 at Z < 0.004 = 0.2 Z_{⊙}. 

Open with DEXTER  
In the text 
Fig. 7. Source frame BH–BH merger rate density Gpc^{−3} yr^{−1} as a function of redshift for our three models: M10 (binary distributions of Sana et al. 2012), I1 (binary distributions of MD17), and I2 (same as I1 but with a metallicitydependent IMF). Upper panel: peak around z ≈ 2 in all the models corresponds to the maximum in the SFR (see text). Lower panel: range of redshifts accessible with the current and future sensitivity of LIGO detectors. We indicate the horizon redshifts for observation of an optimally inclined BH–BH merger with a total mass of ∼ 80 M_{⊙} (roughly the most massive systems in our simulations) for the sensitivities of O1 and O2 observing runs, and for the Advanced LIGO (AdLIGO). The detection distances for NS–NS mergers (Chen et al. 2017) were about d_{NSNS} = 70 Mpc for O1 and O2 sensitivity (Abbott et al. 2018). The local BH–BH merger rate density (within z < 0.02) in our models are about 203 (M10), 89 (I1), and 181 Gpc^{−3} yr^{−1} (I2). The current limits imposed by the existing LIGO detections are 12–219 Gpc^{−3} yr^{−1} (Abbott et al. 2017a). 

Open with DEXTER  
In the text 
Fig. 8. Detectorframe distribution of metallicities of BH–BH, BH–NS, and NS–NS merger progenitors. More than 90% of BH–BH merger detections in all our models are systems formed with metallicites lower than 10% of the solar metallicity, i.e. Z < 0.1 Z_{⊙}. Such a strong dependence of compact binary merger formation with metallicity (see also Fig. 6), especially in the case of BH–BH systems, indicates the significance of SFRs at different metallicites across the Universe for the compact binary merger rates (see discussion in Sect. 6.3). 

Open with DEXTER  
In the text 
Fig. 9. Relative numbers of massive stars ℛ formed within multiple initial mass ranges as a function of powerlaw exponent α_{3} for the highmass IMF component (M > 1 M_{⊙}). See Eq. (6) for the exact definition of ℛ. Upper panel: same SFR, irrespective of the IMF changes. The number of stars is thus calculated with respect to a fixed stellar mass formed. Lower panel: SFR corrected for the IMF changes (see Eq. (5)). The number of stars is thus calculated with respect to a fixed amount of UV radiation at 1500 Å wavelength that is constrained by the observations (Madau & Dickinson 2014). 

Open with DEXTER  
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.