Black hole mergers in the Illustris

The cosmic merger rate of stellar black hole binaries from the Illustris simulation

Abstract

The cosmic merger rate density of black hole binaries (BHBs) can give us an essential clue to constraining the formation channels of BHBs, in light of current and forthcoming gravitational wave detections. Following a Monte Carlo approach, we couple new population-synthesis models of BHBs with the Illustris cosmological simulation, to study the cosmic history of BHB mergers. We explore six population-synthesis models, varying the prescriptions for supernovae, common envelope, and natal kicks. In most considered models, the cosmic BHB merger rate follows the same trend as the cosmic star formation rate. The normalization of the cosmic BHB merger rate strongly depends on the treatment of common envelope and on the distribution of natal kicks. We find that most BHBs merging within LIGO’s instrumental horizon come from relatively metal-poor progenitors ( Z). The total masses of merging BHBs span a large range of values, from to M. In our fiducial model, merging BHBs consistent with GW150914, GW151226 and GW170104 represent , , and per cent of all BHBs merging within the LIGO horizon, respectively. The heavy systems, like GW150914, come from metal-poor progenitors ( Z). Most GW150914-like systems merging in the local Universe appear to have formed at high redshift, with a long delay time. In contrast, GW151226-like systems form and merge all the way through the cosmic history, from progenitors with a broad range of metallicities. Future detections will be crucial to put constraints on common envelope, on natal kicks, and on the BHB mass function.

keywords:
stars: black holes – gravitational waves – methods: numerical – stars: mass-loss – black hole physics

1 Introduction

The first direct detection of gravitational waves (GWs, Abbott et al. 2016b) opens a new perspective for the study of compact object (CO) binaries. Black hole binaries (BHBs) have been predicted and studied for a long time (e.g. Tutukov & Yungelson 1973; Thorne 1987; Schutz 1989; Kulkarni et al. 1993; Sigurdsson & Hernquist 1993; Portegies Zwart & McMillan 2000; Colpi et al. 2003; Belczynski et al. 2004), but the three events observed by LIGO so far (GW150914, GW151226, and GW170104) are the first observational confirmation of their existence.

Moreover, the two black holes (BHs) associated with GW150914 and one of the two BHs associated with GW170104 are surprisingly massive: , (Abbott et al., 2016c, a), and M (Abbott et al., 2017), respectively. If they are the remnants of massive stars, such massive BHs should have formed from relatively metal-poor ( Z) progenitors, which are expected to collapse directly to BHs (e.g. Mapelli et al. 2009; Mapelli et al. 2010, 2013; Belczynski et al. 2010; Spera et al. 2015). Dynamical processes (such as exchanges or runaway collisions in dense star clusters) might also contribute to enhancing the formation of massive BHBs similar to GW150914 and GW170104 (e.g. Ziosi et al. 2014; Chatterjee et al. 2017; Rodriguez et al. 2015, 2016; Mapelli 2016). Alternatively, GW150914 and GW170104 might be the result of primordial BHs born from gravitational collapse in the early Universe (e.g. Bird et al. 2016; Carr et al. 2016; Inomata et al. 2017).

Constraining the formation epoch and the birthplace of BHBs is one of the key points to interpret the nature of GW events associated with BHB mergers. This requires to model the formation and evolution of BHBs in a cosmological context. This task is currently a challenge, because of the huge dynamical range between the scale of cosmological structures (tens of Mpc) and the scale of binary evolution ( few AU). Moreover, since the progenitor’s metallicity appears to be crucial for the BH mass (Spera et al., 2015), any attempt to reconstruct the cosmic formation of BHBs should account for the local and global evolution of metallicity in the proper way.

For these reasons, only few authors attempted to put the formation of BHBs in a cosmological frame. Dominik et al. (2013) plant CO binaries into the cosmic history through a Monte Carlo-based algorithm. They generate a sample of galaxies based on a Press-Schechter like function (Fontana et al., 2006), adopt the average metallicity evolution described by Pei et al. (1999), and finally associate the CO binaries to a given redshift bin based on the cosmic star formation rate (SFR) evolution (Strolger et al., 2004). This gives a BHB merger rate density reaching its maximum at and then slowly decreasing down to .

Similarly, Belczynski et al. (2016b) generate isolated BHBs and then distribute them as a function of redshift, adopting an updated version of the cosmic SFR density and of the average metallicity evolution (Madau & Dickinson, 2014). This approach does not account for the mass-metallicity relation observed in galaxies (Maiolino et al., 2008). The resulting merger rate density peaks at . If only GW150914-like systems are considered, the distribution of the formation times of these systems is markedly bimodal with two peaks, one Gyr ago and the second one Gyr ago.

In contrast, Lamberts et al. (2016) account for the cosmological evolution through a Press-Schechter like formalism (Cole et al., 2008) with a redshift-dependent mass-metallicity relation (Ma et al., 2016). This ensures that the metallicity of a galaxy depends on its mass, consistent with the observations (Maiolino et al., 2008). Lamberts et al. (2016) do not recover the strongly bimodal birth-time distribution of GW150914-like systems reported by Belczynski et al. (2016b). Their predicted BHB merger rate is Gpc yr, significantly larger than inferred from LIGO observations ( Gpc yr, Abbott et al. 2016a). A conceptually similar approach was followed also by Dvorkin et al. (2016) and Elbert et al. (2017).

The formalism adopted by Dominik et al. (2013), Belczynski et al. (2016b), and even Lamberts et al. (2016) cannot give us detailed information on the evolution of the host galaxy of a CO binary. Thus, O’Shaughnessy et al. (2017b) follow a complementary approach: they start from a cosmological simulation and pick up four test galaxies, which they re-simulate at high resolution, by doing a “zoom-in”. Then, they add BHBs to the location of star forming particles in the simulation. They find a significantly higher merger rate per unit mass in dwarf galaxies than in Milky-Way-like galaxies.

Recently, Schneider et al. (2017) characterize the formation and coalescence sites of GW events, by coupling the metallicity-dependent binary population synthesis code SeBa (Portegies Zwart & Verbunt, 1996; Mapelli et al., 2013) with a simulation performed with the GAMESH pipeline (Graziani et al., 2015; Graziani et al., 2017). GAMESH interfaces an N-body simulation with a semi-analytic model for galaxy formation, and a radiative-transfer code. With this approach, Schneider et al. (2017) find that the observed GW events occur most likely in star forming galaxies with stellar mass M.

In this paper, we follow a new approach, complementary to previous work: we draw the cosmic history of the BHB merger rate by coupling up-to-date population synthesis simulations of BHBs with the public Illustris-1 cosmological simulation (Vogelsberger et al., 2014b). The Illustris-1 is the highest resolution hydrodynamical simulation run in the frame of the Illustris project (Vogelsberger et al., 2014a). In the following, we refer to the Illustris-1 simply as Illustris. The Illustris box (length Mpc comoving) is considerably larger than the one adopted by Schneider et al. (2017), ensuring that we are considering a less biased portion of the Universe, even if with lower resolution.

We plant our BHBs in the cosmological simulation through a Monte Carlo approach, based on the metallicity of star particles. Our BHBs were generated by evolving isolated stellar binaries with a new version of the BSE code (Hurley et al., 2002), which includes up-to-date recipes for stellar evolution and stellar winds (Vink et al., 2011; Chen et al., 2015). Moreover, we include the effect of pulsational pair-instability and pair-instability supernovae (Woosley, 2017) which were neglected in most previous studies. This approach allows us to follow the merger history of BHBs, accounting for the evolution of their environment.

2 Methods

To reconstruct the cosmic history of BHB mergers, we couple the Illustris simulation with a large set of population-synthesis simulations of isolated binaries. The main ingredients of our model are the following.

2.1 The BHBs

We simulate the evolution of isolated stellar binaries through an updated version of the BSE code (Hurley et al., 2000, 2002). The changes with respect to the original version of BSE are described in a companion paper by Giacobbo et al. (in prep.). Here we summarize the most important prescriptions. Stellar winds have been updated based on the equations described in Belczynski et al. (2010). Namely, a treatment of stellar winds following Vink et al. (2001) and Vink & de Koter (2005) is included for O-type and Wolf-Rayet stars, respectively. In this model, mass loss by stellar winds depends on metallicity, both in the main sequence (MS) and in later evolutionary stages.

