DOUBLE COMPACT OBJECTS. II. COSMOLOGICAL MERGER RATES

In order to determine the merger rates of DCOs, we need the star formation rate (SFR). We adopt the formula provided by Strolger et al. (2004):

$$\begin{equation} {\rm SFR}=10^9a \left(t^b e^{-t/c} +de^{d(t-t_0)/c} \right) \,{M}_{\odot }\ {\rm yr^{-1} Gpc^{-3}}, \end{equation} \tag{ 1 }$$

where t is the age of the universe (Gyr) as measured in the rest frame, t₀ is the present age of the universe (13.47 Gyr; see Section 4), and the parameters have values a = 0.182, b = 1.26, c = 1.865, and d = 0.071. The SFR described above is expressed in comoving units of length and time.

For redshifts z < 4, we describe the distribution of galaxy masses using a Schechter-type probability density function, calibrated to observations (Fontana et al. 2006):

$$\begin{equation} \Phi (M_{{\rm gal,z}})=\Phi ^*(z)\ln (10)a^{1+\alpha (z)}e^{-a}, \end{equation} \tag{ 2 }$$

where Φ*(z) = 0.0035(1 + z)^−2.2, $a=M_{\rm gal} \cdot 10^{-M_{\rm z}}$ (M_z = 11.16 + 0.17z − 0.07z²), and α(z) = −1.18 − 0.082z. A galaxy mass is drawn from this distribution in solar units (M_☉) and in the range 7 < log (M_gal) < 12. Beyond redshift z = 4, we assume no further evolution in galaxy mass, fixing the mass distribution to the value at z = 4. This assumption reflects the lack of information on galaxy mass distribution at high redshift.

We assume the average oxygen-to-hydrogen number ratio (F_OH = log (10¹²O/H)) in a typical galaxy to be given by

$$\begin{equation} \log (F_{\rm OH})=s+1.847 \log (M_{\rm gal})- 0.08026 (\log (M_{\rm gal}))^2. \end{equation} \tag{ 3 }$$

As suggested by Erb et al. (2006) and Young & Fryer (2007), the functional form of this mass–metallicity relation is redshift independent, with only the normalization factor s varying with redshift. We describe the redshift dependence of galaxy metallicity using the average metallicity relation from Pei et al. (1999):

$$\begin{equation} Z \propto \left\lbrace \begin{array}{@{}l r}10^{-a_2 z} & z < 3.2\\ 10^{-b_1-b_2 z} & 3.2 \le z < 5\\ 10^{-c_1-c_2 z} & z \ge 3.2\\ \end{array}, \right. \end{equation} \tag{ 4 }$$

which implies the evolution of s with redshift:

$$\begin{equation} s \,{\propto} \left\lbrace \begin{array}{@{}l l}-a_2z \,{-}\,1.492 & z < 3.2 \\ -b_2z \,{-}\,3.2(a_2\,{-}\,b_2)\,{-}\,1.492 & 3.2 \le z < 5 \\ -c_2z \,{-}\,5(b_2\,{-}\,c_2) \,{-}\,3.2(a_2\,{-}\,b_2)\,{-}\,1.492 & z \ge 3.2 \\ \end{array} \right. \!\!\!\!. \end{equation} \tag{ 5 }$$

We assume that the oxygen abundance (used in F_OH) correlates linearly with the average abundance of elements heavier than helium (encapsulated in the metallicity measure, Z).

