THE RATE OF BINARY BLACK HOLE MERGERS INFERRED FROM
ADVANCED LIGO OBSERVATIONS SURROUNDING GW150914
LIGO Scienti
fi
c Collaboration and Virgo Collaboration
(
See the Supplement, Abbott et al.
2016g
, for the full list of authors.
)
Received 2016 February 12; revised 2016 September 19; accepted 2016 September 19; published 2016 November 30
ABSTRACT
A transient gravitational-wave signal, GW150914, was identi
fi
ed in the twin Advanced LIGO detectors on 2015
September 2015 at 09:50:45 UTC. To asse
ss the implications of this discovery,
the detectors remained in operation with
unchanged con
fi
gurations over a period of 39 days around the time of t
he 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 yr
61
, yielding a
p
-value for GW150914 of
<
́
-
210
7
. Parameter estimation follo
w-up on this trigger identi
fi
es its source as a binary black hole
(
BBH
)
merger
with component masses
(
)(
)
=
-
+
-
+
mm
M
,36,29
12
4
5
4
4
at redshift
=
-
+
z
0.09
0.04
0.03
(
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 yr
31
(
comoving frame
)
. Incorporating all search triggers that
pass a much lower threshold while accounting for the uncerta
inty in the astrophysical origin of each trigger, we estimate
a higher rate, ranging from
–
--
13 600 Gpc yr
31
depending on assumptions about the BBH mass distribution. All
together, our various rate estimat
es fall in the conservative range
–
--
2
600 Gpc yr
31
.
Key words:
black holes
–
gravitational waves
–
stars: massive
Supporting material:
data behind
fi
gures
1. INTRODUCTION
The
fi
rst 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 con
fi
gurations 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
<
́
--
4.9 10 yr
61
, yielding a
p
-value for GW150914 of
<
́
-
210
7
(
Abbott et al.
2016c
)
.
GW150914 is consistent with a GW signal from the merger of
two black holes with masses
(
)(
)
=
-
+
-
+
mm
M
,36,29
12
4
5
4
4
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 supple-
mentary 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
--
Gpc yr
31
. 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
fi
rst 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
14 Gpc yr
100
12
39
31
(
see Section
2.2
)
.
Merger rates inferred from a single highly signi
fi
cant trigger
are sensitive to the choice of threshold. Less signi
fi
cant search
triggers eliminated under the strict FAR threshold can also
provide information about the merger rate. For example,
thresholded at the signi
fi
cance of the second-most-signi
fi
cant
trigger
(
designated LVT151012
)
, our search FAR is
-
0
.43 yr
1
,
yielding a
p
-value for this trigger of 0.02. This trigger cannot
con
fi
dently be claimed as a detection on the basis of such a
p
-
The Astrophysical Journal Letters,
833:L1
(
8pp
)
, 2016 December 10
doi:10.3847
/
2041-8205
/
833
/
1
/
L1
© 2016. The American Astronomical Society. All rights reserved.
1
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
(
)(
)
=
-
+
-
+
mm
M
,23,13
12
6
18
5
4
at redshift
=
-
+
z
0.21
0.09
0.09
(
Abbott et al.
2016c
)
. Based on two different
implementations of a matched-
fi
lter search, we
fi
nd 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 signi
fi
cant 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
con
fi
dent 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 600 Gpc yr
31
.
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
fi
rst
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
600 Gpc yr
31
.
137
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
--
Gpc yr
31
)
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.
2. RATE INFERENCE
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 speci
fi
c 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-
fi
lter 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
()
()
()
=L +L
dN
dx
px
px
,1
1
1
0
0
where the subscripts
“
1
”
and
“
0
”
refer to the astrophysical and
terrestrial origin,
Λ
1
and
Λ
0
are the Poisson mean numbers of
137
Following submission but before acceptance of this Letter we identi
fi
ed 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 corresponding 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 yr
31
. The correction does not affect the astrophysical interpreta-
tion appearing in Abbott et al.
(
2016d
)
or Abbott et al.
(
2016a
)
.
2
The Astrophysical Journal Letters,
833:L1
(
8pp
)
, 2016 December 10
Abbott et al.
triggers of astrophysical and terrestrial type, and
p
1
and
p
0
are
the
(
normalized
)
density of triggers of astrophysical and
terrestrial origin over detection 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 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 dif
fi
cult 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
Supplement
.
The parameter
Λ
1
is the mean number of signals of
astrophysical origin above the chosen threshold; it is not the
mean number of signals con
fi
dently 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
()
L= á ñ
RVT
,2
1
where
á
ñ
VT
is the time- and population-averaged spacetime
volume to which the detector is
sensitive at the chosen search
threshold, de
fi
ned in Equation
(
15
)
. Because the astrophysical
rate enters the likelihood only in the combination
áñ
R
VT
,
which represents a dimensionless count, we
fi
rst discuss
estimation of
Λ
in this section, and then discuss the
relationship between the posterior on
Λ
andontherate
R
in
Section
2.2
.
The likelihood for a trigger set with detection statistics
{
∣}
=¼
xj
M
1, ,
j
is
(
Loredo & Wasserman
1995
; Farr
et al.
2015
)
({ ∣
}∣
)
[()
()] [
]()
=¼ LL
= L +L
-L-L
=
⎪
⎪
⎪
⎪
⎧
⎨
⎩
⎫
⎬
⎭
xj
M
px
px
1, ,
,
exp
. 3
j
j
M
jj
10
1
1
1
0
0
10
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
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 summar-
ized 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,
138
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
()
()
LLμ
LL
p
,
11
.4
10
10
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
)
:
(∣{∣
})
[()
()]
[]
()
LL = ¼
μL+L
́-L-L
LL
=
⎪
⎪
⎪
⎪
⎧
⎨
⎩
⎫
⎬
⎭
pxjM
px
px
,1,,
exp
1
.5
j
j
M
jj
10
1
1
1
0
0
10
10
Posterior distributions for
Λ
0
and
Λ
1
were obtained using a
Markov Chain Monte Carlo; details are given in the
Supplement
, along with the resulting expected counts
Λ
1
.
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 prob-
ability that an event at detection statistic
x
comes from an
astrophysical source is given by
(
Guglielmetti et al.
2009
; Farr
et al.
2015
)
(∣
)
()
()
()
()
LL=
L
L+L
Px
px
px
px
,.6
101
1
1
1
1
0
0
Marginalizing over the posterior for the expected counts gives
(∣{ ∣
})
(∣
)
(∣{∣
})
()
ò
=¼
ºLL LL
́LL =¼
Px xj
M
ddPx
pxjM
1, ,
,
,1,,, 7
j
j
1
011 01
10
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 dif
fi
cult 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
0
, is analytically extended into this
region. We estimate that the probability of astrophysical origin
for GW150914 is larger than 1
–
10
−
6
.
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
(
)(
)
=
-
+
-
+
mm
M
,36,29
12
4
5
4
4
at redshift
-
+
0
.09
0.04
0.03
(
GW150914
)
and
(
)(
)
=
-
+
-
+
mm
M
,23,13
12
6
18
5
4
at redshift
-
+
0
.21
0.09
0.0
9
(
LVT151012; Abbott et al.
2016c
,
2016e
)
. Following Kim
138
Actually three, as this article goes to press
(
Abbott et al.
2016f
; The LIGO
Scienti
fi
c Collaboration et al.
2016
)
.
3
The Astrophysical Journal Letters,
833:L1
(
8pp
)
, 2016 December 10
Abbott et al.
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:
()
()
()
()
=L +L +L
dN
dx
px
p x
px
,8
1
1
2
2
0
0
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 S
/
Ns 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 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
()
()
()
=L +L
dN
dx
px
p x
,9
i
i
i
0
0
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,
{
∣}
=¼
xj
M
1, ,
j
, exceeding our detection statistic
threshold is similar to Equation
(
3
)
, but we now account for the
distinct classi
fi
cation of GW150914 and LVT151012 based on
PE:
({ ∣
}∣
)
[
()
()][
()
()]
[()
()
()]
[]
()
=¼ LLL
=L +L L +L
́ L +L +L
́-L-L-L
=
⎪
⎪
⎪
⎪
⎧
⎨
⎩
⎫
⎬
⎭
xj
M
px
px
px
px
px
p x
px
1, ,
, ,
exp
.
10
j
j
M
jjj
12 0
1
1
10
0
12
2
20
0
2
3
1
1
2
2
0
0
12 0
The
fi
rst 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
()
L= á ñ
RVT
,11
ii i
where
á
ñ
VT
i
is the time- and population-averaged spacetime
volume to which the detector is sensitive for event class
i
,
de
fi
ned in Equation
(
15
)
under the population assumption in
Equation
(
16
)
.
We impose a prior for the total astrophysical and terrestrial
expected counts:
()
()
LLLμ
L+L L
p
,,
11
.12
12 0
120
This prior is chosen to match the prior in Equation
(
4
)
.Itis
adequate for the two events that we are analyzing here, but
should be modi
fi
ed 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
)
:
(∣{})
({}∣)()()
LLL
μLLLLLL
px
xp
,,
,,
,,
13
j
j
12 0
12 0 12 0
We again use Markov Chain Monte Carlo samplers to obtain
resulting expected counts for
Λ
1
,
Λ
2
,and
LºL+L
12
.These
parameters represent the Poisson mean number of events of type
1
(
GW150914-like
)
,type2
(
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 signi
fi
cance. 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
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 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
()
L= á ñ
RVT
,14
ii i
where
() ( )
( )
ò
qqq
áñ=
+
VT
T
dz d
dV
dz
z
sfz
1
1
,15
i
c
i
(
see the
Supplement
)
. Here,
R
i
is the spacetime rate density in
the comoving frame,
()
q
fz
0
,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
)
.
139
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.
139
Throughout this Letter, we use the
“
TT
+
lowP
+
lensing
+
ext
”
cosmologi-
cal parameters from Table 4 of Planck Collaboration et al.
(
2016
)
.
4
The Astrophysical Journal Letters,
833:L1
(
8pp
)
, 2016 December 10
Abbott et al.
The Kim et al.
(
2003
)
assumption is that the population
follows the observed sources:
()
(
)
( )
qdqq
=-
s
,16
ii
where
δ
is the Dirac delta function and
θ
i
are the parameters of
source type
i
. Because of the
fi
nite 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-
signi
fi
cance 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
fi
xed redshift distribution.
140
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
fi
t to the observed
data
(
Abbott et al.
2016e
)
and we therefore expect they will
accurately represent our true detection ef
fi
ciency for sources of
this type.
We search this modi
fi
ed 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
á
ñ=
-
+
VT
0.082
Gpc years
100
0.032
0.053
3
. 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 Section
2.1
)
, thus the
posterior on the associated rate,
R
100
becomes
(∣
)
[]()
μáñ -áñ
pR
RVT
RVT
GW150914
exp
,
17
100
100
100
100
100
from which we infer
=
-
+
--
R
14 Gpc yr
100
12
39
31
.
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 signi
fi
cant 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
1
,
R
2
,
and
º+
R
RR
12
from our estimates of
á
ñ
VT
1,2
and the
posteriors on the expected counts from Section
2.1
.Results
areshowninTable
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 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
fi
rst source distribution we take to have spins aligned
with the orbital angular momentum, with magnitude uniform in
()
-
am
0.99
0.99
1,2
and masses
fl
at in
()
m
log
1
and
(
)
m
log
2
,
()
( )
q
~
s
mm
11
,18
12
with
m
1
,
m
2
5
M
and
+
mm
M
100
12
. 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
fl
at
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
–
100
M
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
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
140
The source- and observer-frame masses are related by a redshift
factor,
()
+=
zM
M
1
source
observer
.
5
The Astrophysical Journal Letters,
833:L1
(
8pp
)
, 2016 December 10
Abbott et al.
below
∼
5
M
. 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,
141
()
()
~
-
pm
m
,19
1
1
2.35
with the smaller mass distributed uniformly in
q
≡
m
2
/
m
1
and
with
m
1
,
m
2
5
M
and
+
mm
M
100
12
. 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
Λ
1
to
rates via Equation
(
14
)
, with the
á
ñ
VT
for the astrophysical
distributions given in the
Supplement
.We
fi
nd
=
R
flat
-
+
--
61 Gpc yr
48
124
31
and
=
-
+
--
R
200 Gpc yr
pl
160
400
31
. 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
.
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 assumptions about the
BBH distribution from Section
3
imply rates that are signi
fi
cantly
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 yr
31
for GW150914-like systems and
--
4
20 Gpc yr
31
for LVT151012-like systems; compared to
=
-
+
--
R
17 Gpc yr
1
13
39
31
and
=
-
+
--
R
62 Gpc yr
2
54
165
31
,itisclear
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 con
fi
dent 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 signi
fi
cant events we can expect to see in future
observations.
For a Poisson mean occurrence number
Λ
in an experiment
with sensitive time-volume
á
ñ
VT
0
using a high FAR
(
low
signi
fi
cance
)
threshold, the number of triggers with FARs
smaller than one per century in subsequent experiments with
sensitive time-volume
á
ñ
¢
VT
will follow a Poisson distribution
with mean
()
L¢ = L
áñ¢
áñ
VT
VT
0.49
.
20
0
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
signi
fi
cant events in future observations. The 90% credible
interval for the expected number of highly signi
fi
cant events
lies above one once
á
ñ
¢
VT
is approximately 1.5 times
á
ñ
VT
0
.As
Table 1
Rates of BBH Mergers Estimated under Various Assumptions
Mass Distribution
(
)
--
R
Gpc yr
31
/
pycbcgstlal
Combined
GW150914
-
+
1
6
13
38
-
+
1
7
14
39
-
+
1
7
13
39
LVT151012
-
+
61
53
152
-
+
62
55
16
4
-
+
62
54
165
Both
-
+
82
61
155
-
+
84
64
172
-
+
83
63
168
Astrophysical
Flat in log mass
-
+
63
49
121
-
+
60
48
12
2
-
+
61
48
12
4
Power Law
(
−
2.35
)
-
+
2
00
160
39
0
-
+
2
00
160
41
0
-
+
2
00
160
40
0
Note.
See Section
2.2
. All results are reported as a posterior median and 90%
symmetric credible interval.
Figure 1.
Posterior density on the rate of GW150914-like BBH inspirals,
R
1
(
green
)
; LVT151012-like BBH inspirals,
R
2
(
red
)
; and the inferred total rate,
=+
R
RR
12
(
blue
)
. The median and 90% credible levels are given in Table
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 data used to create this
fi
gure are available.
)
141
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 binary would follow the IMF; the
initial mass
–
fi
nal mass relation for massive stars is complicated and nonlinear
(
Fryer & Kalogera
2001
; Dominik et al.
2012
; Fryer 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.
6
The Astrophysical Journal Letters,
833:L1
(
8pp
)
, 2016 December 10
Abbott et al.
a point of reference, we show the expected value of
á
ñ
VT
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 coincident duty cycle as during
the
fi
rst 39 days of O1. We
fi
nd estimates of
á
ñáñ
VT
VT
O2
0
between 7 and 20, and
á
ñáñ
VT
VT
O3
0
between 30 and 70. We
also show
á
ñ
VT
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.
Conditional on the count of loud events,
Λ
′
, we can compute
the probability of having more than
n
high-signi
fi
cance events
in a subsequent observation:
(∣) []
!
()
å
>L¢= -L¢
L¢
=+
¥
pN n
k
exp
.
21
kn
k
1
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-signi
fi
cance events in a
subsequent observation with sensitivity
á
ñ
¢
VT
given our current
observations:
(∣{} )
(∣(
))(∣{})()
ò
>áñ¢
=L>L¢Láñ¢L
pN n x
VT
dpN n
VT p x
,
,.22
j
j
111
The right panel of Figure
3
shows this probability for various
values of
n
and
á
ñ
¢
VT
.
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 signi
fi
cant; only
∼
8% of the
astrophysical distribution of s
ources appearing in our search
with a threshold at FARs of one per century will be more
signi
fi
cant than GW150914. However, it is not so signi
fi
cant
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 obser
ving 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 acknowl-
edge the Italian Istituto Nazionale di Fisica Nucleare
(
INFN
)
,
the French Centre National de la Recherche Scienti
fi
que
(
CNRS
)
and the Foundation for Fundamental Research on
Matter supported by the Netherlands Organisation for Scienti
fi
c
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 Scienti
fi
c 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
Scienti
fi
c 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
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
fl
at in
()
(
)
-
mm
log
log
12
, as in Equation
(
18
)
(
green;
“
Flat
”
)
, are exactly GW150914-like or LVT151012-like as described in
Section
2
(
blue;
“
Reference
”
)
, or are distributed as in Equation
(
19
)(
red;
“
Power Law
”
)
. The
pycbc
results are shown as solid lines, and the
gstlal
results are shown as dotted lines. Though the searches differ in their models of
the astrophysical and terrestrial triggers, the rates inferred from each search are
very similar. The posterior median rates and symmetric 90% credible level
(
CL
)
intervals are given in Table
1
. Comparing 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 population. In
spite of this seemingly large variation, all three rate posterior distributions are
consistent within our statistical uncertainties.
(
The data used to create this
fi
gure are available.
)
7
The Astrophysical Journal Letters,
833:L1
(
8pp
)
, 2016 December 10
Abbott et al.
(
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
.
REFERENCES
Aasi, J., Abadie, J., Abbott, B. P., et al. 2013,
PhRvD
,
87, 022002
Abadie, J., Abbott, B. P., Abbott, R., et al. 2010,
CQGra
,
27, 173001
Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016a,
ApJL
,
818, L22
Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016b, arXiv:
1602.03845
Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016c,
PhRvD
,
93, 122003
Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016d,
PhRvL
,
116, 061102
Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016e,
PhRvL
,
116, 241102
Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016f,
PhRvL
,
116, 241103
Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016g,
ApJS
, 227, 14
Allen, B., Anderson, W. G., Brady, P. R., Brown, D. A., & Creighton, J. D. E.
2012,
PhRvD
,
85, 122006
Chen, H.-Y., & Holz, D. E. 2014, arXiv:
1409.0522
Dominik, M., Belczynski, K., Fryer, C., et al. 2012,
ApJ
,
759, 52
Farr, W. M., Gair, J. R., Mandel, I., & Cutler, C. 2015,
PhRvD
,
91, 023005
Farr, W. M., Sravan, N., Cantrell, A., et al. 2011,
ApJ
,
741, 103
Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012,
ApJ
,
749, 91
Fryer, C. L., & Kalogera, V. 2001,
ApJ
,
554, 548
Gladman, B., Kavelaars, J. J., Nicholson, P. D., Loredo, T. J., & Burns, J. A.
1998,
AJ
,
116, 2042
Guglielmetti, F., Fischer, R., & Dose, V. 2009,
MNRAS
,
396, 165
Hogg, D. W. 1999, arXiv:
astro-ph
/
9905116
Kelly, B. C., Vestergaard, M., & Fan, X. 2009,
ApJ
,
692, 1388
Kim, C., Kalogera, V., & Lorimer, D. R. 2003,
ApJ
,
584, 985
Kreidberg, L., Bailyn, C. D., Farr, W. M., & Kalogera, V. 2012,
ApJ
,
757, 36
Littenberg, T. B., Baker, J. G., Buonanno, A., & Kelly, B. J. 2013,
PhRvD
,
87,
104003
Loredo, T. J., & Wasserman, I. M. 1995,
ApJS
,
96, 261
Loredo, T. J., & Wasserman, I. M. 1998a,
ApJ
,
502, 75
Loredo, T. J., & Wasserman, I. M. 1998b,
ApJ
,
502, 108
Messick, C., Blackburn, K., Brady, P., et al. 2016, arXiv:
1604.04324
Özel, F., Psaltis, D., Narayan, R., & McClintock, J. E. 2010,
ApJ
,
725, 1918
Petit, J.-M., Kavelaars, J. J., Gladman, B., & Loredo, T. 2008, in Size
Distribution of Multikilometer Transneptunian Objects, ed. M. A. Barucci
et al.
(
Tucson, AZ: Univ. Arizona Press
)
,
71
Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016,
A&A
,
594, A13
Pürrer, M. 2016,
PhRvD
,
93, 064041
Salpeter, E. E. 1955,
ApJ
,
121, 161
Schutz, B. F. 2011,
CQGra
,
28, 125023
Spera, M., Mapelli, M., & Bressan, A. 2015,
MNRAS
,
451, 4086
Taracchini, A., Buonanno, A., Pan, Y., et al. 2014,
PhRvD
,
89, 061502
The LIGO Scienti
fi
c Collaboration, The Virgo Collaboration, Abbott, B. P.,
et al. 2016,
PhRvX
,
6, 041015
Usman, S. A., Kehl, M. S., Nitz, A. H., et al. 2016,
CQGra
,
33, 215004
Veitch, J., Raymond, V., Farr, B., et al. 2015,
PhRvD
,
91, 042003
Figure 3.
Left panel: the median value and 90% credible interval for the expected number of highly signi
fi
cant events
(
FARs
<
1
/
century
)
as a function of surveyed
time-volume in an observation
(
shown as a multiple of
á
ñ
VT
0
)
. The expected range of values of
áñ
VT
for the observations in O2 and O3 are shown as vertical bands.
Right panel: the probability of observing
N
>
0
(
blue
)
,
N
>
10
(
green
)
,
N
>
35
(
red
)
, and
N
>
70
(
purple
)
highly signi
fi
cant events, as a function of surveyed time-
volume. The vertical line and bands show, from left to right, the expected sensitive time-volume for each of the O1
(
dashed line
)
, O2, and O3 observations.
(
The data used to create this
fi
gure are available.
)
8
The Astrophysical Journal Letters,
833:L1
(
8pp
)
, 2016 December 10
Abbott et al.