With respect to Belczynski et al. (2010), there is one crucial update: we take into account the dependence of the mass loss on the electron-scattering Eddington ratio (Gräfener & Hamann, 2008; Vink et al., 2011; Vink, 2016). Following Chen et al. (2015), the mass loss scales as (where is the star metallicity) with if the electron-scattering Eddington ratio of a star is , and if . This ensures that the dependence of mass loss on metallicity almost vanishes if the star is radiation-pressure dominated. With this relatively small change, we obtain a mass spectrum of BHs similar to the one published by Spera et al. (2015) and based on the PARSEC stellar evolution tracks (Bressan et al., 2012; Tang et al., 2014; Chen et al., 2015).

Our new version of BSE also includes new fitting formulas for the core radii, as described in Hall & Tout (2014). This is a crucial ingredient for the study of BHBs, because the fate of a common envelope phase depends on the core radius.

Furthermore, we included in BSE new recipes for core-collapse supernovae (SNe). In particular, we implemented both the rapid (R) and the delayed (D) models for SN explosion presented by Fryer et al. (2012). In Appendix A, we detail the prescriptions for the CO mass in the rapid and in the delayed SN model. Finally, we added a formalism to account for pair-instability and pulsational pair-instability SNe, following Spera & Mapelli (2017) (see also Belczynski et al. 2016c; Woosley 2017).

With the new code, we ran six sets of population-synthesis simulations. The details of the six sets are given in Table 1. In particular, the simulation set labelled as ‘R’ adopts the rapid SN model, while all the others (labelled as ‘D’) adopt the delayed model for core-collapse SNe. Pulsational pair instability and pair instability SNe are included in all runs.

For the common envelope (CE) phase, we use the same formalism as described in Hurley et al. (2002), which depends on two free parameters, and . According to this formalism, is the fraction of binding energy converted into kinetic energy of the envelope, while describes the geometry of the envelope. In the formalism by Hurley et al. (2002), and always appear as their product . In simulations D, R, DHG, and DK we use , . The latter choice of is quite well motivated for massive stars (e.g. Xu & Li 2010; Loveridge et al. 2011). In simulation D0.02 we assume , , while in simulation D1.5 we use and .

The treatment of Hertzsprung gap (HG) donors was found to be critical in previous studies (e.g. Dominik et al. 2012). A HG star lacks a steep density gradient between the core and envelope. Thus, its response to a CE should be similar to that of a MS star (Ivanova & Taam, 2004). In the standard version of BSE, MS donors entering a CE phase are forced to merge with the accretor, while HG donors are allowed to survive the CE phase. In our run DHG, we adopted the default setting of BSE allowing HG donors to survive a CE phase. In all the other simulations we modified BSE, by imposing that a HG donor merges with its companion if they enter a CE phase.

Finally, the natal kick of the CO is another essential ingredient, because it can unbind a binary. There are no conclusive observational constraints on the natal kick of BHs, even if some recent studies indicate that high-velocity kicks are possible (Repetto et al., 2012; Repetto et al., 2017; O’Shaughnessy et al., 2017a). Thus, we draw the natal kicks from the Maxwellian distribution described in Hobbs et al. (2005), with dispersion km s. This distribution was obtained from the proper motions of 233 isolated Galactic pulsars.

In model DK, we assume that BH kicks follow the same distribution as neutron star (NS) kicks. In all the other models, we scale the velocities drawn from this distribution by the amount of fallback, according to:

(1)

where is the natal kick for the BH, is the natal kick for a NS (drawn from the distribution proposed by Hobbs et al. 2005), and (ranging from 0 to 1) is the amount of fallback on the proto-NS (Fryer et al., 2012; Spera et al., 2015). The definition of depends on the adopted core-collapse SN prescription. In Appendix A, we detail the values of for the two considered core-collapse SN models.

For each of the six simulation sets described in Table 1, we simulate 12 sub-sets with metallicity 0.0004, 0.0008, 0.0012, 0.0016, 0.002, 0.004, 0.006, 0.008, 0.012, 0.016, and 0.02. Throughout the paper, we define solar metallicity as Z. Thus, the 12 sub-sets correspond to metallicity , 0.02, 0.04, 0.06, 0.08, 0.1, 0.2, 0.3, 0.4, 0.6, 0.8, and 1.0 Z. In each sub-set we simulate stellar binaries. Thus, each of the six sets of simulations is composed of massive binaries. The mass of the primary () is randomly drawn from a Kroupa initial mass function (Kroupa, 2001) ranging1 from 5 to 150 M, and the mass of the secondary () is sampled according to the distribution (where ) in a range . The orbital period and the eccentricity are randomly extracted from the distribution , with , and , with , as suggested by Sana et al. (2012).

Name SN HG Kick
D Delayed 1.0 0.1 new F12
R Rapid 1.0 0.1 new F12
DHG Delayed 1.0 0.1 BSE F12
DK Delayed 1.0 0.1 new H05
D0.02 Delayed 0.2 0.1 new F12
D1.5 Delayed 3.0 0.5 new F12

Column 1: model name; column 2: SN model (delayed and rapid from Fryer et al. 2012); column 3: value of ; column 4: value of ; column 5: treatment for HG stars (‘BSE’ means same treatment as in BSE, ‘new’ means that we force all CE binaries with a HG donor to merge); column 6: model for the SN kick. H05 means that we use the distribution from Hobbs et al. (2005). F12 means that we rescale the natal kicks by the fallback, as described in Fryer et al. (2012). See also equation 1 and the text for details.

Table 1: Properties of the population-synthesis simulations.

2.2 The Illustris

The Illustris simulation covers a comoving volume of , and has an initial dark matter and baryonic matter mass resolution of and M, respectively (Vogelsberger et al., 2014b, a). At redshift zero the softening length is pc, while the smallest hydrodynamical cells have a length of 48 pc. The large size of the Illustris’ box ensures that we are modelling an unbiased portion of the Universe, satisfying the cosmological principle. The main drawback is that the population of dwarf galaxies is heavily under-resolved. In Appendix B, we estimate the impact of resolution on our main results, by comparing the Illustris-1 with the lower-resolution Illustris-3 simulation. Moreover, in a companion paper (Schneider et al., 2017), we follow a complementary approach: we combine our population-synthesis models with the GAMESH simulation, which has a box of only , but much higher resolution, in order to quantify the contribution of dwarf galaxies to the BHB merger rate.

The Illustris was run with the moving mesh code AREPO, to solve the inviscid Euler equations (Springel, 2010). The Illustris includes a treatment for sub-grid physics (cooling, star formation, SNe, super-massive BH formation, accretion and merger, AGN feedback, etc), as described in Vogelsberger et al. (2013). The model of sub-grid physics adopted in the Illustris is known to produce a mass-metallicity relation (Genel et al., 2014; Genel, 2016) which is sensibly steeper than the observed one (see the discussion in Vogelsberger et al. 2013 and Torrey et al. 2014). Moreover, the simulated mass-metallicity relation does not show the observed turnover at high stellar mass ( M). In Appendix B we estimate that the impact of these differences between simulated and observed mass-metallicity relation on the BHB merger rate is per cent.

As for the cosmology, the Illustris adopts WMAP-9 results for the cosmological parameters (Hinshaw et al., 2013), that is , , , and km s Mpc, with .