In this paper, we employ two distinct scenarios for metallicity evolution with redshift in order to investigate the uncertainties of the chemical evolution of the universe. The construction of these scenarios consists of several steps. (1) We utilize two normalizations of Equation (3). In the first, provided by Pei et al. (1999), the coefficients are a₂ = 0.5, b₁ = 0.8, b₂ = 0.25, c₁ = 0.2, and c₂ = 0.4. This grants a rate of average metallicity evolution, which we label slow. The second, provided by Young & Fryer (2007), uses a₂ = 0.12, b₁ = −0.704, b₂ = 0.34, c₁ = 0.0, and c₂ = 0.1992. It is based on ultraviolet Galaxy Evolution Explorer, Sloan Digital Sky Survey (SDSS), and Spitzer infrared and Super Kamiokande neutrino observations (Hopkins & Beacom 2006). This normalization produces a faster rate of chemical evolution and we label the results fast. At this point, for each galaxy mass value at a given redshift (Equation (2)), we have two metallicity values (Equation (3)). (2) We then combine these (slow and fast) metallicities into a single value that is an average of the two; we label this profile as initial. However, this profile yields an unrealistically high number of galaxies with extrasolar (up to 3 Z_☉; Z_☉ = 2% of stellar mass) metallicities at redshift z ∼ 0. (3) In order to be consistent with observational data, we scale down the profile so that it agrees with the observed metallicities of galaxies in the local universe (at z ∼ 0). We explore two such "extreme" scalings resulting in a pair of final metallicity evolution profiles. In the first, we divide the metallicity values from the initial profile by a factor of 1.7. This produces a median value of metallicity of 1.5 Z_☉ at z ∼ 0 (see Figure 1), which corresponds to 8.9 in the "12+log(O/H)" formalism. This calibration was designed to match the upper 1σ scatter of metallicities according to Yuan et al. (2013) (see their Figure 2, top-right panel). We label this profile as high-end, as it is the upper limit on metallicity at z ∼ 0. In the second, we utilize SDSS observations (Panter et al. 2008), from which we infer that one half of the star forming mass of the galaxies at z ∼ 0 has 20% solar metallicity, while the other half has 80% solar metallicity. This yields a median metallicity value of 0.8 Z_☉ and requires the division of the initial profile by a factor of three. We label this profile as low-end.

We distinguish three stellar populations:

$$\begin{equation} \begin{array}{l l}F_{\rm OH,gal} < 10^{-4} & {\rm Population III} \\ 10^{-4} \le F_{\rm OH,gal} \le 10^{-1} & {\rm Population II} \\ F_{\rm OH,gal} > 10^{-1} & {\rm Population I} \end{array}. \end{equation} \tag{ 6 }$$

We choose F_{OH, gal} = 10⁻⁴ as the delineation point between Population II and III stars. A lower abundance of metals provides insufficient cooling in the collapse of gas clouds and thus significantly alters the star formation for Population III stars (e.g., Mackey et al. 2003; Smith et al. 2009). The border point between Population II and I stars is dictated by observations in the Milky Way (e.g., Binney & Merrifield 1998; Beers & Christlieb 2005).

We assume that the binary fraction is 50%: for each single star, there exists one binary. We additionally assume that all the stars within each galaxy share the same metallicity value. The use of average metallicity seems to be appropriate since we draw a large (10⁴) number of galaxies (Equation (2)) via Monte Carlo simulations.

To calculate the evolution of the stellar populations, we utilize the recently updated (Belczynski et al. 2010a, 2012; Dominik et al. 2012) Startrack population synthesis code (Belczynski et al. 2002, 2008). This code can evolve isolated binary stars that are interacting in quasi-hydrostatic equilibrium from the zero-age main sequence (ZAMS), through mass transfer, to the formation of compact objects and to the ultimate merger of the binary components. The code makes use of an extensive set of formulae and prescriptions that adequately approximate more detailed binary evolution calculations (see Hurley et al. 2000).

StarTrack allows one to investigate stable and unstable mass transfers between the binary components. Stable transfer calculations have been calibrated on massive binaries that are relevant to DCO formation (e.g., Tauris & Savonije 1999; Wellstein et al. 2001). It is yet unknown exactly how conservative the stable mass transfer is. Dewi & Pols (2003) suggest that in massive binaries the fraction of the envelope of the donor accreted by its companion ranges between 40% and 70%. In our calculations, we fix this value to be 50% or, in mathematical terms:

$$\begin{equation} \dot{M_{\rm acc}}=f_{\rm a}\dot{M_{\rm don}}, \end{equation} \tag{ 7 }$$

where $\dot{M_{\rm acc}}$ is the accretion rate, $\dot{M_{\rm don}}$ is the mass transfer rate from the donor, and f_a is the fraction of the rate transferred (here equal to 0.5). The remaining mass is expelled to infinity. The dynamically unstable mass transfers (CE) is calculated according to the energy balance formula (Webbink 1984), with the envelope binding energy parameter λ adopted from Xu & Li (2010).

