of 15
Draft version September 21, 2016
Preprint typeset using L
A
T
E
X style AASTeX6 v. 1.0
THE RATE OF BINARY BLACK HOLE MERGERS INFERRED FROM ADVANCED LIGO OBSERVATIONS
SURROUNDING GW150914
B. P. Abbott,
1
Deceased, May 2015.
Deceased, March 2015.
(LIGO Scientific Collaboration and Virgo Collaboration)
(and The LIGO Scientific Collaboration and Virgo Collaboration)
1
LIGO, California Institute of Technology, Pasadena, CA 91125, USA
ABSTRACT
A transient gravitational-wave signal, GW150914, was identified in the twin Advanced LIGO detectors
on September 14, 2015 at 09:50:45 UTC. To assess the implications of this discovery, the detectors
remained in operation with unchanged configurations over a period of 39 d around the time of the
signal. At the detection statistic threshold corresponding to that observed for GW150914, our search
of the 16 days of simultaneous two-detector observational data is estimated to have a false alarm
rate (FAR) of
<
4
.
9
×
10
6
yr
1
, yielding a
p
-value for GW150914 of
<
2
×
10
7
. Parameter estimation
followup on this trigger identifies its source as a binary black hole (BBH) merger with component
masses (
m
1
,m
2
) =
(
36
+5
4
,
29
+4
4
)
M
at redshift
z
= 0
.
09
+0
.
03
0
.
04
(median and 90% credible range). Here
we report on the constraints these observations place on the rate of BBH coalescences. Considering
only GW150914, assuming that all BBHs in the Universe have the same masses and spins as this
event, imposing a search FAR threshold of 1 per 100 years, and assuming that the BBH merger rate is
constant in the comoving frame, we infer a 90% credible range of merger rates between 2–53 Gpc
3
yr
1
(comoving frame). Incorporating all search triggers that pass a much lower threshold while accounting
for the uncertainty in the astrophysical origin of each trigger, we estimate a higher rate, ranging from
13–600 Gpc
3
yr
1
depending on assumptions about the BBH mass distribution. All together, our
various rate estimates fall in the conservative range 2–600 Gpc
3
yr
1
.
1.
INTRODUCTION
The first detection of a gravitational wave (GW) sig-
nal in the twin Advanced LIGO detectors on September
14, 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 d 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
<
4
.
9
×
10
6
yr
1
, yielding a
p
-value
for GW150914 of
<
2
×
10
7
(Abbott et al. 2016c).
GW150914 is consistent with a gravitational-wave sig-
nal from the merger of two black holes with masses
(
m
1
,m
2
) =
(
36
+5
4
,
29
+4
4
)
M
at redshift
z
= 0
.
09
+0
.
03
0
.
04
(Abbott et al. 2016e). Here and throughout, we report
posterior medians and 90% symmetric credible intervals.
In this paper, we discuss inferences on the rate of binary
black hole (BBH) mergers from this detection and the
surrounding data. This Letter is accompanied by Ab-
bott et al. (2016) (hereafter the Supplement) containing
supplementary information on our methods and compu-
tations.
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–300 Gpc
3
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 con-
straining electromagnetic observations produce a wide
range of rate estimates.
Observations of GWs can
tightly constrain this rate with minimal modeling as-
sumptions, and thus provide useful input on the astro-
physics 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 con-
strain the model space of BBH merger rates.
It is possible to obtain a rough estimate of the BBH co-
arXiv:1602.03842v3 [astro-ph.HE] 20 Sep 2016
2
alescence 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 de-
pend 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 de-
tected 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 in-
ferred posterior median rate and 90% credible range is
R
100
= 14
+39
12
Gpc
3
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 signifi-
cance of the second-most-significant trigger (designated
LVT151012), our search FAR is 0
.
43 yr
1
, yielding a
p
-
value for this trigger of 0
.
02. This trigger cannot confi-
dently be claimed as a detection on the basis of such a
p
-value, but neither is it obviously consistent with a ter-
restrial origin, i.e. a result of either instrumental or envi-
ronmental 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
+18
6
,
13
+4
5
)
M
at redshift
z
= 0
.
21
+0
.
09
0
.
09
(Abbott et al. 2016c). Based on two different imple-
mentations of a matched-filter search, we find posterior
probabilities 0
.
84 and 0
.
91 that LVT151012 is of astro-
physical origin (see Section 2.1). This is the only trig-
ger 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 un-
certain 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 (GRBs). A subsequent paper used similar mod-
els 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) addresses 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 terres-
trial (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 con-
stant in comoving volume and source-frame time, and
making various assumptions about the mass distribu-
tion of merging BBH systems as described in Sections
2.2 and 3, we derive merger rates that lie in the range
13–600 Gpc
3
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 giving rates corre-
spond to two different search algorithms (called
pycbc
and
gstlal
, described in the Supplement) with different
models of the astrophysical and terrestrial trigger distri-
butions. Because the rate posteriors from the different
searches are essentially identical (see Figures 1 and 2),
the third column giving rates provides a combined esti-
mate that results from an average of the posterior densi-
ties from each search. Including the rate estimate with a
strict threshold that considers only the GW150914 trig-
ger as described in Section 2.2 all our rate estimates lie
in the conservative range 2–600 Gpc
3
yr
1
.
1
All our rate estimates are consistent within their sta-
tistical uncertainties, and these estimates are also con-
sistent with the broad range of rate predictions re-
viewed in Abadie et al. (2010) with only the low end
(
<
1 Gpc
3
yr
1
) of rate predictions being excluded.
The astrophysical implications of the GW150914 detec-
tion and these inferred rates are further discussed in
1
Following submission but before acceptance of this paper we
identified a mistake in our calculation of the sensitive spacetime
volume for the “Flat” and “Power Law” BBH populations (see
Section 3) that reduced those volumes and increased the corre-
sponding rates by a factor of approximately two. Since the upper
limit of this rate range is driven by the rate estimates for the
“Power Law” population, the range given here increased when
the mistake was corrected. Previous versions of this paper posted
to the arXiv, Abbott et al. (2016d), Abbott et al. (2016a), and
others used the mistaken rate range 2–400 Gpc
3
yr
1
. The cor-
rection does not affect the astrophysical interpretation appearing
in Abbott et al. (2016d) or Abbott et al. (2016a).
3
Abbott et al. (2016a).
Table 1
. Rates of BBH mergers estimated under various
assumptions. See Section 2.2. All results are reported as a
posterior median and 90% symmetric credible interval.
Mass Distribution
R/
(
Gpc
3
yr
1
)
pycbc gstlal
Combined
GW150914
16
+38
13
17
+39
14
17
+39
13
LVT151012
61
+152
53
62
+164
55
62
+165
54
Both
82
+155
61
84
+172
64
83
+168
63
Astrophysical
Flat in log mass
63
+121
49
60
+122
48
61
+124
48
Power Law (-2.35) 200
+390
160
200
+410
160
200
+400
160
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 for-
mation and evolution in the Universe.
2.
RATE INFERENCE
A rate estimate requires counting the number of sig-
nals 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 trig-
ger 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 (Us-
man et al. 2015) and the
gstlal
pipeline (Messick et al.
2016) perform matched-filter searches for CBC signals
using aligned-spin templates (Taracchini et al. 2014;
P ̈urrer 2015) when searching the BBH parts of the CBC
parameter space. In these searches, single-detector trig-
gers are recorded at maxima of the signal-to-noise ra-
tio (SNR) 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 astrophys-
ical and those whose origin is terrestrial. Terrestrial
triggers are the result of either instrumental or environ-
mental 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 ap-
pear above threshold in a Poisson process with number
density in detection space
d
N
d
x
= Λ
1
p
1
(
x
) + Λ
0
p
0
(
x
)
,
(1)
where the subscripts “1” and “0” refer to the astrophys-
ical and terrestrial origin, Λ
1
and Λ
0
are the Poisson
mean numbers of triggers of astrophysical and terres-
trial type, and
p
1
and
p
0
are the (normalized) density of
triggers of astrophysical and terrestrial origin over de-
tection space. We estimate the densities,
p
0
and
p
1
, of
triggers of the two types empirically as described in the
Supplement and in Abbott et al. (2016c). Here we ig-
nore the time of arrival of the triggers in our data set,
averaging the rates of each type of trigger and the sensi-
tivity of the detector to astrophysical signals over time.
We do this because it is difficult to estimate
p
0
and
p
1
over short times and because we see no evidence of time
variation in
p
0
and
p
1
; for more details see the Supple-
ment.
The parameter Λ
1
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, Λ
1
is related
to the astrophysical rate of BBH coalescences
R
by
Λ
1
=
R
V T
,
(2)
where
V T
is the time- and population-averaged space-
time volume to which the detector is sensitive at the
chosen search threshold, defined in Eq. (15). Because
the astrophysical rate enters the likelihood only in the
combination
R
V T
, 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 statis-
tics
{
x
j
|
j
= 1
,...,M
}
is (Loredo & Wasserman 1995;
4
Farr et al. 2015)
L
(
{
x
j
|
j
= 1
,...,M
}|
Λ
1
,
Λ
0
)
=
M
j
=1
1
p
1
(
x
j
) + Λ
0
p
0
(
x
j
)]
exp [
Λ
1
Λ
0
]
.
(3)
See the Supplement for a derivation of this likeli-
hood function for our Poisson mixture model. In each
pipeline, the
shape
of the astrophysical trigger distribu-
tion
p
1
(
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
summarised by a trigger with detection statistic
x
. For
example, to obtain information about the source asso-
ciated 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
2
, 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, we here adopt this simpler method.
In the future, as detections accumulate, we expect to
transition to a method of analysis that incorporates es-
timation of population parameters. Here we deal with
the uncertainty in the astrophysical population by esti-
mating rates under several different assumptions about
the population; see Section 2.2.
We impose a prior on the Λ parameters of:
p
1
,
Λ
0
)
1
Λ
1
1
Λ
0
.
(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 Eq.
(3) and the prior from Eq. (4):
p
1
,
Λ
0
|{
x
j
|
j
= 1
,...,M
}
)
M
j
=1
1
p
1
(
x
j
) + Λ
0
p
0
(
x
j
)]
×
exp [
Λ
1
Λ
0
]
1
Λ
1
Λ
0
.
(5)
Posterior distributions for Λ
0
and Λ
1
were obtained us-
ing a Markov-chain Monte Carlo; details are given in the
Supplement, along with the resulting expected counts
Λ
1
.
2
Actually three, as this article goes to press (Abbott et al.
2016; The LIGO Scientific Collaboration et al. 2016).
Using the posterior on the Λ
0
and Λ
1
, 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)
P
1
(
x
|
Λ
0
,
Λ
1
) =
Λ
1
p
1
(
x
)
Λ
1
p
1
(
x
) + Λ
0
p
0
(
x
)
.
(6)
Marginalizing over the posterior for the expected counts
gives
P
1
(
x
|{
x
j
|
j
= 1
,...,M
}
)
0
1
P
1
(
x
|
Λ
0
,
Λ
1
)
×
p
1
,
Λ
0
|{
x
j
|
j
= 1
,...,M
}
)
,
(7)
which is the posterior probability that an event at de-
tection 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 ex-
ploring 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 estima-
tion in this region, so the probability estimate is sensi-
tive to how the background density,
p
0
, is analytically
extended into this region. We estimate that the prob-
ability of astrophysical origin for GW150914 is larger
than 1
10
6
.
Under
the
assumption
that
GW150914
and
LVT151012
are
astrophysical,
posterior
distri-
butions for system parameters can be derived
(Veitch et al. 2015).
Both triggers are con-
sistent with BBH merger sources with masses
(
m
1
,m
2
) =
(
36
+5
4
,
29
+4
4
)
M
at redshift 0
.
09
+0
.
03
0
.
04
(GW150914) and (
m
1
,m
2
) =
(
23
+18
6
,
13
+4
5
)
M
at
redshift 0
.
21
+0
.
09
0
.
09
(LVT151012) (Abbott et al. 2016e,c).
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, trig-
gers 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 pro-
5
cess with three terms:
d
N
d
x
= Λ
1
p
1
(
x
) + Λ
2
p
2
(
x
) + Λ
0
p
0
(
x
)
,
(8)
where Λ
1
, Λ
2
, and Λ
0
are the mean number of triggers
of each type in the data set, and
p
1
,
p
2
, and
p
0
are
the probability densities for triggers of each type over
the detection statistic. The shape of the distribution of
SNRs is independent of the event properties, so
p
1
=
p
2
(Schutz 2011; Chen & Holz 2014); we cannot distinguish
BBH classes based only on their detection statistic dis-
tributions, but rather require PE.
When an event’s parameters are known to come from
a certain class
i
, under the astrophysical origin assump-
tion, then the trigger rate becomes
d
N
i
d
x
= Λ
i
p
i
(
x
) + Λ
0
p
0
(
x
)
,
(9)
i.e. we permit the event to belong to either its astro-
physical 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
,...,M
}
, exceeding our de-
tection statistic threshold is similar to Eq. (3), but we
now account for the distinct classification of GW150914
and LVT151012 based on PE:
L
(
{
x
j
|
j
= 1
,...,M
}|
Λ
1
,
Λ
2
,
Λ
0
)
= [Λ
1
p
1
(
x
1
) + Λ
0
p
0
(
x
1
)] [Λ
2
p
2
(
x
2
) + Λ
0
p
0
(
x
2
)]
×
M
j
=3
1
p
1
(
x
j
) + Λ
2
p
2
(
x
j
) + Λ
0
p
0
(
x
j
)]
×
exp [
Λ
1
Λ
2
Λ
0
]
.
(10)
The first two terms in this product are the rates of the
form of Eq. (9) for the GW150914 and LVT151012 trig-
gers, whose class, if not terrestrial, is known; the re-
maining 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 correspond-
ing events by
Λ
i
=
R
i
V T
i
,
(11)
where
V T
i
is the time- and population-averaged space-
time volume to which the detector is sensitive for event
class
i
, defined in Eq. (15) under the population assump-
tion in Eq. (16).
We impose a prior for the total astrophysical and ter-
restrial expected counts:
p
1
,
Λ
2
,
Λ
0
)
1
Λ
1
+ Λ
2
1
Λ
0
(12)
This prior is chosen to match the prior in Eq. (4). It is
adequate for the two events that we are analysing here,
but should be modified if a large number of events are
being analysed 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 ex-
pected counts given the trigger set is proportional to the
product of likelihood, Eq. (10), and prior, Eq. (12):
p
1
,
Λ
2
,
Λ
0
|{
x
j
}
)
∝L
(
{
x
j
}|
Λ
1
,
Λ
2
,
Λ
0
)
p
1
,
Λ
2
,
Λ
0
) (13)
We again use Markov-chain Monte Carlo samplers to
obtain resulting expected counts for Λ
1
, Λ
2
, and Λ
Λ
1
+ Λ
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 es-
timates, which are given in the Supplement, are consis-
tent 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 ex-
pected 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 func-
tion, 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 de-
scribe 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
1
(
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 calcu-
lations (Loredo & Wasserman 1995, 1998a,b; Gladman
et al. 1998), where information about the source prop-
erties is contained in each trigger. Under these assump-
tions, a count at a chosen threshold Λ
i
is related to an
astrophysical rate
R
i
by
Λ
i
=
R
i
V T
i
,
(14)
where
V T
i
=
T
d
z
d
θ
d
V
c
d
z
1
1 +
z
s
i
(
θ
)
f
(
z,θ
)
(15)
(see Supplement). Here
R
i
is the space-time rate den-
sity in the comoving frame, 0
f
(
z,θ
)
1 is the se-
lection function,
T
is the total observation time in the
observer frame, and
V
c
(
z
) is the comoving volume con-
6
tained within a sphere out to redshift
z
(Hogg 1999).
3
In other words, the posterior on
R
i
is obtained by sub-
stituting Eq. (14) into Eq. (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 popula-
tion follows the observed sources:
s
i
(
θ
) =
δ
(
θ
θ
i
)
,
(16)
where
δ
is the Dirac delta function and
θ
i
are the param-
eters of source type
i
. Because of the finite SNR of the
events, we do not know these parameters perfectly; we
marginalize over our imperfect knowledge by integrat-
ing over the PE posterior for the intrinsic, source-frame
parameters from the followup on each trigger. This ef-
fectively replaces the Dirac delta in Eq. (16) by the PE
posterior distribution in the integral in Eq. (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 prob-
ability (FAP), and a correspondingly restrictive selec-
tion 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 Eq. (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. (2016e,c), random orientations and sky locations,
and a fixed redshift distribution.
4
The waveforms cor-
respond 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 ̈urrer 2015); in nature we would never expect
perfect spin alignment, but nevertheless the EOB wave-
forms are a good fit to the observed data (Abbott et al.
2016e) and we therefore expect they will accurately rep-
resent 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 injec-
tions appropriately, we can estimate the integral in
Eq. (15), and, accounting for the effect on the recov-
3
Throughout this paper, we use the “TT+lowP+lensing+ext”
cosmological parameters from Table 4 of Planck Collaboration
et al. (2015).
4
The source- and observer-frame masses are related by a red-
shift factor, (1 +
z
)
M
source
=
M
observer
.
ered luminosity distance from a 10% amplitude cali-
bration uncertainty (Abbott et al. 2016b), we obtain
V T
100
= 0
.
082
+0
.
053
0
.
032
Gpc
3
yr. 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
'
0 (see Sec-
tion 2.1), thus the posterior on the associated rate,
R
100
becomes:
p
(
R
100
|
GW150914)
R
100
V T
100
exp [
R
100
V T
100
]
,
(17)
from which we infer
R
100
= 14
+39
12
Gpc
3
yr
1
.
2.2.2.
Rates Incorporating All Triggers
As discussed in Section 2.1, there is useful informa-
tion 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 trig-
gers of terrestrial origin. As before, we perform a Monte
Carlo estimation of the integral in Eq. (15) using pos-
terior distributions from the PE of both the GW150914
and LVT151012 described in Abbott et al. (2016e,c), 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
1
,
R
2
, and
R
R
1
+
R
2
from our estimates of
V T
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
1
and
R
2
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.
3.
SENSITIVITY TO ASTROPHYSICAL MASS
DISTRIBUTION
The assumptions in the Kim et al. (2003) method
about the distribution of intrinsic BBH population pa-
rameters 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 sec-
tion on two additional estimates of the rate using dif-
ferent source distributions
s
(
θ
) that bracket possible as-
trophysical scenarios.
7
10
0
10
1
10
2
10
3
R
(
Gpc
3
yr
1
)
0
.
0
0
.
1
0
.
2
0
.
3
0
.
4
0
.
5
0
.
6
Rp
(
R
)
R
R
1
R
2
Figure 1
. The posterior density on the rate of GW150914-
like BBH inspirals,
R
1
(green), LVT151012-like BBH inspi-
rals,
R
2
(red), and the inferred total rate,
R
=
R
1
+
R
2
(blue). The median and 90% credible levels are given in Ta-
ble 1. Solid lines give the rate inferred from the
pycbc
trigger
set, while dashed lines give the rate inferred from the
gstlal
trigger set.
The first source distribution we take to have spins
aligned with the orbital angular momentum, with mag-
nitude uniform in
0
.
99
(
a/m
)
1
,
2
0
.
99 and masses
flat in log (
m
1
) and log (
m
2
),
s
(
θ
)
1
m
1
1
m
2
,
(18)
with
m
1
,m
2
5 M
and
m
1
+
m
2
100 M
. The spin
distribution of merging BHs is very uncertain, but is un-
likely to be concentrated at
a
= 0, so a uniform distribu-
tion 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; Fryer et al. 2012; Dominik et al. 2012;
Spera et al. 2015). Coalescences of higher-mass BHs
from the 5 M
–100 M
range produce higher signal-to-
noise ratios 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
for the lower mass limit be-
cause 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 (
̈
Ozel et al. 2010; Farr et al. 2011); but see
Kreidberg et al. (2012) for an alternative explanation for
the dearth of BH mass estimates below
5 M
. Using
an injection campaign as described above, and incorpo-
rating 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,
5
p
(
m
1
)
m
2
.
35
1
,
(19)
with the smaller mass distributed uniformly in
q
m
2
/m
1
, and with
m
1
,m
2
5 M
and
m
1
+
m
2
100 M
. The results of using this distribution in an in-
jection campaign are given in the Supplement. This dis-
tribution likely produces more low-mass BHs than the
true astrophysical distribution, and therefore the sen-
sitive time-volume is probably smaller than would be
obtained with the true distribution; the estimated rate
is correspondingly higher (Fryer & Kalogera 2001; Fryer
et al. 2012; Dominik 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 as-
trophysical and a terrestrial population, as in the anal-
ysis at the beginning of Section 2.1 (see Eq. (5)). We
relate expected counts Λ
1
to rates via Eq. (14), with
the
V T
for the astrophysical distributions given in the
Supplement. We find
R
flat
= 61
+124
48
Gpc
3
yr
1
and
R
pl
= 200
+400
160
Gpc
3
yr
1
. Posteriors on the rates, to-
gether 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.
4.
DISCUSSION
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 as-
sumptions 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 Gpc
3
yr
1
for GW150914-
like systems and 420 Gpc
3
yr
1
for LVT151012-like sys-
tems; compared to
R
1
= 17
+39
13
Gpc
3
yr
1
and
R
2
=
62
+165
54
Gpc
3
yr
1
, it is clear that the sensitive time-
volume reach of Advanced LIGO, even from only 16 days
5
The power chosen here is the same as the Salpeter initial
mass function (Salpeter 1955), but this should not be understood
to suggest that the distribution of the more massive BH in a bi-
nary would follow the IMF; the initial mass–final mass relation
for massive stars is complicated and nonlinear (Fryer & Kalogera
2001; Fryer et al. 2012; Dominik et al. 2012; Spera et al. 2015).
Instead, as described in the text, this distribution is designed to
provide a reasonable lower-limit for the sensitive time-volume and
upper limit for the rate.
8
10
0
10
1
10
2
10
3
R
(
Gpc
3
yr
1
)
0
.
0
0
.
1
0
.
2
0
.
3
0
.
4
0
.
5
0
.
6
Rp
(
R
)
Flat
Reference
Power Law
Figure 2
. Sensitivity of the inferred BBH coalescence rate to
the assumed astrophysical distribution of BBH masses. The
curves represent the posterior assuming that BBH masses
are flat in log (
m
1
)–log (
m
2
), as in Eq. (18) (green; “Flat”),
are exactly GW150914-like or LVT151012-like as described
in Section 2 (blue; “Reference”), or are distributed as in
Eq. (19) (red; “Power Law”). The
pycbc
results are shown
in solid lines and the
gstlal
results are shown in dotted
lines. Though the searches differ in their models of the as-
trophysical and terrestrial triggers, the rates inferred from
each search are very similar. The posterior median rates and
symmetric 90% CL intervals are given in Table 1. Compar-
ing to the total rate computed using the assumptions in Kim
et al. (2003) in Section 2 we see that the rate can change by
a factor of a few depending on the assumed BBH popula-
tion. In spite of this seemingly-large variation, all three rate
posterior distributions are consistent within our statistical
uncertainties.
of coincident observations, is vastly larger than that of
any previous gravitational-wave observations.
The search thresholds used in our analysis are much
smaller than those required to produce a confident detec-
tion. We estimate that a fraction 0
.
49 of the events ex-
ceeding 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 ex-
periment with sensitive time-volume
V T
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
V T
will fol-
low a Poisson distribution with mean
Λ
= 0
.
49Λ
V T
V T
0
.
(20)
We plot the median value for Λ
obtained in this way, as
well as the 90% credible interval, as a function of sur-
veyed time-volume in the left panel of Figure 3. There
is, unsurprisingly, a wide range of reasonable possibili-
ties 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
V T
is approximately 1.5 times
V T
0
. As a point of
reference, we show the expected value of
V T
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
, assuming an observation of 6
months for O2 and 9 months for O3 with the same co-
incident duty cycle as during the first 39 d of O1. We
find estimates of
V T
O2
/
V T
0
between 7 and 20, and
V T
O3
/
V T
0
between 30 and 70. We also show
V T
for the recently-completed O1 observing run, which sur-
veyed approximately three times the spacetime volume
discussed here. A paper describing rate estimates us-
ing this methodology from the full O1 BBH search is in
preparation.
Conditional on the count of loud events, Λ
, we can
compute the probability of having more than
n
high-
significance events in a subsequent observation:
p
(
N > n
|
Λ
) = exp [
Λ
]
k
=
n
+1
Λ
k
k
!
.
(21)
Applying Eq. (20), and integrating over our posterior on
Λ from the analysis in Section 2.1, we obtain the poste-
rior probability of more than
n
high-significance events
in a subsequent observation with sensitivity
V T
given
our current observations:
p
(
N > n
|{
x
j
}
,
V T
) =
1
p
(
N > n
|
Λ
1
,
V T
))
p
1
|{
x
j
}
)
.
(22)
The right panel of Figure 3 shows this probability for
various values of
n
and
V T
.
The rates presented here are consistent with the the-
oretical 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 im-
plications 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 uni-
formly in comoving volume and source time. As we accu-
mulate more BBH sources with ongoing Advanced LIGO
observing runs, we will eventually be able to test this as-
sumption. Similarly, as we accumulate more sources and
observation time, we will learn more about the mass dis-
tribution of BBH systems. This is only the beginning.