Draft version September 21, 2016
Preprint typeset using L
A
T
E
X style AASTeX6 v. 1.0
SUPPLEMENT: 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
Supplemental information for a Letter reporting the rate of binary black hole (BBH) coalescences in-
ferred from 16 days of coincident Advanced LIGO observations surrounding the transient gravitational
wave (GW) signal GW150914. In that work we reported various rate estimates whose 90% credible
intervals (CIs) fell in the range 2–600 Gpc
−
3
yr
−
1
. Here we give details of our method and com-
putations, including information about our search pipelines, a derivation of our likelihood function
for the analysis, a description of the astrophysical search trigger distribution expected from merging
BBHs, details on our computational methods, a description of the effects and our model for calibration
uncertainty, and an analytic method of estimating our detector sensitivity that is calibrated to our
measurements.
The first detection of a gravitational wave (GW) sig-
nal from a merging binary black hole (BBH) system is
described in Abbott et al. (2016d). Abbott et al. (2016g)
reports on inference of the local BBH merger rate from
surrounding Advanced LIGO observations. This Sup-
plement provides supporting material and methodolog-
ical details for Abbott et al. (2016g), hereafter referred
to as the Letter.
1.
SEARCH PIPELINES
Both the
pycbc
and
gstlal
pipelines are based on
matched filtering against a bank of template waveforms.
See Abbott et al. (2016c) for a detailed description of the
pipelines in operation around the time of GW150914;
here we provide an abbreviated description.
In the
pycbc
pipeline, the single-detector signal-to-
noise ratio (SNR) is re-weighted by a chi-squared fac-
tor (Allen 2005) to account for template-data mismatch
(Babak et al. 2013); the re-weighted single-detector
SNRs are combined in quadrature to produce a detec-
tion statistic for search triggers.
The
gstlal
pipeline’s detection statistic, however, is
based on a likelihood ratio (Cannon et al. 2013, 2015)
constructed from the single-detector SNRs and a signal-
consistency statistic. An analytic estimate of the distri-
bution of astrophysical signals in multiple-detector SNR
and signal consistency statistic space is compared to a
measured distribution of single-detector triggers without
a coincident counterpart in the other detector to form a
multiple-detector likelihood ratio.
Both pipelines rely on an empirical estimate of the
search background, making the assumption that trig-
gers of terrestrial origin occur independently in the two
detectors. The background estimate is built from ob-
servations of single-detector triggers over a long time
(
gstlal
) or through searching over a data stream with
one detector’s output shifted in time relative to the
other’s by an interval that is longer than the light travel
time between detectors, ensuring that no coincident as-
trophysical signals remain in the data (
pycbc
). For both
pipelines it is not possible to produce an instantaneous
background estimate at a particular time; this drives our
choice of likelihood function as described in Section 2.
The
gstlal
pipeline natively determines the functions
p
0
(
x
) and
p
1
(
x
) for its detection statistic
x
. For this
analysis a threshold of
x
min
= 5 was applied, which
is sufficiently low that the trigger density is dominated
by terrestrial triggers near threshold. There were
M
=
15 848 triggers observed above this threshold in the 17
days of observation time analyzed by
gstlal
.
For
pycbc
, the quantity
x
′
is the re-weighted SNR de-
tection statistic.
1
We set a threshold
x
′
min
= 8, above
1
When quoting pipeline-specific values we distinguish
pycbc
quantities with a prime.
arXiv:1606.03939v2 [astro-ph.HE] 20 Sep 2016
2
8
.
0
8
.
5
9
.
0
9
.
5
10
.
0
x
′
10
−
3
10
−
2
10
−
1
10
0
10
1
p
(
x
′
)
p
0
p
1
Figure 1
. Inferred terrestrial (
p
0
; blue) and astrophysical
(
p
1
; green) trigger densities for the
pycbc
pipeline as de-
scribed in Section 1.
which
M
′
= 270 triggers remain in the search. We use a
histogram of triggers collected from time-shifted data to
estimate the terrestrial trigger density,
p
0
(
x
′
), and a his-
togram of the recovered triggers from the injection sets
described in Section 2.2 of the Letter to estimate the as-
trophysical trigger density,
p
1
(
x
′
). These estimates are
shown in Figure 1. The uncertainty in the distribution of
triggers from this estimation procedure is much smaller
than the uncertainty in overall rate from the finite num-
ber statistics (see, for example, Figure 5). The empirical
estimate is necessary to properly account for the inter-
action of the various single- and double-interferometer
thresholds in the
pycbc
search (Abbott et al. 2016c).
At high SNR, where these thresholds are irrelevant, the
astrophysical triggers follow an approximately flat-space
volumetric density (see Section 3) of
p
1
(
x
′
)
'
3
x
′
3
min
x
′
4
,
(1)
but they deviate from this at smaller SNR due to thresh-
old effects in the search.
For the
pycbc
pipeline, a detection statistic
x
′
≥
10
.
1
corresponds to an estimated search false alarm rate
(FAR) of one per century.
2.
DERIVATION OF POISSON MIXTURE MODEL
LIKELIHOOD
In this section we derive the likelihood function in
Eq. (3) of the Letter. Consider first a search of the
type described in Section 1 over
N
T
intervals of time
of width
δ
i
,
{
i
= 1
,...,N
T
}
. Triggers above some fixed
threshold occur with an instantaneous rate in time and
detection statistic
x
given by the sum of the terrestrial
and astrophysical rates:
d
N
d
t
d
x
(
t,x
) =
R
0
(
t
)
p
0
(
x
;
t
) +
R
1
(
t
)
V
(
t
)
p
1
(
x
;
t
)
,
(2)
where
R
0
(
t
) is the instantaneous rate (number per unit
time) of terrestrial triggers,
R
1
(
t
) is the instantaneous
rate density (number per unit time per unit comoving
volume) of astrophysical triggers,
p
0
is the instantaneous
density in detection statistic of terrestrial triggers,
p
1
is
the instantaneous density in detection statistic of astro-
physical triggers, and
V
(
t
) is the instantaneous sensitive
comoving redshifted volume (Abbott et al. 2016a, see
also Eq. (15) of the Letter) of the detectors to the as-
sumed source population. The astrophysical rate
R
1
is
to any reasonable approximation constant over our ob-
servations so we will drop the time dependence of this
term from here on.
2
Note that
R
0
and
R
1
have different
units in this expression; the former is a rate (per time),
while the latter is a
rate density
(per time-volume). The
density
p
1
is independent of source parameters as de-
scribed in Section 3. Let
d
N
d
t
≡
∫
d
x
d
N
d
t
d
x
=
R
0
(
t
) +
R
1
V
(
t
)
.
(3)
If the search intervals
δ
i
are sufficiently short, they
will contain at most one trigger and the time-dependent
terms in Eq. (2) will be approximately constant. Then
the likelihood for a set of times and detection statistics
of triggers,
{
(
t
j
,x
j
)
|
j
= 1
,...,M
}
, is a product over in-
tervals containing a trigger (indexed by
j
) and intervals
that do not contain a trigger (indexed by
k
) of the cor-
responding Poisson likelihoods
L
=
M
∏
j
=1
d
N
d
t
d
x
(
t
j
,x
j
) exp
[
−
δ
j
d
N
d
t
(
t
j
)
]
×
{
N
T
−
M
∏
k
=1
exp
[
−
δ
k
d
N
d
t
(
t
k
)
]
}
(4)
(cf. Farr et al. (2015, Eq. (21)) or Loredo & Wasserman
(1995, Eq. (2.8))).
3
Now let the width of the observation
intervals
δ
i
go to zero uniformly as the number of inter-
vals goes to infinity. Then the products of exponentials
in Eq. (4) become an exponential of an integral, and we
have
L
=
M
∏
j
=1
[
d
N
d
t
d
x
(
t
j
,x
j
)
]
exp [
−
N
]
,
(5)
2
The astrophysical rate can, in principle, also depend on red-
shift, but in this paper we assume that the BBH coalescence rate
is constant in the comoving frame.
3
There is a typo in Eq. (2.8) of Loredo & Wasserman (1995).
The second term in the final bracket is missing a factor of
δt
.
3
where
N
=
∫
d
t
d
N
d
t
(6)
is the expected number of triggers of both types in the
total observation time
T
.
As discussed in Section 1, in our search we observe
that
R
0
remains approximately constant and that
p
0
re-
tains its shape over the observation time discussed here;
this assumption is used in our search background esti-
mation procedure (Abbott et al. 2016c). The astrophys-
ical distribution of triggers is universal (Section 3) and
also time-independent. Finally, the detector sensitivity
is observed to be stable over our 16 days of coincident
observations, so
V
(
t
)
'
const (Abbott et al. 2016b). We
therefore choose to simply ignore the time dimension in
our trigger set. This generates an estimate of the rate
that is sub-optimal (i.e. has larger uncertainty) but con-
sistent with using the full data set to the extent that the
detector sensitivity varies in time; since this variation is
small, the loss of information about the rate will be cor-
respondingly small. We
do
capture any variation in the
sensitivity with time in our Monte-Carlo procedure for
estimating
〈
V T
〉
that is described in Section 2.2 of the
Letter.
If we ignore the trigger time, then the appropriate
likelihood to use is a marginalization of Eq. (5) over the
t
j
. Let
̄
L≡
∫
∏
j
d
t
j
L
=
∏
j
[Λ
0
p
0
(
x
j
) + Λ
1
p
1
(
x
j
)] exp [
−
Λ
0
−
Λ
1
]
,
(7)
where
Λ
0
p
0
(
x
) =
∫
d
tR
0
(
t
)
p
0
(
x
;
t
)
,
(8)
and
Λ
1
p
1
(
x
) =
∫
d
tR
1
V
(
t
)
p
1
(
x
;
t
)
,
(9)
with
∫
d
xp
0
(
x
) =
∫
d
xp
1
(
x
) = 1
.
(10)
If we assume that
R
1
is constant in (comoving) time,
and measure
p
1
(
x
) by accumulating recovered injections
throughout the run as we have done, then this expres-
sion reduces to the likelihood in Eq. (3) of the Letter. A
similar argument with an additional population of trig-
gers produces Eq. (10) of the Letter.
2.1.
The Expected Number of Background Triggers
The procedure for estimating
p
0
(
x
) in the
pycbc
pipeline also provides an estimate of the mean number
of background events per experiment Λ
0
(Abbott et al.
200
220
240
260
280
300
320
340
Λ
0
0
2
4
6
8
10
12
Λ
1
Figure 2
. The two-dimensional posterior on terrestrial and
astrophysical trigger expected counts (Λ
0
and Λ
1
in Eq. (5)
of the Letter) for the
pycbc
search. Contours are drawn at
the 10%, 20%, . . . , 90%, and 99% credible levels. There is
no meaningful correlation between the two variables. The
Poisson uncertainty in the terrestrial count is
∼
√
270, or
16, which is also very nearly the Poisson uncertainty in the
total count. Because this uncertainty is much larger than
the astrophysical count, changes in the astrophysical count
do not force the terrestrial count to adjust in a meaningful
way and the variables are uncorrelated in the posterior.
2016c). The procedure for estimating
p
0
used in the
gstlal
pipeline, however, does not naturally provide an
estimate of Λ
0
; instead
gstlal
estimates Λ
0
by fitting
the observed number of triggers to a Poisson distribu-
tion. We have chosen to leave Λ
0
as a free parameter
in our canonical analysis with a broad prior and infer
it from the observed data, rather than using the
pycbc
background estimate to constrain the prior, which would
result in a much narrower posterior on Λ
0
. This is equiv-
alent to the
gstlal
procedure for Λ
0
estimation in the
absence of signals; the presence of a small number of
signals in our data here do not substantially change the
Λ
0
estimate due to the overwhelming number of back-
ground triggers in the data set.
Using a broad prior on Λ
0
is
conservative
in the sense
that it will broaden the posterior on Λ
1
from which we
infer rates. However, because there are so many more
triggers in both searches of terrestrial origin than astro-
physical there is little correlation between Λ
0
and Λ
1
,
and so there is little difference between the posterior we
obtain on Λ
1
and the posterior we would have obtained
had we implemented the tight prior on Λ
0
. Figure 2
shows the two-dimensional posterior we obtain from Eq.
(5) of the Letter on Λ
0
and Λ
1
.
We have checked that using a
δ
-function prior
p
(Λ
0
) =
δ
(Λ
0
−
270)
(11)
4
10
0
10
1
10
2
R
(
Gpc
−
3
yr
−
1
)
0
.
0
0
.
1
0
.
2
0
.
3
0
.
4
0
.
5
0
.
6
Rp
(
R
)
R
R
fixed
Figure 3
. The posterior on the population-based rate ob-
tained from our canonical analysis (blue) and an analysis
where the expected background count, Λ
0
, is fixed to the
value measured by the
pycbc
pipeline, Λ
0
= 270 (green).
There is no meaningful change in the rate posterior between
the two analyses.
in the
pycbc
analysis that is the result of the pipeline Λ
0
estimate from timeslides
4
(Abbott et al. 2016c) and us-
ing a looser prior that is the result of a
gstlal
estimate
on a single set of time-slid data produces no meaningful
change in our results. Figure 3 shows our canonical rate
posterior inferred with the
pycbc
Λ
0
prior in Eq. (11)
and our canonical broad prior.
3.
UNIVERSAL ASTROPHYSICAL TRIGGER
DISTRIBUTION
Both the
pycbc
and
gstlal
pipelines rely on the SNR
as part of their detection statistic,
x
. The SNR of an
astrophysical trigger is a function of the detector noise
at the time of detection and the parameters of the trig-
ger. Schutz (2011) and Chen & Holz (2014) demonstrate
that the distribution of the expected SNR
〈
ρ
〉
in a sim-
ple model of a detection pipeline that simply thresholds
on SNR,
ρ
≥
ρ
th
, with sources in the local universe is
universal
, that is, independent of the source properties.
It follows
p
(
〈
ρ
〉
) =
3
ρ
3
th
〈
ρ
〉
4
.
(12)
This result follows from the fact that the expected value
of the SNR in a matched-filter search for compact binary
coalescence (CBC) signals scales inversely with trans-
4
While the statistical uncertainty on the pipeline Λ
0
estimate is
not precisely zero,
σ
Λ
0
/
Λ
0
.
10
−
3
, it is so small that a
δ
-function
prior is appropriate.
verse comoving distance (Hogg 1999):
〈
ρ
〉
=
A
(
m
1
,m
2
,~a
1
,~a
2
,S
(
f
)
,z
)
B
(angles)
D
M
,
(13)
where
A
is an amplitude factor that depends on the
intrinsic properties (source-frame masses and spins) of
the source, the detector sensitivity expressed as a noise
power spectral density
S
(
f
) as a function of observer
frequency and redshift
z
, and
B
is an angular factor
depending on the location of the source in the sky and
the relative orientations of binary orbit and detector.
The redshift enters
A
only through shifting the source
waveform to lower frequency at higher redshift, changing
A
because the sensitivity varies with observer frequency
f
. For the redshifts to which we are sensitive to BBH
in this observation period this effect on
A
is small.
If we assume that the distribution of source parame-
ters is constant over the range of distances to which we
are sensitive, and ignore the small redshift-dependent
sensitivity correction mentioned above, then the distri-
bution of SNR will be governed entirely by the distri-
bution of distances of the sources, which, in the local
universe is approximately
p
(
D
M
)
∝
D
2
M
,
(14)
yielding the distribution of SNR given in Eq. (12).
Both the
pycbc
and
gstlal
pipelines use goodness-of-
fit statistics in addition to SNR and employ a more com-
plicated system of thresholds than this simple model,
but the empirical distribution of detection statistics re-
mains, to an approximation suitable for our purposes,
independent of the source parameters. Figure 4 shows
the distribution of recovered detection statistics for the
various injection campaigns with varying source distri-
bution used to estimate sensitive time-volumes in the
pycbc
pipeline. In each injection campaign
O
(1000)
signals were recovered. For loud signals, the detec-
tion statistic is proportional to SNR in this pipeline,
and the distribution is not sensitive to the complicated
thresholding in the pipeline, so we recover Eq. (12); for
quiet signals the interaction of various single-detector
thresolds in the pipeline causes the distribution to de-
viate from this analytic approximation, but it remains
independent of the distribution of sources. Note that
the
empirical
distribution of detection statistics, not the
analytic one, forms the basis for
p
1
, the foreground dis-
tribution used in this rate estimation work.
To quantify the deviations from universality, we have
preformed two-sample Kolmogorov-Smirnov (KS) tests
between all six pairings of the sets of detections statistics
recovered in the injection campaigns described in Sec-
tions 2 and 3 of the Letter and featured in Figure 4. The
most extreme KS
p
-value occurred with the comparison
between the injection set with BBH masses drawn flat in
5
10
1
10
2
x
10
−
3
10
−
2
10
−
1
10
0
p
(
x
)
Analytic
GW150914
LVT151012
Flat in Log
Power Law
All
Figure 4
. The distribution of detection statistics in the
pycbc
pipeline for the signals recovered in the injection cam-
paigns used to estimate sensitive time-volumes for various
BBH population assumptions (see Sections 2 and 3 of the
Letter). The solid line gives the analytic approximation to
the distribution from Eq. (12), which agrees well with the
recovered statistics for loud signals; for quieter signals the
interaction of various thresholds in the pipeline causes the
distribution to deviate from the analytic approximation, but
it remains independent of the source distribution.
log
m
and the one with masses drawn from a power law
(both described in Section 3 of the Letter); this test gave
a
p
-value of 0
.
013. Given that we have performed six
identical comparisons we cannot reject the null hypoth-
esis that the empirical distributions used for rate esti-
mation from the
pycbc
pipeline are identical even at the
relatively weak significance
α
= 0
.
05. Certainly any dif-
ferences in detection statistic distribution attributable
to the BBH population are far too small to matter with
the few astrophysical signals in our data set (compared
with
O
(1000) recovered injections in each campaign).
Because the distribution of detection statistics is, to
a very good approximation,
universal
, we cannot learn
anything about the source population from the detec-
tion statistic alone; we must instead resort to parameter
estimation (PE) followup (Veitch et al. 2015; Abbott
et al. 2016e) of triggers to determine their parameters.
The parameters of the waveform template that produced
the trigger can be used to guess the parameters of the
source that generated that trigger, but the bias and un-
certainty in this estimate are very large compared to
the PE estimate. We therefore ignore the parameters of
the waveform template that generated the trigger in the
assignment of triggers to BBH classes.
4.
COUNT POSTERIOR
We impose a prior on the Λ parameters of:
p
(Λ
1
,
Λ
0
)
∝
1
√
Λ
1
1
√
Λ
0
.
(15)
The posterior on expected counts is proportional to
the product of the likelihood from Eq. (3) of the Letter
and the prior from Eq. (15):
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
.
(16)
For estimation of the Poisson rate parameter in a simple
Poisson model, the Jeffreys prior is 1
/
√
Λ. With this
prior, the posterior mean on Λ is
N
+1
/
2 for
N
observed
counts. With a prior proportional to 1
/
Λ the mean is
N
for
N >
0, but the posterior is improper when
N
= 0.
For a flat prior, the mean is
N
+1. Though the behaviour
of the mean is not identical with our mixture model
posterior, it is similar; because we find
〈
Λ
1
〉
1
/
2, the
choice of prior among these three reasonable options has
little influence on our results here.
For the
pycbc
data set we find the posterior me-
dian and 90% credible range Λ
1
= 3
.
2
+4
.
9
−
2
.
4
above our
threshold. For the
gstlal
set we find the posterior me-
dian and 90% credible range Λ
1
= 4
.
8
+7
.
9
−
3
.
8
. Though
we have only one event (GW150914) at exceptionally
high significance, and one other at marginal significance
(LVT151012), the counting analysis shows these to be
consistent with the possible presence of several more
events of astrophysical origin at lower detection statistic
in both pipelines.
The thresholds applied to the
pycbc
and
gstlal
trig-
gers for this analysis are
not
equivalent to each other in
terms of either SNR or false alarm rate; instead, both
thresholds have been chosen so that the rate of trig-
gers of terrestrial origin (Λ
0
p
0
) dominates near thresh-
old. Since the threshold is set at
different
values for
each pipeline, we do not expect the counts to be the
same between pipelines.
The estimated astrophysical and terrestrial trigger
rate densities (Eq. (1) of the Letter) for
pycbc
are plot-
ted in Figure 5. We select triggers from a subset of
the search parameter space (i.e. our bank of template
waveforms) that contains GW150914 as well as the mass
range considered for possible alternative populations of
BBH binaries in Section 3 of the Letter. There are
M
′
= 270 two-detector coincident triggers in this range
in the
pycbc
search (Abbott et al. 2016c). Figure 5 also
shows an estimate of the density of triggers that com-
prise our data set which agrees well with our inference
of the trigger rate.