Tidal interactions and their influence on eccentricity, the semi-major axis, and rotation is also evaluated. The calculations are done with the standard equilibrium-tide, weak-friction approximation (Zahn 1977, 1989), using the formalism of Hut (1981). However, the code does not allow one to investigate the influence that rotation of the components has on the internal structure of the components.

Stellar winds are taken into account as a function of the metallicity and evolutionary stage of the star. This piece of physics is especially important as it has a significant impact on the masses of remnant objects, which are the centerpiece of this study. In short, the wind mass loss rates are divided into categories specific to the evolutionary stage of the star: O/B–type, red giant, asymptotic red giant, Wolf–Rayet stars, and luminous blue variable (LBV) stars. The magnitude of the rates increases with metallicity of the star except for the LBV phase. In this stage, the winds are set to be of the order of 10⁻⁴ M_☉ yr⁻¹. This value was calibrated to account for the highest mass BHs in the Milky Way ∼15 M_☉ (Cyg X-1 and GRS 1915). A detailed description of wind mass loss rates can be found in Belczynski et al. (2010a).

Besides stellar winds, the code also calculates changes of the angular momentum arising from gravitational radiation and magnetic braking. The latter is adopted from Ivanova & Taam (2003).

Additionally, the code utilizes the convection-driven, neutrino-enhanced SN engines (Fryer et al. 2012) to determine the properties of the remnant objects (neutron stars (NSs) and BHs).

For each metallicity value in each model, we evolve 2 × 10⁶ binaries, assuming that each component is created at the same time. Each binary system is initialized by four parameters that are assumed to be independent. These are: primary mass, M₁ (initially more massive component), mass ratio, q = M₂/M₁, where M₂ is the mass of the secondary component (initially less massive), the semi-major axis, a, of the orbit, and the eccentricity, e. The mass of the primary component is randomly chosen from the initial mass function adopted from Kroupa et al. (1993) and Kroupa & Weidner (2003):

$$\begin{equation} \Psi (M_1) \propto \left\lbrace \begin{array}{@{}l c}M_1^{-1.3} & \quad 0.08 \ {M}_{\odot }\le M_1 < 0.5 \ {M}_{\odot }\\ M_1^{-2.2} & \quad 0.5 \ {M}_{\odot }\le M_1 < 1.0 \ {M}_{\odot }\\ M_1^{-\alpha } & \quad 1.0 \ {M}_{\odot }\le M_1 < 150 \ {M}_{\odot }, \\ \end{array} \right. \end{equation} \tag{ 8 }$$

where α = 2.7 is our standard choice for field populations. The choice of the upper initial mass function (IMF) cutoff (150 M_☉) is justified by observations of massive stars in the Milky Way (Figer 2005; Oey & Clarke 2005). Stars are generated from within an initial mass range, with the limits based on the targeted stellar population. For example, studies of single NSs require their evolution within the range 8–20 M_☉, while for single BHs the lower limit is 20 M_☉. Binary evolution may broaden these ranges due to mass transfer episodes and we therefore set the minimum mass of the primary to 5 M_☉. We assume a flat mass ratio distribution, Φ(q) = 1, over the range q = 0–1, in agreement with recent observations (Kobulnicky & Fryer 2007). Given a value of the primary mass and the mass ratio, we obtain the mass of the secondary from M₂ = qM₁. However, for the same reasons as for the primary, we do not consider binaries where the mass of the secondary is below 3 M_☉. The distribution of initial binary separations is assumed to be flat in log (a) (Abt 1983) and so ∝1/a, with a ranging from values such that at ZAMS the primary fills no more than 50% of its Roche lobe to 10⁵ R_☉. For the initial eccentricity, we adopt a thermal equilibrium distribution (e.g., Heggie 1975; Duquennoy & Mayor 1991): Ξ(e) = 2e, with e ranging from 0 to 1. We find that the adopted parameters are in accordance with the most recent observations of O-star populations (Sana et al. 2013).

