The properties of merging black holes and neutron stars across cosmic time

TL;DR

This study uses population-synthesis and cosmological simulations to analyze the properties of merging black holes and neutron stars across cosmic history, revealing mild dependence of mass on redshift and differences in black hole mass distributions over time.

Contribution

It provides a comprehensive characterization of merging compact objects across cosmic time, highlighting the influence of progenitor metallicity and delay times on mass distributions.

Findings

Black hole and neutron star merger properties are largely independent of redshift.

Black hole mass distribution is dominated by metal-poor progenitors across cosmic time.

Recent mergers show an excess of heavy black holes (20-35 M_sun) due to longer delay times.

Abstract

The next generation ground-based gravitational wave interferometers will possibly observe mergers of binary black holes (BBHs) and binary neutron stars (BNSs) to redshift $z\gtrsim{}10$ and $z\gtrsim{}2$, respectively. Here, we characterize the properties of merging BBHs, BNSs and neutron star-black hole binaries across cosmic time, by means of population-synthesis simulations combined with the Illustris cosmological simulation. We find that the mass of merging compact objects does not depend (or depends very mildly) on the merger redshift. Even the mass distribution of black holes depends only mildly on redshift, because BBHs originating from metal-poor progenitors ($Z\leq{}4\times{}10^{-3}$) dominate the entire population of merging BBHs across cosmic time. For a common-envelope efficiency $\alpha{}\ge{}3$, the main difference between the mass distribution of BBHs merging in the last Gyr and that of BBHs merging more than 11 Gyr ago is that there is an excess of heavy merging black holes ($20-35$ M$_\odot$) in the last Gyr. This excess is explained by the longer delay time of massive BBHs.

Figures

Tables

Table 1. Table 1: Properties of the population-synthesis simulations.

Run	$\alpha$	${\sigma }_{\mathrm{ECSN}}$	${\sigma }_{\mathrm{CCSN}}$
		[km s-1]	[km s-1]
$\alpha 5$	5.0	15	265
CC15$\alpha 5$	5.0	15	15

Table 2. Table 2: Percentage of systems that merge in the last Gyr and have t delay > 1 subscript 𝑡 delay 1 t_{\rm delay}>1 Gyr or > 10 absent 10 >10 Gyr.

Run	Binary type	$t_{\mathrm{delay}}>1$ Gyr	$t_{\mathrm{delay}}>10$ Gyr
$\alpha 5$	BBH	0.97	0.49
CC15$\alpha 5$	BBH	0.98	0.52
$\alpha 5$	NSBH	0.78	0.25
CC15$\alpha 5$	NSBH	0.87	0.28
$\alpha 5$	BNS	0.17	0.02
CC15$\alpha 5$	BNS	0.44	0.03

Table 3. Table 3: Properties of the supplementary population-synthesis simulations.

Run	$\alpha$	${\sigma }_{\mathrm{ECSN}}$	${\sigma }_{\mathrm{CCSN}}$	Fallback	CCSN
		[km s-1]	[km s-1]		model
$\alpha 1$	1.0	15	265	yes	rapid
$\alpha 3$	3.0	15	15	yes	rapid
K	1.0	265	265	no	delayed

Equations

$$n_{\mathrm{CO},\mathrm{i}}=N_{\mathrm{BSE},\mathrm{i}}\frac{M_{\mathrm{Ill}}}{M_{\mathrm{BSE}}}{f}_{\mathrm{corr}}{f}_{\mathrm{bin}},$$

$$t_{\mathrm{merge}}=t_{\mathrm{form}}-t_{\mathrm{delay}},$$

$$t_{\mathrm{form}}=\frac{1}{H_0}{\int }_0^{z_{\mathrm{Ill}}}\frac{1}{(1+z){\left[\Omega {}_M(1+z{)}^3+\Omega {}_{\Lambda }\right]}^{1/2}}\mathrm{d}z,$$

Full text

1 Introduction

The first two observing runs of the advanced LIGO (Aasi et al., 2015) and Virgo interferometers (Acernese et al., 2015) have led to the detection of ten binary black hole (BBH) mergers (Abbott et al., 2016b, c, a, 2017a, 2017d; Abbott, 2018a) and one binary neutron star (BNS) merger (Abbott et al., 2017b, c). From these detections, we can attempt to reconstruct the properties of BBHs merging in the local Universe: the black holes (BHs) detected thus far are consistent with a power law mass distribution with index $1.6^{+1.5}_{-1.7}$ (at 90% credibility) and with maximum mass $\sim{}45$ M⊙ (e.g. Abbott 2018b). The third observing run of LIGO-Virgo (O3), which will start in few months, is expected to bring tens of new detections, improving significantly our knowledge of the mass distribution of BBH mergers in the local Universe.

Third-generation ground-based gravitational wave (GW) detectors are now being planned (Punturo et al., 2010). The European Einstein Telescope and its American peer Cosmic Explorer will detect BBH mergers across the entire Universe and will probe BNS mergers up to redshift z\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}{}2, which is the peak of cosmic star formation rate (Madau & Dickinson, 2014). Thus, it is of crucial importance to investigate the redshift evolution of merging compact binaries, in preparation for O3 and especially for third-generation ground-based GW detectors. In particular, the predicted dependence of BH mass on the metallicity of the stellar progenitors (Heger et al., 2003; Mapelli et al., 2009; Mapelli et al., 2010, 2013; Belczynski et al., 2010a; Fryer et al., 2012; Spera et al., 2015) might suggest that more massive BHs form in the higher redshift Universe, where the average metallicity was lower.

