Double compact objects II: Cosmological merger rates
The development of advanced gravitational wave (GW) observatories, such as Advanced LIGO and Advanced Virgo, provides impetus to refine theoretical predictions for what these instruments might detect. In particular, with the range increasing by an order of magnitude, the search for GW sources is extending beyond the “local” Universe and out to cosmological distances. Double compact objects (neutron star-neutron star (NS-NS), black hole-neutron star (BH-NS) and black hole-black hole (BH-BH) systems) are considered to be the most promising gravitational wave sources. In addition, NS-NS and/or BH-NS systems are thought to be the progenitors of gamma ray bursts (GRBs), and may also be associated with kilonovae. In this paper we present the merger event rates of these objects as a function of cosmological redshift. We provide the results for four cases, each one investigating a different important evolution parameter of binary stars. Each case is also presented for two metallicity evolution scenarios. We find that (i) in most cases NS-NS systems dominate the merger rates in the local Universe, while BH-BH mergers dominate at high redshift; (ii) BH-NS mergers are less frequent than other sources per unit volume, for all time; and (iii) natal kicks may alter the observable properties of populations in a significant way, allowing the underlying models of binary evolution and compact object formation to be easily distinguished. This is the second paper in a series of three. The third paper will focus on calculating the detection rates of mergers by gravitational wave telescopes.
Among the potential sources of GW, the merger of double compact objects (DCOs) is considered the most promising one for the first detection. The next generation of gravitational wave observatories (i.e., Advanced LIGO, Advanced Virgo, KAGRA) will probe the Universe in search for DCO signatures at unprecedented distances, reaching cosmological scales (). In this paper we present predictions for DCO merger rates from isolated (i.e., field population) DCO progenitors as a function of cosmological redshift.
The distribution of binary coalescence as a function of redshift has been investigated by several authors. An important initial work was the investigation of the redshift distribution of GRBs (e.g. Totani (1999)). Preliminary work on the importance of GW measurements of chirp mass distributions was done by Bulik et al. (2004), while initial studies of the GW confusion background have been presented in Regimbau & Hughes (2009).
In the first paper in this series (Dominik et al. (2012), first in the series) we investigated the sensitivity of DCO formation to major uncertainties of binary evolution (regarding mostly supernovae and common envelope episodes (CE)). We presented several models to bracket the current uncertainty in the phenomena deciding the fate of DCO systems. Building on this work, in the current study we present a set of four evolutionary models. In addition to a standard (reference) model, we have added models investigating a range of Hertzsprung gap CE donors, supernova (SN) explosion engines, and BH natal kicks (see Section 3.2 and Table 1). Additionally, for each model we have performed the evolutionary calculations for 11 metallicity values, allowing us to cover the abundance of metals in Population I and II stars (see Sections 2.4 and 3).
To account for the varied chemical composition of the Universe, we perform the cosmological calculations for two scenarios of metallicity evolution, that we will call “low–end” and “high–end” , respectively. These yield distinct rates of average metallicity growth, allowing us to “bracket” the associated uncertainties (see Section 2.3, Fig. 1 and Fig. 2).
In this study we investigate field stellar populations only. However, recent studies (e.g., Morscher et al. (2012)) suggest that mergers in globular clusters may add a significant contribution to the overall coalescence rates in the Universe. In this sense, our results can be taken as conservative lower limits.
We present the intrinsic merger rate densities and observer frame merger rates of all three types of DCOs in Figures 3, 4, 5, 6, and 7. Figs. 8 and 9 show the BH-BH merger rate densities as a function of the total masses of the systems. The results acquired in this study are available online at www.syntheticuniverse.org.
2 Stellar populations
In this section we describe the properties of stellar populations, and their evolution with redshift. The formalism is mostly adopted from Belczynski et al. (2010c).
2.1 Star Formation History
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):
where is the age of the Universe (Gyr) as measured in the rest frame, is the present age of the Universe ( Gyr, see Section 4) and the parameters have values: , , and . The SFR described above is expressed in comoving units of length and time.
2.2 Galaxy Mass Distribution
For redshifts we describe the distribution of galaxy masses using a Schechter-type probability density function, calibrated to observations (Fontana et al., 2006):
where , (), and . A galaxy mass is drawn from this distribution in solar units () and in the range . Beyond redshift we assume no further evolution in galaxy mass, fixing the mass distribution to the value at . This assumption reflects the lack of information on galaxy mass distribution at high redshift.
2.3 Galaxy Metallicity
We assume the average oxygen to hydrogen number ratio () in a typical galaxy to be given by
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 varying with redshift. We describe the redshift dependence of galaxy metallicity using the average metallicity relation from Pei et al. (1999):
which implies the evolution of with redshift:
We assume that the oxygen abundance (used in ) correlates linearly with the average abundance of elements heavier than Helium (encapsulated in the metallicity measure, ).
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 Eq. 3. In the first, provided by Pei et al. (1999), the coefficients are: , , , , . This grants a rate of average metallicity evolution, which we label slow. The second, provided by Young & Fryer (2007), uses , , , , . It is based on ultraviolet-GALEX, SDSS, infrared-Spitzer and neutrino-Super Kamiokande 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 (Eq. 2) we have two metallicity values (Eq. 3). (2.) We then combine these (slow and fast) metallicities into a single value being 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; of stellar mass) metallicities at redshift . (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 ). 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 . This grants a median value of metallicity of at (see Fig. 1), which corresponds to in the ”12+log(O/H)” formalism. This calibration was designed to match the upper scatter of metallicities according to Yuan et al. (2013) (see their Fig. 2, top-right panel). We label this profile as high–end, as it is the upper limit on metallicity at . 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 has solar metallicity, while the other half has solar metallicity. This yields a median metallicity value of and requires the division of the initial profile by a factor of . We label this profile as low–end.
2.4 Galaxy Stellar Populations
We distinguish three stellar populations:
We choose 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 : 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 () number of galaxies (Eq. 2) via Monte Carlo simulations.
3 Binary Star Modeling
We present our calculations for a set of 4 models, each differing in major input physics (see Table 1 and the subsequent sections). For each model we use a grid of 11 metallicity values (, (solar,), , , , , , , , , ) in order to accurately account for the average metallicity evolution of the stellar populations with redshift.
3.1 The StarTrack code
To calculate the evolution of the stellar populations we utilize the recently updated (Belczynski et al., 2010a, 2012b; 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 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 and . In our calculations we fix this value to be or in mathematical terms:
where is the accretion rate, is the mass transfer rate from the donor and is the fraction of the rate transferred (here equal to ). The remaining mass is expelled to infinity. The dynamically unstable mass transfers (common envelope) 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 to investigate the influence the rotation of the components has on their internal structure.
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 yr. This value was calibrated to account for the highest mass black holes in the Milky Way (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 utilizes the convection driven, neutrino–enhanced supernovae engines (Fryer et al., 2012) to determine the properties of the remnant objects (neutron stars and black holes).
For each metallicity value in each model we evolve binaries, assuming that each component is created at the same time. Each binary system is initialized by four parameters which are assumed to be independent. These are: primary mass, (initially more massive component), mass ratio, , where is the mass of the secondary component (initially less massive), the semi-major axis, , of the orbit, and the eccentricity, . The mass of the primary component is randomly chosen from the initial mass function adopted from Kroupa et al. (1993) and Kroupa & Weidner (2003):
where is our standard choice for field populations. The choice of the upper IMF cutoff () 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 neutron stars require their evolution within the range –, while for single BHs the lower limit is . Binary evolution may broaden these ranges due to mass transfer episodes, and we therefore set the minimum mass of the primary to . We assume a flat mass ratio distribution, , over the range –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 . However, for the same reasons as for the primary, we don’t consider binaries where the mass of the secondary is below . The distribution of initial binary separations is assumed to be flat in (Abt, 1983), and so , with ranging from values such that at ZAMS the primary fills no more than 50% of its Roche lobe to . For the initial eccentricity we adopt a thermal equilibrium distribution (e.g. Heggie, 1975; Duquennoy & Mayor, 1991): , with ranging from to . We find that the adopted parameters are in accordance with the most recent observations of O-star populations (Sana et al., 2013).
3.2 The model suite
3.2.1 The Standard Model
In this subsection we define a reference model for this paper. This model is identical with 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) common envelope (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 (–).
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 (P. Eggleton, 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 a HG star, regardless of the energy balance (labeled as “Submodel B”) or allow 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 ( yr; Belczynski & Dominik, 2012), exceeding even the empirically estimated rates based on IC10 X-1 and NGC 300 X-1 ( yr; Bulik et al., 2011; Belczynski et al., 2012a). 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 leads to the formation of DCOs.
Observations suggest (Hobbs et al., 2005) that neutron stars formed in supernovae receive natal kicks, with velocities drawn from a Maxwellian distribution with . We employ these findings in our simulations, and extend them so that black hole 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, we modify the Maxwellian kicks by the amount of matter falling back on the newly formed compact object:
where is the final magnitude of the natal kick, is the velocity drawn from a Maxwellian kick distribution, and is the fallback factor describing the fraction of the ejecta returning to the object. The values of 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” one, 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 supernova 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 –. For low mass stars () the result is a very strong (high velocity kick) supernova, which generates a NS. For higher mass stars a BH is formed through a direct collapse (failed supernova). This engine, incorporated into binary evolution, successfully reproduces the mass gap (Belczynski et al., 2012b) observed in Galactic X-ray binaries (Bailyn et al., 1998; Özel et al., 2010).
3.2.2 Variations on the standard model
The uncertainties in the CE and the SN engine argue for exploring a range of input physics beyond that in the standard model described in the previous subsection. In this subsection we present three additional models which we have found to encapsulate the full range of possible binary evolutions. All subsequent models use the same input physics as the reference model, except for the parameters described below.
Optimistic Common Envelope. In this model we allow HG stars to be CE donors (see Section 3.2.1). When the donor initiates the CE phase the energy balance determines the outcome. This model is identical to the “standard model – submodel A” from our previous paper in this series (Dominik et al., 2012).
Delayed SN. This model utilizes the “Delayed” supernova engine instead of the Rapid one. The Delayed is also a convection driven, neutrino enhanced engine, but is sourced from the standing accretion shock instability (SASI), and can produce an explosion as late as after bounce. The Delayed engine produces a continuous mass spectrum of compact objects, from NSs, through light BHs, to massive BHs (see Belczynski et al. (2012b)). This model is identical to the “Variation 10 – submodel B” model from our previous paper in this series (Dominik et al., 2012).
High BH kicks. In this model the BHs receive full natal kicks. The newly formed BH acquires a velocity drawn from a Maxwellian distribution (see Section 3.2.1) regardless of the fallback factor (see Eq. 9). This model is identical to the “Variation 8 – submodel B” model in our previous paper in this series (Dominik et al., 2012).
4 Cosmology calculations
We utilize a flat cosmology with , , , and . The relationship between redshift and time is given by:
where Gyr is the Hubble time (e.g. Hogg, 1999) and .
The comoving volume element is given by:
where is the speed of light in vacuum, is the solid angle, and is the comoving distance given by :
There are a series of steps to calculate the rates of events, as we now describe. We employ time as our reference coordinate and start by creating time bins across the entire history of Universe, each bin Myrs wide, from Gyrs (birth) to Gyrs (today). At the center of each bin we evaluate the star formation rate according to Eq. 1. For the redshift value corresponding to the center of a given time bin we generate a Monte Carlo sample of galaxies (a number sufficient to produce a smooth distribution) with masses drawn from the distribution given in Eq. 2. For each time bin we obtain a total mass of galaxies . For each galaxy we then estimate its average metallicity using Eq. 3. We assume that all stars within a given galaxy have identical metallicity as obtained from Eq. 3. Since we draw a large number of galaxies in each time bin, and each galaxy has its own mass, and therefore is described by its own average metallicity, we end up with a distribution of metallicity in each time bin. This also yields a total mass of galaxies with a specific metallicity () within each time bin. We then define the fraction of the total galaxy mass capable of forming stellar population with a specific metallicity by
However, because we use a finite number of metallicity points in our simulations (see Section 3) we need to extrapolate our results in order to account for the continuous spectrum given by Eq. 3. Therefore, the metallicity points are extended into bins delineated by the average value of neighbouring points. For example, given the set of points , the value now corresponds to a bin that extends from to . The border points and extend to lower and higher values, respectively, to cover the rest of the spectrum.
Population synthesis provides us with a representative sample of DCOs. The formation of a single DCO within a time bin corresponds to a fraction, , of the total formation rate:
where is the total mass in our population synthesis calculations (see Section 3). Repeating this calculation of for each metallicity yields a total formation rate, , within a given time bin.
We now need to know the delay time until merger, , for each DCO formed. The delay time is defined as the interval between the formation of the progenitors of a DCO and the coalescence of two compact objects. For each DCO originating from a specific metallicity we randomly choose a birth point, (ZAMS), within each time bin. We then propagate the system forward in time toward its merger using the delay time:
As long as we consider DCOs with , and as long as the width of the time bins throughout the time line is constant, the formation rate () of a DCO from a given bin translates into a merger rate in a later bin, propagated forward in time by . Repeating the above calculations for every time bin yields a total density of rest frame merger events, , in units of Gpc. In other words
where sums over each representative DCOs.
We now provide results from our four models, presenting the intrinsic merger rate densities and the observer frame merger rates, given by
with given by Eq. 11, integrated over the solid angle (hence the factor of ). In the case of the standard/reference model (details in Section 3.2.1) we explain the general redshift behavior of all three types of DCOs (NS-NS, BH-NS, and BH-BH) and compare the reference model for two scenarios of metallicity evolution (high–end and low–end). For our three variations (Optimistic CE, Delayed SN, and High BH kicks, Section 3.2.2) we investigate deviations from the reference model, again incorporating our different metallicity evolution scenarios.
5.1 Standard Model
NS-NS. As shown in Fig. 3 the intrinsic merger rate densities of double neutron star systems peaks at redshift ( yrGpc). As a general rule, the merger rates of all types of DCO are directly related to the star formation rate. However, for a given SFR value the formation efficiency of different DCO 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 Fig. 4). When combined with the peak of the SFR at (average high–end metallicity is , see Fig. 2), high metallicity NS-NS formation efficiency is enhanced, thus creating the profile shown on Fig. 3. What is characteristic for this profile is the ”hump” that arises at , approaching from high redshifts. As can be seen in Fig. 4, this shape is dominated by mergers originating from environments. The reason for this increase in merger rate densities, when transiting from 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 – boundary, allowing for the survival of a larger number of NS-NS progenitors.
By comparison, for the low–end metallicity profile the NS-NS systems dominate the rates only up to (Fig. 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.
BH-NS. For the high–end metallicity evolution, the rest frame merger rate densities for BH-NS systems shown in Fig. 3 peaks at a value of at redshift . However, the merger rate efficiency drops for low () and high () redshifts. This is because of properties of the progenitor masses. For metallicities the bulk of the progenitors masses are in the range – for the primary component, and – 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 lower than 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 are: – for the primary and – for the secondary. Given that in this mass range the Initial Mass Function (IMF) scales as (where 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 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 (), which leaves a narrow mass range: – for the primary and – for the secondary. There are fewer progenitors in these mass ranges when compared to 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 Fig. 4).
BH-BH. For these systems the intrinsic merger rate has a peak-plateau at a rate of – at –, for the high–end case. The low metallicity galaxies abundant at high redshifts are efficient black hole factories (see Fig. 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 to 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 and a BH pair. These systems originate from the extremely low metallicity environments (). We find that such systems merge up up until redshift and 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 a – BH-BH, for – a typical system for the average metallicity acquired in our study: t=0 Myr. The components start with masses and for the primary and secondary, respectively and an orbital separation . t=6.7 Myr. The primary, after becoming a 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 (the secondary component has ). The orbital separation prior to RLOF was and after. t=7.0 Myr. The primary explodes as a supernova, forming a BH. The orbital separation after the explosion was t=8.7 Myr. The secondary () initiates a CE phase and becomes a Wolf-Rayet star with as a result of the outcome (the primary gained ) The orbital separation prior to the CE was and after. t=9.4 Myr. The secondary undergoes a SN explosion and becomes a BH. The orbital separation prior to the explosion was and after. t=26 Myr. The coalescence of a – BH-BH system occurs. This example is illustrated by a diagram in Fig. 10.
On a side note, the formation of the most massive BH-BH systems on close orbits may be questionable. The progenitors of the aforementioned - BH-BH system are massive stars (– at ZAMS). A recent theoretical study by Yusof et al. (2013) suggests that such objects (–) 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 Hubble time.
In the observer frame, BH-BH systems begin to dominate the merger rates at . For the low–end case this happens closer to .
5.2 Optimistic CE
In this model we relax one of the conditions on CE survivability. Specifically, Hertzsprung gap 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, high–end case, the intrinsic merger rate of NS-NS systems peaks at higher redshift () and at higher values ( Gpc yr). The shift in the peak towards 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 towards .
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 . In the latter case they are always sub-dominant compared to BH-BH systems.
BH-NS. The binaries forming these systems usually undergo two CE events in their lifetime, 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 mentioned 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 towards lower redshift (, high–end) when contrasted with the reference model (see Fig. 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 Figs. 8 and 9). For the low–end case this results in a peak-plateau between redshifts . 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 in the high–end case. For the low–end case these DCOs are always primary mergers.
5.3 Delayed SN
In this model we change the supernova explosion engine with respect to the standard model. The standard model uses the Rapid engine, which yields a gap between – 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., 2012b). 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 on Fig. 8 and Fig. 9. The minimal total mass for this model is . Such a system is composed of two BHs of each ( being the delineation between upper NS and lower BH mass). For other models the minimal BH mass is , thus yielding a minimal total mass . However, the supernova engine effects do not play a significant role on the merger rates of any type of DCOs.
5.4 High BH kicks
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 Eq. 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 Fig. 3. In 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 would collapse directly into a black hole after the SN explosion; with no asymmetric ejecta, they do not receive a kick (, Eq. 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 on the bottom panel of Figs. 8 and 9.
6 Summary & Discussion
We have performed a series of cosmological calculations for four populations of DCOs. Each population was generated with different input physics for describing binary evolution and compact object formation. The first model (standard) utilizes the current state-of-the-art description of physical mechanisms governing DCOs. In particular, it uses a Rapid explosion engine, which yields results accurately describing the mass distribution of X-ray binaries (see Section 3.2.1 and references therein). Another major improvement in the model is the realistic treatment of the common envelope parameter , which now depends on the evolutionary stage, radius, mass, metallicity, etc. of the donor star. The three subsequent models explore alternative outcomes of binary evolution, and the resulting properties of remnants. The mechanisms investigated in these models are: the sensitivity of the CE outcome to the type of donor, the Delayed SN explosion mechanism, and the natal kick survivability of DCOs containing BHs (see Section 3.2.2). Additionally, for each model we have created a grid of 11 metallicities to account for the chemical evolution throughout the lifetime of the Universe. We present both the intrinsic and the observer frame merger rates as a function of redshift.
The variation in the rates of our different binary systems as a function of redshift depends upon metallicity, as well as common envelope and supernova physics. In this paper we have studied how these impact the rates for different types of DCOs. Here we review our main findings.
We find that NS-NS systems merge most efficiently at low redshifts (; see Figs. 3 and 5), where metallicities become relatively high (). However, in the case of the Optimistic CE model the merger rate densities peak at higher redshifts (–). This results from relaxing the condition for the termination of binaries initiating a CE with a Hertzsprung gap donor. This optimistic CE treatment enriches the merging population with systems with short merger times. As a result the overall number of NS-NS systems increases and, due to the shorter merger times, these systems coalesce earlier (see Sections 3.2.2 and 5.2).
BH-NS systems merge most infrequently in all but one of the models. The exception is the Full BH Kicks model, where full natal kicks are applied to BH remnants. The kicks eliminate binaries containing BHs from the populations by disrupting them. However, this doesn’t affect BH-NS systems, as strongly as BH-BH systems since they contain only one BH. In general, the low merger rates of BH-NS systems arise from their unique mixed nature. Forming two different compact objects in a single binary generally requires the masses of the progenitors to be significantly separated. This plays an important role at first contact between the components, since if the mass ratio of the binary is larger than – the otherwise stable mass transfer through Roche lobe overflow may become a CE event. These episodes often cause a premature merger and eliminate further binary evolution. Another important factor in making the BH-NS systems small in numbers is that the progenitors don’t have a large range of masses to draw from. The upper limit is set by the binary containing enough mass to form a BH-BH system instead, while the lower limit is set by not having enough mass and instead forming a NS-NS system.
For BH-BH systems, the highest merging efficiency occurs earlier in the Universe when compared with other DCOs (–). This arises from the fact that these systems form most efficiently at the lowest metallicities. For any of the two scenarios of metallicity evolution, the Optimistic CE model blurs this trend. In this case the population is enriched by BHs, which originated from high metallicity environments (see Section 5.2). Another interesting case is the model with High BH Kicks, where BH-BH systems are efficiently disrupted by natal kicks throughout the lifetime of the Universe. This is clearly visible on the bottom panel of Figs. 3 and 5. The kicks affect high mass systems the most. As a consequence of the full natal kicks, the formation and merger rates for BH-BH systems in low metallicity galaxies (higher redshifts) are reduced significantly, and this effect is even more dramatic for high metallicity environments (lower redshifts; see Figs. 8 and 9). The High BH kicks model produces a difference between the merger rates in the observer frame of BH-BH and NS-NS systems that is roughly 100 times larger than within the standard model. This may be a promising avenue for determining the magnitude of the natal kicks imparted to BHs during their formation.
Since (only) NS-NS systems have been observed, we can use observed rates to put constraints on our models. The NS-NS merger rates in each of our models, at , fit within the observational limits for NS-NS systems in the Milky Way: – (Kim et al., 2006), using the galaxy density . Petrillo et al. (2013) used the observed rate of short GRBs to calculate the merger rates of NS-NS and BH-NS systems, since these systems are thought to be the progenitors of short GRBs. The resulting merger rates of DCOs (NS-NS + BH-NS) in the local Universe ranges between and . At our models find a NS-NS merger rate of , with a BH-NS rate lower by a factor of . However, the authors of the aforementioned study state that their results are sensitive primarily to the poorly constrained beaming angle of the colimated emission from short GRBs. They used a beaming angle of deg, while to match our rate the beaming angle would have to be deg (see Fig. 3 therein). In our previous study (Dominik et al., 2012), we found one model that would reproduce the merger rates of NS-NS + BH-NS from Petrillo et al. (2013) ( at ). It is the model described by fully conservative mass transfer episodes and optimistic CE description (labeled ”Variation 12 – submodel A”).
Additional constraints may be provided by observing the potential electromagnetic signatures, other than GRBs, of DCO mergers. One example is the optical/radio afterglow of the GRB, which can be detected even if the GRB itself is not seen (an “orphan afterglow”). Another possibility is a “kilonova”, resulting from the ejection of matter from a neutron star. Since this matter may be enriched in heavy elements through the r–process, the resulting radioactive decay may generate observable light, thereby providing a promising electromagnetic counterpart to the gravitational wave emission (Metzger & Berger, 2012; Piran et al., 2013; Barnes & Kasen, 2013).
- Abt (1983) Abt, H. A. 1983, ARA&A, 21, 343
- Bailyn et al. (1998) Bailyn, C. D., Jain, R. K., Coppi, P., & Orosz, J. A. 1998, ApJ, 499, 367
- Barnes & Kasen (2013) Barnes, J. & Kasen, D. 2013, ArXiv e-prints
- Beers & Christlieb (2005) Beers, T. C. & Christlieb, N. 2005, ARA&A, 43, 531
- Belczynski et al. (2010a) Belczynski, K., Bulik, T., Fryer, C. L., et al. 2010a, ApJ, 714, 1217
- Belczynski et al. (2012a) Belczynski, K., Bulik, T., Mandel, I., et al. 2012a, ArXiv e-prints
- Belczynski & Dominik (2012) Belczynski, K. & Dominik, M. 2012, ArXiv e-prints
- Belczynski et al. (2010b) Belczynski, K., Dominik, M., Bulik, T., et al. 2010b, ApJ, 715, L138
- Belczynski et al. (2010c) Belczynski, K., Holz, D. E., Fryer, C. L., et al. 2010c, ApJ, 708, 117
- Belczynski et al. (2002) Belczynski, K., Kalogera, V., & Bulik, T. 2002, ApJ, 572, 407
- Belczynski et al. (2008) Belczynski, K., Kalogera, V., Rasio, F. A., et al. 2008, ApJS, 174, 223
- Belczynski et al. (2007) Belczynski, K., Taam, R. E., Kalogera, V., Rasio, F. A., & Bulik, T. 2007, ApJ, 662, 504
- Belczynski et al. (2012b) Belczynski, K., Wiktorowicz, G., Fryer, C. L., Holz, D. E., & Kalogera, V. 2012b, ApJ, 757, 91
- Binney & Merrifield (1998) Binney, J. & Merrifield, M. 1998, Galactic Astronomy
- Bulik et al. (2011) Bulik, T., Belczynski, K., & Prestwich, A. 2011, ApJ, 730, 140
- Bulik et al. (2004) Bulik, T., Belczyński, K., & Rudak, B. 2004, A&A, 415, 407
- Dewi & Pols (2003) Dewi, J. D. M. & Pols, O. R. 2003, MNRAS, 344, 629
- Dominik et al. (2012) Dominik, M., Belczynski, K., Fryer, C., et al. 2012, ApJ, 759, 52
- Duquennoy & Mayor (1991) Duquennoy, A. & Mayor, M. 1991, A&A, 248, 485
- Erb et al. (2006) Erb, D. K., Shapley, A. E., Pettini, M., et al. 2006, ApJ, 644, 813
- Figer (2005) Figer, D. F. 2005, Nature, 434, 192
- Fontana et al. (2006) Fontana, A., Salimbeni, S., Grazian, A., et al. 2006, A&A, 459, 745
- Fryer et al. (2012) Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, ApJ, 749, 91
- Gerosa et al. (2013) Gerosa, D., Kesden, M., Berti, E., O’Shaughnessy, R., & Sperhake, U. 2013, Phys. Rev. D, 87, 104028
- Heggie (1975) Heggie, D. C. 1975, MNRAS, 173, 729
- Hobbs et al. (2005) Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974
- Hogg (1999) Hogg, D. W. 1999, ArXiv Astrophysics e-prints
- Hopkins & Beacom (2006) Hopkins, A. M. & Beacom, J. F. 2006, ApJ, 651, 142
- Hurley et al. (2000) Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543
- Hut (1981) Hut, P. 1981, A&A, 99, 126
- Ivanova & Taam (2003) Ivanova, N. & Taam, R. E. 2003, ApJ, 599, 516
- Kim et al. (2006) Kim, C., Kalogera, V., & Lorimer, D. R. 2006, ArXiv Astrophysics e-prints
- Kobulnicky & Fryer (2007) Kobulnicky, H. A. & Fryer, C. L. 2007, ApJ, 670, 747
- Kreidberg et al. (2012) Kreidberg, L., Bailyn, C. D., Farr, W. M., & Kalogera, V. 2012, ApJ, 757, 36
- Kroupa et al. (1993) Kroupa, P., Tout, C. A., & Gilmore, G. 1993, MNRAS, 262, 545
- Kroupa & Weidner (2003) Kroupa, P. & Weidner, C. 2003, ApJ, 598, 1076
- Mackey et al. (2003) Mackey, J., Bromm, V., & Hernquist, L. 2003, ApJ, 586, 1
- Mandel (2010) Mandel, I. 2010, Phys. Rev. D, 81, 084029
- Metzger & Berger (2012) Metzger, B. D. & Berger, E. 2012, ApJ, 746, 48
- Mirabel & Rodrigues (2003) Mirabel, I. F. & Rodrigues, I. 2003, Science, 300, 1119
- Morscher et al. (2012) Morscher, M., Umbreit, S., Farr, W. M., & Rasio, F. A. 2012, ArXiv e-prints
- Oey & Clarke (2005) Oey, M. S. & Clarke, C. J. 2005, ApJ, 620, L43
- O’Shaughnessy (2012) O’Shaughnessy, R. 2012, ArXiv e-prints
- Özel et al. (2010) Özel, F., Psaltis, D., Narayan, R., & McClintock, J. E. 2010, ApJ, 725, 1918
- Panter et al. (2008) Panter, B., Jimenez, R., Heavens, A. F., & Charlot, S. 2008, MNRAS, 391, 1117
- Pei et al. (1999) Pei, Y. C., Fall, S. M., & Hauser, M. G. 1999, ApJ, 522, 604
- Petrillo et al. (2013) Petrillo, E., Dietz, A., & Cavaglià, M. 2013, ApJ, 767, 140
- Piran et al. (2013) Piran, T., Nakar, E., & Rosswog, S. 2013, MNRAS, 430, 2121
- Regimbau & Hughes (2009) Regimbau, T. & Hughes, S. A. 2009, Phys. Rev. D, 79, 062002
- Sana et al. (2013) Sana, H., de Koter, A., de Mink, S. E., et al. 2013, A&A, 550, A107
- Smith et al. (2009) Smith, B. D., Turk, M. J., Sigurdsson, S., O’Shea, B. W., & Norman, M. L. 2009, ApJ, 691, 441
- Strolger et al. (2004) Strolger, L.-G., Riess, A. G., Dahlen, T., et al. 2004, ApJ, 613, 200
- Tauris & Savonije (1999) Tauris, T. M. & Savonije, G. J. 1999, A&A, 350, 928
- Totani (1999) Totani, T. 1999, ApJ, 511, 41
- Webbink (1984) Webbink, R. F. 1984, ApJ, 277, 355
- Wellstein et al. (2001) Wellstein, S., Langer, N., & Braun, H. 2001, A&A, 369, 939
- Xu & Li (2010) Xu, X.-J. & Li, X.-D. 2010, ApJ, 716, 114
- Young & Fryer (2007) Young, P. A. & Fryer, C. L. 2007, ApJ, 670, 584
- Yuan et al. (2013) Yuan, T.-T., Kewley, L. J., & Richard, J. 2013, ApJ, 763, 9
- Yusof et al. (2013) Yusof, N., Hirschi, R., Meynet, G., et al. 2013, MNRAS, 433, 1114
- Zahn (1977) Zahn, J.-P. 1977, A&A, 57, 383
- Zahn (1989) —. 1989, A&A, 220, 112
|Standard||Nanjing/physical, BH kicks: decreased, SN: Rapid|
|HG CE donors not allowed|
|Optimistic CE||HG CE donors allowed|
|Delayed SN||Delayed supernova engine|
|High BH kicks||Full kicks of BHs|