The Cosmic Merger Rate Density Evolution of Compact Binaries Formed in Young Star Clusters and in Isolated Binaries

We derive the cosmic MRD of CBs as

$$\begin{eqnarray}\begin{array}{rcl}{ \mathcal R }(z) & = & \displaystyle \frac{{\text{}}d}{{\text{}}{{dt}}_{\mathrm{lb}}(z)}{\displaystyle \int }_{{z}_{\max }}^{z}\psi (z^{\prime} )\displaystyle \frac{{\text{}}{{dt}}_{{lb}}(z^{\prime} )}{{\text{}}{dz}^{\prime} }\,{\text{}}{dz}^{\prime} \\ & & \times {\displaystyle \int }_{{Z}_{\min }(z^{\prime} )}^{{Z}_{\max }(z^{\prime} )}\eta (Z){ \mathcal F }(z^{\prime} ,z,Z){\text{}}{dZ},\end{array}\end{eqnarray} \tag{ 1 }$$

where ${t}_{\mathrm{lb}}(z)$ is the look-back time at redshift z, $\psi (z^{\prime} )$ is the cosmic SFR density at redshift $z^{\prime}$, ${Z}_{\min }(z^{\prime} )$ and ${Z}_{\max }(z^{\prime} )$ are the minimum and maximum metallicities of stars formed at redshift $z^{\prime}$, $\eta (Z)$ is the merger efficiency at metallicity Z, and ${ \mathcal F }(z^{\prime} ,z,Z)$ is the fraction of CBs that form at redshift $z^{\prime}$ from stars with metallicity Z and merge at redshift z, normalized to all CBs that form from stars with metallicity Z. To calculate the look-back time we take the cosmological parameters (H₀, ${{\rm{\Omega }}}_{{\rm{M}}}$, and ${{\rm{\Omega }}}_{{\rm{\Lambda }}}$) from Ade et al. (2016). The maximum considered redshift in Equation (1) is ${z}_{\max }=15$, which we assume to be the epoch of formation of the first stars.

The cosmic SFR density $\psi (z)$ is given by the following fitting formula (Madau & Fragos 2017)

$$\begin{eqnarray}&&\psi (z)=0.01\,\displaystyle \frac{{\left(1+z\right)}^{2.6}}{1+{[(1+z)/3.2]}^{6.2}}\,{M}_{\odot }\,{\mathrm{Mpc}}^{-3}\,{\mathrm{yr}}^{-1}.\end{eqnarray} \tag{ 2 }$$

To estimate the uncertainty on $\psi (0)$, we assume that the errors follow a log-normal distribution with mean $\mathrm{log}\psi (0)=-2$ and standard deviation ${\sigma }_{\mathrm{log}\psi }=0.2$ (taking into account the typical 1σ error bars on single data points, see Figure 9 of Madau & Dickinson 2014).

We define the merger efficiency $\eta (Z)$ as

$$\begin{eqnarray}&&\eta (Z)=\displaystyle \frac{{{ \mathcal N }}_{\mathrm{TOT}}(Z)}{{M}_{* }(Z)},\end{eqnarray} \tag{ 3 }$$

where ${{ \mathcal N }}_{\mathrm{TOT}}(Z)$ is the total number of CBs (BBHs, BHNSs, or BNSs) that have delay time (i.e., the time elapsed from the formation of the binary star to the merger of the two compact objects) ${t}_{\mathrm{del}}\leqslant 14\,\mathrm{Gyr}$ born from stars with metallicity Z in our population-synthesis simulations, and ${M}_{* }(Z)$ is the total initial stellar mass (corresponding to the zero-age main-sequence mass) simulated with metallicity Z. Thus, the merger efficiency is the number of mergers occurring in a population of initial stellar mass M_* and metallicity Z, integrated over a Hubble time (see, e.g., Giacobbo & Mapelli 2018; Klencki et al. 2018).

In Equation (1), the values of $\eta (Z)$ and ${ \mathcal F }(z^{\prime} ,z,Z)$ are estimated from catalogs of CB mergers obtained with population synthesis and with dynamical simulations, as detailed in the next sections. The catalogs contain information on the masses of the two compact objects, the delay time and the metallicity of the progenitor stars. In practice, since we have 6 (3) catalogs corresponding to 6 (3) different metallicities for isolated (dynamical) binary compact objects, the values of $\eta (Z)$ are linearly interpolated between the available metallicities (Figure 1).