Several studies have investigated the cosmic evolution of the merger rate density (e.g. O’Shaughnessy et al. 2010; Dominik et al. 2013, 2015; Belczynski et al. 2016; Dvorkin et al. 2016; de Mink & Mandel 2016; Mandel & de Mink 2016; Elbert et al. 2018; Mapelli et al. 2017; Mapelli & Giacobbo 2018; Eldridge et al. 2019; Fishbach et al. 2018; Rodriguez & Loeb 2018), suggesting that the merger rate of BBHs, neutron star-black hole binaries (NSBHs) and BNSs increases with redshift, reaching a peak at $z\sim{}2-4$. Other works have tried to characterize the host galaxies of merging compact objects (e.g. O’Shaughnessy et al. 2010; Lamberts et al. 2016; O’Shaughnessy et al. 2017; Schneider et al. 2017; Cao et al. 2018; Mapelli et al. 2018; Lamberts et al. 2018; Marassi et al. 2019; Artale et al. 2019). However, none of the previous studies focuses on the evolution of the population of merging compact objects across cosmic time. Is the typical mass of two BHs merging in the local Universe the same of two BHs merging at redshift $z\sim{}7$?

This paper studies the evolution of merging compact objects (BBHs, NSBHs and BNSs) across cosmic time. We investigate the redshift evolution of the mass spectrum of merging compact objects, of the metallicity of their progenitors and of the delay time (i.e. the time elapsed from the formation of the stellar progenitors and the merger of two compact objects). We use population-synthesis simulations (Giacobbo et al., 2018; Giacobbo & Mapelli, 2018) combined with the outputs of the Illustris cosmological simulation (Vogelsberger et al., 2014b, a) through a Monte Carlo procedure (Mapelli et al., 2017; Mapelli & Giacobbo, 2018). This methodology associates merging compact objects to a given galaxy (O’Shaughnessy et al., 2010; Kelley et al., 2010; Lamberts et al., 2016; O’Shaughnessy et al., 2017; Schneider et al., 2017) based on its star formation rate and on the metallicity of its stellar particles. Remarkably, we find that the mass distribution of merging compact objects depends only mildly on merger redshift.

2 Methods

We present the results of binary population-synthesis simulations run with the code mobse convolved with the outputs of the Illustris-1 cosmological simulation by means of a simple Monte Carlo algorithm. The methodology has already been described in Mapelli et al. (2017) and Mapelli & Giacobbo (2018). Thus, herebelow we briefly summarize the main steps of our methodology and we refer to the aforementioned papers for more details.

2.1 Population-synthesis simulations with mobse

mobse (Giacobbo et al., 2018; Giacobbo & Mapelli, 2018) is an upgraded version of the bse code (Hurley et al., 2000, 2002). The main updates concern mass loss by stellar winds of massive hot stars, the mass of compact remnants and the magnitude of natal kicks.

Mass loss of massive hot stars $\dot{M}$ is assumed to depend on both the metallicity $Z$ and the Eddington ratio $\Gamma_{e}=L_{\ast}{}/L_{\rm Edd}$ (where $L_{\ast}$ and $L_{\rm Edd}$ are the stellar luminosity and its Eddington value, respectively). In particular, we describe mass loss as $\dot{M}\propto{}Z^{\beta{}}$, where $\beta=0.85$ if $\Gamma{}_{e}<2/3$, $\beta=2.45-2.4\,{}\Gamma_{e}$ if $1>\Gamma{}_{e}\geq{}2/3$ and $\beta{}=0.05$ if $\Gamma_{e}\geq{}1$ (Chen et al. 2015, see also Vink et al. 2001; Vink & de Koter 2005; Gräfener & Hamann 2008; Vink et al. 2011).

The mass of compact objects which form via core-collapse supernova (SN) is derived following Fryer et al. (2012). In particular, we assume that the final mass of the compact object depends on the final mass of the Carbon-Oxygen core and on the final total mass of the star. If the final mass of the Carbon-Oxygen core is $m_{\rm CO}\geq{}11$ M⊙ the star is assumed to directly collapse into a black hole (BH) without SN explosion. In this paper, we adopt the rapid model presented in Fryer et al. (2012), which enforces the possible mass gap between the more massive neutron stars (NSs, $m_{\rm NS}\leq{}2$ M⊙) and the lighter BHs ($m_{\rm BH}\geq{}5$ M⊙, see Özel et al. 2010; Farr et al. 2011; Kreidberg et al. 2012; Littenberg et al. 2015). The outcomes of electron-capture SNe are also included in mobse, as described in Giacobbo & Mapelli (2019).

