of 11
SUPPLEMENT:
THE RATE OF BINARY BLACK HOLE MERGERS INFERRED FROM ADVANCED LIGO
OBSERVATIONS SURROUNDING GW150914
(
2016, ApJL, 833, L1
)
LIGO Scienti
fi
c Collaboration and Virgo Collaboration
(
See the end matter for the full list of authors.
)
Received 2016 May 27; accepted 2016 September 22; published 2016 November 30
ABSTRACT
This article provides supplemental information for a Letter reporting the rate of
(
BBH
)
coalescences inferred from 16
days of coincident Advanced LIGO observations surrounding the transient
(
GW
)
signal GW150914. In that work we
reported various rate estimates whose 90% con
fi
dence intervals fell in the range 2
600
Gpc
3
yr
1
. Here we give
details on our method and computations, 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 for estimating our detector sensitivity, which is calibrated to our measurements.
Key words:
gravitational waves
stars: black holes
Supporting material:
data behind
fi
gures
The
fi
rst detection of a gravitational-wave
(
GW
)
signal from
a merging binary black hole
(
BBH
)
system was described in
Abbott et al.
(
2016d
)
. Abbott et al.
(
2016g
)
reported on
inference of the local BBH merger rate from surrounding
Advanced LIGO observations. This Supplement provides
supporting material and methodological 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
fi
ltering 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
(
S
/
N
)
is re-weighted by a chi-squared factor
(
Allen
2005
)
to account for template-data mismatch
(
Babak et al.
2013
)
; the
re-weighted single-detector S
/
Ns are combined in quadrature
to produce a detection 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 S
/
Ns and a signal-consistency
statistic. An analytic estimate of the distribution of astro-
physical signals in multiple-detector S
/
N and signal-consis-
tency 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 triggers of terrestrial
origin occur independently in the two detectors. The back-
ground estimate is built from observations 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 astrophysical
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
px
0
()
and
p
1
(
x
)
for its detection statistic
x
. For this analysis a
threshold of
=
x
5
min
was applied, which is suf
fi
ciently low
that the trigger density is dominated by terrestrial triggers near
a 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 S
/
N detection
statistic.
137
We set a threshold
¢
=
x
8
min
, above 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,
¢
px
0
()
, and a histogram of the recovered triggers from the
injection sets described in Section 2.2 of the
Letter
to estimate the
astrophysical trigger density,
¢
px
1
()
. These estimates are shown in
Figure
1
. The uncertainty in the distribution of triggers from this
estimation procedure is much sma
ller than the uncertainty in the
overall rate from the
fi
nite number statistics
(
see, for example,
Figure
5
)
. The empirical estimate is necessary to properly account
for the interaction of the various s
ingle- and double-interferometer
thresholds in the
pycbc
search
(
Abbott et al.
2016c
)
.Athigh
S
/
N, where these thresholds are irrelevant, the astrophysical
triggers follow an approximately
fl
at-space volumetric density
(
see Section
3
)
of
¢
¢
¢
px
x
x
3
,1
1
min
3
4
()
()
but they deviate from this at smaller S
/
N due to threshold
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 Equation
(
3
)
of the
Letter
. Consider a search of the type described in
The Astrophysical Journal Supplement Series,
227:14
(
11pp
)
, 2016 December
doi:10.3847
/
0067-0049
/
227
/
2
/
14
© 2016. The American Astronomical Society. All rights reserved.
137
When quoting pipeline-speci
fi
c values we distinguish
pycbc
quantities
with a prime.
1
Section
1
over
N
T
intervals of time, of width
δ
i
,
iN
1, ,
T
{
}
.
Triggers above some
fi
xed threshold occur with an instanta-
neous rate in time and detection statistic
x
given by the sum of
the terrestrial and astrophysical rates:
=+
dN
dtdx
t x
R tp x t
R tV tp x t
,; ;,2
0
0
1
1
(
)
() (
)
() () (
)
( )
where
R
t
0
()
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 the
detection statistic of terrestrial triggers,
p
1
is the instantaneous
density in the detection statistic of astrophysical triggers, and
V
(
t
)
is the instantaneous sensitive comoving redshifted volume
(
Abbott et al.
2016a
; see also Equation
(
15
)
of the
Letter
)
of the
detectors to the assumed source population. The astrophysical
rate
R
1
is to any reasonable approximation constant over our
observations so we will drop the time dependence of this term
from here on.
138
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 described in Section
3
. Let
ò
º=+
dN
dt
dx
dN
dtdx
Rt RVt
.3
01
()
()
( )
If the search intervals
δ
i
are suf
fi
ciently short, they will
contain at most one trigger and the time-dependent terms in
Equation
(
2
)
will be approximately constant. Then the
likelihood for a set of times and detection statistics of triggers,
tx j
M
,1,,
jj
{
()∣
}
, is a product over intervals containing a
trigger
(
indexed by
j
)
and intervals that do not contain a trigger
(
indexed by
k
)
of the corresponding Poisson likelihoods
d
d
=-
́-
=
=
-
dN
dtdx
tx
dN
dt
t
dN
dt
t
,exp
exp
4
j
M
jj
j
j
k
NM
kk
1
1
T
()
()
()
()
(
see Farr et al.
2015
, Equation
(
21
)
, or Loredo & Wasser-
man
1995
Equation
(
2.8
))
.
139
Now let the width of the
observation intervals
δ
i
go to zero uniformly as the number of
intervals goes to in
fi
nity. Then the products of exponentials in
Equation
(
4
)
become an exponential of an integral, and we have
=-
=
dN
dtdx
tx
N
,exp ,
5
j
M
jj
1
() []
()
where
ò
=
Ndt
dN
dt
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
retains its shape
over the observation time discussed here; this assumption is
used in our search background estimation procedure
(
Abbott
et al.
2016c
)
. The astrophysical 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 consistent 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 correspondingly small. We
do
capture
any variation in the sensitivity with time in our Monte Carlo
procedure for estimating
á
ñ
VT
, which 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 Equation
(
5
)
over the
t
j
. Let

