THE RATE OF BINARY BLACK HOLE MERGERS INFERRED FROM ADVANCED LIGO OBSERVATIONS SURROUNDING GW150914

The first detection of a gravitational-wave (GW) signal in the twin Advanced LIGO detectors on 2015 September 2015, 09:50:45 UTC was reported in Abbott et al. (2016d). This transient signal is designated GW150914. To assess the implications of this discovery, the detectors remained in operation with unchanged configurations over a period of 39 days around the time of the signal. At the detection statistic threshold corresponding to that observed for GW150914, the false-alarm rate (FAR) of the search of the available 16 days of coincident data is estimated to be $\lt 4.9\times {10}^{-6}\,{\mathrm{yr}}^{-1}$, yielding a p-value for GW150914 of $\lt 2\times {10}^{-7}$ (Abbott et al. 2016c). GW150914 is consistent with a GW signal from the merger of two black holes with masses $({m}_{1},{m}_{2})=({36}_{-4}^{+5},{29}_{-4}^{+4})\,{M}_{\odot }$ at redshift $z={0.09}_{-0.04}^{+0.03}$ (Abbott et al. 2016e). Here and throughout, we report posterior medians and 90% symmetric credible intervals. In this Letter, we discuss inferences on the rate of binary black hole (BBH) mergers from this detection and the surrounding data. This Letter is accompanied by Abbott et al. (2016g, hereafter the Supplement) containing supplementary information on our methods and computations.