Finally, pair instability and pulsational pair instability SNe are implemented in mobse following Spera & Mapelli (2017), as described in Giacobbo et al. (2018). If the Helium core of a star grows larger than $\sim{}64$ M⊙, the star is completely destroyed by pair instability and leaves no compact object. In contrast, if the star has a Helium core mass $32\leq{}m_{\rm He}/m_{\odot}\leq{}64$, it undergoes pulsational pair instability, which leads to significant mass loss and to a smaller BH mass.

The stellar wind model combined with these prescriptions for SNe produces a mass spectrum of compact objects which depends on the metallicity of the progenitor star (see e.g. Figure 4 of Giacobbo et al. 2018). BHs with mass up to $\sim{}65$ M⊙ are allowed to form at low metallicity.

Natal kicks are another critical ingredient of binary population-synthesis codes. In mobse we draw natal kicks from a Maxwellian distribution with one-dimensional root-mean-square velocity $\sigma{}$. For electron-capture SNe, we assume $\sigma{}_{\rm ECSN}=15$ km s*-1* (see the discussion in Giacobbo & Mapelli 2019). For core-collapse SNe, we adopt $\sigma{}_{\rm CCSN}=265$ km s*-1* in run $\alpha{}5$ and $\sigma{}_{\rm CCSN}=15$ km s*-1* in run CC15$\alpha{}5$. The former value is inferred by Hobbs et al. (2005) from the analysis of the proper motion of $\sim{}230$ Galactic pulsars. The latter value was introduced by Giacobbo & Mapelli (2018) to account for the small kicks associated with ultra-stripped SNe (Tauris et al., 2015; Tauris et al., 2017) and for the possible evidence of small kicks in Galactic BNSs (Beniamini & Piran, 2016).

Moreover, we take into account the effect of fallback by reducing the kick velocity as $\tilde{v}_{\rm kick}=(1-f_{\rm fb})\,{}v_{\rm kick}$, where $v_{\rm kick}$ is the kick extracted from the Maxwellian distribution, while $f_{\rm fb}$ is the fraction of mass which falls back onto the proto-neutron star (see Fryer et al. 2012). If a BH forms via direct collapse $f_{\rm fb}=1$. Thus, BHs formed via direct collapse undergo no kick, apart from the Blaauw kick (Blaauw, 1961; Boersma, 1961) associated with symmetric neutrino loss111In the case of a direct collapse, we assume that the gravitational remnant mass is $m_{\rm rem,\,{}grav}=0.9\,{}m_{\rm f}$, where $m_{\rm f}$ is the stellar mass at the onset of collapse, (Fryer et al., 2012). Thus, even the Blaauw kick is generally negligible for direct-collapse BHs..

The main processes of binary evolution (wind mass transfer, Roche-lobe mass transfer, common envelope and tidal evolution) are implemented in mobse as described in bse (Hurley et al., 2002). In particular, our treatment of common envelope (CE) depends on two parameters: $\alpha{}$ (describing the efficiency of energy transfer) and $\lambda{}$ (describing the geometry of the envelope and the importance of recombinations). In the current paper, $\lambda{}$ is defined as in Claeys et al. (2014), to account for the contribution of recombinations, while $\alpha$ is a constant. We adopt $\alpha{}=5$. The main change in the description of CE with respect to bse consists in the treatment of Hertzsprung gap (HG) stars. In the standard version of bse, HG donors entering a CE phase are allowed to survive the CE phase. In the version of mobse we use here, HG donors are forced to merge with their companion if they enter a CE. In fact, mobse simulations in which HG donors are allowed to survive a CE phase produce a local BBH merger rate $R_{\rm BBH}\sim{}600-800$ Gpc*-3* yr*-1*, which is not consistent with LIGO-Virgo results (Mapelli et al., 2017).

Decay by gravitational-wave emission is described as in Peters (1964). In contrast to bse, we account for gravitational-wave decay for all compact-object binaries, even with semi-major axis $a>10$ R⊙.

In this paper, we consider two different runs performed with mobse and already presented in Giacobbo & Mapelli (2018): run $\alpha{}5$ and run CC15$\alpha{}5$ (see Table 1). In both runs, we assume CE efficiency $\alpha{}=5$. We have chosen $\alpha{}=5$ as a fiducial value, because Figure 15 of Giacobbo & Mapelli (2018) shows that the BNS merger rate inferred from GW170817 (Abbott, 2018a) is difficult to match if we assume a lower value of $\alpha{}$.

Runs $\alpha{}5$ and CC15$\alpha 5$ differ only by the magnitude of the natal kick of core-collapse SNe: $\sigma{}_{\rm CCSN}=265$ and $15$ km s*-1* in run $\alpha{}5$ and run CC$15\alpha{}5$, respectively. In Giacobbo & Mapelli (2018), we have shown that the local merger rate of BNSs derived from run CC15$\alpha{}5$ is consistent with the one inferred from GW170817, while the local BNS merger rate derived from run $\alpha{}5$ is only marginally consistent with GW170817. The local merger rate of NSBHs and BBHs is consistent with the numbers derived by the LIGO-Virgo collaborations in both run $\alpha{}5$ and CC15$\alpha{}5$ (see Figure 15 of Giacobbo & Mapelli 2018). In the Appendix we discuss three supplementary runs where we consider different values of $\alpha{}$ ($=1,$ 3) and different assumptions for the natal kick.