3.2.1. The Standard Model

In this subsection, we define a reference model for this paper. This model is identical to the "Standard model–submodel B" in the previous paper in this series (Dominik et al. 2012).

The list of major parameters describing the input physics of binary evolution in this model begins with the Nanjing λ (Xu & Li 2010) CE coefficient used in the energy balance prescription (Webbink 1984). This λ value depends on the evolutionary stage of the donor, its mass at ZAMS and the mass of its envelope, and its radius. In addition, all of these quantities depend on metallicity, which in our simulations changes within a broad range (Z = 10⁻⁴–0.03).

However, before calculating the aforementioned energy balance to determine the outcome of the CE, we check the evolutionary type of the donor star. For example, main sequence (MS) stars do not have clear core-envelope division, as the helium core is still in the process of being developed. Donors on the HG behave similarly, although it remains unclear if such a division can appear on the HG or not until later stages, like the core helium burning stage (P. Eggleton 2010, private communication). In our previous work, we investigated two possibilities of the CE outcome associated with the type of the donor star: an automatic (premature) merger if the donor is an HG star, regardless of the energy balance (labeled as "submodel B") or allowing the CE energy balance to unfold ("submodel A").

The case in which we allow for potential survival of systems with HG donors results in very high Advanced LIGO/VIRGO detection rates (∼8000 yr⁻¹; Belczynski & Dominik 2012), exceeding even the empirically estimated rates based on IC10 X-1 and NGC 300 X-1 (∼2000 yr⁻¹; Bulik et al. 2011; Belczynski et al. 2013). Therefore, we only show one model with this generous assumption on CE physics, which leads to the most optimistic of our predictions. This model (Optimistic CE) will be tested (and probably quickly eliminated) by the upper limits from the Advanced LIGO/VIRGO engineering runs. For all the other models, including our reference model, we make the conservative assumption that none of the HG donor CE phases lead to the formation of DCOs.

Observations suggest (Hobbs et al. 2005) that NSs formed in SNe receive natal kicks, with velocities drawn from a Maxwellian distribution with σ = 265 km s⁻¹. We employ these findings in our simulations and extend them so that BH natal kicks match this distribution as well. However, it is possible that some matter ejected during the explosion will not reach the escape velocity and will thus fall back on the remnant object, potentially stalling the initial kick. To account for this process, we modify the Maxwellian kicks by the amount of matter falling back onto the newly formed compact object:

$$\begin{equation} V_{\rm k}=V_{\rm max}(1-f_{\rm fb}), \end{equation} \tag{ 9 }$$

where V_k is the final magnitude of the natal kick, V_max is the velocity drawn from a Maxwellian kick distribution, and f_fb is the fallback factor describing the fraction of the ejecta returning to the object. The values of f_fb range between 0–1, with 0 indicating no fallback/full kick and 1 representing total fallback/no kick (a "silent supernova," e.g., Mirabel & Rodrigues 2003). We label this the "constant velocity" formalism. An alternative approach is the "constant momentum" formalism, where the kick velocity is inversely proportional to the mass of the remnant object. In general, constant velocity provides larger natal kicks on average than constant momentum, resulting in more frequent disruptions of binaries, especially for systems with BHs. Therefore, we choose the "constant velocity" formalism over the "constant momentum" as it provides a more conservative limit on the number and therefore merger rates of systems containing BHs.

This model also utilizes the "Rapid" convection-driven, neutrino-enhanced SN engine (Fryer et al. 2012). It allows for a successful explosion without the need for the artificial injection of energy into the exploding star. In this scenario, the explosion starts from the Rayleigh–Taylor instability and occurs within the first 0.1–0.2 s. For low-mass stars (M_zams ≲ 25 M_☉), the result is a very strong (high-velocity kick) SN, which generates an NS. For higher mass stars, a BH is formed through a direct collapse (a failed SN). This engine, incorporated into binary evolution, successfully reproduces the mass gap (Belczynski et al. 2012) observed in Galactic X-ray binaries (Bailyn et al. 1998; Özel et al. 2010).