Previous estimates of the BBH merger rate based on population modeling are reviewed in Abadie et al. (2010). The range of rates given there spans more than three orders of magnitude, from 0.1 to 300 ${\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$. The rate of BBH mergers is a crucial output from BBH population models, but theoretical uncertainty in the evolution of massive stellar binaries and a lack of constraining electromagnetic observations produce a wide range of rate estimates. Observations of GWs can tightly constrain this rate with minimal modeling assumptions, and thus provide useful input on the astrophysics of massive stellar binaries. In the absence of detections until GW150914, the most constraining rate upper limits from GW observations, as detailed in Aasi et al. (2013), lie above the model predictions. Here, for the first time, we report on GW observations that constrain the model space of BBH merger rates.

It is possible to obtain a rough estimate of the BBH coalescence rate from the GW150914 detection by setting a low search FAR threshold that eliminates other search triggers (Abbott et al. 2016c). The inferred rate will depend on the detector sensitivity to the BBH population, which strongly depends on BBH masses. However, our single detection leaves a large uncertainty in the mass distribution of merging BBH systems. Kim et al. (2003) faced a similar situation in deriving binary neutron star merger rates from the small sample of Galactic double neutron star systems. They argued that a good rate estimate follows from an approach assuming each detected system belongs to its own class, deriving merger rates for each class independently, and then adding the rates over classes to infer the overall merger rate. If we follow Kim et al. (2003), assume that all BBH mergers in the universe have the same source-frame masses and spins as GW150914, and set a nominal threshold on the search FAR of one per century—eliminating all triggers but the one associated with GW150914—then the inferred posterior median rate and 90% credible range is ${R}_{100}={14}_{-12}^{+39}\,{\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$ (see Section 2.2).

Merger rates inferred from a single highly significant trigger are sensitive to the choice of threshold. Less significant search triggers eliminated under the strict FAR threshold can also provide information about the merger rate. For example, thresholded at the significance of the second-most-significant trigger (designated LVT151012), our search FAR is $0.43\,{\mathrm{yr}}^{-1}$, yielding a p-value for this trigger of 0.02. This trigger cannot confidently be claimed as a detection on the basis of such a p-value, but neither is it obviously consistent with a terrestrial origin, i.e., a result of either instrumental or environmental effects in the detector. Under the assumption that this trigger is astrophysical in origin, parameter estimation (PE; Veitch et al. 2015) indicates that its source is also a BBH merger with source-frame masses $({m}_{1},{m}_{2})=({23}_{-6}^{+18},{13}_{-5}^{+4})\,{M}_{\odot }$ at redshift $z={0.21}_{-0.09}^{+0.09}$ (Abbott et al. 2016c). Based on two different implementations of a matched-filter search, we find posterior probabilities 0.84 and 0.91 that LVT151012 is of astrophysical origin (see Section 2.1). This is the only trigger besides GW150914 that has probability greater than 50% of being of astrophysical origin. Farr et al. (2015) presented a method by which a set of triggers of uncertain origin like this can be used to produce a rate estimate that is more accurate than that produced by considering only highly significant events.

The mixture model of Farr et al. (2015) used here is similar to other models used to estimate rates in astrophysical contexts. Loredo & Wasserman (1995, 1998b) used a similar foreground/background mixture model to infer the rate and distribution of gamma-ray bursts. A subsequent paper used similar models in a cosmological context, as we do here (Loredo & Wasserman 1998a). Guglielmetti et al. (2009) used the same sort of formalism to model ROSAT images, and it has also found use in analysis of surveys of trans-Neptunian objects (Gladman et al. 1998; Petit et al. 2008). Kelly et al. (2009) address selection effects in the presence of population models, an issue that also appears in this work. In contrast to previous analyses, here we operate in the background-dominated regime, setting a search threshold where the FAR is relatively high so that we can be confident that triggers of terrestrial (as opposed to astrophysical) origin dominate near threshold (see Section 2.1).

Incorporating our uncertainty about the astrophysical origin of all search triggers that could represent BBH signals (Abbott et al. 2016c) using the Farr et al. (2015) method, assuming that the BBH merger rate is constant in comoving volume and source-frame time, and making various assumptions about the mass distribution of merging BBH systems as described in Sections 2.2 and 3, we derive merger rates that lie in the range $13\mbox{--}600\,{\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$.

Our rate estimates are summarized in Table 1; see Section 2.2 for more information. Each row of Table 1 represents a different assumption about the BBH mass distribution. The first two columns give rates that correspond to two different search algorithms (called pycbc and gstlal, described in the Supplement) with different models of the astrophysical and terrestrial trigger distributions. Because the rate posteriors from the different searches are essentially identical (see Figures 1 and 2), the third column gives rates that provide a combined estimate that results from an average of the posterior densities from each search. Including the rate estimate with a strict threshold that considers only the GW150914 trigger as described in Section 2.2, all our rate estimates lie in the conservative range $2\mbox{--}600\,{\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$.¹³⁷

Zoom In Zoom Out Reset image size

Table 1.  Rates of BBH Mergers Estimated under Various Assumptions

Mass Distribution	$R/({\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1})$
	pycbc	gstlal	Combined
GW150914	${16}_{-13}^{+38}$	${17}_{-14}^{+39}$	${17}_{-13}^{+39}$
LVT151012	${61}_{-53}^{+152}$	${62}_{-55}^{+164}$	${62}_{-54}^{+165}$
Both	${82}_{-61}^{+155}$	${84}_{-64}^{+172}$	${83}_{-63}^{+168}$
Astrophysical
Flat in log mass	${63}_{-49}^{+121}$	${60}_{-48}^{+122}$	${61}_{-48}^{+124}$
Power Law (−2.35)	${200}_{-160}^{+390}$	${200}_{-160}^{+410}$	${200}_{-160}^{+400}$

Note.  See Section 2.2. All results are reported as a posterior median and 90% symmetric credible interval.

Download table as:  ASCIITypeset image

All our rate estimates are consistent within their statistical uncertainties, and these estimates are also consistent with the broad range of rate predictions reviewed in Abadie et al. (2010) with only the low end (<1 ${\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$) of rate predictions being excluded. The astrophysical implications of the GW150914 detection and these inferred rates are further discussed in Abbott et al. (2016a).

This Letter presents the results of our rate inference. For methodological and other details of the analysis, see the Supplement.

The results presented here depend on assumptions about the masses, spins, and cosmological distribution of sources. As GW detectors acquire additional data and their sensitivities improve, we will be able to test these assumptions and deepen our understanding of BBH formation and evolution in the universe.

A rate estimate requires counting the number of signals in an experiment and then estimating the sensitivity to a population of sources to transform the count into an astrophysical rate. Individually, the count of signals and the sensitivity will depend on specific detection and trigger generation thresholds imposed by the pipeline, but the estimated rates should not depend strongly on such thresholds. We consider various methods of counting signals, employ two distinct search pipelines, and obtain a range of broadly consistent rate estimates.

2.1. Counting Signals

Two independent pipelines searched the coincident data for signals matching a compact binary coalescence (CBC; Abbott et al. 2016c), each producing a set of coincident search triggers. Both the pycbc pipeline (Usman et al. 2016) and the gstlal pipeline (Messick et al. 2016) perform matched-filter searches for CBC signals using aligned-spin templates (Taracchini et al. 2014; Pürrer 2016) when searching the BBH parts of the CBC parameter space. In these searches, single-detector triggers are recorded at maxima of the signal-to-noise ratio (S/N) time series for each template (Allen et al. 2012); coincident search triggers are formed when pairs of triggers, one from each detector, occur in the same template with a time difference of ±15 ms or less. Our data set here consists of the set of coincident triggers returned by each search over the 16 days of coincident observations. See the Supplement for more information about the generation of triggers.

The Farr et al. (2015) framework considers two classes of coincident triggers: those whose origin is astrophysical and those whose origin is terrestrial. Terrestrial triggers are the result of either instrumental or environmental effects in the detector. The two types of sources produce triggers with different densities in the space of detection statistics, which we denote as x. We consider all triggers above a threshold chosen so that triggers of terrestrial origin dominate at the threshold. Triggers appear above threshold in a Poisson process with number density in detection space

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dx}}={{\rm{\Lambda }}}_{1}{p}_{1}(x)+{{\rm{\Lambda }}}_{0}{p}_{0}(x),\end{eqnarray} \tag{ 1 }$$

where the subscripts "1" and "0" refer to the astrophysical and terrestrial origin, Λ₁ and Λ₀ are the Poisson mean numbers of triggers of astrophysical and terrestrial type, and p₁ and p₀ are the (normalized) density of triggers of astrophysical and terrestrial origin over detection space. We estimate the densities, p₀ and p₁, of triggers of the two types empirically as described in the Supplement and in Abbott et al. (2016c). Here, we ignore the time of arrival of the triggers in our data set, averaging the rates of each type of trigger and the sensitivity of the detector to astrophysical signals over time. We do this because it is difficult to estimate p₀ and p₁ over short times and because we see no evidence of time variation in p₀ and p₁; for more details see the Supplement.

The parameter Λ₁ is the mean number of signals of astrophysical origin above the chosen threshold; it is not the mean number of signals confidently detected (see Section 4). Under the assumptions we make here of a rate that is constant in the comoving frame, Λ₁ is related to the astrophysical rate of BBH coalescences R by

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{1}=R\langle {VT}\rangle ,\end{eqnarray} \tag{ 2 }$$

where $\langle {VT}\rangle$ is the time- and population-averaged spacetime volume to which the detector is sensitive at the chosen search threshold, defined in Equation (15). Because the astrophysical rate enters the likelihood only in the combination $R\langle {VT}\rangle$, which represents a dimensionless count, we first discuss estimation of Λ in this section, and then discuss the relationship between the posterior on Λ and on the rate R in Section 2.2.

The likelihood for a trigger set with detection statistics $\{{x}_{j}| j=1,\ldots ,M\}$ is (Loredo & Wasserman 1995; Farr et al. 2015)

$$\begin{eqnarray}&&{ \mathcal L }(\{{x}_{j}| j=1,\ldots ,M\}| {{\rm{\Lambda }}}_{1},{{\rm{\Lambda }}}_{0})\\ &&\quad =\left\{\displaystyle \prod _{j=1}^{M}[{{\rm{\Lambda }}}_{1}{p}_{1}({x}_{j})+{{\rm{\Lambda }}}_{0}{p}_{0}({x}_{j})]\right\}\exp [-{{\rm{\Lambda }}}_{1}-{{\rm{\Lambda }}}_{0}].\end{eqnarray} \tag{ 3 }$$

See the Supplement for a derivation of this likelihood function for our Poisson mixture model. In each pipeline, the shape of the astrophysical trigger distribution p₁(x) is, to a very great extent, universal (Schutz 2011; Chen & Holz 2014); that is, it does not depend on the properties of the source (see Supplement). It is this remarkable property that motivates this approach to our analysis. In principle, the data from the LIGO detector contain much more information than can be summarized by a trigger with detection statistic x. For example, to obtain information about the source associated to a trigger (if any) we can follow it up with a separate PE analysis (Veitch et al. 2015). Unfortunately, with only two likely astrophysical sources,¹³⁸ the amount of information available about the distribution of source properties is minimal. Since we cannot eliminate the dominant astrophysical systematic of uncertainty about the distribution of source parameters through a more detailed analysis, here we adopt this simpler method. In the future, as detections accumulate, we expect to transition to a method of analysis that incorporates estimation of population parameters. Here, we deal with the uncertainty in the astrophysical population by estimating rates under several different assumptions about the population; see Section 2.2.

We impose a prior on the Λ parameters of

$$\begin{eqnarray}&&p({{\rm{\Lambda }}}_{1},{{\rm{\Lambda }}}_{0})\propto \displaystyle \frac{1}{\sqrt{{{\rm{\Lambda }}}_{1}}}\displaystyle \frac{1}{\sqrt{{{\rm{\Lambda }}}_{0}}}.\end{eqnarray} \tag{ 4 }$$

See Section 4 of the Supplement for a discussion of our choice of prior. The posterior on expected counts is proportional to the product of the likelihood from Equation (3) and the prior from Equation (4):

$$\begin{eqnarray}&&p({{\rm{\Lambda }}}_{1},{{\rm{\Lambda }}}_{0}| \{{x}_{j}| j=1,\ldots ,M\})\\ &&\quad \propto \left\{\displaystyle \prod _{j=1}^{M}[{{\rm{\Lambda }}}_{1}{p}_{1}({x}_{j})+{{\rm{\Lambda }}}_{0}{p}_{0}({x}_{j})]\right\}\\ &&\qquad \times \exp [-{{\rm{\Lambda }}}_{1}-{{\rm{\Lambda }}}_{0}]\displaystyle \frac{1}{\sqrt{{{\rm{\Lambda }}}_{1}{{\rm{\Lambda }}}_{0}}}.\end{eqnarray} \tag{ 5 }$$

Posterior distributions for Λ₀ and Λ₁ were obtained using a Markov Chain Monte Carlo; details are given in the Supplement, along with the resulting expected counts Λ₁.

Using the posterior on the Λ₀ and Λ₁, we can compute the posterior probability that each particular trigger comes from an astrophysical versus terrestrial source. The conditional probability that an event at detection statistic x comes from an astrophysical source is given by (Guglielmetti et al. 2009; Farr et al. 2015)

$$\begin{eqnarray}&&{P}_{1}(x| {{\rm{\Lambda }}}_{0},{{\rm{\Lambda }}}_{1})=\displaystyle \frac{{{\rm{\Lambda }}}_{1}{p}_{1}(x)}{{{\rm{\Lambda }}}_{1}{p}_{1}(x)+{{\rm{\Lambda }}}_{0}{p}_{0}(x)}.\end{eqnarray} \tag{ 6 }$$

Marginalizing over the posterior for the expected counts gives

$$\begin{eqnarray}&&{P}_{1}(x| \{{x}_{j}| j=1,\ldots ,M\})\\ &&\quad \equiv \displaystyle \int d{{\rm{\Lambda }}}_{0}\,d{{\rm{\Lambda }}}_{1}\,{P}_{1}(x| {{\rm{\Lambda }}}_{0},{{\rm{\Lambda }}}_{1})\\ &&\qquad \times p({{\rm{\Lambda }}}_{1},{{\rm{\Lambda }}}_{0}| \{{x}_{j}| j=1,\ldots ,M\}),\end{eqnarray} \tag{ 7 }$$

which is the posterior probability that an event at detection statistic x is astrophysical in origin given the observed event set (and associated count inference). In particular, we calculate the posterior probability that LVT151012 is of astrophysical origin to be 0.84 with the gstlal pipeline and 0.91 with the pycbc pipeline. These probabilities, while not high enough to claim LVT151012 as a second detection, are large enough to motivate exploring a second class of BBHs in the Kim et al. (2003) prescription.

It is more difficult to estimate the posterior probability that GW150914 is of astrophysical origin because there are no samples from the empirical background estimation in this region, so the probability estimate is sensitive to how the background density, p₀, is analytically extended into this region. We estimate that the probability of astrophysical origin for GW150914 is larger than 1–10⁻⁶.

Under the assumption that GW150914 and LVT151012 are astrophysical, posterior distributions for system parameters can be derived (Veitch et al. 2015). Both triggers are consistent with BBH merger sources with masses $({m}_{1},{m}_{2})=({36}_{-4}^{+5},{29}_{-4}^{+4})\,{M}_{\odot }$ at redshift ${0.09}_{-0.04}^{+0.03}$ (GW150914) and $({m}_{1},{m}_{2})=({23}_{-6}^{+18},{13}_{-5}^{+4})\,{M}_{\odot }$ at redshift ${0.21}_{-0.09}^{+0.09}$ (LVT151012; Abbott et al. 2016c, 2016e). Following Kim et al. (2003), we consider the second event, if astrophysical, to be a separate class of BBH from GW150914.

We can incorporate a second class of BBH merger into our mixture model in a straightforward way. Let there be two classes of BBH mergers: type 1, which are GW150914-like, and type 2, which are LVT151012-like. Our trigger set then consists of triggers of type 1, triggers of type 2, and triggers of terrestrial origin, denoted as type 0. The distribution of triggers over detection statistic, x, now follows an inhomogeneous Poisson process with three terms:

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dx}}={{\rm{\Lambda }}}_{1}{p}_{1}(x)+{{\rm{\Lambda }}}_{2}{p}_{2}(x)+{{\rm{\Lambda }}}_{0}{p}_{0}(x),\end{eqnarray} \tag{ 8 }$$

where Λ₁, Λ₂, and Λ₀ are the mean number of triggers of each type in the data set and p₁, p₂, and p₀ are the probability densities for triggers of each type over the detection statistic. The shape of the distribution of S/Ns is independent of the event properties, so p₁ = p₂ (Schutz 2011; Chen & Holz 2014); we cannot distinguish BBH classes based only on their detection statistic distributions, but rather require PE.

When an event's parameters are known to come from a certain class i, under the astrophysical origin assumption, then the trigger rate becomes

$$\begin{eqnarray}&&\displaystyle \frac{{{dN}}_{i}}{{dx}}={{\rm{\Lambda }}}_{i}{p}_{i}(x)+{{\rm{\Lambda }}}_{0}{p}_{0}(x),\end{eqnarray} \tag{ 9 }$$

i.e., we permit the event to belong to either its astrophysical class or to an terrestrial source, but not to the other astrophysical class. The Poisson likelihood for the set of M triggers, $\{{x}_{j}| j=1,\ldots ,M\}$, exceeding our detection statistic threshold is similar to Equation (3), but we now account for the distinct classification of GW150914 and LVT151012 based on PE:

$$\begin{eqnarray}&&{ \mathcal L }(\{{x}_{j}| j=1,\ldots ,M\}| {{\rm{\Lambda }}}_{1},{{\rm{\Lambda }}}_{2},{{\rm{\Lambda }}}_{0})\\ &&\quad =[{{\rm{\Lambda }}}_{1}{p}_{1}({x}_{1})+{{\rm{\Lambda }}}_{0}{p}_{0}({x}_{1})][{{\rm{\Lambda }}}_{2}{p}_{2}({x}_{2})+{{\rm{\Lambda }}}_{0}{p}_{0}({x}_{2})]\\ &&\qquad \times \left\{\displaystyle \prod _{j=3}^{M}[{{\rm{\Lambda }}}_{1}{p}_{1}({x}_{j})+{{\rm{\Lambda }}}_{2}{p}_{2}({x}_{j})+{{\rm{\Lambda }}}_{0}{p}_{0}({x}_{j})]\right\}\\ &&\qquad \times \exp [-{{\rm{\Lambda }}}_{1}-{{\rm{\Lambda }}}_{2}-{{\rm{\Lambda }}}_{0}].\end{eqnarray} \tag{ 10 }$$

The first two terms in this product are the rates of the form of Equation (9) for the GW150914 and LVT151012 triggers, whose class, if not terrestrial, is known; the remaining terms in the product over coincident triggers represent the other events, whose class is not known.

As above, the expected counts of type 1 and 2 triggers are related to the astrophysical rates of the corresponding events by

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{i}={R}_{i}{\langle {VT}\rangle }_{i},\end{eqnarray} \tag{ 11 }$$

where ${\langle {VT}\rangle }_{i}$ is the time- and population-averaged spacetime volume to which the detector is sensitive for event class i, defined in Equation (15) under the population assumption in Equation (16).

We impose a prior for the total astrophysical and terrestrial expected counts:

$$\begin{eqnarray}&&p({{\rm{\Lambda }}}_{1},{{\rm{\Lambda }}}_{2},{{\rm{\Lambda }}}_{0})\propto \displaystyle \frac{1}{\sqrt{{{\rm{\Lambda }}}_{1}+{{\rm{\Lambda }}}_{2}}}\displaystyle \frac{1}{\sqrt{{{\rm{\Lambda }}}_{0}}}.\end{eqnarray} \tag{ 12 }$$

This prior is chosen to match the prior in Equation (4). It is adequate for the two events that we are analyzing here, but should be modified if a large number of events are being analyzed in this formalism; with N categories of foreground event under this prior, the expected number of total counts becomes $N+1/2$. The posterior on expected counts given the trigger set is proportional to the product of likelihood, Equation (10), and prior, Equation (12):

$$\begin{eqnarray}&&p({{\rm{\Lambda }}}_{1},{{\rm{\Lambda }}}_{2},{{\rm{\Lambda }}}_{0}| \{{x}_{j}\})\\ &&\quad \propto { \mathcal L }(\{{x}_{j}\}| {{\rm{\Lambda }}}_{1},{{\rm{\Lambda }}}_{2},{{\rm{\Lambda }}}_{0})p({{\rm{\Lambda }}}_{1},{{\rm{\Lambda }}}_{2},{{\rm{\Lambda }}}_{0})\end{eqnarray} \tag{ 13 }$$

We again use Markov Chain Monte Carlo samplers to obtain resulting expected counts for Λ₁, Λ₂, and ${\rm{\Lambda }}\equiv {{\rm{\Lambda }}}_{1}+{{\rm{\Lambda }}}_{2}$. These parameters represent the Poisson mean number of events of type 1 (GW150914-like), type 2 (LVT151012-like), and both types over the observation, above a very low detection statistic threshold. The estimates, which are given in the Supplement, are consistent with one event of astrophysical origin (GW150914) at very high probability, a further trigger (LVT151012) with high probability, and possibly several more of each type in the set of triggers at lower significance. In the next subsection, we will describe how to turn these expected counts of events into astrophysical rates.

2.2. Rates

The crucial element in the step from expected counts to rates is to determine the sensitivity of the search. Search sensitivity is described by the selection function, which gives, as a function of source parameters, the probability of generating a trigger above the chosen threshold. Here we assume that events are uniformly distributed in comoving volume and source time and describe the distribution of the other parameters (masses, spins, orientation angles, etc., here denoted by θ) for events of type i by a distribution function s_i(θ). Because the shape of the distribution p₁(x) is universal (Schutz 2011; Chen & Holz 2014), the source population enters the likelihood only through the search sensitivity; this situation differs from previous astrophysical rate calculations (Loredo & Wasserman 1995, 1998a, 1998b; Gladman et al. 1998), where information about the source properties is contained in each trigger. Under these assumptions, a count at a chosen threshold Λ_i is related to an astrophysical rate R_i by

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{i}={R}_{i}{\langle {VT}\rangle }_{i},\end{eqnarray} \tag{ 14 }$$

where

$$\begin{eqnarray}&&{\langle {VT}\rangle }_{i}=T\int {dz}\,d\theta \,\displaystyle \frac{{{dV}}_{c}}{{dz}}\displaystyle \frac{1}{1+z}{s}_{i}(\theta )f(z,\theta )\end{eqnarray} \tag{ 15 }$$

(see the Supplement). Here, R_i is the spacetime rate density in the comoving frame, $0\leqslant f(z,\theta )\leqslant 1$ is the selection function, T is the total observation time in the observer frame, and V_c(z) is the comoving volume contained within a sphere out to redshift z (Hogg 1999).¹³⁹ In other words, the posterior on R_i is obtained by substituting Equation (14) into Equation (13). We need to know (or assume) s_i, the population distribution for events of type i, before we can turn expected counts into rates.

The Kim et al. (2003) assumption is that the population follows the observed sources:

$$\begin{eqnarray}&&{s}_{i}(\theta )=\delta (\theta -{\theta }_{i}),\end{eqnarray} \tag{ 16 }$$

where δ is the Dirac delta function and θ_i are the parameters of source type i. Because of the finite S/N of the events, we do not know these parameters perfectly; we marginalize over our imperfect knowledge by integrating over the PE posterior for the intrinsic, source-frame parameters from the follow-up on each trigger. This effectively replaces the Dirac delta in Equation (16) by the PE posterior distribution in the integral in Equation (15).

2.2.1. Rate Using GW150914 Only

It is possible to obtain a rough estimate of the BBH coalescence rate from only the GW150914 trigger using a high-significance threshold on the false-alarm probability (FAP), and a correspondingly restrictive selection function to estimate sensitive time-volumes, so that only this one event is above threshold. We impose a nominal one-per-century threshold on the search FAR. The GW150914 trigger is the only one that exceeds this threshold. We estimate the integral in Equation (14) via a Monte Carlo procedure. We add simulated BBH signals to the detector data streams with source-frame masses and spins sampled from the posterior distributions from the PE of the GW150914 trigger described in Abbott et al. (2016c, 2016e), random orientations and sky locations, and a fixed redshift distribution.¹⁴⁰ The waveforms correspond to BBH systems with spins aligned with the orbital angular momentum and are generated using the effective one body (EOB) formalism (Taracchini et al. 2014; Pürrer 2016); in nature we would never expect perfect spin alignment, but nevertheless the EOB waveforms are a good fit to the observed data (Abbott et al. 2016e) and we therefore expect they will accurately represent our true detection efficiency for sources of this type.

We search this modified data stream and record all injections found above the threshold for inclusion in our trigger sets. By weighting the recovered injections appropriately, we can estimate the integral in Equation (15), and, accounting for the effect on the recovered luminosity distance from a 10% amplitude calibration uncertainty (Abbott et al. 2016b), we obtain ${\langle {VT}\rangle }_{100}={0.082}_{-0.032}^{+0.053}\,{\mathrm{Gpc}}^{3}\,\mathrm{years}$. See the Supplement for details. Systematic uncertainties in the waveforms used for the injections and search are estimated to induce an uncertainty in the sensitive volume calculation that is much smaller than the calibration uncertainty (Aasi et al. 2013; Littenberg et al. 2013).

With such a high threshold, any trigger is virtually certain to be astrophysical in origin, so p₀ ≃ 0 (see Section 2.1), thus the posterior on the associated rate, R₁₀₀ becomes

$$\begin{eqnarray}&&p({R}_{100}| \mathrm{GW}150914)\\ &&\quad \propto \sqrt{{R}_{100}{\langle {VT}\rangle }_{100}}\exp [-{R}_{100}{\langle {VT}\rangle }_{100}],\end{eqnarray} \tag{ 17 }$$

from which we infer ${R}_{100}={14}_{-12}^{+39}\,{\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$.

2.2.2. Rates Incorporating All Triggers

As discussed in Section 2.1, there is useful information about the merger rate from triggers with FAR less significant than one per century. Following Farr et al. (2015), we set a lower acceptance threshold such that the trigger density at threshold is dominated by triggers of terrestrial origin. As before, we perform a Monte Carlo estimation of the integral in Equation (15) using posterior distributions from the PE of both the GW150914 and LVT151012 described in Abbott et al. (2016c, 2016e), but with the lower thresholds used in Section 2.1; the results are given in the Supplement.

Figure 1 shows the posterior we infer on the rates R₁, R₂, and $R\equiv {R}_{1}+{R}_{2}$ from our estimates of $\langle {VT}\rangle {}_{\mathrm{1,2}}$ and the posteriors on the expected counts from Section 2.1. Results are shown in Table 1 in the rows GW150914, LVT151012, and both. Because the two pipelines give rate estimates that are in excellent agreement with each other, we report a combined rate that gives the median and 90% symmetric credible range for a posterior that is the average of the posterior derived from each pipeline independently. Here, R₁ and R₂ are the contributions to the rate from systems of each class, and R should be interpreted as the total rate of BBH mergers in the local universe.

The assumptions in the Kim et al. (2003) method about the distribution of intrinsic BBH population parameters are strong and almost certainly unrealistic. To test the sensitivity of our rate estimate to assumptions about black hole (BH) masses, we report in this section on two additional estimates of the rate using different source distributions s(θ) that bracket possible astrophysical scenarios.

The first source distribution we take to have spins aligned with the orbital angular momentum, with magnitude uniform in $-0.99\leqslant {(a/m)}_{\mathrm{1,2}}\leqslant 0.99$ and masses flat in $\mathrm{log}({m}_{1})$ and $\mathrm{log}({m}_{2})$,

$$\begin{eqnarray}&&s(\theta )\sim \displaystyle \frac{1}{{m}_{1}}\displaystyle \frac{1}{{m}_{2}},\end{eqnarray} \tag{ 18 }$$

with m₁, m₂ ≥ 5 ${M}_{\odot }$ and ${m}_{1}+{m}_{2}\leqslant 100\,{M}_{\odot }$. The spin distribution of merging BHs is very uncertain, but is unlikely to be concentrated at a = 0, so a uniform distribution of a is a reasonable choice. The leading-order term in the GW amplitude depends only on the masses of the system, so our results are not particularly sensitive to the choice of spin distribution. This flat distribution in mass probably weights more heavily toward high-mass BHs than the true astrophysical distribution (Fryer & Kalogera 2001; Dominik et al. 2012; Fryer et al. 2012; Spera et al. 2015). Coalescences of higher-mass BHs from the 5 ${M}_{\odot }$–100 ${M}_{\odot }$ range produce higher S/Ns in the detectors at the same distance, so this time-volume estimate is probably higher than that for the true astrophysical distribution; the corresponding rate estimate is therefore probably lower than the true BBH rate. We choose 5 ${M}_{\odot }$ for the lower mass limit because it encompasses the inferred mass range from PE on LVT151012 and because there are indications of a mass gap between the heaviest neutron stars and the lightest BHs (Özel et al. 2010; Farr et al. 2011); however, see Kreidberg et al. (2012) for an alternative explanation for the dearth of BH mass estimates below ∼5 ${M}_{\odot }$. Using an injection campaign as described above, and incorporating calibration uncertainty, we estimate the sensitive time-volume for this population; the results are given in the Supplement.

The second source distribution we take to have the same spin distribution with masses following a power-law on the larger BH mass,¹⁴¹

$$\begin{eqnarray}&&p({m}_{1})\sim {m}_{1}^{-2.35},\end{eqnarray} \tag{ 19 }$$

with the smaller mass distributed uniformly in q ≡ m₂/m₁ and with m₁, m₂ ≥ 5 ${M}_{\odot }$ and ${m}_{1}+{m}_{2}\leqslant 100\,{M}_{\odot }$. The results of using this distribution in an injection campaign are given in the Supplement. This distribution likely produces more low-mass BHs than the true astrophysical distribution, and therefore the sensitive time-volume is probably smaller than would be obtained with the true distribution; the estimated rate is correspondingly higher (Fryer & Kalogera 2001; Dominik et al. 2012; Fryer et al. 2012; Spera et al. 2015).

We use the same astrophysical and terrestrial trigger densities as described in Section 2.1; we now consider all triggers to belong to only two populations, an astrophysical and a terrestrial population, as in the analysis at the beginning of Section 2.1 (see Equation (5)). We relate expected counts Λ₁ to rates via Equation (14), with the $\langle {VT}\rangle$ for the astrophysical distributions given in the Supplement. We find ${R}_{\mathrm{flat}}\,={61}_{-48}^{+124}\,{\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$ and ${R}_{\mathrm{pl}}={200}_{-160}^{+400}\,{\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$. Posteriors on the rates, together with the reference BBH coalescence rate R from Section 2, appear in Figure 2. A summary of the various inferred rates appears in Table 1.

Zoom In Zoom Out Reset image size

In the absence of clear detections, previous LIGO–Virgo observing runs have yielded merger rate upper limits (Aasi et al. 2013). Even the most optimistic assumptions about the BBH distribution from Section 3 imply rates that are significantly below the rate upper limits for the same distribution of masses implied by the results of Aasi et al. (2013). For the rate estimates from Section 2.2, the corresponding upper limits from Aasi et al. (2013) are $140\,{\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$ for GW150914-like systems and $420\,{\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$ for LVT151012-like systems; compared to ${R}_{1}={17}_{-13}^{+39}\,{\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$ and ${R}_{2}={62}_{-54}^{+165}\,{\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$, it is clear that the sensitive time-volume reach of Advanced LIGO, even from only 16 days of coincident observations, is vastly larger than that of any previous GW observations.

The search thresholds used in our analysis are much smaller than those required to produce a confident detection. We estimate that a fraction 0.49 of the events exceeding our search threshold in the pycbc pipeline would also exceed a one-per-century FAR threshold. One may wonder, then, how many of these significant events we can expect to see in future observations.

For a Poisson mean occurrence number Λ in an experiment with sensitive time-volume ${\langle {VT}\rangle }_{0}$ using a high FAR (low significance) threshold, the number of triggers with FARs smaller than one per century in subsequent experiments with sensitive time-volume $\langle {VT}\rangle ^{\prime}$ will follow a Poisson distribution with mean

$$\begin{eqnarray}&&{\rm{\Lambda }}^{\prime} =0.49{\rm{\Lambda }}\displaystyle \frac{\langle {VT}\rangle ^{\prime} }{{\langle {VT}\rangle }_{0}}.\end{eqnarray} \tag{ 20 }$$

We plot the median value for Λ' obtained in this way, as well as the 90% credible interval, as a function of surveyed time-volume in the left panel of Figure 3. There is, unsurprisingly, a wide range of reasonable possibilities for the number of highly significant events in future observations. The 90% credible interval for the expected number of highly significant events lies above one once $\langle {VT}\rangle ^{\prime}$ is approximately 1.5 times ${\langle {VT}\rangle }_{0}$. As a point of reference, we show the expected value of $\langle {VT}\rangle$ for the second and third planned observing runs, O2 and O3. These volumes are calculated as in Abbott et al. (2016a), for an equal-mass binary with non-spinning components and total mass 60 ${M}_{\odot }$, assuming an observation of 6 months for O2 and 9 months for O3 with the same coincident duty cycle as during the first 39 days of O1. We find estimates of $\langle {VT}\rangle {}_{{\rm{O}}2}/{\langle {VT}\rangle }_{0}$ between 7 and 20, and $\langle {VT}\rangle {}_{{\rm{O}}3}/{\langle {VT}\rangle }_{0}$ between 30 and 70. We also show $\langle {VT}\rangle$ for the recently completed O1 observing run, which surveyed approximately three times the spacetime volume discussed here. A paper describing rate estimates using this methodology from the full O1 BBH search is in preparation.

Zoom In Zoom Out Reset image size

Conditional on the count of loud events, Λ', we can compute the probability of having more than n high-significance events in a subsequent observation:

$$\begin{eqnarray}&&p(N\gt n| {\rm{\Lambda }}^{\prime} )=\exp [-{\rm{\Lambda }}^{\prime} ]\displaystyle \sum _{k=n+1}^{\infty }\displaystyle \frac{{\rm{\Lambda }}{^{\prime} }^{k}}{k!}.\end{eqnarray} \tag{ 21 }$$

Applying Equation (20), and integrating over our posterior on Λ from the analysis in Section 2.1, we obtain the posterior probability of more than n high-significance events in a subsequent observation with sensitivity $\langle {VT}\rangle ^{\prime}$ given our current observations:

$$\begin{eqnarray}&&p(N\gt n| \{{x}_{j}\},\langle {VT}\rangle ^{\prime} )\\ &&\quad =\displaystyle \int d{{\rm{\Lambda }}}_{1}\,p(N\gt n| {\rm{\Lambda }}^{\prime} ({{\rm{\Lambda }}}_{1},\langle {VT}\rangle ^{\prime} ))p({{\rm{\Lambda }}}_{1}| \{{x}_{j}\}).\end{eqnarray} \tag{ 22 }$$

The right panel of Figure 3 shows this probability for various values of n and $\langle {VT}\rangle ^{\prime}$.

The rates presented here are consistent with the theoretical expectations detailed in Abadie et al. (2010), but rule out the lowest theoretically allowed rates. See Abbott et al. (2016a) for a detailed discussion of the implications of our rate estimates for models of the binary BH population.

GW150914 is unusually significant; only ∼8% of the astrophysical distribution of sources appearing in our search with a threshold at FARs of one per century will be more significant than GW150914. However, it is not so significant as to call into question the assumption used here that BBH coalescences are distributed uniformly in comoving volume and source time. As we accumulate more BBH sources with ongoing Advanced LIGO observing runs, we will eventually be able to test this assumption. Similarly, as we accumulate more sources and observation time, we will learn more about the mass distribution of BBH systems. This is only the beginning.

The authors gratefully acknowledge the support of the United States National Science Foundation (NSF) for the construction and operation of the LIGO Laboratory and Advanced LIGO as well as the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO600 detector. Additional support for Advanced LIGO was provided by the Australian Research Council. The authors gratefully acknowledge the Italian Istituto Nazionale di Fisica Nucleare (INFN), the French Centre National de la Recherche Scientifique (CNRS) and the Foundation for Fundamental Research on Matter supported by the Netherlands Organisation for Scientific Research, for the construction and operation of the Virgo detector and the creation and support of the EGO consortium. The authors also gratefully acknowledge research support from these agencies as well as by the Council of Scientific and Industrial Research of India, Department of Science and Technology, India, Science & Engineering Research Board (SERB), India, Ministry of Human Resource Development, India, the Spanish Ministerio de Economía y Competitividad, the Conselleria d'Economia i Competitivitat and Conselleria d'Educació Cultura i Universitats of the Govern de les Illes Balears, the National Science Centre of Poland, the European Commission, the Royal Society, the Scottish Funding Council, the Scottish Universities Physics Alliance, the Hungarian Scientific Research Fund (OTKA), the Lyon Institute of Origins (LIO), the National Research Foundation of Korea, Industry Canada and the Province of Ontario through the Ministry of Economic Development and Innovation, the Natural Science and Engineering Research Council Canada, Canadian Institute for Advanced Research, the Brazilian Ministry of Science, Technology, and Innovation, Russian Foundation for Basic Research, the Leverhulme Trust, the Research Corporation, Ministry of Science and Technology (MOST), Taiwan and the Kavli Foundation. The authors gratefully acknowledge the support of the NSF, STFC, MPS, INFN, CNRS, and the State of Niedersachsen/Germany for provision of computational resources. This article has been assigned the document number LIGO-P1500217.