For each run, we have simulated 12 sub-sets of binaries with metallicity $Z=0.0002,$ 0.0004, 0.0008, 0.0012, 0.0016, 0.002, 0.004, 0.006, 0.008, 0.012, 0.016, and 0.02. In each sub-set we have simulated $10^{7}$ stellar binaries. Thus, each of the two runs is composed of $1.2\times{}10^{8}$ massive binaries.

For each binary, the mass of the primary is randomly drawn from a Kroupa initial mass function (Kroupa, 2001) ranging from 5 to 150 M⊙, and the mass of the secondary is sampled according to the distribution $\mathcal{F}(q)\propto{}q^{-0.1}$ (where $q$ is the ratio between mass of the secondary and mass of the primary) in a range $[0.1\,{}-1]\,{}m_{\rm p}$. The orbital period $P$ and the eccentricity $e$ are randomly extracted from the distribution $\mathcal{F}(P)\propto{}(\log_{10}{P})^{-0.55}$, with $0.15\leq{}\log_{10}{(P/{\rm day})}\leq{}5.5$, and $\mathcal{F}(e)\propto{}e^{-0.4}$, with $0\leq{}e<1$ (Sana et al., 2012).

2.2 Coupling with the Illustris cosmological simulation

The Illustris-1 is the highest resolution hydrodynamical simulation run in the frame of the Illustris project (Vogelsberger et al., 2014b, a; Nelson et al., 2015). In the following, we refer to it simply as the Illustris. It covers a comoving volume of $(106.5\,{}{\rm Mpc})^{3}$, and has an initial dark matter and baryonic matter mass resolution of $6.26\times{}10^{6}$ and $1.26\times{}10^{6}$ M⊙, respectively (Vogelsberger et al., 2014b, a). 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).

As for the cosmology, the Illustris adopts WMAP-9 results for the cosmological parameters (Hinshaw et al., 2013), that is $\Omega{}_{M}=0.2726$, $\Omega{}_{\Lambda}=0.7274$, $\Omega{}_{b}=0.0456$, and $H_{0}=100\,{}h$ km s*-1* Mpc*-1*, with $h=0.704$.

We combine the catalogues of merging BNSs, NSBHs and BBHs with the Illustris simulations through a Monte Carlo procedure, as first described in Mapelli et al. (2017). In particular, we extract all new-born star particles from the Illustris and we associate to each of these star particles a number $n_{\rm CO,\,{}i}$ of merging compact-object binaries (where the index $i=$ BBH, NSBH or BNS) based on the initial mass $M_{\rm Ill}$ and on the metallicity $Z_{\rm Ill}$ of each Illustris star particle. In particular,

[TABLE]

where $M_{\rm BSE}$ is the total initial stellar mass of the population simulated with mobse whose metallicity is closer to $Z_{\rm Ill}$, while $N_{\rm BSE,\,{}i}$ is the number of merging compact-object binaries which form in the stellar population of mass $M_{\rm BSE}$. The term $f_{\rm corr}=0.285$ is a correction factor, accounting for the fact that we actually simulate only primaries with zero-age main sequence mass $m_{\rm ZAMS}\geq{}5$ M⊙, neglecting lower mass stars. The term $f_{\rm bin}$ accounts for the fact that we simulate only binary systems. Here we assume that 50 % of stars are in binaries, thus $f_{\rm bin}=0.5$. We note that this is a simplifying assumption, because $f_{\rm bin}$ should also depend on stellar mass (O’Shaughnessy et al., 2010). This dependence might be important when calculating the merger rate density, while it does not affect significantly the properties of merging binaries.

Finally, we randomly draw $n_{\rm CO,\,{}i}$ merging compact objects from the simulated population with mass $M_{\rm BSE}$ and we calculate the look-back time of the merger of each compact object binary as

[TABLE]

where $t_{\rm delay}$ is the delay time (i.e. the time elapsed between the formation of the stellar binary and the merger of the two compact objects) and $t_{\rm form}$ is the look back time at which the Illustris’ particle has formed, calculated as

[TABLE]

where the cosmological parameters are set to WMAP-9 values (for consistency with the Illustris) and $z_{\rm Ill}$ is the formation redshift of the Illustris’ particle. We have a minus sign in equation 2 because both $t_{\rm merge}$ and $t_{\rm form}$ are look-back times.

3 Results

3.1 Binary black holes (BBHs)

Figure 1 shows the mass distribution of BBHs merging in the local Universe ($\leq{}1$ Gyr ago) and at higher redshift. In particular, we consider BBHs merging $\leq{}1$, $3-4$, $6-7$, $9-10$, $11-12$ and $12-13$ Gyr ago (corresponding to redshift $z\leq{}0.08$, $0.26-0.37$, $0.64-0.82$, $1.35-1.78$, $2.43-3.65$ and $3.65-7.15$, respectively). From this Figure, it is apparent that the mass of BHs in merging BBHs evolves only mildly with redshift.