Zoom In Zoom Out Reset image size

The value of ${ \mathcal F }(z^{\prime} ,z,Z)$ depends on the metallicity Z of stars that form at redshift $z^{\prime}$. To derive the average metallicity evolution as a function of redshift we use the following fitting formula:

$$\begin{eqnarray}&&\mu (z)=\mathrm{log}\left(\displaystyle \frac{Z(z)}{{Z}_{\odot }}\right)=\mathrm{log}a\,+\,b\,z,\end{eqnarray} \tag{ 4 }$$

where $a=1.04\pm \,0.14$ and $b=-0.24\pm 0.14$. In the above equation, the slope b comes from De Cia et al. (2018), who provide a fit to the metallicity evolution of a large sample of damped Lyα systems with redshift between 0 and 5. The original fit by De Cia et al. (2018) yields a metallicity $Z(z=0)=0.66$ ${\text{}}{Z}_{\odot }$, which is low compared to the average stellar metallicity measured at redshift zero (see, e.g., the discussion in Madau & Dickinson 2014). Hence, in Equation (4), we have rescaled the fitting formula provided by De Cia et al. (2018) to yield $Z(z=0)=(1.04\,\pm 0.14)\,{\text{}}{Z}_{\odot }$, where ${Z}_{\odot }=0.019$, consistent with the average metallicity of galaxies at $z\sim 0$ from the Sloan Digital Sky Survey (Gallazzi et al. 2008). The value of $a=1.04\pm 0.14$ adopted in Equation (4) is the result of this rescaling. The quoted uncertainties on both a and b are at 1σ, assuming (as done in the original papers by Gallazzi et al. 2008 and De Cia et al. 2018) that the observational values follow a Gaussian distribution.

We model the distribution of stellar metallicities $\mathrm{log}(Z/{\text{}}{Z}_{\odot })$ at a given redshift as a normal distribution with mean value $\mu (z)$ from Equation (4) and standard deviation⁶ ${\sigma }_{Z}=0.20$

$$\begin{eqnarray}&&p(z^{\prime} ,Z)=\displaystyle \frac{1}{\sqrt{2\pi \,{\sigma }_{Z}^{2}}}\,\exp \left\{-\,\displaystyle \frac{{\left[\mathrm{log}(Z/{\text{}}{Z}_{\odot })-\mu (z^{\prime} )\right]}^{2}}{2\,{\sigma }_{Z}^{2}}\right\}.\end{eqnarray} \tag{ 5 }$$

Based on our definition, ${ \mathcal F }(z^{\prime} ,z,Z)$ and $p(z^{\prime} ,Z)$ are connected by the following relation:

$$\begin{eqnarray}&&{ \mathcal F }(z^{\prime} ,z,Z)=\displaystyle \frac{{ \mathcal N }(z,Z)}{{{ \mathcal N }}_{\mathrm{TOT}}(Z)}\,p(z^{\prime} ,Z),\end{eqnarray} \tag{ 6 }$$

where ${ \mathcal N }(z,Z)$ is the number of CBs that form from stars with metallicity Z and merge at redshift z, while ${{ \mathcal N }}_{\mathrm{TOT}}(Z)$ is the total number of CBs that merge within a Hubble time and form from stars with metallicity Z (as already detailed above).

We performed 10³ realizations of Equation (1) per each considered model, in order to estimate the impact of observational uncertainties on the MRD. At each realization, we randomly draw the normalization value of the SFR density (Equation (2)), the intercept and the slope of the average metallicity (Equation (4)) from three Gaussian distributions with mean (standard deviation) equal to $\mathrm{log}\psi (0)=-2$ (${\sigma }_{\mathrm{log}\psi }=0.2$), a = 1.04 (${\sigma }_{a}=0.14$), and $b=-0.24$ (${\sigma }_{b}\,=0.14$), respectively. The value of the intercept and that of the slope are drawn separately, assuming no correlation. This procedure is implemented in the new python script cosmo ${ \mathcal R }$ate, which allows us to calculate up to 10³ models per day on a single core.

The catalogs of isolated binaries have been generated with our population-synthesis code mobse (Mapelli et al. 2017; Giacobbo et al. 2018; Giacobbo & Mapelli 2018; Mapelli & Giacobbo 2018). In mobse, the mass loss of massive hot stars is described as $\dot{M}\propto {Z}^{\beta }$, where β is defined as in Giacobbo et al. (2018):