The list of major physical parameters used in this and subsequent models is given in Table 1. More details on the physics described above can be found in Dominik et al. (2012).

3.2.2. Variations on the Standard Model

NS–NS. As shown in Figure 3, the intrinsic merger rate densities of double NS systems peak at redshift z ≈ 1 (∼200 yr⁻¹ Gpc⁻³). As a general rule, the merger rates of all types of DCO are directly related to the SFR. However, for a given SFR value, the formation efficiency of different DCOs may vary. In other words, the proportions of NS–NS, BH–NS, and BH–BH systems may differ beyond the regime set by the IMF (e.g., Dominik et al. 2012). For example, NS–NS systems are on average efficiently created in high-metallicity environments (see Figure 4). When combined with the peak of the SFR at z ∼ 2 (the average high-end metallicity is ∼0.4 Z_☉; see Figure 2), high-metallicity NS–NS formation efficiency is enhanced, thus creating the profile shown in Figure 3. What is characteristic for this profile is the "hump" that arises at z ∼ 1.6, approaching from high redshifts. As can be seen in Figure 4, this shape is dominated by mergers originating from 0.75 Z_☉ environments. The reason for this increase in merger rate densities, when transiting from 0.5 Z_☉ environments (higher redshifts) to higher metallicity ones (lower redshifts), is a consequence of the applied CE approach. Within the framework of the Nanjing CE treatment adapted for the Startrack code, the binding energy of the CE decreases at the 0.5–0.75 Z_☉ boundary, allowing for the survival of a larger number of NS–NS progenitors.

In comparison, for the low-end metallicity profile, the NS–NS systems dominate the rates only up to z ≈ 0.5 (Figure 5). This is a consequence of the adopted metallicity evolution scenario. Specifically, galaxies of a given metallicity are shifted to lower redshift when compared with the high–end scenario, causing a shift of the NS–NS systems also to lower redshifts.

As shown in Figure 6, in the observer frame the systems dominate the merger rates up to redshift z ≈ 2.4. However, decreased metallicity for the low-end case shifts this point to z ≈ 0.6 (Figure 7).

BH–NS. For the high-end metallicity evolution, the rest-frame merger rate densities for BH–NS systems shown in Figure 3 peaks at a value of ∼50 Gpc⁻³ yr⁻¹ at redshift z ∼ 3. However, the merger rate efficiency drops for low (z ∼ 0) and high (z ∼ 6) redshifts. This is because of properties of the progenitor masses. For metallicities ∼Z_☉, the bulk of the progenitors masses are in the range 45–60 M_☉ for the primary component and 22–32 M_☉ for the secondary. Pairs of progenitors outside these ranges are unlikely. The upper mass limit delineates between BH–NS and BH–BH systems; crossing it results in the formation of the latter systems instead of the former. The lower mass limit is set by a similar phenomenon, only this time through BH–NS/NS–NS formation. Progenitors of these systems for metallicities a factor of ∼10 lower than Z_☉ must have lower masses on average, primarily because of the decreased stellar wind mass losses. Otherwise, the binary would retain enough mass to form a BH–BH system or go through a terminal CE event. Therefore, the mass ranges for BH–NS progenitors for Z ∼ 10% Z_☉ are 20–50 M_☉ for the primary and 12–25 M_☉ for the secondary. Given that in this mass range the IMF scales as M^−2.7 (where M is the mass of the progenitor), there are more BH–NS progenitors available at moderately low metallicity than at higher values. This, in turn, translates into increased merger rates arising from these environments. Decreasing the metallicity to ∼1% Z_☉ decreases the masses of the progenitors even further due to the same wind effects. However, in this case the BH progenitors are closely approaching their lower mass limit (∼20 M_☉), which leaves a narrow mass range: 20–25 M_☉ for the primary and 18–22 M_☉ for the secondary. There are fewer progenitors in these mass ranges when compared with the previous case and therefore we find a lower merger rate. Overall, the BH–NS merger rates peak originates from systems being created at moderate metallicities (see Figure 4).