This might seem a bit surprising but comes from two facts. First, metal-rich binaries are much less efficient than metal-poor ones in producing merging BBHs (Belczynski et al., 2010b; Giacobbo et al., 2018; Klencki et al., 2018). In particular, the number of mergers per unit stellar mass is approximately $3-4$ orders of magnitude higher at low metallicity than at solar metallicity, according to the population-synthesis runs adopted this paper (see e.g. Figure 14 of Giacobbo & Mapelli 2018). Thus, the contribution of metal-rich stars to the merger rate is so small that most merging BBHs come from metal-poor stars even in the local Universe. Second, metallicity in the Universe is quite patchy: even at low redshift we can find a number of galaxies with very low metallicity, while even at high redshift the most massive galaxies have already reached a significantly high metallicity (see e.g. Maiolino & Mannucci 2019 for a recent review).

The main difference between BBHs merging in the local Universe and BBHs merging $>12$ Gyr ago concerns BBHs with mass $\sim{}20-35$ M⊙. BHs with mass $>20$ M⊙ are more common in BBHs that merge $\leq{}1$ Gyr ago than in BBHs merging $\geq{}12$ Gyr ago. This might be surprising because in our model more massive BHs are produced by metal-poor stars. On the other hand, this can be explained with a different delay time: the most massive BBHs tend to have longer delay time than light BBHs. This is apparent from Figure 2, where we show the delay time distribution of all simulated merging BBHs (regardless of their merger redshift) distinguishing between massive BHs ($m_{\rm BH}\geq{}20$ M⊙) and light BHs ($m_{\rm BH}<20$ M⊙): there is a clear dearth of massive BHs with t_{\rm delay}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}{}1 Gyr. Thus, even if they form preferentially in the early Universe, massive BHs tend to merge locally with a long delay time.

It is important to note that the distributions of delay times shown in Figure 2 refer to all BBHs merging within the cosmological simulation, regardless of when their progenitor stars formed: this is conceptually different from the delay time distribution of a coeval stellar population. The latter scales approximately as $t^{-1}$ (O’Shaughnessy et al., 2010; Dominik et al., 2012; Belczynski et al., 2016), with a moderate dependence on the assumed population-synthesis parameters (see e.g. Figure 6 of Giacobbo et al. 2018).

Finally, Figure 1 does not show any significant difference between the two fiducial population-synthesis models we have considered ($\alpha{}5$ and CC15$\alpha{}5$). In the Appendix, we discuss the impact on BBH masses of the CE parameter $\alpha{}$ and of the natal kick assumption. Our supplementary runs confirm that BBH mass changes very mildly with redshift. In particular, we find that if $\alpha{}=1$ even the slight dependence of BBH mass on redshift disappears.

The distribution of delay times of BBHs grouped by different merger time (Figure 3) confirms our analysis on BBH masses: most BBHs merging in the last Gyr form in the early Universe and have an extremely long delay time. Table 2 shows that $>97$% of all BBHs merging in the last Gyr have a delay time $>1$ Gyr, and $\sim{}50$% have a delay time $>10$ Gyr (i.e. they formed $>10$ Gyr ago). This result is in good agreement with previous work (Dominik et al., 2013, 2015; Mapelli et al., 2017).

The contribution of metal-rich binaries with short delay time to the local BBH merger rate is very minor (see also Mapelli et al. 2017; Mapelli & Giacobbo 2018; Mapelli et al. 2018). Thus, the distribution of BH masses we can reconstruct from GW detections of BBHs merging in the local Universe reflects the distribution of BH masses formed several Gyr ago, rather than the local distribution of BH masses. This might be the key to explain the difference between the mass distribution of BHs in local X-ray binaries (Özel et al., 2010; Farr et al., 2011) and the mass distribution of BHs in GW events.

If we look at the metallicity of stellar progenitors (Figure 4), we find a significant shift between the peak metallicity of the progenitors of BBHs which merge $12-13$ Gyr ago ($Z_{\rm peak}\sim{}2-3\times{}10^{-4}$) and that of BBHs merging $\leq{}1$ Gyr ago ($Z_{\rm peak}\sim{}10^{-3}$). At all redshifts, merging BBHs with metallicity \lower 2.15277pt\hbox{;\buildrel>\over{\sim};}{}4-6\times{}10^{-3} are extremely rare (see also Mapelli et al. 2017; Mapelli & Giacobbo 2018).

Stars with metallicity $Z\ll{}10^{-5}$ (i.e. population III stars) give a small contribution even to BBHs merging at high redshift, consistent with previous work (Hartwig et al., 2016; Belczynski et al., 2017). However, we warn that neither our population synthesis models include a specific treatment for population III stars nor the Illustris simulation contains any specific sub-grid model to describe them, and thus we cannot specifically quantify their impact in the current study.

3.2 Neutron star - black hole binaries (NSBHs)

Even the mass of BHs in merging NSBHs does not change significantly with time (see the left-hand panel of Figure 5). Again, this originates from the fact that the merger efficiency of NSBHs born from metal-poor progenitors is orders of magnitude higher than the merger efficiency of NSBHs born from metal-rich progenitors (see Giacobbo & Mapelli 2018).

Light BHs are more frequent than massive BHs in NSBHs: most BHs in merging NSBHs have mass $\leq{}10$ M⊙, regardless of the merger time. The dearth of massive BHs is particularly strong in the case of run CC$15\alpha{}5$.