ò
º
= L +L
-L -L
dt
px
px
exp
,
7
j
j
j
jj
0
0
1
1
01
̄
[]
[()
()][
] ()
where
ò
L=
px
dtRtpxt
;,
8
0
0
0
0
()
() ( )
()
and
ò
L=
px
dtRVtpxt
;,
9
1
1
1
1
()
() ( )
()
with
òò
==
dx p x
dx p x
1.
10
01
()
()
( )
Figure 1.
Inferred terrestrial
(
p
0
; blue
)
and astrophysical
(
p
1
; green
)
trigger
densities for the
pycbc
pipeline as described in Section
1
.
(
The data used to create this
fi
gure are available.
)
138
The astrophysical rate can, in principle, also depend on redshift, but in this
paper we assume that the BBH coalescence rate is constant in the comoving
frame.
139
There is a typo in Equation
(
2.8
)
of Loredo & Wasserman
(
1995
)
. The
second term in the
fi
nal bracket is missing a factor of
δ
t
.
2
The Astrophysical Journal Supplement Series,
227:14
(
11pp
)
, 2016 December
Abbott et al.
If we assume that
R
1
is constant in
(
comoving
)
time, and
measure
px
1
()
by accumulating recovered injections through-
out the run as we have done, then this expression reduces to the
likelihood in Equation
(
3
)
of the
Letter
. A similar argument
with an additional population of triggers produces Equation
(
10
)
of the
Letter
.
2.1. The Expected Number of Background Triggers
The procedure for estimating
px
0
()
in the
pycbc
pipeline
also provides an estimate of the mean number of background
events per experiment
Λ
0
(
Abbott et al.
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
fi
tting the observed number of triggers to a
Poisson distribution. 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 equivalent to the
gstlal
procedure for
Λ
0
estimation in the absence of signals;
the presence of a small number of signals in our data here does
not substantially change the
Λ
0
estimate due to the over-
whelming number of background 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 searches
of terrestrial origin than there are in those of astrophysical
origin there is little correlation between
Λ
0
and
Λ
1
, 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 Equation
(
5
)
of the
Letter
on
Λ
0
and
Λ
1
.
We have checked that using a
δ
-function prior
d
L=L-
p
270
11
00
() (
)
()
in the
pycbc
analysis that is the result of the pipeline
Λ
0
estimate from timeslides
140
(
Abbott et al.
2016c
)
and using a
looser prior that is the result of a
gstlal
estimate on a single
set of time-slid data produce no meaningful change in our
results. Figure
3
shows our canonical rate posterior inferred
with the
pycbc
Λ
0
prior in Equation
(
11
)
and our canonical
broad prior.
3. UNIVERSAL ASTROPHYSICAL
TRIGGER DISTRIBUTION
Both the
pycbc
and
gstlal
pipelines rely on the S
/
Nas
part of their detection statistic,
x
. The S
/
N of an astrophysical
trigger is a function of the detector noise at the time of
detection and the parameters of the trigger. Schutz
(
2011
)
and
Chen & Holz
(
2014
)
demonstrate that the distribution of the
expected S
/
N
r
á
ñ
in a simple model of a detection pipeline that
simply thresholds on S
/
N,
r
r
th
, with sources in the local
universe is
universal
, that is, independent of the source
properties. It follows
r
r
r
áñ=
áñ
p
3
.12
th
3
4
()
()
This result follows from the fact that the expected value of the
S
/
N in a matched-
fi
lter search for compact binary coalescence
(
CBC
)
signals scales inversely with transverse comoving
distance
(
Hogg
1999
)
:
r
áñ=
aa
Am m
S f zB
D
, , , ,
,
angles
,13
M
1212
(())()
()
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 the 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 parameters is
constant over the range of distances to which we are sensitive,
and ignore the small redshift-dependent sensitivity correction
mentioned above, then the distribution of S
/
N will be governed
entirely by the distribution of distances of the sources, which,
in the local universe, is approximately
μ
pD
D
,14
M
M
2
()
()
yielding the distribution of S
/
N given in Equation
(
12
)
.
Figure 2.
Two-dimensional posterior on terrestrial and astrophysical trigger
expected counts
(
Λ
0
and
Λ
1
in Equation
(
5
)
of the
Letter
)
for the
pycbc
search. Contours are drawn at the 10%, 20%,
K
, 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
near 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.
(
The data used to create this
fi
gure are available.
)
140
While the statistical uncertainty on the pipeline
Λ
0
estimate is not precisely
zero,
s
L
L
-
10
0
3
0
, it is so small that a
δ
-function prior is appropriate.
3
The Astrophysical Journal Supplement Series,
227:14
(
11pp
)
, 2016 December
Abbott et al.
Both the
pycbc
and
gstlal
pipelines use goodness-of-
fi
t
statistics in addition to S
/
N and employ a more complicated
system of thresholds than this simple model, but the empirical
distribution of detection statistics remains, to an approximation
suitable for our purposes, independent of the source para-
meters. Figure
4
shows the distribution of recovered detection
statistics for the various injection campaigns with varying
source distribution used to estimate sensitive time-volumes in
the
pycbc
pipeline. In each injection campaign
1000
(
)
signals were recovered. For loud signals, the detection statistic
is proportional to S
/
N in this pipeline, and the distribution is
not sensitive to the complicated thresholding in the pipeline, so
we recover Equation
(
12
)
; for quiet signals the interaction of
various single-detector thresholds in the pipeline causes the
distribution to deviate 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 distribution 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 Sections
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
fl
at in
m
log
and the one with masses drawn from a power
law
(
both described in Section
3
of the
Letter
)
;thistestgavea
p
-
value of 0.013. Given that we have performed six identical
comparisons we cannot reject the null hypothesis that the
empirical distributions used for rate estimation from the
pycbc
pipeline are identical even at the relatively weak signi
fi
cance
α
=
0.05. Certainly any differences in detection statistic distribu-
tion attributable to the BBH population are far too small to matter
with the few astrophysical signals in our data set
(
compared with
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 detection statistic alone; we
must instead resort to parameter estimation
(
PE
)
follow-up
(
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 uncertainty 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:
LLμ
LL
p
,
11
.15
10
10
()
()
The posterior on expected counts is proportional to the
product of the likelihood from Equation
(
3
)
of the
Letter
and
the prior from Equation
(
15
)
:
LL = ¼
μL+L
́-L-L
LL
=
pxjM
px
px
,1,,
exp
1
.16
j
j
M
jj
10
1
1
1
0
0
10
10
(∣{∣
})
[()
()]
[]
()
For estimation of the Poisson rate parameter in a simple
Poisson model, the Jeffreys prior is
L
1
. With this prior, the
Figure 3.
Posterior on the population-based rate obtained from our canonical
analysis
(
blue
)
and an analysis where the expected background count,
Λ
0
,is
fi
xed to the value measured by the
pycbc
pipeline,
L=
27
0
0
(
green
)
. There is
no meaningful change in the rate posterior between the two analyses.
(
The data used to create this
fi
gure are available.
)
Figure 4.
Distribution of detection statistics in the
pycbc
pipeline for the
signals recovered in the injection campaigns 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
Equation
(
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.
(
The data used to create this
fi
gure are available.
)
4
The Astrophysical Journal Supplement Series,
227:14
(
11pp
)
, 2016 December
Abbott et al.
posterior mean on
Λ
is
+
N
1
2
for
N
observed counts. With a
prior proportional to
L
1
the mean is
N
for
N
>
0, but the
posterior is improper when
N
=
0. For a
fl
at prior, the mean is
N
+
1. Though the behavior of the mean is not identical with
our mixture model posterior, it is similar; because we
fi
nd
á
1
2
1
, the choice of prior among these three reasonable
options has little in
fl
uence on our results here.
For the
pycbc
data set we
fi
nd the posterior median and
90% credible range
L=
-
+
3.2
1
2.4
4.9
above our threshold. For the
gstlal
set we
fi
nd the posterior median and 90% credible
range
L=
-
+
4.8
1
3.8
7.9
. Though we have only one event
(
GW150914
)
at exceptionally high signi
fi
cance, and one other
at marginal signi
fi
cance
(
LVT151012
)
, the counting analysis
shows these to be consistent with the possible presence of
several more events of astrophysical origin at a lower detection
statistic in both pipelines.
The thresholds applied to the
pycbc
and
gstlal
triggers
for this analysis are not equivalent to each other in terms of
either S
/
N or FAR; instead, both thresholds have been chosen
so that the rate of triggers of terrestrial origin
(
L
p
0
0
)
dominates
near the threshold. 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 te
rrestrial trigger rate densities
(
Equation
(
1
)
of the
Letter
)
for
pycbc
areplottedinFigure
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 tr
iggers in this range in the
pycbc
search
(
Abbott et al.
2016c
)
. Figure
5
also shows an
estimate of the density of triggers that comprise our data set,
which agrees well with our inference of the trigger rate.
Based on the probability of astr
ophysical origin inferred for
LVT151012 from the two-component mixture model in
Equation
(
16
)
and shown in Figure
6
, we introduce a third class
of signals and use a three-compon
ent mixture model with expected
counts
Λ
0
(
terrestrial
)
,
Λ
1
(
GW150914-like
)
,and
Λ
2
(
LVT151012-
like
)
to infer rates in Sections
2.1
and 2.2 of the
Letter
.
We use the Stan and
emcee
Markov Chain Monte Carlo
samplers
(
Foreman-Mackey et al.
2013
; Stan Development
Team
2015a
,
2015b
)
to draw samples from the posterior in
Equation
(
5
)
of the
Letter
for the two pipelines. We have
assessed the convergence and mixing of our chains using
empirical estimates of the autocorrelation length in each
parameter
(
Sokal
1996
)
, the Gelman-Rubin
R
convergence
statistic
(
Gelman & Rubin
1992
)
, and through visual inspection
of chain plots. By all measures, the chains appear to be well-
converged to the posterior distribution.
Table
1
contains the full results on expected counts and
associated sensitive time-volumes for both pipelines.
5. CALIBRATION UNCERTAINTY
The LIGO detectors are subject to uncertainty in their
calibration, in both the measured amplitude and the phase of
the GW strain. Abbott et al.
(
2016b
)
discussed the methods
used to calibrate the strain output of the detector during the 16
days of coincident observations discussed here. Abbott et al.
(
2016b
)
estimated that the reported strain is accurate to within
10% in amplitude and 10 degrees in phase between 20 Hz and
1 kHz throughout the observations.
The S
/
Ns reported by our searches are quadratically
sensitive to calibration errors because they are maximized
over arrival time, waveform phase, and a template bank of
waveforms
(
Allen
1996
; Brown & LIGO Scienti
fi
c Collabora-
tion
2004
)
. Abbott et al.
(
2016c
)
demonstrated that the other
search pipeline outputs are also not affected to a signi
fi
cant
degree by the calibration uncertainty present during our
observing run. Therefore, we ignore the effects of calibration
on the pipeline detection statistics
x
and
x
that we use here to
estimate rates from the
pycbc
and
gstlal
pipelines.
The amplitude calibration uncertainty in the detector results,
at leading order, in a corresponding uncertainty between the
luminosity distances of sources measured from real detector
outputs
(
Abbott et al.
2016e
)
and the luminosity distances used
to produce injected waveforms used to estimate sensitive time-
volumes in this work. A 10% uncertainty in
d
L
at these
redshifts corresponds to an approximately 30% uncertainty in
volume. We model this uncertainty by treating
á
ñ
VT
as a
Figure 5.
Inferred number density of astrophysical
(
green
)
, terrestrial
(
blue
)
, and all
(
red
)
triggers as a function of
¢
x
for the
pycbc
search
(
see Equation
(
1
)
of the
Letter
)
, using the models for each population described in Section
2.1
of the
Letter
. The solid lines give the posterior median and the shaded regions give the
symmetric 90% credible interval from the posterior in Equation
(
5
)
of the
Letter
. We also show a binned estimate of the trigger number density from the search
(
black
)
;
bars indicate the 68% con
fi
dence Poisson uncertainty on the number of triggers in the vertical-direction and bin width in the horizontal-direction.
(
The data used to create this
fi
gure are available.
)
5
The Astrophysical Journal Supplement Series,
227:14
(
11pp
)
, 2016 December
Abbott et al.