BINARY COMPACT OBJECT COALESCENCE RATES: THE ROLE OF ELLIPTICAL GALAXIES

R. O'Shaughnessy; V. Kalogera; Krzysztof Belczynski

doi:10.1088/0004-637X/716/1/615

1. INTRODUCTION

Over the last decade the question of the Galactic inspiral rate of binaries with two compact objects (neutron stars (NSs) or black holes (BHs)) has attracted attention primarily because of the development and planning of gravitational-wave interferometric detectors both on the Earth and in space (e.g., LIGO and GEO600, described in Abbott et al. 2004; Virgo, described at the Virgo project Web site www.virgo.infn.it and in Acernese et al. 2006; and LISA, at http://lisa.nasa.gov). These rate estimates have been widely used in the assessment of gravitational inspiral detectability, given assumed instrument sensitivities. A number of different groups have calculated inspiral rates using population synthesis calculations, most commonly with Monte Carlo methods (Fryer et al. 1998, 1999; Portegies Zwart & Yungelson 1998; Brown & Bethe 1994; Belczynski et al. 2002; Voss & Tauris 2003). Such studies consider the complete formation history of double compact objects through long sequences of binary evolution phases, terminated by a gravitational-wave-driven inspiral toward merger. Though less critical for studies of the Milky Way, contributions to the present-day merger rate per unit volume from early star formation have most often been ignored in merger rate calculations. Since long inspiral delays before merger are not sufficiently uncommon, even for double NSs, merger detection rate calculations should account for the full time history of star formation, including star formation in the early universe (see, for example, de Freitas Pacheco et al. 2006). Additionally, our current understanding of single and binary star evolution is incomplete. The many uncertainties have been parameterized and the resulting parameter space explored in some studies to determine the range of plausible results (see, e.g., Belczynski et al. 2008b), both by focusing narrowly on physically motivated regions of parameter space as in Belczynski et al. (2008b) and by broadly exploring all plausible population synthesis simulations (see, e.g., O'Shaughnessy et al. 2008a, henceforth denoted PS-GRB, and references therein). Finally, by assigning equal prior likelihood to any population synthesis model, PS-GRB estimated the relative likelihood of any BH–NS or NS–NS merger rate.⁶

The population of binary BHs, however, behaves qualitatively differently from BH–NS or NS–NS binaries: (1) most of those that merge do so only after a few to several Gyr delay after their birth as binary stars and (2) those BH–BH binaries with merger times from 1to10 Gyr have masses quite different from the typical assumption of a 10 M_☉ BH (see Appendix A). The rarity, long delay times, high masses, and thus cosmologically significant detection horizons of binary BHs have already been discussed in the literature (see, e.g., Kulczycki et al. 2006, and references therein). However, population synthesis simulations extensive enough to contain a statistically meaningful binary BH population come at a significant computational cost. Even with adequate population synthesis simulations, given the long-lived nature of the progenitors, BH–BH merger and detection rates depend critically on the earliest and least certain star formation rates (SFRs). Certain exotic but not exceptionally uncommon binary evolution scenarios could lead to a "high-rate tail"—a scenario with rare massive mergers, for example, that could plausibly produce a merger detection in the near future. In this paper, we strive to understand what realistic scenarios could lead to a high-rate tail that LIGO could detect or, conversely, constrain.

Given the technical challenges and significant uncertainties involved, only one paper (O'Shaughnessy et al. 2005) has previously estimated the relative likelihood of different BH–BH detection rate estimates expected from isolated binary star evolution, only for the Milky Way. In the present study, we calculate the expected detection rate for LIGO while for the first time simultaneously including (1) all past star formation, particularly the overwhelming importance of elliptical galaxies, (2) a large database of results that accounts for the dominant model uncertainties, and (3) a careful treatment of the mass distribution of BH–BH binaries. Additionally, unlike previous analyses presenting distributions of merger rates, we have post-facto varied the birth binary fraction of stars from 100% to 15%; this addition allows us to better compare our results with the fiducial population synthesis model and BH–BH merger rates presented in Belczynski et al. (2007). Finally, whereas previous comparisons relied essentially on confidence intervals (e.g., PS-GRB; O'Shaughnessy et al. 2008b, henceforth denoted PSC2), in this paper we generate true posterior distributions: each model is weighted by its relative conditional likelihood given observations of merging NS–NS binaries in the Milky Way.

To fully assess the total probability for a LIGO detection, rather than limit attention to BH–BH mergers we explore the rate at which all LIGO-detectable binaries (BH–BH, BH–NS, and NS–NS; henceforth collectively denoted as compact binary coalescences or CBCs) merge through gravitational-wave emission. As the methods used to determine LIGO detection rates from population synthesis simulations are already extensively discussed in the literature (see, e.g., Bulik et al. 2004; PS-GRB), our presentation only reviews those tools, emphasizing unique features of our present analysis (i.e., our methods to estimate BH–BH masses). Specifically, in Section 2 we briefly review the ingredients that enter into an estimate of the gravitational-wave detection rate for a short-ranged network (z_max ≪ 1). While the local universe is emphasized in the text, networks of advanced ground-based interferometers can detect optimally oriented BH–BH binaries with component masses M ≃ 10–15 M_☉ at cosmologically significant distances z ⩾ 0.1. Merger rates on the past light cone of a detector become inhomogeneous (versus redshift) at this scale, simply because the SFR increases dramatically near z ≃ 1–2 during the epoch of galaxy assembly. At these distances, the detection rate must be integrated over the mass distribution, the full networked orientation-dependent sensitivity, and redshift-dependent merger rate. This integral becomes increasingly sensitive to rare, high-mass events that can be seen in the ancient universe. Therefore, we generalize the short-range method described in Section 2 to include all sources on the past light cone, to arbitrary redshift and with accurate orientation-dependent sensitivity. In Section 3, we review the model for star formation in the universe adopted here. This experience is applied in Section 4 (for BH–BH binaries) and Section 5 (including NS–NS and BH–NS binaries too), where we review our population synthesis calculations and explain what features of those calculations influence our predictions for the relative likelihood of a LIGO detection. These sections also describe Figure 7, which contains our predictions regarding initial and advanced LIGO detection rates. Finally, in Section 5.2 we show how rate constraints from the observed sample of Milky Way NS–NS binaries affect the predictions for LIGO detection rates.

To summarize, by fully simulating the past history of the local universe, this paper develops models for the present-day detection rate of short-range (z ≪ 1) and long-range (z not much less than 1) gravitational-wave detectors. Our results are rate distributions, where each distribution includes some normalization uncertainties (SFR and fraction of stars born in binaries), certain population synthesis model parameters, and our simulation Monte Carlo uncertainty. Further, our Bayesian approach to model constraints surveys many key uncertainties that should be included, whatever the model family involved, when attempting to interpret upper limits or detections from pulsar populations and gravitational-wave observatories. We believe our distributions represent the best predictions of a concrete, conservative model family that includes both elliptical and spiral star formation yet is also consistent with initial-LIGO upper limits on CBC merger rates.

2. GRAVITATIONAL-WAVE DETECTION

This section primarily reviews the final stage in any calculation of a LIGO detection rate: the connection between, on the one hand, the event rate per comoving volume ${\cal R}$ and mass distribution p(m₁, m₂) of merging binaries and, on the other hand, the detection rate R_D (Equation (9)). This relation depends critically on the range to which LIGO could observe each merger. A thorough discussion of the LIGO range requires careful review of data analysis strategies and interferometer network geometry, and remains substantially beyond the scope of this paper. For details see Abbott et al. (2008, 2009), and references therein. For example, with modern multi-mass-region searches using inspiral templates, one essential effect is a merger detection threshold in signal-to-noise ratio that changes with total binary mass, as the well-understood portions of detection templates grow shorter. In this paper, we adopt the customary set of approximations (described below) and ignore the errors they introduce.

The distance to which LIGO is sensitive can be calculated using the sensitivity of the detector (strictly, the strain noise spectral distribution S_h), the masses of the inspiralling bodies, and the emitted waveform. (For simplicity, we will assume the BHs initially nonspinning.) Specifically, for a source with masses m₁ and m₂, merging at a luminosity distance D (or equivalently at redshift z, or comoving distance r) with a binary inclination ι and a relative orientation to the detector given by θ, ϕ (the orientation on the sky) and ψ (a phase angle in the plane of the detector), the signal to noise seen by a matched-filter search in the LIGO data stream can be expressed relative to the strain noise spectrum of a single interferometer S_h and the Fourier transform $\tilde{h}(f)$ of gravitational-wave strain h(t) received at the detector by an arbitrarily oriented and located source by

$\begin{flequation} \rho ^2=4\frac{w(\theta,\phi,\psi,\iota)^2}{D^2}\!\! \int _0^{\infty } \!\! df \frac{|r \tilde{h}(m_1(1+z),m_2(1+z),f)|^2}{S_h(f)} \end{flequation} \tag{ 1 }$

using units with G = c = 1, the known dependence of planar-symmetry gravitational waves on angle, and a standard one-sided Fourier convention for $\tilde{h}$ and S_h versus frequency f. The function w(θ, ϕ, ψ, ι) takes values between 0 and 1 and completely encompasses the detector- and source-orientation-dependent sensitivity for the most common sources, those that are dominated by l = m = 2 quadrupole emission⁷

$\begin{eqnarray} w(\theta,\phi,\psi,\iota)&=& \sqrt{F_+^2 (1 + \cos ^{2} \iota)^2/4+ F_{\times }^2\cos ^{2} \iota } \end{eqnarray} \tag{ 2 }$

$\begin{eqnarray} F_\mathrm{+} &=& \frac{1}{2}(1+\cos ^{2} \theta) \cos 2 \phi \cos 2 \psi \nonumber \\ &&+ \cos \theta \sin 2 \phi \sin 2 \psi \, \end{eqnarray} \tag{ 3 }$

$\begin{eqnarray} F_\mathrm{\times } &=&{-} \frac{1}{2}(1+\cos ^{2} \theta) \sin 2 \phi \cos 2 \psi \nonumber \\ &&+ \cos \theta \cos 2 \phi \sin 2 \psi. \end{eqnarray} \tag{ 4 }$

To a rough approximation that depends on the data analysis strategy used and the amount of non-Gaussian noise present in the detector, a single LIGO interferometer can detect the gravitational-wave signature of a merging binary if ρ>ρ_c (see, e.g., Abbott et al. 2006, 2008, and references therein; henceforth we adopt ρ_c = 8). For this reason, the previous expression for the signal to noise (Equation (1)) is often re-expressed as

$\begin{eqnarray} \rho &=& \rho _cw(\theta,\phi,\psi,\iota) \frac{ {D_{\rm H}}(m_1(1+z),m_2(1+z))}{D} \end{eqnarray} \tag{ 5 }$

which, by comparison with Equation (1), implicitly defines the horizon distance D_H—the maximum luminosity distance to which the detectors are sensitive—as a function of the redshifted masses m₁(1 + z) and m₂(1 + z) of the binary (see, e.g., Anderson et al. 2001, Finn & Chernoff 1993, and Flanagan & Hughes 1998 for a brief review of the theory underlying the direction-dependent LIGO sensitivity as well as for expressions that approximate D_H(m₁, m₂) for low-mass CBC binaries⁸).

More commonly used is the volume-averaged distanceD_v, chosen so the volume of a sphere of radius D_v agrees with the average volume enclosed in a (Euclidean) detection surface:

$\begin{eqnarray} {D_{\rm v}}^3&=& \frac{1}{4\pi } {D_{\rm H}}^3 \int d\Omega \frac{d\psi }{\pi } \frac{d \cos \iota }{2} w^3 \nonumber \\ &\simeq & ({D_{\rm H}}/2.26)^3. \end{eqnarray} \tag{ 6 }$

The volume-averaged distance is a meaningful measure of sensitivity only when the horizon range is much smaller than the Hubble scale.

A network of detectors can coherently add signals from each interferometer to increase the signal-to-noise ρ associated with each signal and by implication its reach. To a first approximation, ignoring small differences in sensitivity due to both their intrinsic differences and their orientations,⁹ the combined signal-to-noise ratio of all interferometers is higher by roughly $\sqrt{\sum _k L_k^2/L_1^2}$ , where L_k is the length of the kth interferometer. For example, this factor is approximately $g_N\simeq \sqrt{1+1+1/4}\simeq \sqrt{2.25}$ for the initial LIGO network, consisting of two 4 km and one 2 km interferometer (see Cutler & Flanagan 1994 for details on realistic multidetector beampatterns). However, in part because a network performs more trials of the same data streams (e.g., for each sky position) and is sensitive to both polarizations simultaneously, for the same false alarm rate the network detection threshold ρ_c,net is generally greater than the threshold ρ_c for the individual detectors (see, for example, Cutler & Flanagan 1994 and Buonanno et al. 2003b for a discussion in the case of Gaussian noise). One can therefore speak of a single-interferometer (D_H) and network (D_H⁽ⁿ⁾ ≃ g_ND_H(ρ_c/ρ_c,net)) horizon distance. Unfortunately, only expert analysis of real data can determine the sensitivity of real gravitational-wave networks, not the least because noise statistics (and therefore ρ_c,net) are highly detector and even search dependent. To avoid ambiguity and misleading approximations, we provide results for single, idealized interferometers assuming a detection threshold of ρ_c = 8 (e.g., appropriate to Gaussian noise). This reference sensitivity agrees with the customary sensitivity measure used in the LIGO coalescence search, a threshold of ρ_c = 8 for a single detector as part of a network.¹⁰ For sources in the local universe, without loss of generality the reader can scale our results to any realistic network and search.

For sources of sufficiently low mass (M_tot ≲ 20 M_☉), LIGO is primarily sensitive to less-relativistic and therefore better-understood phases of the spiral-in, allowing us to approximate the waveform by its Newtonian limit (Peters 1964). After considerable algebra (see, e.g., Flanagan & Hughes 1998; Cutler & Flanagan 1994), the signal to noise due to an optimally oriented inspiral at comoving distance r (corresponding to a redshift z(r)) can be related to the strain sensitivity S_h by

$\begin{eqnarray} \rho ^2 &=& \frac{5 w^2 [ {{\cal M}_c}(1+z)]^{5/3}}{6 \pi ^{4/3} D^2} \int _0^{\infty } df \, \frac{{ f}^{-7/3}}{S_h(f)}, \end{eqnarray} \tag{ 7 }$

where ${{\cal M}_c}\equiv (m_1 m_2)^{3/5}/(m_1+m_2)^{1/5}$ is the "chirp mass" of the binary. Based on this expression, individual interferometers in the initial,¹¹ and advanced¹² LIGO networks can see low-mass binaries (m₁ + m₂ ≲ 40 M_☉)¹³ out to a horizon (volume-averaged) distance

$\begin{equation} D_{\rm H(v)} \simeq {\cal C}_{\rm H(v)} ({{\cal M}_c}/1.2\,M_\odot)^{5/6} \end{equation} \tag{ 8 }$

where ${\cal C}_{\rm H(v)}=31 (14) \, {\rm Mpc}$ and ${\cal C}_{\rm H(v)}=445 (197) \, {\rm Mpc}$ , respectively. For clarity, the distance ${\cal C}_v$ will henceforth be denoted D_bns, the volume-averaged distance to which a binary NS inspiral can be detected by a single interferometer as part of a network.

In terms of this horizon distance D_H, the intrinsic present-day rate of events per unit comoving volume ${\cal R}(0)$ , and the chirp mass distribution of merging binaries $p({{\cal M}_c})$ , the average rate of events with ρ>ρ_c occurring in the nearby universe can be expressed as

$\begin{eqnarray} R_D = {\cal R}(0)&&\int _{\frac{w {D_{\rm H}}}{D}>1} r^2 dr d\Omega \frac{d\psi }{\pi } \frac{d\cos \iota }{2} p({{\cal M}_c}) d {{\cal M}_c} \end{eqnarray} \tag{ 9 }$

$\begin{eqnarray} \hspace{-1pc}= {\cal R}(0) && \left[\int d {{\cal M}_c}\frac{4\pi }{3}{D_{\rm v}}({{\cal M}_c})^3 p({{\cal M}_c}) \right]\,, \end{eqnarray} \tag{ 10 }$

where the integral is taken over all detectable combinations of masses, orientation, and location. For low-mass binary mergers (M < 40 M_☉), for which Equation (8) applies, the chirp mass-weighted average simplifies to

$\begin{eqnarray} R_D&\simeq & {\cal R}(0) \frac{4\pi ({\cal C}_v)^3}{3} \int p({{\cal M}_c}) ({{\cal M}_c}/1.2\,M_\odot)^{15/6} d {{\cal M}_c}{,}\qquad \end{eqnarray} \tag{ 11 }$

an expression that has been applied extensively in the literature to translate event rates per unit volume into LIGO detection rates (see, e.g., O'Shaughnessy et al. 2008b, and references therein).

2.1. Rates for Advanced Detectors

As the product of the present-day merger rate and an effective detection volume, the detection rate estimate Equation (9) assumes a locally homogeneous and isotropic universe. Advanced detectors can see back to sky-position-dependent epochs with noticeably different star formation history (z>0.1). For advanced detectors, a more generic expression must be used, which integrates over the mass distribution, the redshift-dependent merger rate, and the orientation-dependent reach w.

For the purposes of this paper, we continue to assume the maximum luminosity distance to which a single interferometer can identify a merger in noise is well described by the inspiral phase (Equation (7)). To normalize our estimate, we adopt a single-interferometer advanced LIGO range of 445 Mpc. Based on the optimal range for a given chirp mass, each simulation's the chirp mass distributions, and the single-interferometer beampattern w (Equation (2)), we determine the fraction of all mergers at redshift z that could be detected, P_ok(z):

$\begin{eqnarray} P_{\rm ok}(z)&=& \int _{D(z)<w {D_{\rm H}}} \frac{d\Omega }{4\pi } \frac{d\psi }{\pi } \frac{d\cos \iota }{2} d {{\cal M}_c}p({{\cal M}_c}). \end{eqnarray} \tag{ 12 }$

The overall single-interferometer detection rate is therefore the sum over the past light cone of the (redshifted) rate of mergers on it, times the fraction P_ok of mergers that could be detected:

$\begin{eqnarray} R_D &=& \int dz \frac{dV}{dz} \frac{{\cal R}(t)}{1+z} P_{\rm ok}(z). \end{eqnarray} \tag{ 13 }$

For a maximum reach z ≪ 1, these two expressions are equivalent to Equation (7).

3. MULTICOMPONENT STAR FORMATION HISTORY

The LIGO detection rate depends on the mass distribution of merging binaries $p({{\cal M}_c})$ and on the average rate of mergers per unit comoving volume, ${\cal R}(t)$ . This rate, in turn, generally depends on the net contribution from multiple star-forming populations, a contribution that convolves the SFR in each component (the subject of the present section) with the rate at which each components' star-forming gas yields CBC mergers (the subject of Section 4). More specifically, the merger rate density ${\cal R}(t)$ , a quantity required to calculate the detection rate, is obtained from (1) the SFR in each component dρ_C/dt where C indexes the different star-forming types, C = e for ellipticals and s for spirals, defined as the mass per unit comoving volume and time that forms as stars; (2) the mass efficiency λ_C at which each type of CBC binaries form, defined as the total number of binaries that survive isolated evolution and form CBC binaries per unit star-forming mass; and (3) the probability distribution dP_mC/dt for merger events to occur after a delay time t after star formation, where dP_mC is the fraction of mergers occurring between t and t+dt after the progenitor binary forms at t = 0:

$\begin{eqnarray} {\cal R}(t) &=& \sum _{C}{\cal R}_{C}(t) \end{eqnarray} \tag{ 14 }$

$\begin{eqnarray} {\cal R}_{C}(t) &=& \int _{-T}^t dt_b \lambda _{ C}\frac{d\!P_{m C}}{dt}(t-t_b) \frac{d\rho _{ C}}{dt}(t_b), \end{eqnarray} \tag{ 15 }$

where the integration variable t_b is time at which λdρ binaries per unit volume are born. As first recognized by de Freitas Pacheco et al. (2006) and as demonstrated systematically in PS-GRB, a multicomponent star formation history that includes both elliptical and spiral galaxies must be applied when calculating present-day merger rates of double compact objects. Simply put, even though elliptical galaxies formed stars long ago, dP/dt decays slowly enough (roughly as 1/t; see, e.g., PS-GRB) that the high rate of star formation in the early universe can lead to a significant elliptical-galaxy-hosted merger rate density in the present-day universe. Because elliptical galaxies have distinctly different star-forming conditions (e.g., metallicities) than present-day star-forming galaxies, models for stellar evolution in the universe should not require identical behavior in present-day solar-metallicity galaxies and in young ellipticals. A full multicomponent treatment is necessary. The first multicomponent star formation history was applied to predict BH–NS and NS–NS merger rates in PS-GRB; this paper employs the same framework and preferred star formation history, discussed below.

We adopt the two-component star formation history model presented by Nagamine et al. (2006). This model consists of an early "elliptical" component and a fairly steady "spiral" component, with SFRs given by

$\begin{equation} \dot{\rho } =\dot{\rho }_e + \dot{\rho }_s, \end{equation} \tag{ 16 }$

$\begin{equation} \dot{\rho }_{C}= A_{C} (t/\tau _{C}) e^{-t/\tau _{C}}, \end{equation} \tag{ 17 }$

where cosmological time t is measured starting from the beginning of the universe and where the two components' decay timescales are τ_e,s= 1.5 and 4.5 Gyr, respectively (see Section 2 and Table 2 of Nagamine et al. 2006). These normalization constants A_e,s = 0.22, 0.06 M_☉ yr⁻¹ Mpc⁻³ were chosen by Nagamine et al. (2006) so the integrated amount of elliptical and spiral star formation reproduce the present-day census of baryonic matter in ellipticals and spirals, allowing for a certain fraction of gas recycling between different generations of stars.¹⁴ Figure 1 illustrates the SFRs assumed by the model adopted here.

**Figure 1.** Star formation history of the universe used in this paper vs. time, relative to the present day. Solid line: net star formation history implied by Equation (16). Dashed, dotted line: the star formation history due to elliptical and spiral components.
Download figure:
Standard image High-resolution image

Each component forms stars in its own distinctive conditions, set by comparison with observations of the Milky Way and elliptical galaxies. We assume mass converted into stars in the fairly steady "spiral" component with solar metallicity and with a fixed high-mass initial mass function (IMF) power law (p = −2.7 in the broken-power-law Kroupa IMF; see Kroupa & Weidner 2003). On the other hand, we assume stars born in the "elliptical" component are drawn from a broken-power-law IMF with high-mass index within p ∈ [ − 2.27, − 2.06] and metallicity Z within 0.56 < Z/Z_☉ < 1.5. These elliptical birth conditions agree with observations of both old ellipticals in the local universe (see Li et al. 2006 and references therein) as well as of young starburst clusters (see Fall et al. 2005; Zhang & Fall 1999, and references therein).

Since two independent stellar populations give birth to CBCs, the net detection rate R_D is the sum of the detection rate due to ellipticals and spirals. In particular, the probability distribution function for R_D is necessarily the convolution of the distribution functions for detection rates due to elliptical and spiral galaxies, respectively. Because of the many orders of magnitude uncertainty in merger rates and the need to smooth across relative errors, however, all our probability distribution functions (PDFs) are stored and smoothed in a logarithmic representation, as p(log R_D; see, for example, PS-GRB). For this reason, the convolution takes on an unusual form, involving an integral over the fraction f of the detection rate due to ellipticals (i.e., f = R_D,e/R_D). Using this variable to form the logarithmic convolution, the distribution for the net rate log R_D can be calculated from the distributions for ellipticals and spirals (p_e,s) as

$\begin{eqnarray} p(\log R_D) &=& \int _0^1\frac{df}{f(1-f)\ln (10)} p_e(\log R_D + \log f) \nonumber \\ &&\times p_s(\log R_D + \log (1-f))\\ &\equiv& p_s \oplus p_e. \end{eqnarray} \tag{ 18 }$

Limitations of two-component model. Our two-component model implicitly adopts two extremely strong assumptions: (1) that metallicity in each component does not significantly evolve with redshift and (2) that all star-forming gas in each component is identical. While sufficient for most problems involving binary evolution, these approximations may be particularly ill-suited to determining the progenitors of gravitational-wave detections. Very low metallicity environments could produce BHs with exceptionally high mass; also, the BH binary merger rate increases strongly with metallicity (see, e.g., Belczynski et al. 2010a). Because of the long typical delay between progenitor birth and BH–BH merger, very low metallicity environments in the early universe could contribute significantly to the present-day gravitational-wave detection rate (Belczynski et al. 2010b). Even at present, low-metallicity environments exist often enough to potentially dominate the local merger rate (see, e.g., O'Shaughnessy et al. 2008c). For example, the nearby BH–BH progenitor binary IC 10 X-1 both lies in a low-metallicity environment and suggests a high BH–BH detection rate for initial LIGO (0.5 yr, multiplied by a factor of order unity, strongly dependent on survey selection effects; see Bulik et al. 2008).

Despite strong theoretical and observational evidence for a significant low-metallicity contribution, we have neither the confidence in the necessary astrophysical inputs (i.e., a redshift-dependent metallicity distribution) nor the required database of simulations to include these metallicity effects at this time. Further, unless initial LIGO detects BH–BH mergers, the peak detection rates cannot be significantly larger than what our model family predicts (see Figure 7). For this reason, we continue to adopt the simple and conservative two-component mechanism outlined above.

Relation to blue light density. Previously, merger rate estimates for the Milky Way have been scaled up to the local universe, assuming mergers track blue light (see, e.g., Kopparapu et al. 2008 and references therein). As roughly speaking blue light traces star formation (see, e.g., Gallagheret al. 1984; Kennicutt 1998; O'Shaughnessy et al. 2008c), this extrapolation roughly agrees with our estimate for binary merger rates in spiral-like galaxies, with a suitable choice of normalization (i.e., SFR per unit blue light); see the conclusions for a detailed numerical comparison. O'Shaughnessy et al. (2008c) explains why blue light normalization was not adopted in this paper.

4. CBC POPULATION SYNTHESIS MODELS

As outlined in Sections 2 and 3 (see particularly Equations (9) and (15)), in order to estimate how often CBC binary mergers could be seen with gravitational-wave detectors we must know several features of the CBC population: (1) how many CBC binaries form from any given star-forming mass (characterized by the mass efficiency λ), (2) how long these CBC binaries last between birth as stars and merger through gravitational-wave emission (characterized by dP/dt), and (3) what masses these merging CBC binaries have (characterized by a probability distribution in chirp mass $p({{\cal M}_c})$ ). In the absence of observations, these properties must be obtained from theoretical models of stellar and binary evolution. We study the formation of compact objects with the StarTrack population synthesis code, first developed by Belczynski et al. (2002) and recently significantly extended as described in detail in Belczynski et al. (2008b). Since our understanding of the evolution of single and binary stars is incomplete, this code parameterizes several critical physical processes with a great many parameters (∼30), many of which influence compact object formation dramatically. In this specific study, in addition to the IMF and metallicity (which vary depending on whether a binary is born in an elliptical or spiral galaxy), seven parameters strongly influence compact object merger rates: the supernova kick distribution (modeled as the superposition of two independent Maxwellians, using three parameters: one parameter for the probability of drawing from each Maxwellian and the other to characterize the dispersion of each Maxwellian), the stellar wind strength, the common-envelope energy transfer efficiency, the fraction of angular momentum lost to infinity in phases of non-conservative mass transfer, and the relative distribution of masses in the binary. For this reason, following methods largely outlined in our previous work in PS-GRB and PSC2, we perform a broad parameter study to ensure all plausible population synthesis models have been explored. We include simulations of two different classes of star-forming conditions: "spiral" conditions, with Z = Z_☉ and a high-mass IMF slope of p = −2.7, and "elliptical" conditions, with a much flatter IMF slope p ≃ −2.3 and a range of allowed metallicities 0.56 < Z/Z_☉ < 1.5. Though we have performed thousands of simulations of these two types of conditions in the past (e.g., PS-GRB), to obtain better statistics on rare BH–BH events, we have constructed a new set of archives of N_S ≡ 282 and N_E ≡ 206 simulations of these two type of conditions, respectively Within each set, we select parameters randomly according to the distributions described above (see Appendix B for more details).

The model database described above is statistically very similar to the database PS-GRB previously used to estimate NS–NS and BH–NS merger rates, with one critical exception: the simulations used here are extensive enough to have a statistically significant number of merging BH–BH binaries. Merging BH–BH binaries are a small fraction of the total BH–BH population and an exceedingly rare consequence of isolated binary evolution. For example, the number (n) of BH–BH binaries in each of our simulations was always by construction 3000 and could be as high as few ×10⁵ (see Figure 11). On the contrary, the number (m) of potentially merging BH–BH binaries whose delay between formation and merger is less than 13.5 Gyr is typically 30, multiplied by a factor of order unity; despite our efforts, a spiral-galaxy simulation produced only one merging BH–BH binary (Figure 11). Further, among those few simulated binaries the most massive few can be seen to the greatest distance; thus, only a fraction of these m binaries significantly impact the average chirp mass $\langle {{\cal M}_c}^{15/6}\rangle$ (Equation (9)) calculated therefore our estimate of the detection rate. Given the small number statistics upon which we build our BH–BH detection rate predictions, we very carefully quantify and propagate uncertainties into our estimated BH–BH (and also BH–NS, NS–NS) detection rates.

4.1. Model Uncertainty and Recent Literature

Our simulations do not explore some known model uncertainties that could significantly increase the predicted merger rate. As noted previously, our star-forming conditions involve only near-solar metallicity (Bulik et al. 2008; Belczynski et al. 2010b, 2009). Additionally, we assumed that the common-envelope evolution of a Hertzsprung gap donor led to stellar merger, not a compact binary. However, as discussed in Belczynski et al. (2007), the alternative and still plausible option will increase the BH–BH merger rate by × 500. Finally, we assume binaries form through isolated evolution alone, not allowing for any enhancement in young globular proto-clusters (see, e.g., O'Shaughnessy et al. 2007b; Sadowski et al. 2008). As mentioned previously, because these factors should increase the detection rate, our results are the predictions of a conservative model family.

Our simulations also do not self-consistently explore all model uncertainties corresponding that are comparable or smaller than the statistical simulation errors described extensively below. For example, as the StarTrack code has evolved over the extended assembly of our archive and this paper, our simulations differ somewhat from other contemporary papers that employ it. Notably, unlike Belczynski et al. (2007), we adopt the full Bondi accretion rate during common-envelope evolution; by trapping more mass in the binary, this generally leads to larger detection rates.¹⁵ Also, the StarTrack code does not yet fully explore all uncertainties in single star evolution, such as uncertainties in stellar radii (Fryer et al. 1999) and due to rotational effects (see, e.g., de Mink et al. 2009; Maeder & Meynet 2001, and references therein). Nor did we explore modifications to the recipes used by StarTrack for single and binary evolution. For example, StarTrack models orbital decay in common-envelope evolution with a single factor αλ, though in principle λ differs between and can be calculated for different donor star configurations. For a supernova, the initial (pre-fallback) compact object remnant mass is estimated by tabulated estimates of the core mass (see Belczynski et al. 2008b). Compact remnants (post-fallback) more massive than 2.5 M_☉ are assumed to form BHs.

Finally, all of our results are sensitive to the total number of massive progenitor stars and our estimate for the high-mass IMF. While our computationally limited exploration cannot give an idealized, fully marginalized theoretical prediction, our approach is extremely useful for the reverse problem: how this particular concrete, conservative model family can be tested against future gravitational-wave observations.

4.2. Properties of All CBC Binaries

Though the calculation was applied only to BH–NS and NS–NS mergers, PS-GRB described how to calculate two of the three essential ingredients needed to calculate the CBC detection rate (Equations (9) and (15)): (1) the number of CBC merger events per unit mass in progenitors (the mass efficiency λ_C,α) and (2) the probability that given a CBC progenitor, a merger occurs between t and t+dt since the binary's birth as two stars (the delay-time distribution dP_c,α/dt). For example, as in PS-GRB we estimate the mass efficiency λ for forming a merging binary of type K(=BH–BH, BH–NS, and NS–NS) from the number n of binary progenitors of K with

$\begin{eqnarray} \lambda &=& \frac{n}{N} \frac{f_{\rm cut}}{\left< M \right>}, \end{eqnarray} \tag{ 20 }$

where N is the total number of binaries simulated, from which the n progenitors of K were drawn; <M> is the average mass of all possible binary progenitors; and f_cut is a correction factor accounting for the great many very low mass binaries (i.e., with primary mass m₁ < m_c = 4 M_☉) not included in our simulations at all. Expressions for both <M> and f_cut in terms of population synthesis model parameters are provided in Equations (1) and (2) of O'Shaughnessy et al. (2007a). We also estimate the delay-time distribution as before, by smoothing (in log t) the n simulated delays t between binary formation and merger (see particularly the Appendix of PS-GRB). Finally, though PS-GRB do not mention masses, we use precisely the same logarithmic smoothing technique to estimate $d\!P/d {{\cal M}_c}$ from a set of binaries (see Appendix A).

The procedures described above lead to very reliable results for NS–NS and BH–NS binaries, because many binaries n and even merging binaries m are present in each simulation. However, most simulations have only a few BH–BH binaries whose delays between birth and merger are less than the age of the universe. With so few merging binaries per simulation, statistical uncertainties would severely limit our ability to determine the physically relevant portion of delay-time (dP/dt) and chirp mass ( $p({{\cal M}_c})d {{\cal M}_c}$ ) distributions based on those binaries alone. However, since gravitational-wave decay depends sensitively on the post-supernova orbital parameters of a newly born BH–BH binary, the population of merging BH–BH binaries should be very similar to a population with marginally wider orbits but often dramatically longer decay timescales. For this reason, in this paper we will improve our statistics for dP/dt and $p({{\cal M}_c})$ by including nearly merging binaries with delay times t between birth and merger less than a cutoff T = 10⁵ Myr, as justified by our studies in Appendix A.

4.3. Varying the Binary Fraction

In our population synthesis simulations, we systematically varied almost all parameters that could significantly impact the present-day merger rate. One parameter left unchanged in past studies, however, was the binary fraction f_b, defined as the fraction of stellar systems that are binaries. Without loss of generality, our simulations assume all stars form in binaries. The merger rates per unit mass implied by a population with a lower binary fraction f_b can be related to the merger rates we calculate using

$\begin{equation} {\cal R}(t|f_b) = {\cal R}(t|f_b=1) \frac{f_b(1+{\left<q\right>})}{1+f_b{\left<q\right>}}, \end{equation} \tag{ 21 }$

where <q> = <m₂/m₁> is the average mass ratio of our initial stellar population and which is varied between population synthesis models. In particular, because the effect of this parameter trivially influences our rate predictions, we can analytically incorporate the influence of any prior assumptions regarding the binary fraction f_b.

The true initial binary fraction of stellar populations is not known; it is believed to be between f_b,min = 15% and f_b,max = 100% (Duquennoy & Mayor 1991). For this paper, we will assume f_b could equally likely take on any value in this range. For any given population synthesis model, the binary fraction and mass ratio distribution multiplicatively influence the observed detection rate by a factor

$\begin{equation} X\equiv \frac{f_b(1+{\left<q\right>})}{1+f_b{\left<q\right>}}\;, \end{equation} \tag{ 22 }$

which is distributed between X_min and X = 1 (at f_b = 1) according to

$\begin{equation} p(X)dX = \frac{1+{\left<q\right>}}{[{\left<q\right>}(X-1)-1]^2}\frac{dX}{1-f_{b,{\rm min}}}. \end{equation} \tag{ 23 }$

4.4. Uncertainties and Smoothing: Sampling Errors

Roughly speaking, the relatively small numbers of merging binary BHs in any simulation (m ≲ 100) limit our ability to predict the overall detection rate log R_D more accurately than $1/\sqrt{m}\ln 10$ , because of inaccuracies in reproducing both the mass distribution $d\!P/d {{\cal M}_c}$ and the delay-time distribution dP/dt. Additionally, the relatively small number of simulations of elliptical and spiral galaxies N_E,S ≃ 300–400 provides little chance of discovering high-rate models that are more rare than roughly 1/N. To allow for these two errors, we must convolve each individual simulation's merger rate with kernels reflecting the Poisson uncertainty introduced into each measurement by the limited number of samples. For simplicity, we approximate the relative error in log N as normally distributed with standard deviation $1/\sqrt{N}\ln 10$ . Translating this approximation into a PDF, instead of employing a Poisson posterior distribution we adopt a simple Gaussian kernel to describe the relative error in rate given m and N samples:

$\begin{eqnarray} K_{\rm sim}&=& K_o(\log R_D, 1/\sqrt{m}\ln 10), \end{eqnarray} \tag{ 24 }$

$\begin{eqnarray} \hspace{-0.3pc}K_{\rm samp}&=& K_o(\log R_D, 1/\sqrt{N}\ln 10), \end{eqnarray} \tag{ 25 }$

$\begin{eqnarray} \hspace{-2pc}K_o(x,s) &=& \frac{1}{\sqrt{2\pi s^2}}\exp (-x^2/2s^2). \end{eqnarray} \tag{ 26 }$

As noted above, a third uncertainty is the binary fraction, which we incorporate with the following kernel:

$\begin{eqnarray} K_b(z, {\left<q\right>}) &\equiv & \frac{10^z \ln 10}{1-f_{b,{\rm min}}} \frac{1+{\left<q\right>}}{[{\left<q\right>}(10^z-1)-1]^2} \end{eqnarray} \tag{ 27 }$

for log 0.15 < z < 0. These terms capture the most significant simulation-related sources of error. We have also explored and included several other potential sources of error, including (1) calibration uncertainty in the detector or mismatch in the waveform model (small, typically 10%, multiplied by a factor of order unity),¹⁶ (2) uncertainty in the overall star formation history (a factor roughly 2) and in the relative proportions of elliptical and spiral galaxies (roughly 10%), and (3) additional sampling errors introduced by the sensitivity of the chirp mass average to a few high-mass binaries ( $0.1/\sqrt{n_{\rm eff}}/\ln 10$ , multiplied by a factor of order unity).

4.5. Model Prediction for CBC Rate Densities

To review, using the StarTrack population synthesis code, we have explored a representative sample of stellar evolution produced by elliptical galaxies (N_E = 206) and spiral galaxies (N_S = 282).¹⁷ From these simulations, using the tools indicated in Section 4 we carefully extracted their likely properties (e.g., dP/dt), which in turn we convolved with the star formation history of each of the two major components of the universe (Section 3) to generate a preferred present-day merger rate per unit volume ${\cal R}_{k}(0)$ for each model k, for both elliptical and spiral star-forming conditions.

Though the set of model universe merger rates R_(e,s)k encompasses many of the most significant modeling uncertainties, to be more fully comprehensive and to arrive at a smooth PDF we convolve of the "model universe" merger rates R_s,p and R_e,l for p = 1, ...N_S and l = 1...N_E with the uncertainty kernels and binary fraction kernel described previously assuming equal prior likelihood for each model (i.e., P(E_l) = 1/N_E and P(S_p) = 1/N_S)

$\begin{eqnarray} \bar{K} &=&(K_{\rm samp}*K_{\rm sim}*K_b), \end{eqnarray} \tag{ 28 }$

$\begin{eqnarray} p_e(\log R_{e}) &=& \sum _{l=1}^{N_E} P(E_l) \bar{K} \left(\log \frac{R_{e}}{ R_{e,l}}\right), \end{eqnarray} \tag{ 29 }$

$\begin{eqnarray} p_s(\log R_{s})&=& \sum _{p=1}^{N_S} P(S_p) \bar{K} \left(\log \frac{R_{e}}{ R_{e,p}}\right), \end{eqnarray} \tag{ 30 }$

$\begin{eqnarray} p(\log R) &=& (p_e \oplus p_s). \end{eqnarray} \tag{ 31 }$

Figure 2 shows all of the kernels $\bar{K}$ . In Figure 3, the thin (p(log R)) and dashed (p_e(log R)) curves correspond precisely to the output of this procedure. The bottom right panel also shows our estimate for the total merger rate distribution: a three-fold convolution¹⁸

$\begin{eqnarray} &&p(\log R) = p_{\rm BH-BH}\oplus p_{\rm BH-NS} \oplus p_{\rm NS-NS}. \end{eqnarray} \tag{ 32 }$

**Figure 2.** All of the kernels $\bar{K}$ vs. log X used in this paper, for NS–NS (top panel), BH–NS (center panel), and BH–BH (bottom panel). As extremely many binaries are available in each simulation, the NS–NS simulations' uncertainties are dominated by the uncertainty in the binary fraction *f_b*. Because few BH–BH and occasionally BH–NS binaries are available in each simulation, however, the statistical uncertainty becomes significant.
Download figure:
Standard image High-resolution image

**Figure 2.** All of the kernels $\bar{K}$ vs. log X used in this paper, for NS–NS (top panel), BH–NS (center panel), and BH–BH (bottom panel). As extremely many binaries are available in each simulation, the NS–NS simulations' uncertainties are dominated by the uncertainty in the binary fraction *f_b*. Because few BH–BH and occasionally BH–NS binaries are available in each simulation, however, the statistical uncertainty becomes significant.
Download figure:
Standard image High-resolution image

**Figure 3.** Cumulative probability P(<log R) of various merger rates per unit volume (units Mpc⁻³ Myr⁻¹), without (thin) and with (thick) requiring consistency with the binary pulsar population in the Milky Way. These distributions incorporate uncertainty in the initial binary fraction, the unknown population synthesis model parameters, and statistical errors implied by the limited size of our simulations. For comparison, the dashed lines show the contribution from ellipticals only; the dotted lines show the contribution from spiral galaxies without requiring consistency with the Milky Way.
Download figure:
Standard image High-resolution image

In the above, we assume all simulated population synthesis models are equally plausible a priori. However, observations of single and binary pulsars significantly affect our perspective regarding the relative likelihood of different models. NS–NS merger rates derived from the observed sample are in fact typically higher than what most population models predict for the Milky Way. We adopt standard Bayesian methods to determine the posterior probability P(S_p|C) for a given spiral-galaxy model S_p (p = 1 ... N_S)

$\begin{eqnarray} &&P(S_p|C) \propto P(C|S_p)P(S_p)= P(C|S_p)/N_S \end{eqnarray} \tag{ 33 }$

given the observational constraint C on the present-day spiral-galaxy contribution to the NS–NS merger rate implied by present-day observations of pulsars. We adopt as a constraint the requirement developed in PS-GRB from existing NS-pulsar observations, that the spiral-galaxy merger rate density of double NSs lies between 0.15 and 5.8 Myr⁻¹ Mpc⁻³ (90% confidence). The probability P(C|S_p) that the pth spiral-galaxy model predicts a rate consistent with observations of double NSs is found simply by integrating the PDF P(log R|S_p) for the rate over the constraint interval¹⁹C = [0.15, 5.8] Myr⁻¹ Mpc⁻³. After renormalizing the probabilities P(C|S_P) into P(S_p|C) by requiring ∑_pP(S_p|C) = 1, we arrive at revised predictions for CBC merger rates for BH–BH, BH–NS, and NS–NS binaries as in Equations (29)–(31) but with conditional probabilities:

$\begin{eqnarray} p_s(\log R_{s}|C)&=& \sum _{p=1}^{N_S} P(S_p|C) \bar{K} \left(\log \frac{R_{s}}{ R_{s,p}}\right). \end{eqnarray} \tag{ 34 }$

Convolving the constrained spiral rate density distribution p_s with the unconstrained elliptical density p_e leads to our best estimate for the present-day merger rate per unit volume, shown as the thick solid curve in Figure 3. Our best results favor merger rates between 2 × 10⁻³and0.04 Mpc⁻³ Myr⁻¹ for BH–BH mergers (90% confidence), 0.01and0.28 Mpc⁻³ Myr⁻¹ for BH–NS mergers, and 0.11and1.7 Mpc⁻³ Myr⁻¹ for NS–NS mergers.

Our previous papers have similarly constrained binary evolution in the Milky Way (PSC2) and in the population of spiral galaxies (PS-GRB). Both interpret observations as much tighter constraints than we do here, because these papers did not allow for error (e.g., in each simulation's prediction and, for the predictions per unit volume, in the star formation history of the universe). Our Bayesian constraints are significantly less restrictive. Nonetheless, as seen in Figure 4 the top spiral-galaxy models (90%) have a kick velocity distribution roughly consistent with pulsar observations, similar to that recovered in PSC2 (their Figure 5). Our predictions for the posterior NS–NS merger rate are conservative for another, technical reason: we implement the binary fraction as an uncertainty rather than a parameter. We therefore cannot rule out the lowest low binary fractions we allow (f_b = 0.15) in other spiral galaxies, even though a comparison between binary pulsars in the Milky Way and our model database disfavors such low values.

**Figure 4.** Constrained population synthesis parameters for spiral galaxies: compare to the top right and bottom left panels of Figure 5 of PSC2. In each panel, dark points are included when the corresponding model's posterior probability P(*S_p*|C) is among the top 90%; gray points are shown for all other 282 spiral-galaxy simulations.
Download figure:
Standard image High-resolution image

4.6. Discussion

Two-component universe rarely gives low rates. The probability we assign to very low rates depends sensitively on our implicit assumption that elliptical- and spiral-galaxy star-forming conditions are uncorrelated. The convolution p_e ⊕ p_s (Equation (18)) assigns little probability below the sum of the median elliptical and median spiral merger rate (see, e.g., Figure 6). On the contrary, if all galaxies have the same undetermined star-forming conditions, then the relevant merger rate distribution

$\begin{eqnarray} &&p(\log R) = \sum _p P(S_p)\bar{K}\left(\log \frac{R}{R_{e,p}+R_{s,p}}\right) \end{eqnarray} \tag{ 35 }$

assigns more probability to the lowest rates (see, e.g., the bottom panel of Figure 6). Adopting this implicit prior—identical star-forming conditions in elliptical and spiral galaxies—leads to minimum plausible total BH–BH, BH–NS, and NS–NS merger rates ${\cal R}_{5\%}$ roughly 2.4, 7, and 3.3 times larger than $\mbox{max}({\cal R}_{e,5\%}, {\cal R}_{s,5\%})$ , where ${\cal R}_{e,5\%}$ is the 5% probability elliptical-galaxy merger rate and similarly; compare, e.g., the top panel of Figure 6 with Figure 3.

Ellipticals have different star formation? As in de Freitas Pacheco et al. (2006), our fiducial calculations suggest that elliptical galaxies gave birth to a significant fraction of all present-day compact object mergers, particularly BH–BH mergers (Figure 3). Our results inevitably depend on the amount of and conditions for star formation in elliptical versus spiral galaxies. On the one hand, we assume most ancient star formation has occurred in ellipticals. On the other hand, the Salpeter-like IMFs (p ≃ −2.3) adopted for elliptical galaxies produce compact object binaries much more efficiently than the much steeper high-mass slope (p = −2.7) adopted for spiral galaxies.²⁰ To quantify the impact of each factor separately, in Figure 5 we compare our fiducial BH–BH and NS–NS merger rate distributions (thick solid curves) with predictions based on only spiral or elliptical star-forming conditions. Specifically, each new distribution still follows from Equation (15) (to get ${\cal R}={\cal R}_s$ from our set of spiral-galaxy simulations) and Equation (30) (to create a PDF, based on ${\cal R}$ , an error kernel, and probability P(S)). These comparisons allow us to determine the relative impact that three key factors have on our predicted results: (1) significantly increasing ancient star formation, (2) requiring ancient star formation produce high-mass stars more efficiently than present-day star formation, and (3) adopting independent binary evolution models for elliptical and spiral galaxies. First, both in distribution and on a model-by-model basis, the total compact object merger rate due to nearly steady star formation (blue curves, which adopt the "spiral" SFR $\dot{\rho }=\dot{\rho }_s$ ) is less than the amount due to all past star formation history (red curves, which adopt $\dot{\rho }=\dot{\rho }_e+\dot{\rho }_s$ ) only by a simulation-dependent factor between 1.15and1.3 for NS–NS and between 1.5and2 for BH–BH binaries.²¹ Second, on a model-by-model²² and distribution basis, provided the same amount of star formation our elliptical-galaxy models produce merging binaries roughly 1 (NS–NS) or 1–10 (BH–BH) times more frequently than spiral-only star-forming conditions (both multiplied by factors of order unity), entirely due to their shallower IMF and broader metallicity distribution; compare the red and green curves in Figure 5. Third, though ellipticals produce merging binaries more efficiently per unit mass, fairly few BH–NS and NS–NS binaries remain to merge after such a long delay (dP/dt ∝ 1/t). For BH–NS and NS–NS binaries, elliptical galaxies' present-day merger rate balances two positive factors—very efficient compact binary formation (e.g., ∼3 times higher, based on the IMFs) combined with a high rate of early star formation (∼10 times higher)—against an ∼30 lower merger rate dP/dt via a few Gyr versus ∼100 Myr typical delay (see the discussion in PS-GRB). For BH–NS and NS–NS binaries, this balance roughly leads to comparable merger rates in elliptical and spiral galaxies.²³ To summarize, our fiducial unconstrained predictions differ from previous predictions based on steady, spiral-galaxy star formation by three factors: (1) the addition of ancient star formation, increasing the overall merger rate above steady spiral-only by a factor 1.15–1.3 for NS–NS and 1.5–2 for BH–BH; (2) the conversion of ancient star-forming conditions from spiral like to elliptical, increasing the total present-day merger rate by a factor of order unity for BH–NS and NS–NS binaries; and (3) the assumption that elliptical and spiral star-forming conditions are independent, biasing against low total merger rates.²⁴

**Figure 5.** Relative likelihood of present-day merger rates R for NS–NS binaries (top panel) and BH–BH binaries (bottom panel) given five assumptions: (1) only nearly steady spiral-galaxy star formation produces merging binaries (blue), (2) all galaxies form stars like spiral galaxies (red), (3) all galaxies form stars like elliptical galaxies (green), (4) like (1), but using only models that reproduce the present-day Milky Way merger rate (blue, dotted), and (5) like (2), similarly (red, dotted). For comparison, the thick black curves show the fiducial total BH–BH and NS–NS merger rate (solid curve includes all spiral models; dotted curve only models that reproduce the present-day Milky Way).
Download figure:
Standard image High-resolution image

Ellipticals "dominate" BH-BH rate. Adopting the same binary evolution model but different star formation histories, elliptical galaxies usually produce more merging BH–BH binaries at present than spiral galaxies. In this sense, the elliptical-galaxy BH–BH merger and detection rates "dominate" the total merger rate. However, some elliptical-galaxy models do yield lower merger rates than some spiral-galaxy models; neither ubiquitously dominates. For example, Figure 6 suggests that low but plausible total BH–NS merger rates require more mergers in spirals than ellipticals; high but plausible merger rates require the opposite.

**Figure 6.** Convolutions illustrated. Top panel: the BH–NS merger rates ${\cal R}$ due to ellipticals (dashed), spirals (dotted), and overall (solid) assuming uncorrelated elliptical and spiral star-forming conditions, as in the text and Figure 3. The minimum likely merger rate overall is significantly larger than the minimum likely due to ellipticals or spirals alone. Bottom panel: a Gaussian distribution p(log x)dlog x (dotted); p ⊕ p (solid), analogous to assuming a two-component universe where elliptical- and spiral-galaxy rates are independently drawn from p; and the distribution for y = x + x (dashed), analogous to a two-component universe with identical elliptical and spiral star-forming conditions, both drawn from p.
Download figure:
Standard image High-resolution image

Binary evolution priors. Our predictions depend on the relative likelihood of different binary evolution model parameters. To facilitate comparison with PSC2 and to avoid omitting any vital region of parameter space, we allow and treat as equally plausible a wide parameter range, including extremely large (and empirically implausible) supernova kick magnitudes v_kick ≃ 1000 km s^-1. The reader can reproduce our results or explore alternative assumptions by reweighting the data provided in Table 1, available in the online version of the journal.

Table 1. Parameters of Six Random Binary Evolution Simulations, along with Predicted BH–BH Merger Rates R, Detection Rates R_D, and (for Spiral Simulations) Posterior Probabilities

Type	Z/Z_☉	p	r	w	v₁	v₂	s	αλ	β	f_a	N	N_bbh	m_bbh	log R_bbh	log R*_bbh	log R_D,bbh	P(S_p\|C)
S	1.00	−2.70	0.79	0.07	466.63	111.47	0.43	0.83	1	0.17	3670000	35070	20	−3.01	−2.74	−2.15	0.0006
S	1.00	−2.70	0.55	0.40	343.41	733.15	0.58	0.55	1	0.92	8690000	70975	43	−3.14	−2.95	−2.70	0.0000
S	1.00	−2.70	2.59	0.07	14.72	62.71	0.30	0.36	1	0.86	1640000	21287	10	−2.83	−2.48	−2.37	0.0087
E	1.19	−2.07	2.34	0.86	608.36	559.76	0.10	0.38	1	0.97	1580000	42106	32	−2.31	−1.93	−1.81
E	1.46	−2.08	2.15	0.95	705.84	29.51	0.92	0.57	1	0.75	1500000	29031	13	−2.67	−2.29	−2.18
E	0.74	−2.16	0.10	0.16	473.38	922.95	0.74	0.87	1	0.87	2870000	64357	23	−2.80	−2.42	−2.20

Notes. The first column is a label denoting the type of simulation, elliptical or spiral. The next nine columns provide the parameters passed to our StarTrack binary evolution code: p is the high-mass IMF power law; r is the mass ratio distribution exponent; w is the wind strength; v_1,2 (in km s⁻¹) and s characterize the strength of supernova kicks; and αλ, β, and f_a pertain to mass transfer. The next three columns provide the total number of binaries simulated (N), the number of BH–BH binaries seen (N_bbh), and the number of BH–BH binaries that could merge (m_bbh). The next three columns provide present-day BH–BH merger rate densities (R, R*, in Mpc⁻³ Myr⁻¹) and initial LIGO detection rates (R_D, in yr⁻¹). Finally, for spiral-galaxy simulations, the last column provides the posterior probability P(S_p|C), describing how consistent that model's NS–NS merger rate predictions are with the galactic Milky Way population.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as: Data Typeset image

Binary evolution parameters. Finally, as outlined in Section 4.1, we fix a number of assumptions which can dramatically influence merger and even detection rates. For example, the maximum NS mass m_+NS completely determines the nature of final compact binaries, as it subdivides a roughly m_+NS-independent mass distribution into three regions: NS–NS, BH–NS, and BH–BH. Evidently, the individual component merger rates depend sensitively on m_+NS.

However, as relevant gravitational-wave interferometers cannot easily distinguish between two compact objects of comparable mass 1 M_☉ < m < 3 M_☉, the overall detection rate described later depends only weakly on m_+NS. Similarly, while the NS birth mass distribution can influence the relative proportions of each type, as well as the observable pulsar binary mass distribution, our results should not depend significantly on the amount distribution of masses.

Several other binary evolution parameters not explored here have previously been shown to influence merger rates significantly, albeit much less significantly than the parameters we explored here. First, as noted previously in Section 4.1, we adopt the full Bondi–Hoyle accretion rate onto BHs in common envelope, rather than the smaller rates suggested by recent hydrodynamical simulations (cf. Belczynski et al. 2008a, 2007b). A smaller mass accretion rate, if adopted, would reduce the incidence of accretion induced collapse, marginally changing the proportions of BH–BH, BH–NS, and NS–NS binaries; marginally lower the masses of some merging binaries, increasing the factor by which common-envelope evolution tightens orbits; and generally change merger and detection rates by a small factor compared with the overall range considered (footnote 15). Second, we examined only isotropic kick distributions. Several authors have suggested spin-kick alignment on theoretical and empirical grounds (see, e.g., Kalogera et al. 2008; Kuranov et al. 2009; Wang et al. 2007, and references therein). As the same mass transfer that brings binaries close enough to merge also aligns their spins and orbit before the second supernova, polar kicks are particularly efficient at disrupting binaries (see Postnov & Kuranov 2008). Though polar kicks would dramatically transform our predictions, the observed PSR-NS merger rate in the Milky Way not only corresponds to a high, not low rate²⁵ but also implies model supernova kick magnitudes comparable to observed pulsar proper motion velocities (PSC2, as well as our Figure 4). Third, we adopt a specific model for mass loss in massive stars; by changing the evolution of massive progenitors, alternate stellar wind models can dramatically modify merger rates (cf. Belczynski et al. 2010b). Fourth, we adopt a single, time-independent metallicity for each type of star-forming region. Binary merger rates can depend sensitively on metallicity, particularly through metallicity-dependent stellar winds (Belczynski et al. 2010b). A future publication will explore the latter two issues in more detail.

Advanced detectors and redshift-dependent mass evolution. For simplicity, we have adopted a single time-independent mass distributions for each type of merging binary: BH–BH, BH–NS, and NS–NS. Looking backward to higher redshift, however, the relative proportions of each binary type change, as each possesses a different delay-time distribution dP/dt. In other words, despite using "time-independent" mass distributions, the instantaneous total mass distribution, found by gluing these three distributions together in proportion to their instantaneous merger rates ${\cal R}$ , will vary with time. Networks of advanced detectors and particularly third generation detectors will probe a time-dependent mass distribution. The authors and collaborators will also address this issue in a future paper.

5. CBC DETECTION RATES AND DETECTION PROBABILITY

5.1. CBC Detection Rates for LIGO

Using the StarTrack population synthesis code, we have explored a representative sample of stellar evolution produced by elliptical galaxies (N_E = 206) and spiral galaxies (N_S = 282). From these simulations, using the tools indicated in Section 4 we carefully extracted their likely properties (e.g., dP/dt), which in turn we convolved with the star formation history of each of the two major components of the universe (Section 3). Using our understanding of the LIGO range (Section 2), we can convert these merger rate histories into expected single-interferometer LIGO detection rates. Finally, though the set of model universe detection rates R_D,(e,s)k encompasses many of the most significant modeling uncertainties, to be more fully comprehensive and to arrive at a smooth PDF we plot in Figure 7 the convolution of the "model universe" detection rates R_D,s,p and R_D,e,l for p = 1, ...N_p and l = 1...N_E with the uncertainty kernels and binary fraction kernel described previously assuming equal prior likelihood for each model (i.e., P(E_l) = 1/N_E and P(S_p) = 1/N_S)

$\begin{eqnarray} \bar{K} &=&(K_{\rm samp}*K_{\rm sim}*K_b), \end{eqnarray} \tag{ 36 }$

$\begin{eqnarray} p_e(\log R_{D,e}) &=& \sum _{l=1}^{N_E} P(E_l) \bar{K} \left(\log \frac{R_{D,e}}{ R_{D,e,l}}\right), \end{eqnarray} \tag{ 37 }$

$\begin{eqnarray} p_s(\log R_{D,s})&=& \sum _{p=1}^{N_S} P(S_p) \bar{K} \left(\log \frac{R_{D,e}}{ R_{D,e,p}}\right), \end{eqnarray} \tag{ 38 }$

$\begin{eqnarray} p(\log R_D) &=& (p_e \oplus p_s), \end{eqnarray} \tag{ 39 }$

and subsequently combine the elliptical and spiral detection rate distributions (Equation (18)). (Most of these convolutions can be performed trivially by adding the standard deviations of Gaussians in quadrature; the remaining two convolutions are performed numerically.) As previously we estimate the total detection rate via a three-fold convolution (Equation (32)).²⁶ Though the same basic techniques applied above can be and have been applied to BH–NS and NS–NS binaries previously (PS-GRB), we now incorporate error propagation, a simulation-by-simulation estimate of the relevant chirp mass, and a variable binary fraction; the merger rate distributions they provide therefore are not proportional to our predictions.

Figure 7 shows the distribution of CBC detection rates predicted by Equation (10) for individual initial LIGO interferometers as part of a network (thin curves) as well as the contribution from elliptical galaxies alone (dashed curves).²⁷ As indicated by the close proximity of the elliptical and total merger rate distributions, most BH–BH mergers that LIGO detects should be produced in elliptical galaxies. This elliptical bias arises because BH–BH binaries have long delays between their formation and eventual merger; their present-day merger rate is more easily influenced by ancient star formation. On the other hand, the difference between the dashed and thin curves in the bottom left panel indicates that the NS–NS merger rate must be spiral dominated.

Because the initial LIGO detectors essentially probe different-sized volumes of the local universe for each merger type, their detection rate distributions are identical to those of any similar-scale single- or multiple-interferometer network, mod a constant horizontal offset determined by the relative increase in volume NS–NS binaries can be seen [log V_NS-NS/4π(14 Mpc)³/3]. Similarly, because a single advanced LIGO detector's reach to NS–NS and BH–NS binaries is often not cosmologically significant, weighting over detector and source orientation, our predictions for NS–NS and BH–NS detection rates are identical to those for initial LIGO, mod an offset [= log(197/14)³ = 3.44]. However, the typical range of a single detector to BH–BH binaries does extend past z ≃ 0.25, where the SFR is noticeably higher and where cosmological volume and redshift factors become significant. But because of the long delay between BH–BH merger and progenitor birth, the BH–BH merger rate R(t) generally increases much more slowly with redshift than the star formation history. As a result, the single-interferometer BH–BH detection rate distribution is also comparable to but slightly less than what we expect after scaling up initial-LIGO results to larger range: a detailed calculation suggests a merger detection rate between 3and300 yr⁻¹ (90% confidence), in good agreement with the 20–500 yr⁻¹(D_bns/197 Mpc)³ (90% confidence) expected from cubic rescaling of the initial-LIGO result.²⁸ For networks of advanced detectors, however, cubic scaling breaks down severely; the range is significant enough to require a full cosmological treatment and network-sensitivity-beampattern averaging.

5.2. Comparisons with Pulsar Observations of Galactic NS–NS Binaries

Following the previously discussed procedure to re-evaluate the likelihood of our spiral-galaxy population synthesis models, we arrive at a distribution of initial LIGO detection rates as indicated in Figure 7. Comparing these distributions with the unconstrained results of the previous calculation, we arrive at the following conclusions:

Higher NS–NS detection rate. As seen in Figure 3, inferences drawn from the double pulsar population about the present-day NS–NS merger rate lie on the high end of what our models can presently reproduce. By requiring our models reproduce those observations, our constrained distributions naturally favor higher NS–NS detection rates.

Lower BH–BH detection rate from spirals. Compared to our previous results (Figure 7), the highest binary BH merger and thus detection rates are ruled out in spiral galaxies. Physically, the observational constraint supports those population synthesis models where one particular parameter (αλ, related to common-envelope evolution) is not exceptionally small (αλ>0.1; see for comparison the top left panel of Figure 5 in PSC2). In turn, this parameter αλ correlates strongly with the BH–BH merger rate, with low values of αλ being required to produce the largest BH–BH merger rates. Thus, because observations of NS–NS binaries in the Milky Way favor high NS–NS merger rates, they disfavor the highest rates of BH–BH mergers in spiral galaxies.

Observations of the NS–NS merger rate in spiral galaxies do not constrain ellipticals, BH–NS well. With the exception of NS–NS mergers, Figure 3 demonstrates that most constrained predictions for merger rates (and, though not shown here, LIGO's detection rate) agree strikingly well with our a priori expectations. These consistent predictions should be expected, since observations of spiral galaxies cannot constrain differences (e.g., due to metallicity and IMF) in the binary evolution in elliptical galaxies. In particular, the single most likely source of a gravitational-wave detection in the near future, BH–BH mergers, is due to their long characteristic ages expected to occur primarily in elliptical galaxies; as a result, the present-day BH–BH and total CBC detection rate depends only weakly on our ability to rule out certain population synthesis parameters in spiral galaxies. Furthermore, even in spiral galaxies alone, as seen in PSC2 and references therein, the NS–NS and other CBC merger rates are not tightly correlated. Thus, despite the fairly strong and consistent constraints that double pulsar observations imply for binary evolution parameters, as described in detail in PSC2, our overall expectations for LIGO's detection rate remain largely unchanged.

5.3. Net Probability of CBC Detection

Factoring in Poisson statistics for each model, the net prior probability P_d that a single LIGO interferometer with a single-interferometer range to NS–NS inspiral of D_bns operating for time T can make at least one detection is

$\begin{eqnarray} P_{{\rm detect}}(\ge 1| {D_{\rm bns}},T)&=& \int dR_D p(R_D)(1- e^{-R_D T}),\qquad \end{eqnarray} \tag{ 40 }$

where D_bns implicitly enters as a scale factor for R_D (∝D_v³). More generally, for a multi-interferometer configuration with a time-varying network sensitivity, the appropriate detection probability can be easily calculated from our reference probability by replacing T in the above expression by ∫dt(D_bns(t)/D_bns_ref)³; the latter quantity characterizes the 4-volume to which the search is sensitive. This detection probability can be significantly increased by improvements in the LIGO detector or run schedule that improve the detection 4-volume VT. Assuming a fixed, orientation-averaged range D_bns for any fraction of the whole LIGO network, the total probability of detecting one or more events by that fractional network can be approximated by

$\begin{eqnarray} P_{{\rm detect}}(\ge 1|14\, {\rm Mpc},T)&\simeq & 0.4 + 0.5 \log (T/8 \, {\rm yr}), \end{eqnarray} \tag{ 41a }$

$\begin{eqnarray} P_{{\rm detect}}(\ge 1|197 \, {\rm Mpc},T)\!&\simeq & \! 0.4 + 0.5 \log (T/0.01\, {\rm yr}),\qquad \end{eqnarray} \tag{ 41b }$

for T between one year and $100\, {\rm }$ years or 2 × 10⁻³ and 0.1 yr, respectively. These approximations are shown in Figure 8.

Equivalently, our calculations suggest that a gravitational-wave network must have a four-volume sensitivity of 1.3 × 10⁵ Mpc³ yr (roughly equivalent to a 32 Mpc sphere for one year) in order to have a 50% chance of detecting one or more events, as

$\begin{eqnarray} &&P_{\rm detect}(\ge 1|{D_{\rm bns}}, 1\, {\rm yr})\simeq 0.5 + 1.5 \log {D_{\rm bns}}/32\, {\rm Mpc} \end{eqnarray} \tag{ 42 }$

for D_bns between 20 and 70 Mpc.

5.4. Comparison to Earlier Studies

Our fiducial predictions for compact object merger rates (Figure 3) and gravitational-wave detection rates (Figure 7) reflect our uncertainties about the star formation conditions and history of the universe. For example, results consistent with current observational constraints require compact binary formation in elliptical galaxies and therefore disallow the worst-case scenario: mergers formed only in spirals by our least-efficient spiral-galaxy model. Nonetheless, our broad parameter study encompasses many scenarios and predictions. For 100% binary fraction and for each type of compact binary (NS–NS, BH–NS, and BH–BH), the most optimistic and pessimistic single-component scenarios in our comprehensive table of results (Table 1, full table available in the online version of the journal)²⁹ correspond to merger rates extending between 10⁻⁴and1 Myr⁻¹ Mpc⁻³, multiplied by a factor of order unity.

Detailed quantitative comparisons of our results to others appearing in the literature must account for alternate normalization and star formation history choices, which often vary significantly from study to study. Our constant star formation spiral-galaxy predictions³⁰ in Table 1 are comparable to Table 4 in the review article by Postnov & Yungelson (2006), to Table I in Belczynski et al. (2002), to Tables 2 and 3 in Belczynski et al. (2007), and to Table 5 in Hurley et al. (2002). Our constant star formation elliptical-galaxy³¹ predictions can be compared to Tables 4, 6, and 7 of Voss & Tauris (2003) and to Table 1 in Portegies Zwart & Yungelson (1998). Our results assuming all galaxies have spiral-like star-forming conditions (cf. Figure 5) can be compared with Bulik & Belczyński (2003). Finally, our constant and time-dependent elliptical-galaxy results are directly comparable to de Freitas Pacheco et al. (2006). Even after allowing for different inputs, a model-by-model comparison of our results with others available in the literature illustrate physical differences between different authors' population models and our own (Section 4.1). For example, unlike Voss & Tauris (2003), rather than calculate the envelope binding energy fraction λ, we fix a single factor αλ for all binaries in common envelope, leading to smaller BH–BH merger rates for comparable parameters (see, e.g., the discussion in Postnov & Yungelson 2007). As discussed in Section 4.6, unlike Postnov & Kuranov (2008) we adopt isotropic kicks, avoiding the least favorable kick configurations and merger rates (see their Figure 3). Unlike model C in Belczynski et al. (2007), we assumed that the common-envelope evolution of a Hertzsprung gap donor led to stellar merger, not a compact binary. We apply a similar kick to all NS, rather than assign a different kick based on its prior (binary) evolution history (Pfahl et al. 2002) or its supernova mechanism (i.e., a low kick for electron-capture supernova; see Linden et al. 2009; Kalogera et al. 2008, and references therein). Additionally, as described in Sections 4.1 and 4.6, we make several other specific choices consistent with our previous work but not universally adopted in the literature (e.g., for the Bondi accretion rate in episodes of hypercritical common-envelope accretion; for maximum NS mass). However, a detailed exploration of all these correlated effects would require coordinated simulations performed for the purpose of detailed comparisons; such an undertaking is beyond the scope of this paper. Instead, we are satisfied that our results are broadly consistent with early studies of the Milky Way, particularly given the number and parameter survey breadth of the model database used here compared with earlier studies by other groups. For comparison, in the analysis of our present results we do make detailed comparisons to earlier studies from our group (Section 4.6 and again in our concluding discussion), to make the origin of any apparent inconsistencies fully transparent.

Despite providing a distribution, our fiducial results conservatively do not account for any other channels for forming compact objects besides isolated binary evolution, even though several observational, semianalytic, and numerical studies have demonstrated that interacting stellar environments can efficiently form compact binaries. For example, in young very massive stellar clusters, a high-mass BH subcluster can mass segregate and rapidly form and eject merging BH binaries (see, e.g., Portegies Zwart & McMillan 2000; O'Shaughnessy et al. 2007b; Banerjee et al. 2010, and references therein). A similar process occurs in nuclear star clusters (see Miller & Lauburg 2009). Alternatively, young dense clusters could undergo runaway stellar collisions, eventually forming an "intermediate-mass" BH M ∈ 100–10³ M_☉, each of which can capture stellar-mass compact objects (Mandel et al. 2008) or even form binaries and merge themselves (Fregeau et al. 2006; Mandel et al. 2009). Theory aside, observations suggest pulsar and X-ray binaries are highly overabundant in clusters (see, e.g., Rasio et al. 2000; Pooley et al. 2003; Ransom et al. 2005, and references therein). Further, several authors have argued observations of short gamma-ray bursts' (GRBs') hosts, luminosities, and redshift distribution support a multicomponent origin, as cluster dynamics lead to a different merger delay-time distribution P(>t) (see Hopman et al. 2006; Salvaterra et al. 2008; Guetta & Stella 2009, and references therein). Though present-day massive stellar clusters have little mass, if a sufficient fraction of all star formation occurs through clusters and a sufficient fraction occurs in massive clusters, similar highly efficient dynamical processes could produce many merging compact binaries of all types (see, e.g., Sadowski et al. 2008). To summarize, both theory and several different types of observation suggest that all types of dynamically produced compact binary merger exist. These mergers could occur as and potentially orders of magnitude more frequently than mergers produced from field binaries; the estimates shown in Figure 7 therefore underestimate the total merger rate, summed over all channels. Nonetheless, we anticipate that future observations might distinguish between the relative contributions of each component. Assuming even rough separation, our results determine how and how precisely future measurements of compact object merger rates (through gravitational waves, short GRBs, or galactic observations) will constrain binary evolution.

6. CONCLUSIONS AND DISCUSSION

In this paper, we provide the first estimate of the CBC detection rate for initial and advanced LIGO which accounts for (1) star formation in both elliptical and spiral galaxies as a function of redshift, (2) binary BHs, (3) a range of plausible binary fractions, and most significantly (4) a large range of plausible binary evolution scenarios. Our results strongly indicate that that existing gravitational-wave detectors are at the threshold of detection; that moderate improvements in these detectors, such as the enhanced LIGO upgrade of 2009, could plausibly lead to detecting at least one binary merger in a year of operation; that advanced LIGO detectors should be reasonably expected to detect inspirals over 1–2 yr of operation; and that even the absence of a detection by advanced interferometers would very significantly constrain the set of model parameters that could be consistent with observations.

Our predicted detection rates are somewhat higher than those presented in PSC2 because we include elliptical galaxies: mergers in elliptical galaxies likely dominate our estimates of the BH–BH and BH–NS event rates. However, our previous estimates effectively assumed a high present-day SFR per unit volume: an assumed Milky-Way-equivalent galaxy density of 0.01 Mpc⁻³ and SFR per spiral galaxy of 3.5 M_☉ yr⁻¹ correspond to a present-day SFR three times higher than observed. For this reason, a naive product of the merger rates per Milky Way galaxy drawn from PSC2's Figure 6 leads to preferred merger rates per unit volume that agree well our preferred values. For example, our current results and this naive rescaling of PSC2 results both favor local merger rate densities of 10⁻² Mpc⁻³ Myr⁻¹ for BH–NS binaries and 3 × 10⁻² Mpc⁻³ Myr⁻¹ for binary NSs. Correcting for the star formation history normalization, however, our calculations suggest three times as many mergers per volume for BH–NS and NS–NS than in PSC2.

Similarly, our predicted detection rates are higher than those implied by Figure 7 in PSC2, correcting for the previously mentioned SFR bias (reducing the rates shown in that figure by 3) and rescaling to a single interferometer rather than network range (reducing the rate by ×10).

Limiting attention to elliptical- plus spiral-galaxy models that include only those spiral-galaxy history models that reproduce current observations of double pulsars, we find that a single detector in the initial LIGO network should be sensitive to 2.4 × 10⁻²–0.46 (ρ_c/8)⁻³(D_bns/14 Mpc)³ mergers per year (90% confidence level), where we sum over all merger types to produce an overall detection rate, where ρ_c = 8 is the threshold signal-to-noise ratio of the search and D_bns is the radius of the effective volume inside of which a single detector in the appropriate array could observe the inspiral of two 1.4 M_☉ NSs. We estimate that the LIGO detector network has a probability P_detect(⩾1) ⩾ 0.4 + 0.5log(V/V_c)(T/8 yr) of detecting a merger during each year of operation, where V = 4πD_bns³/3 is the volume to which a multidetector search for binary NS inspiral is sensitive, V_c is the volume inside a 14 Mpc radius sphere, and 1/3 < (V/V_c)(T/ yr) < 8. For example, the probability P_detect of detecting one or more mergers can be approximated by (1) P_detect ≃ 0.4 + 0.5log(T/0.01 yr), assuming D_bns = 197 Mpc and it operates for T yr, for T between 2 days and 0.1 yr); or by (2) P_detect ≃ 0.5 + 1.5log D_bns/32 Mpc, for one year of operation and for D_bns between 20 and 70 Mpc.

In this paper, we have continued to employ approximate waveform and detection models to estimate the sensitivity of ground-based gravitational-wave detectors to moderate-mass (M ≃ 10 M_☉) BH–BH mergers (see, e.g., Flanagan & Hughes 1998; Buonanno et al. 2003b for a review). Advanced detectors, however, will be sensitive at cosmological distances to the later phases of high-mass mergers, a complex strong-field and strongly precessing regime which is only beginning to be thoroughly understood through numerical simulations. Progress in numerical waveform modeling and search pipleline analysis is required to better understand how often these candidates could be identified, especially for spinning systems, as expected in nature.

6.1. Anticipated Constraints from Future Observations

Despite the existence of two competing scenarios for compact binary formation, in many cases the electromagnetic and gravitational-wave signatures of binaries formed dynamically and in isolation can be distinguished. For example, BH binaries formed dynamically should have random spins and therefore (generically) beating gravitational inspiral waveforms, while isolated binary evolution models favor spin–orbit alignment. Once binaries with isolated progenitors are identified, their gravitational-wave signals will encode both the number and properties of binaries formed through isolated evolution, information that in turn can tightly constrain isolated binary progenitor parameters. To use an extreme example, future BH–BH observations could imply a gravitational-wave detection rate R_D consistent with only the highest event rates shown in Figure 7. In the context of the models presented here, such an observation would imply (1) that elliptical galaxies host most BH–BH binaries; (2) that massive stars have a high binary fraction; (3) that a certain phenomenological parameter is very small (αλ < 0.1; Section 5.2); and thus (4) since Milky Way observations of binary pulsars exclude low αλ, that isolated binaries formed in elliptical galaxies evolve differently to those in spirals. However, by analogy with the insights provided from observations of isolated Galactic pulsars, typical detection rates will be compatible with many plausible progenitor models (see PSC2; O'Shaughnessy et al. 2005). Rather than constraining single parameters tightly, the information provided by the first few observations will be encoded only in high-order correlations among many model parameters (see, e.g., Figure 5 in PSC2 and the highly correlated merger rate fit presented in Appendix B in O'Shaughnessy et al. 2005). Finally, what we can learn about binary evolution depends on the merger detections that nature provides. Because correlations between parameters change significantly across the model space, we cannot describe posterior uncertainties once and for all.

If NS mergers occasionally have electromagnetic counterparts (i.e., short GRBs), then well-localized electromagnetic counterparts to binary NS mergers, such as (prompt or orphan) short GRB afterglows (Gehrels et al. 2009; Rau et al. 2009; Bloom et al. 2009) or radio bursts, allow optical follow-up studies of its host galaxy's star-forming and chemical enrichment history. Combined with the position of the burst on the galaxy (kick) and, if within range of gravitational-wave detectors, the binary component masses, this host information can tightly limit isolated binary progenitor models (see, e.g., PS-GRB). For example, these observations could provide unique direct probes on the distribution of delays between binary birth and merger (see, e.g., Nakar 2007; Prochaska et al. 2006; Zheng & Ramirez-Ruiz 2007; Rhoads 2010, and references therein). Better still, offsets from host galaxies might provide an independent constraint on delay times (and on supernova kick velocities; see, e.g., Belczynski et al. 2006; Troja et al. 2008, and references therein). Unfortunately, purely electromagnetic observations of short GRBs cannot distinguish between multiple source populations, such as cluster-formed mergers or even non-merger-powered bursts (see, e.g., Virgili et al. 2009; Chapman et al. 2009). Worse, given plausible kick offsets, purely electromagnetic observations may not isolate the appropriate host (Zemp et al. 2009). Fortunately, if existing space-based and ground-based resources are applied in the advanced LIGO epoch, coincident gravitational-wave and electromagnetic observations can resolve these ambiguities.

Finally, as always with binary evolution, any further observations pertaining to binary stars can improve our understanding and limit model-building freedom, such as SNIa rates, galactic pulsar and X-ray binaries, and even galactic WD binaries (see, e.g., the extensive discussion in Nelemans et al. 2009).

The authors appreciate helpful comments received from Ilya Mandel, Alberto Vecchio, Patrick Brady, Chad Hanna, Steve Fairhurst, Tania Regimbau, Alan Weinstein, and all the members of LIGO inspiral search group. R.O. was supported by the National Science Foundation awards PHY 06-53462 and the Center for Gravitational Wave Physics. The Center for Gravitational Wave Physics is supported by the George A. and Margaret M. Downsbrough Endowment and by the National Science Foundation under cooperative agreement PHY 01-14375. V.K. was supported by NSF grant PHY-0653321. K.B. acknowledges the support from the Polish Ministry of Science and Higher Education (MSHE) grant N N203 302835 (2008–2011).

APPENDIX A: RECONSTRUCTING THE MASS AND TIME DISTRIBUTIONS

Selecting binaries for mass estimation. We prefer to smooth dP/dt over a shorter timescale than the maximum allowed by evolution of isolated BH–BH binaries. On physical grounds, double BH binaries whose delay times are comparable to or smaller than the age of the universe are more likely to be similar to merging BH–BH binaries (the focus of our study) than wide binaries that have never interacted. To quantify the point at which a the transition occurs and therefore the timescale inside of which BH–BH binaries are comparatively similar, we calculate the "average" chirp mass of binaries $\tilde{ {{\cal M}_c}}(t)$ merging within a factor 10 of time t:

$\begin{eqnarray} \bar{ {{\cal M}_c}}(t)&=&\left[ \frac{\sum _{j: 0.1 t <t_j< 10 t } {{\cal M}_c}{}_j^{15/6} }{\sum _{j: t 0.1 <t_j< t 10 } 1} \right]^{6/15}. \end{eqnarray} \tag{ A1 }$

Figure 9 shows the results of this calculation. Based on the delay-time distribution P_m(<t) and average chirp mass of similar binaries $\bar{ {{\cal M}_c}}(t)$ , a sharp distinction exists between binaries with t ≲ 10⁷ Myr and wider binaries. These massive wide binaries have orbits largely unaffected by binary evolution and have a delay-time distribution dP/dt ∝ 1/t largely determined by their initial orbital configuration. On the contrary, tighter and less massive binaries have been significantly perturbed by binary evolution. Within each simulation, they appear to be drawn from a similar chirp mass distribution (allowing for significant sampling-induced fluctuations in $\bar{ {{\cal M}_c}}(t)$ ).

Predictions and smoothing: physical timescales. From the above similarity argument, we suspect the n_eff BH–BH binaries with delay times <10⁵ Myr are similar to the set of binaries with delay times <13.5 Gyr. However, even with those binaries we often have only ∼50 with which to estimate $p({{\cal M}_c})$ ; the chirp mass distribution must be estimated with smoothing. We use the same approach as is used for dP/dt in PS-GRB to estimate the chirp mass distribution:

$\begin{eqnarray} \hat{p}({{\cal M}_c}) &\equiv & \frac{1}{ {{\cal M}_c}n_{\rm eff} \ln 10}\nonumber\\ &&\times \sum _{k=1}^{n_{\rm eff}} \frac{1}{\sqrt{2\pi (s_2)^2}} e^{-(\log {{\cal M}_c}- l_{m,k})^2/2(s_2)^2} \end{eqnarray} \tag{ A2 }$

$\begin{eqnarray} \hspace{-1.5pc}s_2 &\equiv & \frac{[({\rm max}_k l_{m,k}) - ({\rm min}_k l_{m,k})]}{1.25\sqrt{n_{\rm eff}}}, \end{eqnarray} \tag{ A3 }$

where the $l_{m,k}=\log _{10}({{\cal M}_c}/M_\odot)$ for k = 1...n_eff are the logarithms of each binaries' chirp mass for the n_eff binaries with delay times <10⁵ Myr. Similarly, the delay time distribution dp/dt is also estimated using all n_eff binaries, exactly as in PS-GRB; see Figure 10.

**Figure 10.** Detailed view of the cumulative probabilities *P_m*(<t) for four randomly selected spiral (left) and elliptical (right) models vs. log t/ Myr (compare with Figure 9). Solid curves indicate the fit; points indicate the raw cumulative *P_m*(<t). The horizontal dashed lines indicate the smallest and largest possibilities for log P ≈ log 1/n, the smallest resolvable cumulative probability; see Figure 11 and the associated text for a discussion of its range and relation to m (i.e., the largest 1/n occurs in only a handful of simulations, most of which have many merging binaries m ≫ 1 and therefore a well-resolved P(<t)). Our estimate reproduces *P_m*(<t) reliably at t ≃ 100 Myr, the smallest scale on which ρ(t) varies for the redshift range of interest. Though our fits necessarily fail to reliably recover *dP_m*/dt on very short timescales, uncertainty in integrals of *dP_m*/dt over longer timescales average out (i.e., in merger rate estimates ${\cal R}(t)$ using Equation (15)).
Download figure:
Standard image High-resolution image

Absent correlations. We identify very few significant correlations on the timescales of interest between $\log {{\cal M}_c}$ and log t. For example, the correlation coefficient between $\log {{\cal M}_c}_k$ and log t_k for each simulation's n_eff close BH–BH binaries is roughly consistent with sampling error ( $\propto 1/\sqrt{n_{\rm eff}}$ ). For BH–NS and particularly NS–NS, nominally significant correlation coefficients occur. For example, some NS–NS $({{\cal M}_c},t)$ distributions concentrate a fraction of mergers on and exclude a region below a critical $t({{\cal M}_c})$ curve, associated with mass-transfer-limited ultracompact mergers on times t ⩽ 100 Myr. Other distributions show no apparent structure, but through large n_eff and therefore extremely narrow mass distributions yield nominally significant correlation. Our manual inspection of all these outliers suggests no structure significant enough to impact merger rates (i.e., associated both with significant change in ${{\cal M}_c}$ and long star formation timescales). On the contrary, as seen in Figure 9, the joint distribution $p({{\cal M}_c},t)$ varies significantly between simulations, as the distribution in each factor [ ${{\cal M}_c}$ and t] changes significantly. The authors can provide all binary properties [(m₁, m₂, t)] for each simulation on request.

APPENDIX B: SAMPLING LIMITATIONS AND ARCHIVE SELECTION

Both sets of simulations are highly heterogeneous: the number of binaries simulated (N), the number of bound BH–BH binaries remaining at the end of stellar evolution (n), and the number of those binaries which merge within 13.5 Gyr (m) all span several orders of magnitude, as shown in Figure 11. Evidently, only for those simulations with many BH–BH binaries (n ≫ 1000)—particularly with several merging binaries (m ≫ 1)—can we reliably extract relevant physics, such as the present-day merger rate and characteristic chirp mass of BH–BH binaries due to a starburst in the early universe. However, as discussed in PS-GRB, the obvious subsets of simulations for which we can reliably make predictions, such as the set of simulations with n>3000, could have been biased, overrepresenting those simulations which predict more mergers per simulated binary. To mitigate the influence of this bias, in addition to requiring that each simulation have many BH–BH binaries (n>3000), we also (or, equivalently, instead) require that nN>1 × 10⁹. The second cut's prefactor (1 × 10⁹) is chosen so that every simulation that satisfies it has n>3000. The subset of simulations which satisfies both conditions has a relatively unbiased distribution of n/N, the fraction of simulated binaries that end up as bound BH–BH binaries (PS-GRB). These two cuts do not introduce biases in any of the parameters: the set of simulations which satisfy these cuts has similar statistical properties as the entire archive (we omit the many plots needed as proof).

**Figure 11.** Left plot: for spiral (top panels) and elliptical (bottom) archives, a scatter plot of the number of BH–BH binaries n seen in a simulation of N progenitor binaries (each progenitor component having mass greater than 4 M_☉; higher primary masses are drawn from the broken-Kroupa IMF). The two solid lines show the cutoffs n>3000 and nN>1 × 10⁹ imposed to insure data quality and reduce sampling bias. Simulations used in this paper, shown in blue, must lie above and to the right of these cutoffs. Right plot: for the simulations above and to the right of the lines shown to left, scatter plot of the total number of bound BH–BH binaries remaining in the simulation (n) vs. the number m of potentially merging binaries (i.e., with a delay between birth and merger less than 13.5 Gyr). All selected simulations (nN>1 × 10⁹) have at least one merger (m>0).
Download figure:
Standard image High-resolution image

Handling undersampled cases. Finally, with each simulation guaranteed to have 3000 binaries, most simulations have a minimum of ∼1 merging BH–BH binaries (see Figure 11 for details). Some simulations still do contain very few binaries that are merging or even nearly so, largely because certain combinations of physical circumstances rarely lead to the formation of many tight BH binaries. Given that we construct model universes from independent elliptical and spiral simulations, these simulations will lead to an indeterminate, small detection rate in only 1/N_EN_S ≃ 1/300² universe models, multiplied by a factor of order unity. Even considering only elliptical or spiral-galaxy star formation, they influence our predictions O(1/N_E,S) ≃ 1% of the time, multiplied by a factor of order unity. We therefore adopt a few percent as the lowest resolvable probability of any event our simulations can resolve.

Potential biases in multiply selected simulations. Just as the set of simulations was reduced to find simulations with many BH–BH binaries in Section 4, here we require our simulations additionally have enough NS–NS binaries to enable accurate estimates.³² Most population synthesis simulations we performed that produce an adequate number of BH-BH binaries also produce more than enough NS–NS and NS–BH binaries (∼10⁴) to allow us to estimate the three properties $(\lambda,d\!P/dt, p({{\cal M}_c})$ needed to reconstruct merger and detection rates for these NS binaries. However, certain combinations of population synthesis parameters strongly favor merging BH–BH production over NS–NS production. Most of our simulations of these conditions did not produce enough NS–NS binaries to allow accurate estimates. In other words, the set of simulations that have both enough NS–NS and close BH–BH binaries is slightly biased against the highest BH–BH and lowest NS–NS merger rates (see Table 1).

The goal of this paper is to determine the relative likelihood of CBC merger and detection rates given existing observations. In particular, in Section 5.2 we examine the distribution of BH–BH merger rates for (spiral-galaxy) simulations for which the NS–NS merger rates are consistent with Milky Way observations. This subset of simulations is inevitably contained within the biased set of simulations for which NS–NS and BH–BH merger rates can be accurately estimated. That being said, the observational constraints reviewed in Section 5.2 point to a high NS–NS merger rate and will therefore inevitably rule out the few models for which NS–NS merger rates could not be accurately determined.

BINARY COMPACT OBJECT COALESCENCE RATES: THE ROLE OF ELLIPTICAL GALAXIES

Article metrics

Permissions

Author e-mails

Author affiliations

Author notes

Dates

ABSTRACT

1. INTRODUCTION

2. GRAVITATIONAL-WAVE DETECTION

2.1. Rates for Advanced Detectors

3. MULTICOMPONENT STAR FORMATION HISTORY

4. CBC POPULATION SYNTHESIS MODELS

4.1. Model Uncertainty and Recent Literature

4.2. Properties of All CBC Binaries

4.3. Varying the Binary Fraction

4.4. Uncertainties and Smoothing: Sampling Errors

4.5. Model Prediction for CBC Rate Densities

4.6. Discussion

5. CBC DETECTION RATES AND DETECTION PROBABILITY

5.1. CBC Detection Rates for LIGO

5.2. Comparisons with Pulsar Observations of Galactic NS–NS Binaries

5.3. Net Probability of CBC Detection

5.4. Comparison to Earlier Studies

6. CONCLUSIONS AND DISCUSSION

6.1. Anticipated Constraints from Future Observations

APPENDIX A: RECONSTRUCTING THE MASS AND TIME DISTRIBUTIONS

APPENDIX B: SAMPLING LIMITATIONS AND ARCHIVE SELECTION

Footnotes

BINARY COMPACT OBJECT COALESCENCE RATES: THE ROLE OF ELLIPTICAL GALAXIES

Article metrics

Permissions

Share this article

Author e-mails

Author affiliations

Author notes

Dates

ABSTRACT

1. INTRODUCTION

2. GRAVITATIONAL-WAVE DETECTION

2.1. Rates for Advanced Detectors

3. MULTICOMPONENT STAR FORMATION HISTORY

4. CBC POPULATION SYNTHESIS MODELS

4.1. Model Uncertainty and Recent Literature

4.2. Properties of All CBC Binaries

4.3. Varying the Binary Fraction

4.4. Uncertainties and Smoothing: Sampling Errors

4.5. Model Prediction for CBC Rate Densities

4.6. Discussion

5. CBC DETECTION RATES AND DETECTION PROBABILITY

5.1. CBC Detection Rates for LIGO

5.2. Comparisons with Pulsar Observations of Galactic NS–NS Binaries

5.3. Net Probability of CBC Detection

5.4. Comparison to Earlier Studies

6. CONCLUSIONS AND DISCUSSION

6.1. Anticipated Constraints from Future Observations

APPENDIX A: RECONSTRUCTING THE MASS AND TIME DISTRIBUTIONS

APPENDIX B: SAMPLING LIMITATIONS AND ARCHIVE SELECTION

Footnotes