From the right-hand panel of Figure 5, we see that massive neutron stars (NSs, $m_{\rm NS}>1.4$ M⊙) are more common than low-mass NSs in NSBHs. This holds in both runs, although the dearth of light NSs is even stronger in run CC15$\alpha{}5$.

Thus, merging NSBHs tend to have the maximum possible mass ratio between the NS and BH ($m_{\rm NS}/m_{\rm BH}\sim{}0.1-0.4$ in our models which enforce the mass gap between NSs and BHs). This can be explained with the fact that non-conservative mass transfer and CE in close binaries tend to lead to compact objects of similar mass.

The distribution of NSBH delay times (Figure 6) is much steeper than that of BBHs: there is a significantly larger number of merging NSBHs with short delay times with respect to BBHs (Dominik et al., 2013, 2015; Mapelli et al., 2018). Table 2 shows that $\sim{}25-28$% of NSBHs merging in the last Gyr formed $\geq{}10$ Gyr ago, but $\sim{}13-22$% formed $<1$ Gyr ago: even NSBHs which formed in the last Gyr give an important contribution to the population of NSBHs merging in the local Universe.

The metallicity distribution of progenitors of NSBHs follows the same trend as that of BBHs: the progenitors of NSBHs merging $\geq{}12$ Gyr ago tend to be more metal-poor than those of NSBHs merging in the last Gyr. However, the offset between the two populations is less evident than in the case of BBHs (the peak metallicity being $Z_{\rm peak}\sim{}10^{-3}$ and $Z_{\rm peak}\sim{}1-3\times{}10^{-3}$ for NSBHs merging $12-13$ Gyr ago and NSBHs merging in the last Gyr, respectively). The progenitors of merging NSBHs have $Z\leq{}1-2\times{}10^{-2}$. There is no significant difference between run $\alpha{}5$ and CC15$\alpha{}5$.

3.3 Binary neutron stars (BNSs)

Figure 8 shows that even the mass spectrum of BNSs merging $12-13$ Gyr ago does not differ dramatically with respect to the mass spectrum of BNSs merging in the last Gyr. However, BNSs merging $12-13$ Gyr ago have a slight preference for larger masses than BNSs merging $<1$ Gyr ago. This is an effect of metallicity, as we already discussed for the same population-synthesis runs in Giacobbo & Mapelli (2018): metal-poor binaries tend to produce slightly more massive NSs than metal-rich ones.

Differently from what happens to BHs (for which the BH mass dependence on metallicity does not translate into a BH mass dependence on merger redshift, because BBH mergers are orders of magnitude more common in a metal-poor population than in a metal-rich population), in the case of BNSs the NS mass dependence with progenitor’s metallicity translates also into a NS mass dependence with merger redshift, because the number of BNS mergers per unit stellar mass is approximately the same at high and low metallicity.

In contrast to NSBHs, NSs in merging BNSs tend to be light222As we already discussed in Mapelli & Giacobbo (2018), the NS mass distribution we derived adopting the prescriptions in Fryer et al. (2012) tends to underestimate the typical mass of NSs by $\sim{}0.1$ M⊙ with respect to observed Galactic NSs. ($m_{\rm NS}<1.3$ M⊙).

Merging BNSs tend to have much shorter delay time than both BBHs and NSBHs (Figure 9). This difference comes mostly from the metallicity dependence of the merger rate density of NSBHs and BBHs (Dominik et al., 2013, 2015; Mapelli et al., 2018).

Even if the delay time of BNSs tends to be shorter than that of BBHs and NSBHs, we stress that $\sim{}17$% and $\sim{}44$% of all BNSs merging in the last Gyr have a delay time longer than 1 Gyr in runs $\alpha{}5$ and CC15$\alpha{}5$, respectively (Table 2, note that delay times of BNSs are significantly longer in run CC15$\alpha{}5$). This is a fundamental clue to understand why GW170817 is associated with an early-type galaxy: even if star formation is low in NGC 4993 nowadays, $\sim{}17-44$% of all NSs which merge in the local Universe have formed $>1$ Gyr ago and are now locked in non-star forming massive galaxies.

The metallicity of BNS progenitors (Figure 10) peaks at solar or super-solar metallicity: $\sim{}2-4\times{}10^{-2}$ for BNSs merging $0-1$ Gyr ago. Figure 10 basically reflects what is the most common metallicity in the Illustris simulation, because the number of BNS mergers per unit stellar mass does not significantly depend on metallicity. The only significant difference between $\alpha{}5$ and CC$15\alpha{}5$ is the secondary peak at $Z\sim{}10^{-4}$ in the latter simulation. This is an effect of the dearth of BNS mergers from stars with metallicity $Z\sim{}4\times{}10^{-4}-4\times{}10^{-3}$ in the population synthesis simulations (see Figure 14 of Giacobbo & Mapelli 2018).

4 Discussion

The main prediction of this paper is that the mass of merging compact objects does not depend (or depends mildly) on their merger redshift.