BH–BH. For these systems, the intrinsic merger rate has a peak-plateau at a rate of ∼300–400 Gpc⁻³ yr⁻¹ at z ∼ 4–8, for the high-end case. The low-metallicity galaxies abundant at high redshifts are efficient BH factories (see Figures 8 and 9). This also means that adopting the low-end metallicity scenario will allow for more BH–BH systems to form at lower redshifts, when compared with the high-end. Additionally, environments with low amounts of metals favor massive BHs. For example, the most massive BH–BH system acquired in this model consisted of a 62 M_☉ and a 74 M_☉ BH pair. These systems originate from the extremely low metallicity environments (Z = 0.0001). We find that such systems merge up until redshift z ∼ 3 and z ∼ 2 for the high-end and low-end metallicity evolution models, respectively. However, due to statistical uncertainties, these redshift values may be even lower. These massive systems originate through the standard BH–BH formation channel. As an instructive example, we detail the formation scenario of an 8.3–5.8 M_☉ BH–BH, for Z = 0.005—a typical system for the average metallicity acquired in our study: t = 0 Myr. The components start with masses 32 M_☉ and 25 M_☉ for the primary and secondary, respectively, and an orbital separation a = 995 R_☉. t = 6.7 Myr. The primary, after becoming an HG star, expands and initiates a mass transfer through Roche lobe overflow (RLOF). The transfer continues until the primary loses almost all of its hydrogen envelope and becomes a Wolf–Rayet star with 10 M_☉ (the secondary component has 35 M_☉). The orbital separation prior to RLOF was a = 1000 R_☉ and a = 1600 R_☉ after t = 7.0 Myr. The primary explodes as an SN, forming a 7.8 M_☉ BH. The orbital separation after the explosion was a = 1760 R_☉ t = 8.7 Myr. The secondary (34 M_☉) initiates a CE phase and becomes a Wolf–Rayet star with 11 M_☉ as a result of the outcome (the primary gained ∼0.5 M_☉). The orbital separation prior to the CE was a = 1780 R_☉ and a = 2.6 R_☉ after t = 9.4 Myr. The secondary undergoes an SN and becomes a 5.8 M_☉ BH. The orbital separation prior to the explosion was a = 2.8 R_☉ and a = 3 R_☉ after t = 26 Myr. The coalescence of an 8.3 M_☉–5.8 M_☉ BH–BH system occurs. This example is illustrated by a diagram in Figure 10.

Zoom In Zoom Out Reset image size

On a side note, the formation of the most massive BH–BH systems on close orbits may be questionable. The progenitors of the aforementioned 62 M_☉–74 M_☉ BH–BH system are massive stars (140 M_☉–150 M_☉ at ZAMS). A recent theoretical study by Yusof et al. (2013) suggests that such objects (150 M_☉–500 M_☉) will not increase in size significantly during their evolution. Therefore, it is more likely for such binaries to bypass the CE phase and avoid the reduction of orbital separation. This in turn will prevent the resulting BH–BH system from merging within a Hubble time.

In the observer frame, BH–BH systems begin to dominate the merger rates at z ∼ 2. For the low-end case, this happens closer to z ∼ 1.

In this model, we relax one of the conditions on CE survivability. Specifically, HG donors are now allowed to undergo full energy balance calculations. In the Standard model, CEs with HG donors resulted in an immediate merger, terminating further binary evolution. This has been shown to have a significant impact on the number of DCOs, altering the merger rates by orders of magnitude (Belczynski et al. 2007, 2010b). When HG donors survive CE, their numbers naturally increase, as do their merger rates.

NS–NS. When compared with the Standard model, the high-end case, the intrinsic merger rate of NS–NS systems peaks at higher redshift (z ∼ 3) and at higher values (∼1000 Gpc⁻³ yr⁻¹). The shift in the peak toward higher redshifts is associated with the systems having shorter delay times on average, which allows them to merge more quickly after formation. As expected, the decrease of average delay times for NS–NS systems is caused by the new CE condition. In the Standard model, the only surviving binaries were those that did not initiate the CE while the potential donor was an HG star. In order to prevent a rapidly expanding HG star from overfilling its Roche lobe, these binaries had to have a significant initial separation, which resulted in relatively large delay times. In this model, the CE phase with an HG donor is allowed, so initial separation is no longer such a crucial issue. Therefore, binaries with smaller initial separations are able (if they have sufficient orbital energy) to survive and form NS–NS systems. This results in shorter delay times. In the low-end case, the same mechanism causes the peak to shift toward z ∼ 2.