Through the web-based interface (API) made available by the the Illustris project (http://www.illustris-project.org/), we downloaded the stellar particles in each snapshot, including information on their formation time, initial mass, and metallicity. A total of 108 snapshots have stellar particles, from redshift to 0. More details on the Illustris can be found in the presentation papers (Vogelsberger et al., 2014b, a), in the release paper (Nelson et al., 2015) and on the aforementioned website.

2.3 Planting BHBs into a cosmological simulation

We wrote a Monte Carlo code to associate the simulated BHBs to the Illustris, performing the following operations.

For each of the simulation sets listed in Table 1, we extract information on those BHBs merging within a Hubble time. Namely, we store the BH masses and the delay time between the formation of the progenitor stellar binary and the merger of the BHB. We also store information on the total initial stellar mass of each sub-set of simulations with the same metallicity (including binary systems which do not evolve into BHBs).

We read each stellar particle from the Illustris only once, when it first appears in the snapshots. We store information on its initial mass , formation redshift , and metallicity . We then find the metallicity that best matches among the 12 metallicities simulated with BSE2.

We then associate to each Illustris’ particle a number of merging BHBs, randomly extracted from the sub-set with the best-matching metallicity, based on the following algorithm:

(2)

where is the initial stellar mass of the Illustris’ particle and is the initial stellar mass in the BSE sub-set with the selected metallicity. In our calculations, . is the number of merging BHBs within the simulated sub-set of initial stellar mass . is a correction factor, accounting for the fact that we actually simulate only primaries with M, neglecting lower mass stars. Finally, accounts for the fact that we simulate only binary systems, whereas a fraction of stars are single. Here we assume that 50 per cent of stars are in binaries, thus . We note that is only a scale factor and our results can be rescaled to a different a posteriori. We notice that is (by definition) the number of merging BHBs per unit stellar mass at a given metallicity.

With this procedure, we associate to each Illustris’ particle a number of randomly selected merging BHBs whose progenitors have metallicity .

We then estimate the look-back time of the merger () of each BHB in the randomly selected sample as

(3)

where is the time between the formation of the progenitor stellar binary and the merger of the BHB, and is the look back time at which the Illustris’ particle has formed, calculated as

(4)

where the cosmological parameters are set to WMAP-9 values (for consistency with the Illustris) and is the formation redshift of the Illustris’ particle.

According to this definition, is also a look back time: it tells us how far away from us the BHB merged. For our analysis, we consider only BHBs with , i.e. we do not consider BHBs that will merge in the future.

We repeat the same procedure for each of the six simulation sets in Table 1 and we obtain six different models of the cosmic BHB merger evolution.

3 Results

Figure 1: Left axis: cosmic merger rate density of BHBs () in the comoving frame, as a function of the look-back time (bottom axis) and of the redshift (top axis) in our models. Red solid line: D (fiducial model); black dashed line: R; violet dash-dot line: DHG; orange dashed line: DK; blue dotted line: D0.02; green dash-dot line: D1.5. Green shaded area: BHB merger rate inferred from LIGO detections (Abbott et al., 2016a). Right axis: cosmic SFR density from the Illustris (grey thin solid line), as a function of the look-back time (bottom axis) and of the redshift (top axis).

3.1 Merger rate

Figure 1 shows the cosmic BHB merger rate density () in the comoving frame, derived from our simulations. To obtain the BHB merger rate shown in this Figure, we extracted the number of BHB mergers () per time bin (each bin spanned over Myr from to ) and then we did the following simple conversion:

(5)

where Mpc is the size of the Illustris box (in the comoving frame) and is the size of the time bin (10 Myr).

Name
[Gpc yr] [Gpc yr]
D 125 181
R 155 228
DHG 572 772
DK 20 29
D1.5 145 181
D0.02 278 279

Column 1: model name; column 2: present-time BHB merger rate density; column 3: BHB merger rate density at .

Table 2: Comoving BHB merger-rate density at redshift and (corresponding to and 2.43 Gyr, respectively).

From Fig. 1 it is apparent that the overall behaviour of the merger rate is the same for all considered BHB models, with the partial exception of D0.02 (see Table 1 for details about the models). The behaviour of the merger rate density as a function of time depends only on the SFR (given by the Illustris and thus common to all BHB models) and on the delay between the formation time of a stellar binary and the merger time of the BHB born from the stellar binary (which depends on the BHB model).

The shape of the merger rate density in Fig. 1 resembles the one of the cosmic SFR density (e.g. Madau & Dickinson 2014) with a peak at Gyr (i.e. ). The decrease of the merger rate approaching is more gentle than the decrease of the SFR density, because of BHBs that formed at high redshift but merge with a delay of several Gyr (see the next section).

The main difference between the considered BHB models is the normalization of the merger rate density, which depends on the BHB merger efficiency (see Table 2 for details). In particular, the present-time merger rate density ranges from Gpc yr to Gpc yr in models D, D1.5 and R. Model D0.02 results in a factor of two larger rate ( Gpc yr). Finally, the rate is much higher ( Gpc yr) in model DHG and much lower ( Gpc yr) in model DK (see Table 2).

The BHB merger rate density inferred from the first LIGO observations (O1 run) is Gpc yr (Abbott et al., 2016a). While this paper was being reviewed, the inferred rate was updated to Gpc yr, based on the first results of the O2 run (Abbott et al., 2017). Thus, the present-day merger rate density of models D, R, DK, and D1.5 are consistent with observations, while D0.02 is slightly above the observed range and DHG gives a much higher rate. Thus, models in which HG stars can survive a CE phase (DHG) are not consistent with the observed merger rate, unless natal kicks are much higher than assumed.

Figure 2: Metallicity of progenitors of merging BHBs in the fiducial model (D). Upper panel: metallicity versus formation time of the stellar progenitors (). Lower panel: metallicity of the stellar progenitors versus merger time of the BHBs (). Both and are expressed as look-back time. The colour-coded map (in logarithmic scale) indicates the number of merging BHBs per cell. The black lines are isocontours enclosing a number of merging BHBs ranging from to (as indicated by the black labels).
Figure 3: Upper (lower) panel: total mass (chirp mass) of merging BHBs as a function of in the fiducial model. is expressed as look-back time. The colour-coded map (in logarithmic scale) indicates the number of merging BHBs per cell. The black lines are isocontours enclosing a number of merging BHBs ranging from to (as indicated by the black labels).

Natal kicks have a strong impact on the BHB merger rate: is a factor of lower in run DK than in run D, which differ only by the kick prescription. In run D the magnitude of the kick depends on the amount of fallback. In run DK all BHs receive a natal kick, drawn from the same distribution as Galactic single pulsars (Hobbs et al., 2005). The merger rate density of both run D and DK are consistent with current observations, but future detections might be able to discriminate between such models.

Runs D, D0.02 and D1.5 differ by the choice of the and CE parameters. Unlike the other considered effects (SN model, SN kicks and HG treatment), the choice of CE parameters affects not only the normalization but also the shape of the BHB merger rate density as a function of time. Since the SFR history is the same for all models, this difference indicates that models with different CE parameters have also different distributions for the delay time.

The effect of the choice of is much more important at high redshift than at low redshift. At , the difference between runs D and D1.5 is negligible, while the difference between run D and D0.02 is about a factor of two.

Finally, runs R and D have similar merger rates (within a factor of 1.3). This indicates that the choice of the core-collapse SN model (rapid or delayed) does not affect the BHB merger rate significantly. In the following, we will consider run D as our fiducial model.

3.2 Formation time, progenitor’s metallicity and BHB masses

In this section, we discuss the main properties of merging BHBs in the Illustris simulation. We consider only our fiducial run D.

Figure 2 maps the metallicity of the stellar progenitors of merging BHBs. In the upper panel the metallicity is plotted against the formation time of the stellar progenitors, while the lower panel shows the metallicity versus the merger time of BHBs. From the comparison between the two panels, it is apparent that a large fraction of metal-poor systems which formed at high redshift merge at relatively low redshift with a long delay time. For example, BHBs with progenitor metallicity Z merge at redshift in the simulation, but only of them form at redshift . This implies that a significant number of merging BHBs visible in the LIGO instrumental horizon were born in the high-redshift Universe and possibly in a metal-poor environment.

Figure 2 also shows that the stellar progenitors of BHBs have all possible metallicities ranging from up to . The most common metallicity of BHBs merging at low redshift is , i.e. significantly sub-solar. This result comes from a combination of two factors. Firstly, relatively metal-poor stars form efficiently even at , as expected from the mass-metallicity relation. Secondly, many BHBs merging at formed at high-redshift, where low metallicity was more common. Mergers associated with solar or super-solar metallicity are strongly suppressed in our models, because stellar radii are larger at higher metallicity, causing early mergers of massive stars before they become BHBs.

At higher redshift, the percentage of merging BHBs born in metal-poor environments increases, and the contribution of lower metallicities becomes more important. However, we stress that even at Gyr there is a non-negligible fraction of systems with metallicity Z.

Figure 3 shows the behaviour of the total mass (, where and are the mass of the primary and secondary BH, respectively) and of the chirp mass () as a function of the merger time. The distribution of total masses (chirp masses) peaks at (), but we find merging systems with total masses (chirp masses) ranging from to M ( to M). From this Figure it is apparent that there is nearly no dependence of the BHB mass on the merging time. This is primarily a consequence of the broad distribution of delay times (see next Section).

Figure 4: Distribution of the delay time (, estimated as the time elapsed between the formation of a stellar binary and the merger of the BHB formed from this stellar binary) for the simulated BHBs in runs D (red solid line), R (black dashed line), DK (orange dashed line), D1.5 (green dash-dot line) and D0.02 (blue dotted line). Only BHBs merging within the mass dependent instrumental horizon of LIGO are shown. Purple dotted line: . Note that (unlike and ) is not a look-back time. The value on the axis is the number of simulated BHBs per time bin ( Myr).

3.3 BHBs merging within the LIGO horizon

In this section, we focus only on simulated BHBs that merge within the LIGO instrumental horizon, defined as the luminosity distance at which GWs from a face-on, equal-mass, overhead binary would be detected with signal-to-noise ratio of 8 (Abbott et al., 2016c). To account for the dependence of the instrumental horizon on the BHB mass, we use the curve reported in Fig. 4 of Abbott et al. (2016c) for the 2015-2016 LIGO sensitivity (left-hand panel).

In order to extract from our simulations all BHBs merging inside the LIGO instrumental horizon, we check whether the luminosity distance at the time of merger is smaller than the instrumental horizon for a BHB with the same total mass, as given in Fig. 4 of Abbott et al. (2016c).

Figure 5: Red solid line: distribution of formation time for the simulated BHBs in the fiducial model (run D). Blue dashed line: distribution of merger time for the simulated BHBs in the fiducial model. Only BHBs merging within the LIGO instrumental horizon are shown. Bottom axis: and expressed as look-back time. Top axis: and expressed as redshift. The value on the axis is the number of simulated BHBs per time bin ( Myr).

In this section we also compare the properties of our fiducial model with the other runs. Figure 4 shows the delay time distribution for the BHB merging within the LIGO horizon. This distribution depends only on the BSE models and is not affected by the cosmological simulation. In run D, the delay distribution matches the behaviour found in previous studies (Belczynski et al., 2016b; Lamberts et al., 2016), but only for Gyr. For longer delay times, the distribution flattens considerably: it becomes nearly independent of time. This explains why a large fraction of metal-poor systems formed at high redshift merge within the LIGO horizon (see Fig. 2 and Section 3.2). Runs R (adopting the rapid SN model) and DK (assuming large BH kicks) behave exactly the same as run D.

In contrast, the delay time distributions of runs D1.5 and especially D0.02 (which differ from run D only for the CE parameters) are significantly different. They decrease more steeply than at short delay time. This produces a much lower number of BHB mergers with delay time . For Gyr, the number of BHB mergers in runs D1.5 and D0.02 becomes significantly higher than that of run D. This difference in the distribution of delay times explains why the BHB merger rate density in run D is higher (lower) at high (low) redshift than that of runs D1.5 and D0.02.

Figure 5 shows the distribution of formation times and the merger times of BHBs that merge within the LIGO horizon in our fiducial model. In Fig. 5, the merger time peaks at . This results from the convolution between the LIGO horizon (which depends on the BHB mass) and the increase of the merger rate as the redshift increases (see Fig. 1). Interestingly, the estimated redshift of GW150914 and GW151226 is , while is the redshift of the candidate event (LVT151012) and of the most recent detection (GW170104).

The distribution of and in the other models is similar to the one shown in Fig. 5. The main differences arise from the distribution of delay times. For this reason, the distribution of formation times in run D0.02 has a much higher peak at high redshift than that of run D.

Figure 6: Distribution of total masses (upper panel) and chirp masses (lower panel) of the simulated BHBs merging within the LIGO horizon. Red solid line: run D; black dashed line: R; orange dashed line: DK; green dash-dot line: D1.5; blue dotted line: D0.02. The value on the axis is the number of simulated BHBs per mass bin ( M for both the chirp and the total mass).

Figure 6 shows the distribution of total masses and chirp masses of the simulated BHBs that merge within the LIGO horizon. The chirp masses (total masses) range between 3 and 35 M (6 and 82 M).

Runs D1.5 and D0.02 are significantly different from the others, with a peak at lower masses ( M, M) and a dearth of massive systems. Runs with are more similar between each other. For high total masses ( M) there is no appreciable difference between runs D and R (which differ for the SN model). At lower masses, run R (assuming a rapid SN model) has a peak at , which is completely absent if the delayed SN model is assumed. Moreover, in run D and in the other runs assuming a delayed SN models, BHs with mass down to 3 M are allowed to form, while no BHs with mass M form in run R.

Finally, the distribution of BHB masses in run DK is very similar to that of run D, except for one important detail: there is no peak for total BHB masses (chirp masses) (). This comes from the fact that in run DK all BHs receive a strong natal kick (regardless of their mass), while in the other runs the kick is modulated by the fallback.

Figure 6 is significantly affected by the fact that the LIGO horizon depends on the BHB mass. More massive binaries can be observed also if they merge at higher redshift ( if M). Thus, their contribution to Fig 6 is enhanced with respect to that of lighter BHBs.

Future LIGO-Virgo detections will enable us to reconstruct total mass and chirp mass distributions. Hopefully, the observed distributions will be able to discriminate between different models, indicating which one captures the main physics of BHB evolution.

Figure 7: Distribution of metallicity () of the stellar progenitors of the simulated BHBs that merge within the LIGO horizon. Red solid line: run D; black dashed line: R; orange dashed line: DK; green dash-dot line: D1.5; blue dotted line: D0.02. The value on the axis is the number of simulated BHBs per metallicity bin (). The step-like features in the plot correspond to the metallicity groups in the BSE simulations.

Finally, Fig. 7 shows the distribution of metallicity of the stellar progenitors of the simulated BHBs merging within LIGO horizon. In runs D, DK and R, the metallicity of BHB progenitors peaks in the range . The metallicity distribution drops for .

Also in this case, runs with different CE parameters (D1.5 and D0.02) have a different trend. In run D0.02 the metallicity of BHB progenitors peaks at , significantly higher than in runs with . In contrast, in run D1.5 the most metal-poor systems are more efficient in producing merging BHBs than in the other runs. The step-like features clearly visible in Fig. 7 are model artifacts, due to the fact that we simulated BHBs only in 12 metallicity bins (see Section 2).

Name GW150914 LVT151012 GW151226 GW170104
D 6.0% 9.9% 3.4% 12.1%
R 5.3% 8.9% 3.3% 10.8%
DK 9.8% 11.0% 2.5% 12.5%
D1.5 8.4% 13.1% 5.0% 10.4%
D0.02 2.8% 12.9% 4.9% 5.2%

Column 1: model name; columns 2, 3, 4, and 5: percentage of simulated BHBs with mass consistent with GW150914, LVT151012, GW151226 and GW170104, respectively, which merge within the LIGO horizon. This percentage is calculated over all simulated BHBs that merge within the LIGO horizon.

Table 3: Percentage of simulated BHBs merging in the LIGO horizon with mass consistent with the detected events.
Figure 8: Distribution of the formation time (, red solid line) and of the merger time (, blue dashed line) for the simulated BHBs matching the mass of GW150914 (top left), LVT151012 (top right), GW151226 (bottom left) and GW170104 (bottom right) in the fiducial model (run D). Thin lines indicate all simulated systems matching the mass of GW150914, LVT151012, GW151226, and GW170104, while thick lines indicate only systems that merge within LIGO’s horizon. Both and are expressed in look-back time (bottom axis) and redshift (top axis). The value on the axis is the number of simulated BHBs per time bin ( Myr for both the merger time and the formation time).

3.4 GW150914, LVT151012, GW151226, and GW170104-like systems

Finally, we focus on the properties of simulated systems that match the three observed GW events, GW150914 (, , Abbott et al. 2016a), GW151226 (, , Abbott et al. 2016a), and GW170104 (, , Abbott et al. 2017), and the fourth possible signal, LVT151012 (, , Abbott et al. 2016a). From our simulations, we extract all merging BHBs with masses consistent with those of the BHs associated with GW150914, GW151226, GW170104, and possibly LVT151012 (within 90% credible intervals, as given by Abbott et al. 2016a). Table 3 shows the percentage of merging BHBs that have masses consistent with GW150914, LVT151012, GW151226 and GW170104, if we consider only BHBs that merge within the LIGO horizon.

BHB mergers consistent with GW150914, LVT151012, GW151226, and GW170104 exist in all of our models. In particular, we expect %, %, % and % of events, inside the LIGO horizon, to have mass consistent with GW150914, LVT151012, GW151226, and GW170104 respectively, in run D (see Table 3 for the other runs). Interestingly, GW150914-like events appear to be more common than GW151226-like ones within the 2015-2016 LIGO’s horizon. Also, LVT151012-like events seem to be the most common ones within LIGO’s horizon, but this is mainly due to the fact that the mass of LVT151012 is much more uncertain than that of both GW150914 and GW151226, thus more systems fall within the allowed window.

Figure 9: Metallicity of progenitors of systems matching the mass of GW150914 in the fiducial model (D). Upper panel: metallicity versus formation time of the stellar progenitors (). Lower panel: metallicity of the stellar progenitors versus merger time of the BHBs (). Both and are expressed as look-back time. The colour-coded map (in logarithmic scale) indicates the number of merging BHBs per cell. The black lines are isocontours enclosing a number of merging BHBs ranging from to (as indicated by the black labels).

Figure 8 shows the behaviour of the merger time and the formation time for simulated systems like GW150914, LVT151012, GW151226 and GW170104 in our fiducial model. It is apparent that the vast majority of GW150914-like systems form at high redshifts (), with a peak at Gyr (corresponding to in the Illustris’ cosmology). Mergers also peak at high redshift ( Gyr). However, the long delay time between formation and merger for several systems causes to decrease more gently than when approaching .

If we consider only GW150914-like systems that merge within the LIGO horizon, their formation time still peaks at high redshift (), even if the slope of is much less steep. This is consistent with Schneider et al. (2017), who predict that GW150914-like events form preferentially at redshift .

In the nearby Universe, we expect a factor of (and up to ) higher merger rate of GW150914-like systems than their formation rate. In run D, the birth rate of new GW150914-like systems at (the redshift associated with the LIGO detection) is only 0.4 Gpc yr, whereas the expected merger rate of GW150914-like systems at redshift is Gpc yr (see Table 4). This result is perfectly consistent with the estimates from the LIGO-Virgo collaboration ( Gpc yr, Abbott et al. 2016a).

Fig. 9 shows that this effect is linked to the metallicity of GW150914 progenitors. In fact, most GW150914-like systems merging at low redshift have , but the number of progenitor systems forming at low redshift in this metallicity range is a factor of lower than the number of merging GW150914-like systems.

Name GW150914 LVT151012 GW151226 GW170104
[Gpc yr] [Gpc yr] [Gpc yr] [Gpc yr]

D
5.2 20.9 8.3 15.9
R 5.2 20.7 9.3 16.0
DK 1.5 3.5 1.0 2.9
D1.5 7.7 21.8 10.9 12.2
D0.02 5.1 37.4 18.7 10.5
A16

Column 1: model name; columns 2, 3, 4, and 5: merger rate of simulated BHBs with mass and merger redshift consistent with GW150914, LVT151012, GW151226, and GW170104, respectively. The last line shows the rate inferred from LIGO observations (see Table II of Abbott et al. 2016a) for GW150914, LVT151012, and GW151226.

Table 4: Expected merger rate for GW150914, LVT151012, GW151226 and GW170104-like events. To estimate these values we consider only systems merging at redshift consistent with the LIGO detections.

We find a similar trend for GW170104-like systems. In run D, the birth rate of new GW170104-like systems at (the redshift associated with the LIGO detection) is only 0.4 Gpc yr, whereas the expected merger rate of GW170104-like systems at redshift is Gpc yr (see Table 4).

In contrast, the shift between formation time and merger time is rather negligible for systems like GW151226. From our simulations, the expected merger rate of GW151226-like systems at redshift (the redshift associated with the LIGO detection) is Gpc yr (see Table 4), similar to their birth rate ( Gpc yr at ). This merger rate is consistent with the estimates from the LIGO-Virgo collaboration, even if close to the low tail ( Gpc yr, Abbott et al. 2016a).

LVT151012-like events behave in a similar way to GW151226, with a mild offset between their current merger time and their current formation time (Fig. 8). The merger rate of LVT151012-like events is Gpc yr at redshift (i.e. the redshift associated with the possible LIGO signal), significantly higher than that of the other two events. This happens because the predicted BHB merger rate at is higher by a factor of than that at and because the mass of LVT151012 is more uncertain than that of GW151226 and GW150914 (thus, more simulated systems are consistent with it). This merger rate is also consistent with the estimates from the LIGO-Virgo collaboration ( Gpc yr, Abbott et al. 2016a).

Figure 10: Metallicity distribution of simulated GW150914-like (blues dashed line), LVT151012-like (green dash-dot line), GW151226-like (red solid line) and GW170104-like systems (violet dotted line) in the fiducial model (run D). The value on the axis is the number of simulated BHBs per metallicity bin (). The step-like features in the plot correspond to the metallicity groups in the BSE simulations.

Figure 10 shows the metallicity distribution of the stellar progenitors of GW150914-like, GW151226-like, LVT151012-like, and GW170104-like systems in our fiducial model. As expected, GW150914-like systems form only at relatively low metallicity: . Metal-poor stars ( Z) are the most common progenitors also for GW170104-like systems. In contrast, GW151226-like systems form nearly at all metallicities up to Z. LVT151012-like systems also form with a broad range of metallicities (up to Z).

4 Discussion: Comparison with previous work and caveats

The method we followed in this paper ensures that we account for the cosmic SFR and for the mass-metallicity relation, at least within the limitations of state-of-the-art cosmological simulations. This is a crucial point for understanding the cosmic evolution of BHB mergers, since the mass of BHs is expected to depend on the metallicity of the progenitor stars.

The cosmic BHB merger rate we obtain from our models (Fig. 1) approximately follows the same trend as the cosmic SFR, suggesting that the slope of the curve is primarily set by star formation. The peak of the BHB merger rate is at , corresponding to the peak of the cosmic SFR (Madau & Dickinson, 2014). This is in reasonable agreement with previous results (e.g. Dominik et al. 2013; Belczynski et al. 2016b).

The normalization of the BHB merger rate in Fig. 1 strongly depends on the treatment of CE (especially for HG stars) and on the distribution of natal kicks. In particular, models in which HG donors can survive the CE phase disagree with the BHB merger rate estimated from LIGO observations, unless very large natal kicks are assumed for most BHs (in agreement with Belczynski et al. 2016b). Different distributions of the natal kicks result in a factor of different BHB merger rate, given the large uncertainties in the distribution of BH natal kicks. Future observations by LIGO-Virgo will likely allow us to put constraints on the natal kicks of BHs (Vitale et al., 2017a, b; Stevenson et al., 2017; Zevin et al., 2017).

The behaviour of the BHB merger rate differs significantly from the cosmic SFR only if the CE parameters are drastically different () from our fiducial values (). According to the CE formalism (Webbink, 1984), low values of imply that is more difficult to eject the envelope. Thus, if is very low, the closest binaries merge prematurely during a CE phase, before they become BHBs. This explains why the merger rate of D0.02 is much lower than that of D at high redshift (). On the other hand, a very low value of implies that the shrinking of the orbit of the two cores within the CE is more efficient. Thus, looser binaries going through a CE phase might shrink enough to produce tight BHBs, which merge in a Hubble time. This might explain why the merger rate density of D0.02 increases at low redshift with respect to that of run D: in run D0.02, there is a higher number of BHBs with long delay time (Fig. 4), which form at high redshift but merge at low redshift.

In contrast, if is high, the CE is ejected easily and the orbit does not shrink efficiently. Thus, even if most binaries survive the CE phase, a lower number of BHBs become sufficiently close to merge in a Hubble time. These merging BHBs will have, on average, a longer delay time. This explains why less BHBs merge at high redshift () in run D1.5 with respect to run D, while at low redshift () the merger rate of D1.5 is slightly higher than that of D.

Among previous studies, Lamberts et al. (2016) and Belczynski et al. (2016b) contain several results that can be compared with ours quite straightforwardly. To produce their sample of BHBs, Belczynski et al. (2016b) use the startrack code (Belczynski et al., 2008), while Lamberts et al. (2016) use another updated version of the BSE code (Hurley et al., 2002) and focus only on the study of GW150914-like systems. The present-time BHB merger rates we obtain with are consistent (within a factor of 2) with those shown by Belczynski et al. (2016b), who adopt a similar choice for the CE parameters. In contrast, Lamberts et al. (2016) obtain a much larger present-time BHB merger rate ( Gpc yr) than expected from LIGO detections. As to the origin of this difference, we note that Lamberts et al. (2016) do not update the recipes for core radii in BSE, and (most importantly) do not include pulsational pair instability SNe, while this paper and Belczynski et al. (2016b) do.

If we restrict our analysis to systems like GW150914, GW151226, and GW170104, we find that their formation history also mimics the cosmic BHB merger rate density, but with a substantial difference. The distribution of formation times for GW150914-like and GW170104-like binaries peaks at high redshift () and then drops much faster than the total BHB merger rate density, while GW151226-like systems follow the general trend as the BHB merger rate. As a consequence, present-day GW150914-like and GW170104-like events are dominated by systems that formed at high redshift and have a long delay time. The reason is that GW150914-like and GW170104-like systems form mainly from metal-poor progenitors. We find a rather smooth distribution of the birth times for GW150914-like systems, consistently with Lamberts et al. (2016) and at odds with Belczynski et al. (2016b). This is not surprising, since both us and Lamberts et al. (2016) account for the mass-metallicity relation. Our main findings for GW150914, LVT151012 and GW151226 are also consistent with the results of Schneider et al. (2017), who predict that most GW150914 candidates form at high redshift (), whereas both GW151226 and LVT151012 have a broad range of possible formation redshifts.

We now discuss the main caveats of the present work. First, we have not included the impact of stellar dynamics on the evolution of BHBs. This might be a serious issue for our work, because massive stars are known to form especially in dense massive star clusters (see e.g. Weidner & Kroupa 2006), which are active dynamical places. Dynamical exchanges might lead to the formation of additional BHBs, which are likely more massive and more eccentric than average (e.g. Downing et al. 2010; Ziosi et al. 2014; Askar et al. 2017; Rodriguez et al. 2016; Banerjee 2017). Very massive BHs or even intermediate-mass BHs might form from runaway collisions (Portegies Zwart et al., 2004; Giersz et al., 2015; Mapelli, 2016). Furthermore, Kozai-Lidov resonances in triple systems might enhance the merger rate both in dense star clusters (e.g. Kimpson et al. 2016; Antonini & Rasio 2016) and in the field (e.g. Antonini et al. 2017). A study of the cosmic BHB merger rate including the effects of dynamics is still missing, because it poses a serious numerical challenge.

Moreover, we neglect the contribution of very massive stars ( M), because the current version of BSE does not include them. While these very massive objects are presumably very rare, they might significantly contribute to the merger rate of the most massive BHBs (e.g. Mapelli 2016).

Furthermore, our calculations neglect the chemically homogeneous evolutionary formation channel, recently proposed by Mandel & de Mink (2016), which might account for Gpc yr (see also Marchant et al. 2016 and de Mink & Mandel 2016).

In our analysis, we simply assume that population III stars behave as “normal” metal-poor stars (with ). While this choice is partially motivated by the fact that stellar winds are already highly inefficient at (see e.g. figure 6 of Spera et al. 2015), we do not know whether the mass function and the binary fraction of population III stars were significantly different from those of population II stars. On the other hand, we have checked that the contribution of population III stars (defined as Illustris’ star particles with metallicity ) to the cosmic BHB merger rate is negligible, especially within the LIGO horizon (see Fig. 11). We refer to other studies (e.g. Kinugawa et al. 2014, 2016; Belczynski et al. 2016a; Hartwig et al. 2016; Inayoshi et al. 2017) for a more accurate treatment of BHBs from population III stars.

Figure 11: Upper panel: cosmic merger rate density of BHBs () in run D, if we include (red solid line) or neglect (blue dashed line) population III stars (i.e. Illustris’ stars with metallicity ). Lower panel: residuals of the model with population III stars with respect to the model without them.

We also stress that we keep using the original fitting formulas of Hurley et al. (2000) for the photospheric radius and luminosity of stars, while our knowledge of the evolution of massive stars changed significantly in the last decades (e.g. Martins & Palacios 2013; Chieffi & Limongi 2013; Chen et al. 2015). In forthcoming studies, we will account for more recent stellar evolution prescriptions in the frame of the new SEVN population-synthesis code (Spera et al., in preparation).

Finally, the Illustris reproduces reasonably well many observational features (such as the cosmic SFR density and the galaxy luminosity function), but has several limitations which might affect the BHB merger rate. For example, it predicts a too mild decline in the cosmic SFR density at (Pillepich et al., 2017). This might lead to an overestimate of the present-day merger rate by several per cents (see Appendix B). Galaxies in the Illustris follow a mass-metallicity relation which is sensibly steeper than the observed relation (Torrey et al., 2014). As we discuss in Appendix B, this should affect the BHB merger rate density by per cent. Thus, it is important to repeat the same exercise that we did in this paper also with other cosmological simulations and to check for discrepancies3.

5 Summary

We reconstructed the cosmic BHB merger rate by planting BHBs into the Illustris cosmological box. The Illustris’ box (length=106.5 Mpc) is large enough to guarantee that we are modelling an unbiased portion of the Universe, satisfying the cosmological principle.

The population of BHBs is estimated through population synthesis simulations of isolated binaries, performed with an updated version of BSE. In particular, our new version of BSE includes up-to-date prescriptions for stellar winds of massive stars, during and after the MS. We account not only for the dependence of stellar winds on metallicity, but also for the effect of the electron-scattering Eddington limit (Chen et al., 2015). Up-to-date recipes for core collapse SNe, pair-instability and pulsational pair-instability SNe are also included.

We perform six different sets of runs with BSE, changing the SN prescription, the CE efficiency, the treatment of HG stars, and the distribution of natal kicks (Table 1). Each of these six simulations produces a population of merging BHBs, depending on the metallicity of the progenitor stars. We then interface the population of merging BHBs with the Illustris simulation through a Monte Carlo model, to obtain the cosmic BHB merger rate density for each of the six BSE simulation sets.

The cosmic BHB merger rate follows the same trend as the cosmic SFR, with a peak at , in all simulation sets (Fig. 1). In contrast, the normalization of the BHB merger rate strongly depends on the specific BSE simulation set. In particular, the treatment of CE and the distribution of SN kicks appear to be the most important processes.

Models in which a HG donor is allowed to survive the CE phase are not consistent with LIGO observations (unless very high natal kicks are assumed), because they give a too high BHB present-time merger rate density. The choice of the CE parameters (we study values of ranging from 0.02 to 1.5, in addition to our fiducial value ) significantly affects the merger rate at high redshift, while it induces only minor changes (a factor of 2) in the low-redshift merger rate ().

Population-synthesis simulations with large natal kicks (distributed according to Hobbs et al. 2005) result in a present-day BHB merger rate Gpc yr, while population-synthesis simulations with lower kicks (accounting for the amount of fallback) give a merger rate Gpc yr. Both rates are still consistent with the constraints from LIGO observations ( Gpc yr, Abbott et al. 2016a), but forthcoming detections might be able to distinguish between them.

From our simulations, we can also trace the metallicity of the progenitors of merging BHBs. We find that most BHB mergers detectable by current GW interferometers come from relatively metal-poor progenitors, ranging from Z to Z (Figs 2 and 7).

The merging BHBs have chirp masses (total masses) ranging from to M ( to M), with a slight dependence on the SN prescription and on CE parameters (Fig. 6).

If we focus on systems similar to GW150914, GW151226 and GW170104, we can obtain some useful hints on their progenitors. The formation of GW150914-like and GW170104-like systems is more efficient at high redshift and drops in the local Universe (Fig. 8). Most GW150914-like and GW170104-like systems merging in the local Universe appear to have formed at higher redshift with a long delay time. This happens because only genuinely metal-poor stars can produce GW150914-like and GW170104-like systems (with metallicity Z and Z, respectively; see Fig. 10). In contrast, GW151226-like systems form and merge all the way through the cosmic history (from to ). The progenitors of GW151226-like systems can be either metal-rich or metal-poor stars with about the same probability.

In our fiducial model (run D) the percentage of GW150914-like systems merging within the LIGO instrumental horizon is higher ( %) than the percentage of GW151226-like systems ( %, see Table 3). This might suggest that LIGO and Virgo will observe more events like GW150914 than like GW151226.

In conclusion, this study provides several clues about merging BHBs and their progenitors. We show that the BHB merger rate density poses constraints on both the CE process and the natal kicks of BHs. Forthcoming GW detections will allow us to further strengthen these constraints. Moreover, our results suggest that most BHBs merging within the LIGO horizon formed from relatively metal-poor progenitors ( Z). The Illustris simulation also includes information on the properties of BHB host galaxies. In a forthcoming paper, we will explore the properties of the galaxies where BHBs form and merge. The same analysis will be done also for NS binary systems and NS-BH binary systems, providing a clue for the detection of electromagnetic counterparts.

Acknowledgments

We thank the anonymous referee for their critical reading of the manuscript. We thank Alessandro Bressan, Marica Branchesi, Raffaella Schneider, Luca Graziani, Elena D’Onghia and Roberto Maiolino for useful discussions. We warmly thank The Illustris team for making their simulations publicly available. Numerical calculations have been performed through a CINECA-INFN agreement and through a CINECA-INAF agreement, providing access to resources on GALILEO and MARCONI at CINECA. MM and MS acknowledge financial support from the Italian Ministry of Education, University and Research (MIUR) through grant FIRB 2012 RBFR12PM1F, and from INAF through grant PRIN-2014-14. MM acknowledges financial support from the MERAC Foundation.

Appendix A Prescriptions for core-collapse SN and for fallback

We adopt two different models for core-collapse SNe: the rapid (R) and the delayed (D) model. Both models were introduced by Fryer et al. (2012) and they differ by the time-scale over which the explosion occurs: ms after the bounce for the rapid model, ms for the delayed mechanism.

Rapid SN model

For the rapid SN mechanism, a fixed mass of the proto-compact object, , is assumed. In this case, the fallback parameter is defined as

(6)

where is the final mass of the star, is the mass of the proto-NS, is the mass of the Carbon-Oxygen core, , and .

For all models, the amount of mass that falls onto the proto-NS is defined as . The final mass of the compact object is given by .

Delayed SN model

For the delayed SN mechanism, the mass of the proto-compact object is defined as

(7)

The amount of fallback is determined using the following relations

(8)

where and .

Appendix B The impact of resolution and mass-metallicity relation

Figure 12: Left axis: cosmic merger rate density of BHBs () in the comoving frame, as a function of the look-back time (bottom axis) and of the redshift (top axis). Red solid line: model D applied to the Illustris-1 simulation (same as in Fig. 1); violet dash-dot line: model D applied to the Illustris-3 simulation; blue dashed line: model D applied to the Illustris-3 simulation and adopting Maiolino et al. (2008; hereafter M08) fitting formulas for the mass-metallicity relation. Green shaded area: BHB merger rate inferred from LIGO detections (Abbott et al., 2016a). Right axis: cosmic SFR density from the Illustris-1 (grey thin solid line) and from Madau & Dickinson 2014 (grey thin dashed line), as a function of the look-back time (bottom axis) and of the redshift (top axis).
Name
[Gpc yr] [Gpc yr]
Il-1 125 181
Il-3 78 114
Il-3 and M08 96 135

Same as Table 1 but for the check runs. Column 1: model name; column 2: present-time BHB merger rate density; column 3: BHB merger rate density at .

Table 5: Comoving BHB merger-rate density at redshift and .
0.07 11.18 9.04
0.7 11.57 9.04
2.2 12.38 8.99
3.5 12.28 8.69

Column 1: redshift; column 2 and 3: values of the parameters in equation 9 at different redshift.

Table 6: Best fit parameters for the mass-metallicity relation in equation 9 at different redshift. The values of and at come from Mannucci et al. (2009), while the other values come from Maiolino et al. (2008).

In this section, we estimate the impact of the numerical resolution and of the simulated mass-metallicity relation on the main results of this paper. In Fig. 12 and Table 5, we compare the BHB merger rate density we obtain from the Illustris-1 simulation (model Il-1) with the BHB merger rate density we obtain from the Illustris-3 simulation (model Il-3), which has a factor of poorer resolution. The current BHB merger rate density in the Illustris-3 simulation is  % lower than in the Illustris-1 simulation, because of the lower resolution. Galaxies with stellar mass M are under-resolved (they consist of star particles) even in the Illustris-1 simulation. It is reasonable to expect that accounting for these low-mass galaxies might further boost the merger rate.

In a companion paper (Schneider et al., 2017), we follow a complementary approach by coupling our population-synthesis models with the GAMESH pipeline (Graziani et al., 2015; Graziani et al., 2017). GAMESH consists of a box and reaches a much higher resolution: even galaxies with a stellar mass of M are effectively resolved. Schneider et al. (2017) predict that all GW150914-like systems in the GAMESH box form in progenitor galaxies with a stellar mass M (see their Figure 2). Thus, the merger rate of GW150914-like systems is probably underestimated in the current study, due to the Illustris resolution.

Fig. 12 also shows the difference between the cosmic SFR density in the Illustris-1 simulation and the cosmic SFR density from Madau & Dickinson (2014). At redshift the Illustris-1 simulation predicts a higher SFR density by a factor of . This suggests that the current BHB merger rate density might be slightly overestimated in our calculations.

As we discussed in Section 2.2, the model of sub-grid physics adopted in the Illustris produces a mass-metallicity relation which significantly differs from the observed one. In particular, the simulated mass-metallicity relation is sensibly steeper than the observed one and does not show the observed turnover at high stellar mass (Vogelsberger et al., 2013; Torrey et al., 2014). We quantify the impact of these differences between simulated and observed mass-metallicity relation through the following procedure.

In the Illustris-3 simulation, we override the metallicity of a given star particle with the metallicity we expect from the observed mass-metallicity relation. For the observed relation, we adopt the fitting formula by Maiolino et al. (2008) and Mannucci et al. (2009):

(9)

where is the total stellar mass of the host galaxy in solar masses, while and are given in Table 6. For intermediate redshifts between those in Table 6, we obtain the metallicity by linear interpolation. At redshift () we simply use the same coefficients as for (). The new metallicity of each Illustris’ star is randomly extracted from a Gaussian distribution with mean value given by equation 9 (where is the total stellar mass of the sub-halo hosting the Illustris’ star) and standard deviation dex (accounting for metallicity dispersion within galaxies).

In Fig. 12 and Table 5 we show the BHB merger rate density we obtain with this procedure, compared to the BHB merger rate density we derive using the simulated metallicity of each Illustris-3 particle. The maximum difference between the two curves is only per cent. The same procedure can be applied to the Illustris-1 simulation, but is a factor of 60 more computationally expensive. We will show the results for the Illustris-1 in our follow-up paper.

Footnotes

  1. The fitting formulas by Hurley et al. (2000) might be inaccurate for very massive stars. To improve the treatment of massive stars, we impose that the values of the radius of single stars are consistent with PARSEC stellar evolution tracks (Chen et al., 2015) for stars with mass , as discussed in Mapelli (2016).
  2. If (), we associate to the Illustris’s particle a BSE set with (), since the maximum (minimum) metallicity we simulated with BSE is 0.02 (0.0002). This procedure is particularly arbitrary for population III stars, whose binarity properties are not known. However, we show in Section 4 that population III stars do not significantly affect the rate of detectable BHB mergers.
  3. An updated Illustris simulation will be available soon, with an improved algorithm for AGN feedback and galactic winds (Pillepich et al., 2017).

References

  1. Abbott B. P., et al., 2016a, Physical Review X, 6, 041015
  2. Abbott B. P., et al., 2016b, Phys. Rev. Lett., 116, 061102
  3. Abbott B. P., et al., 2016c, \apjl, 818, L22
  4. Abbott B. P., et al., 2017, Physical Review Letters, 118, 221101
  5. Antonini F., Rasio F. A., 2016, \apj, 831, 187
  6. Antonini F., Toonen S., Hamers A. S., 2017, \apj, 841, 77
  7. Askar A., Szkudlarek M., Gondek-Rosińska D., Giersz M., Bulik T., 2017, \mnras, 464, L36
  8. Banerjee S., 2017, \mnras, 467, 524
  9. Belczynski K., Sadowski A., Rasio F. A., 2004, \apj, 611, 1068
  10. Belczynski K., Kalogera V., Rasio F. A., Taam R. E., Zezas A., Bulik T., Maccarone T. J., Ivanova N., 2008, \apjs, 174, 223
  11. Belczynski K., Bulik T., Fryer C. L., Ruiter A., Valsecchi F., Vink J. S., Hurley J. R., 2010, \apj, 714, 1217
  12. Belczynski K., Ryu T., Perna R., Berti E., Tanaka T. L., Bulik T., 2016a, preprint, (arXiv:1612.01524)
  13. Belczynski K., Holz D. E., Bulik T., O’Shaughnessy R., 2016b, \nat, 534, 512
  14. Belczynski K., et al., 2016c, \aap, 594, A97
  15. Bird S., Cholis I., Muñoz J. B., Ali-Haïmoud Y., Kamionkowski M., Kovetz E. D., Raccanelli A., Riess A. G., 2016, Physical Review Letters, 116, 201301
  16. Bressan A., Marigo P., Girardi L., Salasnich B., Dal Cero C., Rubele S., Nanni A., 2012, \mnras, 427, 127
  17. Carr B., Kühnel F., Sandstad M., 2016, \prd, 94, 083504
  18. Chatterjee S., Rodriguez C. L., Rasio F. A., 2017, \apj, 834, 68
  19. Chen Y., Bressan A., Girardi L., Marigo P., Kong X., Lanza A., 2015, \mnras, 452, 1068
  20. Chieffi A., Limongi M., 2013, \apj, 764, 21
  21. Cole S., Helly J., Frenk C. S., Parkinson H., 2008, \mnras, 383, 546
  22. Colpi M., Mapelli M., Possenti A., 2003, \apj, 599, 1260
  23. Dominik M., Belczynski K., Fryer C., Holz D. E., Berti E., Bulik T., Mandel I., O’Shaughnessy R., 2012, \apj, 759, 52
  24. Dominik M., Belczynski K., Fryer C., Holz D. E., Berti E., Bulik T., Mandel I., O’Shaughnessy R., 2013, \apj, 779, 72
  25. Downing J. M. B., Benacquista M. J., Giersz M., Spurzem R., 2010, \mnras, 407, 1946
  26. Dvorkin I., Vangioni E., Silk J., Uzan J.-P., Olive K. A., 2016, \mnras, 461, 3877
  27. Elbert O. D., Bullock J. S., Kaplinghat M., 2017, preprint, (arXiv:1703.02551)
  28. Fontana A., et al., 2006, \aap, 459, 745
  29. Fryer C. L., Belczynski K., Wiktorowicz G., Dominik M., Kalogera V., Holz D. E., 2012, \apj, 749, 91
  30. Genel S., 2016, \apj, 822, 107
  31. Genel S., et al., 2014, \mnras, 445, 175
  32. Giersz M., Leigh N., Hypki A., Lützgendorf N., Askar A., 2015, \mnras, 454, 3150
  33. Gräfener G., Hamann W.-R., 2008, \aap, 482, 945
  34. Graziani L., Salvadori S., Schneider R., Kawata D., de Bennassuti M., Maselli A., 2015, \mnras, 449, 3137
  35. Graziani L., de Bennassuti M., Schneider R., Kawata D., Salvadori S., 2017, \mnras, 469, 1101
  36. Hall P. D., Tout C. A., 2014, \mnras, 444, 3209
  37. Hartwig T., Volonteri M., Bromm V., Klessen R. S., Barausse E., Magg M., Stacy A., 2016, \mnras, 460, L74
  38. Hinshaw G., et al., 2013, \apjs, 208, 19
  39. Hobbs G., Lorimer D. R., Lyne A. G., Kramer M., 2005, \mnras, 360, 974
  40. Hurley J. R., Pols O. R., Tout C. A., 2000, \mnras, 315, 543
  41. Hurley J. R., Tout C. A., Pols O. R., 2002, \mnras, 329, 897
  42. Inayoshi K., Hirai R., Kinugawa T., Hotokezaka K., 2017, \mnras, 468, 5020
  43. Inomata K., Kawasaki M., Mukaida K., Tada Y., Yanagida T. T., 2017, \prd, 95, 123510
  44. Ivanova N., Taam R. E., 2004, \apj, 601, 1058
  45. Kimpson T. O., Spera M., Mapelli M., Ziosi B. M., 2016, \mnras, 463, 2443
  46. Kinugawa T., Inayoshi K., Hotokezaka K., Nakauchi D., Nakamura T., 2014, \mnras, 442, 2963
  47. Kinugawa T., Miyamoto A., Kanda N., Nakamura T., 2016, \mnras, 456, 1093
  48. Kroupa P., 2001, \mnras, 322, 231
  49. Kulkarni S. R., Hut P., McMillan S., 1993, \nat, 364, 421
  50. Lamberts A., Garrison-Kimmel S., Clausen D. R., Hopkins P. F., 2016, \mnras, 463, L31
  51. Loveridge A. J., van der Sluys M. V., Kalogera V., 2011, \apj, 743, 49
  52. Ma X., Hopkins P. F., Faucher-Giguère C.-A., Zolman N., Muratov A. L., Kereš D., Quataert E., 2016, \mnras, 456, 2140
  53. Madau P., Dickinson M., 2014, \araa, 52, 415
  54. Maiolino R., et al., 2008, \aap, 488, 463
  55. Mandel I., de Mink S. E., 2016, \mnras, 458, 2634
  56. Mannucci F., et al., 2009, \mnras, 398, 1915
  57. Mapelli M., 2016, \mnras, 459, 3432
  58. Mapelli M., Colpi M., Zampieri L., 2009, \mnras, 395, L71
  59. Mapelli M., Ripamonti E., Zampieri L., Colpi M., Bressan A., 2010, \mnras, 408, 234
  60. Mapelli M., Zampieri L., Ripamonti E., Bressan A., 2013, \mnras, 429, 2298
  61. Marchant P., Langer N., Podsiadlowski P., Tauris T. M., Moriya T. J., 2016, \aap, 588, A50
  62. Martins F., Palacios A., 2013, \aap, 560, A16
  63. Nelson D., et al., 2015, Astronomy and Computing, 13, 12
  64. O’Shaughnessy R., Gerosa D., Wysocki D., 2017a, Physical Review Letters, 119, 011101
  65. O’Shaughnessy R., Bellovary J. M., Brooks A., Shen S., Governato F., Christensen C. R., 2017b, \mnras, 464, 2831
  66. Pei Y. C., Fall S. M., Hauser M. G., 1999, \apj, 522, 604
  67. Pillepich A., et al., 2017, preprint, (arXiv:1703.02970)
  68. Portegies Zwart S. F., McMillan S. L. W., 2000, \apjl, 528, L17
  69. Portegies Zwart S. F., Verbunt F., 1996, \aap, 309, 179
  70. Portegies Zwart S. F., Baumgardt H., Hut P., Makino J., McMillan S. L. W., 2004, \nat, 428, 724
  71. Repetto S., Davies M. B., Sigurdsson S., 2012, \mnras, 425, 2799
  72. Repetto S., Igoshev A. P., Nelemans G., 2017, \mnras, 467, 298
  73. Rodriguez C. L., Morscher M., Pattabiraman B., Chatterjee S., Haster C.-J., Rasio F. A., 2015, Physical Review Letters, 115, 051101
  74. Rodriguez C. L., Chatterjee S., Rasio F. A., 2016, \prd, 93, 084029
  75. Sana H., et al., 2012, Science, 337, 444
  76. Schneider R., Graziani L., Marassi S., Spera M., Mapelli M., Alparone M., de Bennassuti M., 2017, preprint, (arXiv:1705.06781)
  77. Schutz B. F., 1989, Classical and Quantum Gravity, 6, 1761
  78. Sigurdsson S., Hernquist L., 1993, \nat, 364, 423
  79. Spera M., Mapelli M., 2017, \mnras, 470, 4739
  80. Spera M., Mapelli M., Bressan A., 2015, \mnras, 451, 4086
  81. Springel V., 2010, \mnras, 401, 791
  82. Stevenson S., Berry C. P. L., Mandel I., 2017, preprint, (arXiv:1703.06873)
  83. Strolger L.-G., et al., 2004, \apj, 613, 200
  84. Tang J., Bressan A., Rosenfield P., Slemer A., Marigo P., Girardi L., Bianchi L., 2014, \mnras, 445, 4287
  85. Thorne K. S., 1987, Gravitational radiation.. pp 330–458
  86. Torrey P., Vogelsberger M., Genel S., Sijacki D., Springel V., Hernquist L., 2014, \mnras, 438, 1985
  87. Tutukov A., Yungelson L., 1973, Nauchnye Informatsii, 27, 70
  88. Vink J. S., 2016, preprint, (arXiv:1610.00578)
  89. Vink J. S., de Koter A., 2005, \aap, 442, 587
  90. Vink J. S., de Koter A., Lamers H. J. G. L. M., 2001, \aap, 369, 574
  91. Vink J. S., Muijres L. E., Anthonisse B., de Koter A., Gräfener G., Langer N., 2011, \aap, 531, A132
  92. Vitale S., Lynch R., Sturani R., Graff P., 2017a, Classical and Quantum Gravity, 34, 03LT01
  93. Vitale S., Lynch R., Raymond V., Sturani R., Veitch J., Graff P., 2017b, \prd, 95, 064053
  94. Vogelsberger M., Genel S., Sijacki D., Torrey P., Springel V., Hernquist L., 2013, \mnras, 436, 3031
  95. Vogelsberger M., et al., 2014a, \mnras, 444, 1518
  96. Vogelsberger M., et al., 2014b, \nat, 509, 177
  97. Webbink R. F., 1984, \apj, 277, 355
  98. Weidner C., Kroupa P., 2006, \mnras, 365, 1333
  99. Woosley S. E., 2017, \apj, 836, 244
  100. Xu X.-J., Li X.-D., 2010, \apj, 716, 114
  101. Zevin M., Pankow C., Rodriguez C. L., Sampson L., Chase E., Kalogera V., Rasio F. A., 2017, preprint, (arXiv:1704.07379)
  102. Ziosi B. M., Mapelli M., Branchesi M., Tormen G., 2014, \mnras, 441, 3703
  103. de Mink S. E., Mandel I., 2016, \mnras, 460, 3545
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
128002
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description