This result depends on a number of assumptions. First, it depends on the strong dependence of BBH and NSBH merger efficiency on metallicity predicted by population-synthesis simulations. In particular, mobse predicts that the number of BBH and NSBH mergers originating from a metal-poor ($Z\leq{}2\times{}10^{-3}$) population is $3-4$ orders of magnitude larger than the number of BBH and NSBH mergers we expect from a metal-rich population with the same initial total mass. Thus, the number of merging BBHs and NSBHs which form from metal-poor stars is tremendously larger than the number of merging BBHs and NSBHs which form from metal-rich progenitors, even in the local Universe. Most BBHs/NSBHs that merge in the local Universe are expected to come from metal-poor progenitors which formed in the early Universe and have a long delay time. This means that the population of merging BBHs and NSBHs does not evolve with redshift because we are always looking at merging BBHs and NSBHs which formed from metal-poor stars, no matter what merger redshift we are considering.

Other population-synthesis codes (e.g. Belczynski et al. 2010b; Dominik et al. 2013, 2015; Stevenson et al. 2017; Chakrabarti et al. 2017; Kruckow et al. 2018; Spera et al. 2019) have predicted a similar trend (i.e. BBH mergers are more common from metal-poor progenitors than from metal-rich progenitors), but possibly with different strength. For example, in the simulations performed with sevn (Spera et al., 2019), the number of mergers per unit stellar mass is “only” two orders of magnitude lower at solar metallicity than at low metallicity (with mobse we find a difference of at least three orders of magnitude).

A further crucial ingredient is metallicity evolution across cosmic time, which is very uncertain (see e.g. Madau & Dickinson 2014; Maiolino & Mannucci 2019). 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 (\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}{}10^{10} M⊙). In Mapelli et al. (2017), we estimate that the impact on BBH merger rate of the difference between the Illustris mass-metallicity relation and the observational relation (Maiolino et al., 2008; Mannucci et al., 2009) is of the order of 20 %. However, even if metallicity does not affect significantly the merger rate, it might affect the properties of merging compact objects (e.g. the masses).

Thus, in a follow-up paper we will investigate the impact of different population-synthesis prescriptions and different cosmological evolution models, to properly assess the impact of all of these ingredients on the redshift evolution of compact binaries (Santoliquido et al. in preparation).

Moreover, in this work we neglect the dynamical formation channel of merging compact objects. Dynamics of young star clusters (Ziosi et al., 2014; Mapelli, 2016; Di Carlo et al., 2019) and globular clusters (Portegies Zwart & McMillan, 2000; Downing et al., 2010; Rodriguez et al., 2016; Askar et al., 2017; Hong et al., 2018) is known to favour the merger of more massive BBHs and to speed up their delay time (Di Carlo et al., 2019). Thus, dynamics might significantly affect the mass and delay time of BBHs with respect to isolated binary evolution.

Finally, we have not considered primordial BHs, i.e. BHs formed from gravitational instabilities in the early Universe (Carr & Hawking, 1974). Their very existence and their mass range are highly uncertain (Carr et al., 2016). Any deviation from our results which cannot be explained with uncertainties on population synthesis models, dynamical effects, star formation or metallicity evolution in the Universe might suggest a different formation channel for BHs than the stellar one.

5 Summary

We have investigated the redshift evolution of several properties of merging compact objects (BBHs, NSBHs and BNSs) by means of population synthesis simulations. We used our population-synthesis code mobse, which adopts updated stellar wind models and prescriptions for electron-capture, core-collapse and (pulsational) pair instability SNe (Giacobbo et al., 2018; Giacobbo & Mapelli, 2018, 2019). We have considered two very different prescriptions for natal kicks of core-collapse SNe: compact objects receive a natal kick distributed according to a Maxwellian distribution with one-dimensional root mean square velocity $\sigma_{\rm CCSN}=265$ km s*-1* and $\sigma_{\rm CCSN}=15$ km s*-1* in run $\alpha{}5$ and  CC15$\alpha{}5$, respectively (in both cases the kick of BHs is modulated by fallback). We combined the results of mobse with the outputs of the Illustris cosmological simulations (Vogelsberger et al., 2014b, a; Nelson et al., 2015) by means of a Monte Carlo formalism (Mapelli et al., 2017; Mapelli & Giacobbo, 2018; Mapelli et al., 2018). With this procedure, we can account for redshift evolution of star formation rate and metallicity.

We find that the mass distribution of merging compact objects (even BBHs) depends only mildly on merger redshift (Figures 1, 5, 8). This happens because the merger rate of BBHs and NSBHs depends dramatically on metallicity (Belczynski et al., 2010b; Giacobbo & Mapelli, 2018): the entire population of merging BBHs and NSBHs across cosmic time is dominated by metal-poor progenitors. Even if metal-rich stars should be more common in the local Universe and BBHs formed from metal-rich progenitors tend to be less massive than BBHs formed from metal-poor progenitors, we do not see the signature of BBHs originating from metal-rich stars because they are outnumbered by BBHs originating from metal-poor stars, even in the local Universe.