In the observer frame, the merger rate of NS–NS systems is a few times higher than in the Standard model for both high-end and low-end metallicity evolutions. In the former case, NS–NS systems dominate the merger rates up to z ≲ 0.5. In the latter case, they are always sub-dominant compared with BH–BH systems.

BH–NS. The binaries forming these systems usually undergo two CE events in their lifetimes, due to their relatively high initial mass ratios (for details, see Dominik et al. 2012). The two CEs reduce the initial separation, which makes the relaxed CE condition much less relevant than for the NS–NS case discussed above. The result is that there are no significant changes in the intrinsic merger rate density for BH–NS systems. As in the Standard model, the mergers of BH–NS systems are the rarest of all types of DCOs. This is true for both of the metallicity scenarios.

BH–BH. These systems do not experience two CE events, unlike the BH–NS systems and therefore they do not reduce their initial separations as efficiently. The peak of the intrinsic merger rate density shifts slightly toward lower redshift (z ∼ 4; high-end) when contrasted with the reference model (see Figure 3). This is because of the effect of metallicity on the outcome of the CE phase. The larger the fraction of metals in a star, the bigger its radius (e.g., Hurley et al. 2000). This effect is particularly strong during the HG phase. Therefore, high-metallicity BH progenitors are more likely to initiate CE on the HG. In the Standard model, this is not allowed and such systems are removed from the population. However, here we relax this condition and as a consequence we add more BH–BH systems originating from higher metallicities (see Figures 8 and 9). For the low-end case, this results in a peak-plateau between redshifts 3 < z < 4. This is because of the higher metallicities appearing at lower redshifts when compared with the high-end case.

In the observer frame, the BH–BH systems start to dominate the merger rates at z ≈ 0.5 in the high-end case. For the low-end case, these DCOs are always primary mergers.

In this model, we change the SN explosion engine with respect to the Standard model. The Standard model uses the Rapid engine, which yields a gap between 2–5 M_☉ in the masses of the resulting compact objects. Here, we utilize the Delayed scenario (for details, see Section 3.2.2). The main feature of this engine is that it produces a continuous mass spectrum of remnant objects (Belczynski et al. 2012). As suggested by Kreidberg et al. (2012), the presence of the mass gap feature may be a result of systematic errors arising from misinterpretation of the BH binary light curve analyses. The resulting errors in estimating the inclination of the binary may shift low-mass BHs from the gap. The distinction between the two engines is clearly visible in Figures 8 and 9. The minimal total mass for this model is ∼5 M_☉. Such a system is composed of two BHs of 2.5 M_☉ each (2.5 M_☉ being the delineation between upper NS and lower BH mass). For other models, the minimal BH mass is ∼5 M_☉, thus yielding a minimal total mass ∼10 M_☉. However, the SN engine effects do not play a significant role in the merger rates of any type of DCOs.

Here, we employ full natal kicks (as measured for NSs) just on BHs (see Section 3.2.2). This is performed regardless of the amount of fallback (see Equation (9)). The kicks for NS–NS systems remain unchanged, as does their population with respect to the Standard model.

In this variation, the velocity of the natal kick acquired upon BH formation will disrupt many binaries that would otherwise (in the Standard model) form coalescing BH–NS or BH–BH systems, as is clearly visible in Figure 3. As a consequence, the NS–NS systems will dominate the merger rates in the observer frame.

In addition, the full natal kick will affect the most massive BHs. In the Standard model, stars with masses M_zams > 40 M_☉ would collapse directly into a BH after the SN explosion; with no asymmetric ejecta, they do not receive a kick (f_fb = 1, Equation (9)). However, in this model, these stars receive a maximum velocity kick and thus often disrupt the system. As a consequence, the probability of the formation and eventual merger of the most massive BH–BH systems is lowered significantly, which can be seen in the bottom panel of Figures 8 and 9.