$$\begin{eqnarray}\beta =\left\{\begin{array}{ll}0.85, & \mathrm{if}\,{{\rm{\Gamma }}}_{e}\leqslant 2/3\\ 2.45-2.4{{\rm{\Gamma }}}_{e}, & \mathrm{if}\,2/3\lt {{\rm{\Gamma }}}_{e}\leqslant 1\\ 0.05, & \mathrm{if}\,{{\rm{\Gamma }}}_{e}\gt 1.\end{array}\right.\end{eqnarray} \tag{ 7 }$$

In Equation (7), ${{\rm{\Gamma }}}_{e}$ is the Eddington ratio, i.e., the ratio between the luminosity of the star and its Eddington value.

mobse includes two different prescriptions for core-collapse supernovae (SNe) from Fryer et al. (2012): the rapid and the delayed SN models. The former model assumes that the SN explosion is launched ≲250 ms after the bounce, while the latter has a longer timescale (≳500 ms). In both models, a star is assumed to directly collapse into a black hole (BH) if its final carbon–oxygen mass is $\gtrsim 11\,{M}_{\odot }$. For the simulations described in this paper we adopt the rapid model, which enforces a gap in the mass function of compact objects between 2 and 5 ${M}_{\odot }$. Recipes for electron-capture SNe are included in mobse as described in Giacobbo & Mapelli (2019).

Prescriptions for pair instability and pulsational pair instability are implemented using the fitting formulas derived by Spera & Mapelli (2017). In particular, stars which grow a helium core mass $64\leqslant {m}_{\mathrm{He}}/{M}_{\odot }\leqslant 135$ are completely disrupted by pair instability and leave no compact objects, while stars with $32\leqslant {m}_{\mathrm{He}}/{M}_{\odot }\lt 64$ undergo a set of pulsations, which enhance mass loss and cause the final compact-object mass to be significantly smaller than it would be if we had accounted only for core-collapse SNe.

Natal kicks are randomly drawn from a Maxwellian velocity distribution. In the run presented here, we adopt a one-dimensional rms velocity $\sigma =15\,\mathrm{km}$ s⁻¹ for neutron stars. BH natal kicks are drawn from the same distribution as neutron-star kicks, but reduced by the amount of fallback as ${v}_{\mathrm{KICK}}=(1-{f}_{\mathrm{fb}})v$, where ${f}_{\mathrm{fb}}$ is the fallback parameter described in Fryer et al. (2012) and v is the velocity drawn from the Maxwellian distribution.

Binary evolution processes such as tidal evolution, Roche lobe overflow, common envelope, and GW energy loss are taken into account as described in Hurley et al. (2002). In particular, the treatment of the common envelope is described by the efficiency parameter α. In this work, we assume $\alpha =5$, as suggested by recent studies (Fragos et al. 2019; Giacobbo & Mapelli 2020). Orbital decay and circularization by GW emission are calculated according to Peters (1964).

We have simulated $6\times {10}^{7}$ isolated binaries with mobse, 10⁷ per each metallicity we considered (Z = 0.0002, 0.0008, 0.002, 0.008, 0.016, and 0.02). The mass of the primary star is randomly drawn from a Kroupa (2001) initial mass function, with minimum mass 5 ${M}_{\odot }$ and maximum mass 150 ${M}_{\odot }$. The orbital periods, eccentricities, and mass ratios of binaries are drawn from Sana et al. (2012). In particular, we derive the mass ratio $q={m}_{2}/{m}_{1}$ as ${ \mathcal D }(q)\propto {q}^{-0.1}$ with $q\in [0.1\mbox{--}1]$, the orbital period P from ${ \mathcal D }({\rm{\Pi }})\propto {{\rm{\Pi }}}^{-0.55}$ with ${\rm{\Pi }}={\mathrm{log}}_{10}(P/\mathrm{day})\,\in [0.15\mbox{--}5.5]$, and the eccentricity e from ${ \mathcal D }(e)\propto {e}^{-0.42}\,\mathrm{with}0\leqslant {\rm{e}}\leqslant 1$. These simulations are part of run CC15α5 in Giacobbo & Mapelli (2018).

We derive the catalogs of CB mergers from a set of direct N-body simulations already described in Di Carlo et al. (2020b) and Rastello et al. (2020). These dynamical simulations were ran with the direct N-body code nbody6++gpu (Wang et al. 2015, 2016), coupled with the population-synthesis code mobse, as already described in Di Carlo et al. (2019). In this way, the dynamical simulations include binary population synthesis, performed with the same code as the isolated binary simulations.

The masses of the simulated young star clusters range from 300 to 30,000 ${M}_{\odot }$. In particular, we consider $7.5\times {10}^{4}$ star clusters with mass ${M}_{\mathrm{SC}}\in [300,1000]\,{M}_{\odot }$ ($2.5\times {10}^{4}$ runs per each considered metallicity: Z = 0.0002, 0.002, and 0.02, from Rastello et al. 2020) and 3000 star clusters with mass ${M}_{\mathrm{SC}}\in [1000,{\rm{30,000}}]\,{M}_{\odot }$ (1000 runs per each considered metallicity: Z = 0.0002, 0.002, and 0.02, presented as set A in Di Carlo et al. 2020b). The total mass ${M}_{\mathrm{SC}}$ of a star cluster is drawn from a distribution ${dN}/{{dM}}_{\mathrm{SC}}\propto {M}_{\mathrm{SC}}^{-2}$, consistent with the mass function of young star clusters in the Milky Way (Lada & Lada 2003).

The initial half mass–radius r_h of star clusters is distributed according to the Marks and Kroupa relation (Marks et al. 2012), which relates the total mass of the star cluster ${M}_{\mathrm{SC}}$ with its initial half mass–radius r_h as

$$\begin{eqnarray}&&{r}_{h}={0.10}_{-0.04}^{+0.07}\,\mathrm{pc}\,{\left(\displaystyle \frac{{M}_{\mathrm{SC}}}{{M}_{\odot }}\right)}^{0.13\pm 0.04}.\end{eqnarray} \tag{ 8 }$$

The star clusters are initialized in virial equilibrium.

The initial distribution of stellar positions and velocities in the star clusters have been generated through the mcluster code (Küpper et al. 2011), according to a fractal distribution with fractal dimension D = 1.6 (Goodwin & Whitworth 2004). This ensures that the initial conditions of the simulated star clusters are clumpy and asymmetric as observed embedded star clusters. The mass of the stars is drawn from a Kroupa (2001) initial mass function between 0.1 and 150 ${M}_{\odot }$. The total initial binary fraction is ${f}_{\mathrm{bin}}=0.4$. The mass ratios between secondary and primary stars and the orbital properties of the binary systems (period and eccentricity) are drawn according to Sana et al. (2012), to ensure a fair comparison with the isolated binary simulations. The force integration includes a solar neighborhood-like static external tidal field. In particular, the simulated star clusters are assumed to be on a circular orbit around the center of the Milky Way with a semimajor axis of $8\,\mathrm{kpc}$ (Wang et al. 2016). Each star cluster is evolved until its dissolution or for a maximum time $t=100\,\mathrm{Myr}$.

Only three metallicities (Z = 0.0002, 0.002, and 0.02) were available from young star cluster simulations (Di Carlo et al. 2020b; Rastello et al. 2020). Running a larger metallicity set is computationally prohibitive. Thus, we linearly interpolated the merger efficiency $\eta (Z)$ (Figure 1) in our dynamical simulations to infer the values of $\eta (Z)$ for three additional metallicities (Z = 0.0008, 0.008, 0.016). We assigned to these three interpolated metallicities the available catalogs of dynamical CB mergers with the closest metallicity to the interpolated values.

Figure 1 shows the merger efficiency $\eta (Z)$ from young star clusters and isolated binaries. This quantity gives us an idea of the impact of the progenitor's metallicity on the merger rate in the different scenarios (isolated and dynamical) we considered. The trend of BNS merger efficiency with metallicity is similar in young star clusters and in isolated binaries, but isolated binaries are more efficient in producing BNS mergers. The main reason is that dynamical encounters may perturb the evolution of relatively low-mass binaries (such as BNSs and their progenitors), widening their orbit or even leading to their disruption (e.g., Hills & Fullerton 1980; Ye et al. 2020).

As already noted in several other works (e.g., Dominik et al. 2013; Giacobbo & Mapelli 2018; Klencki et al. 2018; Mapelli et al. 2019), the merger efficiency of BNSs is not significantly affected by the progenitor's metallicity.

The most interesting difference between isolated binaries and young star clusters is the behavior of BHNSs and BBHs at solar metallicity. The merger efficiency at solar metallicity is about a factor of 100 higher for BBHs/BHNSs formed in young star clusters than for BBHs/BHNSs formed in isolated binaries. The vast majority of dynamical BBH/BHNS mergers at solar metallicity originate from dynamical exchanges⁷ (see Di Carlo et al. 2020b for further details). This means that dynamical encounters tend to boost the merger rate of BBHs and BHNSs in the solar metallicity environment.

Figure 2 shows the MRD of BBHs as a function of time when considering young star clusters (i.e., dynamical binaries) and isolated binaries. In either case, we assume that the entire population of mergers forms from a single channel (i.e., either from young star clusters or from isolated binaries). It is more likely that a percentage of all mergers comes from young star clusters and another percentage from isolated binaries. In a follow-up paper (Y. Bouffanais et al. 2020, in preparation), we will try to constrain these percentages based on current LVC results. Here, we just want to compare the differences between the two scenarios.

Zoom In Zoom Out Reset image size

The MRD of BBHs (in both young star clusters and isolated binaries) grows with redshift (an MRD uniform in comoving volume would be a horizontal line in the plot), peaks at $z\sim 1.5\mbox{--}2.5$, and finally drops at $z\gt 2.5$. This trend is mostly determined by the cosmic SFR density, which peaks at $z\sim 2$, convolved with the delay time and the metallicity dependence. These results are fairly consistent with previous papers, which consider different population-synthesis models, metallicity evolution, and SFR evolution with redshift (e.g., Dominik et al. 2013; Belczynski et al. 2016; Mapelli et al. 2017; Mapelli & Giacobbo 2018; Artale et al. 2019; Neijssel et al. 2019; Tang et al. 2020).

At z = 0, the median values of the MRD of BBHs formed dynamically in young star clusters (hereafter, dynamical BBHs) and the one of isolated BBHs are ${R}_{\mathrm{BBH}}\sim 64$ and 50 Gpc⁻³ yr⁻¹, respectively. Both values are consistent with the ones inferred from O1 and O2 (Abbott et al. 2019b). The median merger rate of dynamical BBHs is higher than that of isolated BBHs up to $z\sim 4$ (see Table 1 for more details). This trend can be interpreted by looking at the merger efficiency (Figure 1): around solar metallicity, the dynamical channel is more efficient than the isolated channel. Hence, we expect a higher number of dynamical BBH mergers with short delay time in the local universe, where metallicity is higher. In contrast, the merger efficiency of dynamical BBHs formed from metal-poor stars (Z = 0.002) is a factor of $\sim 2$ lower than that of isolated BBHs with the same metallicity. Hence, isolated binaries are associated with a higher merger rate from very metal-poor systems.

Table 1.  MRD in $({\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1})$ for Five Redshift Intervals

	Redshift Intervals
	$z\in [0,0.1]$		$z\in [0.9,1.0]$		$z\in [1.9,2.0]$		$z\in [2.9,3.0]$		$z\in [3.9,4.0]$
CB	Dynamical	Isolated	Dynamical	Isolated	Dynamical	Isolated	Dynamical	Isolated	Dynamical	Isolated
BBH	${64}_{-20}^{+34}$	${50}_{-37}^{+71}$	${150}_{-52}^{+107}$	${92}_{-73}^{+178}$	${220}_{-77}^{+161}$	${207}_{-160}^{+256}$	${168}_{-71}^{+136}$	${130}_{-91}^{+192}$	${101}_{-51}^{+75}$	${105}_{-83}^{+191}$
BHNS	${41}_{-23}^{+33}$	${49}_{-34}^{+48}$	${114}_{-53}^{+80}$	${152}_{-120}^{+227}$	${168}_{-76}^{+138}$	${406}_{-331}^{+516}$	${142}_{-91}^{+129}$	${395}_{-286}^{+286}$	${99}_{-55}^{+62}$	${225}_{-124}^{+131}$
BNS	${151}_{-38}^{+59}$	${283}_{-75}^{+97}$	${473}_{-126}^{+192}$	${856}_{-249}^{+355}$	${460}_{-130}^{+177}$	${777}_{-228}^{+354}$	${247}_{-68}^{+98}$	${379}_{-113}^{+191}$	${110}_{-31}^{+44}$	${190}_{-63}^{+98}$

Note. We show a comparison between dynamical CBs formed in young star clusters and isolated CBs.

Download table as:  ASCIITypeset image

The MRD of isolated BBHs increases by a factor of $\sim 1.8$ from local universe up to z = 1, and then it grows up faster from redshift z = 1 to redshift z = 2 (Table 1). On the other hand, the MRD of dynamical BBHs increases almost with the same trend from z = 0 to $z\sim 2$ (i.e., without a change of slope at redshift $z\sim 1$). The main reason for the change of slope in the MRD of isolated BBHs is again the stronger dependence of the merger efficiency on metallicity. In the isolated model, most mergers at redshift $z\lt 1$ are due to BBHs that formed at higher redshift in lower metallicity environments ($Z\sim 0.0002$) and have a long delay time (Mapelli et al. 2017, 2018).

The uncertainty on MRD resulting from cosmic SFR and metallicity evolution is large, especially for the isolated scenario. For isolated BBHs, the 50% credible interval spreads over more than one order of magnitude between redshift 0 and 4. The 50% credible interval for the MRD of dynamical BBHs is contained within the credible interval of isolated BBHs. The 50% credible interval is smaller for dynamical BBHs, because the merger efficiency is less sensitive to metallicity in the dynamical scenario than in the isolated one (Figure 1).

Figure 3 shows the MRD evolution of BHNSs. At z = 0, ${{ \mathcal R }}_{\mathrm{BHNS}}={41}_{-23}^{+33}$ and ${49}_{-34}^{+48}$ Gpc⁻³ yr⁻¹ for dynamical and isolated BHNSs, respectively. At redshift z = 2, ${{ \mathcal R }}_{\mathrm{BHNS}}\,={168}_{-76}^{+138}$ and ${406}_{-331}^{+516}$ Gpc⁻³ yr⁻¹ for dynamical and isolated BHNSs, respectively. For most of the cosmic time, the boundaries of the 50% credible intervals of our two models have similar values. The higher boundary of the 50% credible interval for both dynamical and isolated BHNSs is below the upper limit from the LVC (${{ \mathcal R }}_{\mathrm{BHNS}}\lt 610$ Gpc⁻³ yr⁻¹; Abbott et al. 2019a), indicating that our model is consistent with O1 and O2 results. In the case of both BBHs and BHNSs, most of the uncertainty comes from metallicity evolution, because BBHs and BHNSs are extremely sensitive to metallicity variations (as shown in Figure 1).

Zoom In Zoom Out Reset image size

Finally, Figure 4 shows the MRD evolution of dynamical and isolated BNSs. At redshift $z\leqslant 0.1$, the MRD of dynamical BNSs (${{ \mathcal R }}_{\mathrm{BNS}}={151}_{-38}^{+59}$Gpc⁻³ yr⁻¹) is a factor of $\sim 2$ lower than the one of isolated BNSs (${283}_{-75}^{+97}$ Gpc⁻³ yr⁻¹). A similar difference is found at z = 2, where the MRD is ${{ \mathcal R }}_{\mathrm{BNS}}={460}_{-130}^{+177}$ and ${777}_{-228}^{+354}$ Gpc⁻³ yr⁻¹, for dynamical and isolated BNSs respectively. Overall, the MRD of dynamical BNSs is significantly lower than that of isolated BNSs, even if the MRD evolution with redshift is similar. This trend is expected by looking at Figure 1, because the merger efficiency of dynamical BNSs is lower at all metallicities. In young star clusters, the formation of BNSs is slightly suppressed with respect to isolated binaries, because such relatively low-mass binaries tend to be broken or softened (i.e., their orbital separation is increased) by dynamical encounters.

Zoom In Zoom Out Reset image size

The local MRD of isolated BNSs is consistent with that inferred from the LVC, while the local MRD of dynamical BNSs is below the 90% credible interval from the LVC. This suggests that (young) star clusters alone might not be able to explain all the BNS mergers detected by the LVC.

The models presented in this work assume small natal kicks for neutron stars, which are in tension with the proper motions of Galactic young pulsars (Giacobbo & Mapelli 2018). We recently proposed a new model for natal kicks that can reproduce the proper motions of Galactic pulsars and gives a value for the MRD close to the one presented in this study (Giacobbo & Mapelli 2020). As a result, we do not expect significant differences in the MRD between the model adopted in this work and the one proposed by Giacobbo & Mapelli (2020).

The 50% credible interval of simulated BNSs is significantly smaller than that of both BHNSs and BBHs, because BNSs are less sensitive to stellar metallicity (Figure 1). Hence, the uncertainty on BNS merger rate comes mostly from the SFR, for a fixed binary evolution model.

Our local MRDs for dynamical BNSs and BHNSs are higher than the values estimated by Ye et al. (2020) for globular clusters (${{ \mathcal R }}_{\mathrm{BNS}}\sim {{ \mathcal R }}_{\mathrm{BHNS}}\sim 0.02$ Gpc⁻³ yr⁻¹). This is not surprising because globular clusters form mostly at $z\gtrsim 2$, while smaller star clusters, like the ones we simulated, form all the time from high to low redshift and are an important channel of star formation in the local universe.

Figure 5 shows the mass distribution of BBHs, BHNSs, and BNSs merging across cosmic time. We plot together binaries merging at different redshift because we find no significant dependence of the mass distribution on the merger redshift, consistent with Mapelli et al. (2019). The main difference between the mass distribution of dynamical BBHs and that of isolated BBHs is that low-mass BBHs are less numerous in the former than in the latter scenario. Moreover, the maximum mass of merging BHs from isolated binaries is ${m}_{\mathrm{BH},\max }\sim 45\,{M}_{\odot }$, whereas dynamics in young star clusters leads to a significantly larger maximum mass ${m}_{\mathrm{BH},\max }\,\sim 90\,{M}_{\odot }$. Quantitatively, the percentage of isolated BBHs that have a primary mass $\gt 40\,{M}_{\odot }$ is equal to 0.07%, while it is 10.6% for dynamical BBHs. This marked difference in the maximum mass of merging BHs between isolated and dynamical BBHs can be understood as follows (see also Di Carlo et al. 2019, 2020a). The stellar wind and core-collapse SN prescriptions adopted in mobse allow the formation of BHs with mass up to $\sim 65\,{M}_{\odot }$ (Giacobbo et al. 2018), but only BHs with masses up to $\sim 45\,{M}_{\odot }$ are able to merge within a Hubble time in isolated BBHs, because of a subtle interplay between mass transfer and stellar radii. In fact, BHs with masses $\gt 45\,{M}_{\odot }$ form only from stars with zero-age main-sequence mass $\sim 60\mbox{--}80\,{M}_{\odot }$ which retain a large fraction of hydrogen envelope and collapse to a BH directly (Figure 4 of Giacobbo et al. 2018). When such stars are members of a tight binary system, most of the hydrogen envelope is removed by mass transfer (or by common envelope) before the collapse; hence, even if they might end up in a BBH merger, the mass of the final BHs will be smaller than the one we expect from single star evolution. In contrast, if such stars are members of loose binaries (initial orbital separation $a\gtrsim {10}^{4}$ ${\text{}}{R}_{\odot }$), which do not undergo mass transfer, they produce BBHs with individual BH masses $\gt 45\,{M}_{\odot }$, but the orbital separation is too large to lead to coalescence.

Zoom In Zoom Out Reset image size

In young star clusters, instead, BHs with masses $\gt 45\,{M}_{\odot }$ are able to merge within a Hubble time, because (i) if they form from the collapse of single stars, they can acquire companions through dynamical exchanges, and (ii) if they are members of loose binaries, these massive binaries are efficiently hardened by three body encounters (Di Carlo et al. 2019). Moreover, (multiple) stellar mergers can even lead to the formation of BHs with masses $\gt 65\,{M}_{\odot }$, as discussed in Di Carlo et al. (2020a). Such massive BHs are single at birth but can acquire a companion by dynamical exchanges.

Figure 5 shows that dynamical BHNSs can host significantly more massive BHs than isolated BHNSs. Only $9\times {10}^{-4} \%$ of BHs in isolated BHNSs have masses ${m}_{\mathrm{BH}}\gt 20\,{M}_{\odot }$, while 1.6% of BHs in dynamical BHNSs have masses above this value. This is another effect of dynamics, which boosts the formation of massive binaries by dynamical exchanges and facilitates the coalescence of binaries with extreme mass ratio by dynamical hardening (see the discussion in Rastello et al. 2020 for additional details). Finally, we do not find any significant difference between the mass distribution of dynamical BNSs and that of isolated BNSs.⁸