The only significant difference between BBHs merging in the last Gyr and BBHs merging more than 11 Gyr ago is that there is an excess of merging BBHs with mass $>20$ M⊙ in the former population with respect to the latter (Figure 1). This happens because massive BHs ($m_{\rm BH}\sim{}20-35$ M⊙) have preferentially long delay times. Thus, even if they form preferentially in the early Universe, they merge mostly in the local Universe. It is important to note that this difference appears only if the CE efficiency parameter is large ($\alpha\geq{}3$). For smaller values of $\alpha{}$ even the BBH mass distribution does not change appreciably with redshift (see Figure 11).

The mass of BHs in BBHs spans from $\sim{}5$ to $\sim{}45$ M⊙, while the mass of BHs in NSBHs is preferentially low: most BHs in NSBHs have mass $<10$ M⊙, especially in run CC$15\alpha{}5$. The mass of NSs in NSBHs is preferentially large ($>1.4$ M⊙), while the mass of NSs in BNSs is preferentially small ($<1.3$ M⊙).

The delay time distribution of BBHs (Figure 3) is significantly flatter than that of NSBHs (Figure 6) and especially BNSs (Figure 9), in agreement with previous work (Dominik et al., 2013, 2015). This is a consequence of the maximum stellar radius of the progenitors, which is significantly larger in BBHs progenitors, suppressing the formation of very close BBH binaries (even after CE). Thus, this result might depend on the stellar evolution models adopted in mobse.

The typical progenitor’s metallicity of merging BBHs and NSBHs is well below solar and evolves mildly with redshift (Figures 4 and 7). In contrast, the typical progenitor’s metallicity of BNSs is solar or super-solar (Figure 10). Since the number of BNS mergers per unit solar mass does not depend on metallicity significantly (Dominik et al., 2015; Giacobbo & Mapelli, 2018), the metallicity evolution of BNS progenitors traces the average metallicity evolution of the Universe, at least in the cosmological simulation.

The main prediction of this paper is that the mass of merging compact objects does not depend (or depends very mildly) on their merger redshift. It is worth noting that we find no dramatic differences between run $\alpha{}5$ and CC15$\alpha{}5$, despite the very different assumption for natal kicks. In this work, we have neglected the effect of dynamics. Dynamics tends to favour the merger of the most massive BHs (e.g. Mapelli 2016; Di Carlo et al. 2019) and might significantly affect our conclusions. Moreover, the evolution of metallicity in the Universe is highly uncertain (Maiolino & Mannucci, 2019). This sums up to uncertainties on binary evolution (Spera et al., 2019). Thus, it is tremendously important to perform follow-up studies with a different treatment of metallicity and binary evolution.

Acknowledgments

We thank the anonymous referee for their critical reading of the manuscript. We thank Daniel Wysocki, Alessandro Ballone, Alessandro Bressan, Emanuele Ripamonti and Mario Spera 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 acknowledges financial support by the European Research Council for the ERC Consolidator grant DEMOBLACK, under contract no. 770017. MCA acknowledges financial support from the Austrian National Science Foundation through FWF stand-alone grant P31154-N27 ‘Unraveling merging neutron stars and black hole - neutron star binaries with population-synthesis simulations’. NG acknowledges financial support from Fondazione Ing. Aldo Gini and thanks the Institute for Astrophysics and Particle Physics of the University of Innsbruck for hosting him during the preparation of this paper. This work benefited from support by the International Space Science Institute (ISSI), Bern, Switzerland, through its International Team programme ref. no. 393 The Evolution of Rich Stellar Populations & BH Binaries (2017-18).

Appendix A Dependence of BBH mass on CE efficiency and natal kicks

The two simulations presented in the main text (fiducial simulations) assume CE parameter $\alpha{}=5$. We have chosen $\alpha{}=5$ as a fiducial value, because Figure 15 of Giacobbo & Mapelli (2018) shows that the BNS merger rate inferred from GW170817 (Abbott, 2018a) is difficult to match if we assume a lower value of $\alpha{}$. The choice of a different $\alpha{}$ does not affect significantly the masses of BNSs and NSBHs (see e.g. Figures 3 and 8 of Giacobbo & Mapelli 2018), but might be important for BBHs (Figure 12 of Giacobbo & Mapelli 2018). Moreover, the natal kicks of BBHs are quite similar in run $\alpha{}5$ and  CC15$\alpha{}5$, because these runs assume the same dependence on fallback.

Thus, in this Appendix, we consider three additional models with different values of $\alpha{}$ and a different kick prescription, to check their impact on BBH masses. These additional runs are described in Table 3. In particular, runs $\alpha{}1$ and $\alpha{}3$ are the same as presented in Giacobbo & Mapelli (2018): they differ from run $\alpha{}5$ only because $\alpha{}=1,\,{}3$, respectively. Run K was already discussed in Mapelli et al. (2017). In this run, all BHs receive a natal kick randomly drawn from a Maxwellian distribution with $\sigma{}_{\rm CCSN}=265$ km s*-1*. The kick is not reduced by the amount of fallback. Moreover, this run adopts the delayed SN model (Fryer et al., 2012).

Figure 11 shows the distribution of primary BH masses in these three check runs. BBH masses vary very mildly with redshift even in these three runs. If $\alpha{}\leq{}1$, even the slight excess of BHs with mass $\sim{}20-35$ M⊙ merging in the last Gyr disappears: the mass of BHs does not evolve significantly with redshift.