of 19
Properties of the Binary Black Hole Merger GW150914
B. P. Abbott
etal.
*
(LIGO Scientific Collaboration and Virgo Collaboration)
(Received 18 February 2016; revised manuscript received 18 April 2016; published 14 June 2016)
On September 14, 2015, the Laser Interferometer Gravitational-Wave Observatory (LIGO) detected a
gravitational-wave transient (GW150914); we characterize the properties of the source and its parameters.
The data around the time of the event were analyzed coherently across the LIGO network using a suite of
accurate waveform models that describe gravitational waves from a compact binary system in general
relativity. GW150914 was produced by a nearly equal mass binary black hole of masses
36
þ
5
4
M
and
29
þ
4
4
M
; for each parameter we report the median value and the range of the 90% credible interval. The
dimensionless spin magnitude of the more massive black hole is bound to be
<
0
.
7
(at 90% probability).
The luminosity distance to the source is
410
þ
160
180
Mpc, corresponding to a redshift
0
.
09
þ
0
.
03
0
.
04
assuming
standard cosmology. The source location is constrained to an annulus section of
610
deg
2
, primarily
in the southern hemisphere. The binary merges into a black hole of mass
62
þ
4
4
M
and spin
0
.
67
þ
0
.
05
0
.
07
.
This black hole is significantly more massive than any other inferred from electromagnetic observations in
the stellar-mass regime.
DOI:
10.1103/PhysRevLett.116.241102
I. INTRODUCTION
In Ref.
[1]
we reported the detection of gravitational
waves (GWs), observed on September 14, 2015 at 09:50:45
UTC, by the twin instruments of the twin instruments of
the Laser Interferometer Gravitational-Wave Observatory
(LIGO) located at Hanford, Washington, and Livingston,
Louisiana, in the USA
[2,3]
. The transient signal, named
GW150914, was detected with a false-alarm-probability of
<
2
×
10
7
and has been associated with the merger of a
binary system of black holes (BHs).
Here we discuss the properties of this source and its
inferred parameters. The results are based on a complete
analysis of the data surrounding this event. The only
information from the search stage is the time of arrival
of the signal. Crucially, this analysis differs from the search
in the following fundamental ways: it is coherent across the
LIGO network, it uses waveform models that include
the full richness of the physics introduced by BH spins,
and it covers the full multidimensional parameter space
of the considered models with a fine (stochastic) sampling;
we also account for uncertainty in the calibration of the
measured strain. The results of this analysis provide the
parameter values quoted in Ref.
[1]
describing the proper-
ties of the source of GW150914. (Following the publication
of Ref.
[1]
we repeated the analysis after fixing minor errors
in the coordinate transformation between precessing and
nonprecessing binary systems that affects waveforms used
in the analysis. The parameter values reported here, in
Table
I
and the figures, are the updated values and should
be used for future studies. The only value different from
what was reported in Ref.
[1]
is the total energy radiated,
which previously was
3
.
0
þ
0
.
5
0
.
5
M
c
2
and now rounds to
3
.
0
þ
0
.
5
0
.
4
M
c
2
).
In general relativity, two objects in orbit slowly spiral
together due to the loss of energy and angular momentum
through gravitational radiation
[7,8]
. This is in contrast
to Newtonian gravity where bodies can follow closed,
elliptical orbits
[9,10]
. As the binary shrinks, the frequency
and amplitude of the emitted GWs increase. Eventually the
two objects merge. If these objects are BHs, they form a
single perturbed BH that radiates GWs as a superposition of
quasinormal ringdown modes. Typically, one mode domi-
nates soon after merger, and an exponentially damped
oscillation at constant frequency can be observed as the BH
settles to its final state
[11,12]
.
An isolated BH is described by only its mass and spin,
since we expect the electric charge of astrophysical BHs
to be negligible
[13
16]
. Merging (BBHs) are therefore
relatively simple systems. The two BHs are described by
eight intrinsic parameters: the masses
m
1
;
2
(defined as the
gravitational masses of the BHs in isolation) and spins
S
1
;
2
(magnitude and orientation) of the individual BHs. For a
BH of mass
m
, the spin can be at most
Gm
2
=c
; hence,
it is conventional to quote the dimensionless spin magni-
tude
a
¼
c
j
S
j
=
ð
Gm
2
Þ
1
. Nine additional parameters are
needed to fully describe the binary: the location (luminosity
distance
D
L
, right ascension
α
, and declination
δ
); orien-
tation (the binary
s orbital inclination
ι
and polarization
ψ
);
*
Full author list given at the end of the article.
Published by the American Physical Society under the terms of
the
Creative Commons Attribution 3.0 License
. Further distri-
bution of this work must maintain attribution to the author(s) and
the published article
s title, journal citation, and DOI.
PRL
116,
241102 (2016)
PHYSICAL REVIEW LETTERS
week ending
17 JUNE 2016
0031-9007
=
16
=
116(24)
=
241102(19)
241102-1
Published by the American Physical Society
time
t
c
and phase
φ
c
of coalescence, and two parameters
describing eccentricity (the magnitude
e
and the argument
of periapsis) of the system.
Radiation reaction is efficient in circularizing orbits
[17]
before the signal enters the sensitive frequency band of the
instruments (
20
1000
Hz). In our analysis, we assume
circular orbits (we therefore do not include the eccentricity
parameters), and we find no evidence for residual eccen-
tricity, see the Discussion. Under the approximation of a
circular orbit, the binary emits GWs primarily at twice the
orbital frequency
[18]
.
The gravitational waveform observed for GW150914
comprises
10
cycles during the inspiral phase from 30 Hz,
followed by the merger and ringdown. The properties of the
binary affect the phase and amplitude evolution of the
observed GWs, allowing us to measure the source param-
eters. Here we briefly summarize these signatures, and
provide an insight into our ability to characterize the
properties of GW150914 before we present the details of
the Results; for methodological studies, we refer the reader
to Refs.
[4,19
24]
and references therein.
In general relativity, gravitational radiation is fully
described by two independent, and time-dependent polar-
izations,
h
þ
and
h
×
. Each instrument
k
measures the strain
h
k
¼
F
ðþÞ
k
h
þ
þ
F
ð
×
Þ
k
h
×
;
ð
1
Þ
a linear combination of the polarizations weighted by the
antenna beam patterns
F
ðþ
;
×
Þ
k
ð
α
;
δ
;
ψ
Þ
, which depend on the
source location in the sky and the polarization of the waves
[25,26]
. During the inspiral and at the leading order, the
GW polarizations can be expressed as
h
þ
ð
t
Þ¼
A
GW
ð
t
Þð
1
þ
cos
2
ι
Þ
cos
φ
GW
ð
t
Þ
;
ð
2a
Þ
h
×
ð
t
Þ¼
2
A
GW
ð
t
Þ
cos
ι
sin
φ
GW
ð
t
Þ
;
ð
2b
Þ
where
A
GW
ð
t
Þ
and
φ
GW
ð
t
Þ
are the GWamplitude and phase,
respectively. For a binary viewed face-on (cos
ι
¼
1
),
TABLE I. Summary of the parameters that characterize GW150914. For model parameters we report the median value as well as the
range of the symmetric 90% credible interval
[4]
; where useful, we also quote 90% credible bounds. For the logarithm of the Bayes
factor for a signal compared to Gaussian noise we report the mean and its 90% standard error from 4 parallel runs with a nested sampling
algorithm
[5]
. The source redshift and source-frame masses assume standard cosmology
[6]
. The spin-aligned EOBNR and precessing
IMRPhenom waveform models are described in the text. Results for the effective precession spin parameter
χ
p
used in the IMRPhenom
model are not shown as we effectively recover the prior; we constrain
χ
p
<
0
.
71
at 90% probability, see left panel of Fig.
5
. The Overall
results are computed by averaging the posteriors for the two models. For the Overall results we quote both the 90% credible interval or
bound and an estimate for the 90% range of systematic error on this determined from the variance between waveform models. The sky
location associated with GW150914 is presented in Fig.
4
and discussed in the text.
EOBNR
IMRPhenom
Overall
Detector-frame total mass
M=M
70
.
3
þ
5
.
3
4
.
8
70
.
9
þ
4
.
0
3
.
9
70
.
6
þ
4
.
6

0
.
5
4
.
5

1
.
3
Detector-frame chirp mass
M
=M
30
.
2
þ
2
.
5
1
.
9
30
.
6
þ
1
.
8
1
.
8
30
.
4
þ
2
.
1

0
.
2
1
.
9

0
.
5
Detector-frame primary mass
m
1
=M
39
.
4
þ
5
.
5
4
.
9
38
.
5
þ
5
.
6
3
.
6
38
.
9
þ
5
.
6

0
.
6
4
.
3

0
.
4
Detector-frame secondary mass
m
2
=M
30
.
9
þ
4
.
8
4
.
4
32
.
2
þ
3
.
6
4
.
8
31
.
6
þ
4
.
2

0
.
1
4
.
7

0
.
9
Detector-frame final mass
M
f
=M
67
.
1
þ
4
.
6
4
.
4
67
.
6
þ
3
.
6
3
.
5
67
.
4
þ
4
.
1

0
.
4
4
.
0

1
.
2
Source-frame total mass
M
source
=M
65
.
0
þ
5
.
0
4
.
4
65
.
0
þ
4
.
0
3
.
6
65
.
0
þ
4
.
5

0
.
8
4
.
0

0
.
7
Source-frame chirp mass
M
source
=M
27
.
9
þ
2
.
3
1
.
8
28
.
1
þ
1
.
7
1
.
6
28
.
0
þ
2
.
0

0
.
3
1
.
7

0
.
3
Source-frame primary mass
m
source
1
=M
36
.
3
þ
5
.
3
4
.
5
35
.
3
þ
5
.
2
3
.
4
35
.
8
þ
5
.
3

0
.
9
3
.
9

0
.
1
Source-frame secondary mass
m
source
2
=M
28
.
6
þ
4
.
4
4
.
2
29
.
6
þ
3
.
3
4
.
3
29
.
1
þ
3
.
8

0
.
1
4
.
3

0
.
7
Source-frame final mass
M
source
f
=M
62
.
0
þ
4
.
4
4
.
0
62
.
0
þ
3
.
7
3
.
3
62
.
0
þ
4
.
1

0
.
7
3
.
7

0
.
6
Mass ratio
q
0
.
79
þ
0
.
18
0
.
19
0
.
84
þ
0
.
14
0
.
20
0
.
82
þ
0
.
17

0
.
01
0
.
20

0
.
03
Effective inspiral spin parameter
χ
eff
0
.
09
þ
0
.
19
0
.
17
0
.
05
þ
0
.
13
0
.
15
0
.
07
þ
0
.
16

0
.
01
0
.
17

0
.
05
Dimensionless primary spin magnitude
a
1
0
.
32
þ
0
.
45
0
.
28
0
.
32
þ
0
.
53
0
.
29
0
.
32
þ
0
.
49

0
.
06
0
.
29

0
.
01
Dimensionless secondary spin magnitude
a
2
0
.
57
þ
0
.
40
0
.
51
0
.
34
þ
0
.
54
0
.
31
0
.
44
þ
0
.
50

0
.
08
0
.
40

0
.
02
Final spin
a
f
0
.
67
þ
0
.
06
0
.
08
0
.
66
þ
0
.
04
0
.
06
0
.
67
þ
0
.
05

0
.
01
0
.
07

0
.
02
Luminosity distance
D
L
=
Mpc
390
þ
170
180
440
þ
150
180
410
þ
160

20
180

40
Source redshift
z
0
.
083
þ
0
.
033
0
.
036
0
.
093
þ
0
.
029
0
.
036
0
.
088
þ
0
.
032

0
.
005
0
.
037

0
.
008
Upper bound on primary spin magnitude
a
1
0.65
0.74
0
.
69

0
.
08
Upper bound on secondary spin magnitude
a
2
0.93
0.78
0
.
89

0
.
13
Lower bound on mass ratio
q
0.64
0.68
0
.
66

0
.
03
Log Bayes factor ln
B
s=n
288
.
7

0
.
2
290
.
3

0
.
1

PRL
116,
241102 (2016)
PHYSICAL REVIEW LETTERS
week ending
17 JUNE 2016
241102-2
GWs are circularly polarized, whereas for a binary observed
edge-on (cos
ι
¼
0
), GWs are linearly polarized.
During the inspiral, the phase evolution
φ
GW
ð
t
;
m
1
;
2
;
S
1
;
2
Þ
can be computed using post-Newtonian (PN) theory, which
is a perturbative expansion in powers of the orbital velocity
v=c
[27]
. For GW150914,
v=c
is in the range
0
.
2
0
.
5
in
the LIGO sensitivity band. At the leading order, the phase
evolution is driven by a particular combination of the two
masses, commonly called the chirp mass
[28]
,
M
¼
ð
m
1
m
2
Þ
3
=
5
M
1
=
5
c
3
G

5
96
π
8
=
3
f
11
=
3
_
f

3
=
5
;
ð
3
Þ
where
f
is the GW frequency,
_
f
is its time derivative, and
M
¼
m
1
þ
m
2
is the total mass. Additional parameters
enter at each of the following PN orders. First, the mass
ratio,
q
¼
m
2
=m
1
1
, and the BH spin components
parallel to the orbital angular momentum vector
L
affect
the phase evolution. The full degrees of freedom of the
spins enter at higher orders. Thus, from the inspiral, we
expect to measure the chirp mass with highest accuracy and
only place weak constraints on the mass ratio and (the
components parallel to
L
of) the spins of the BHs
[21,29]
.
Spins are responsible for an additional characteristic
effect: if misaligned with respect to
L
, they cause the
binary
s orbital plane to precess around the almost-constant
direction of the total angular momentum of the binary,
J
¼
L
þ
S
1
þ
S
2
. This leaves characteristic amplitude and
phase modulations in the observed strain
[30,31]
,as
ψ
and
ι
become time dependent. The size of these modulations
depends crucially on the viewing angle of the source.
As the BHs get closer to each other and their velocities
increase, the accuracy of the PN expansion degrades, and
eventually the full solution of Einstein
s equations is
needed to accurately describe the binary evolution. This
is accomplished using numerical relativity (NR) which,
after the initial breakthrough
[32
34]
, has been improved
continuously to achieve the sophistication of modeling
needed for our purposes. The details of the ringdown are
primarily governed by the mass and spin of the final BH. In
particular, the final mass and spin determine the (constant)
frequency and decay time of the BH
s ringdown to its
final state
[35]
. The late stage of the coalescence allows
us to measure the total mass which, combined with the
measurement of the chirp mass and mass ratio from the
early inspiral, yields estimates of the individual component
masses for the binary.
The observed frequency of the signal is redshifted by a
factor of
ð
1
þ
z
Þ
, where
z
is the cosmological redshift.
There is no intrinsic mass or length scale in vacuum general
relativity, and the dimensionless quantity that incorporates
frequency is
fGm=c
3
. Consequently, a redshifting of
frequency is indistinguishable from a rescaling of the
masses by the same factor
[20,36,37]
. We therefore
measure redshifted masses
m
, which are related to source
frame masses by
m
¼ð
1
þ
z
Þ
m
source
. However, the GW
amplitude
A
GW
, Eq.
(2)
, also scales linearly with the mass
and is inversely proportional to the comoving distance in
an expanding universe. Therefore, the amplitude scales
inversely with the luminosity distance,
A
GW
1
=D
L
, and
from the GW signal alone we can directly measure the
luminosity distance, but not the redshift.
The observed time delay, and the need for the registered
signal at the two sites to be consistent in amplitude and
phase, allow us to localize the source to a ring on the sky
[38,39]
. Where there is no precession, changing the
viewing angle of the system simply changes the observed
waveform by an overall amplitude and phase. Furthermore,
the two polarizations are the same up to overall amplitude
and phase. Thus, for systems with minimal precession, the
distance, binary orientation, phase at coalescence, and sky
location of the source change the overall amplitude and
phase of the source in each detector, but they do not change
the signal morphology. Phase and amplitude consistency
allow us to untangle some of the geometry of the source. If
the binary is precessing, the GW amplitude and phase have
a complicated dependency on the orientation of the binary,
which provides additional information.
Our ability to characterize GW150914 as the signature of
a binary system of compact objects, as we have outlined
above, is dependent on the finite signal-to-noise ratio
(SNR) of the signal and the specific properties of the
underlying source. These properties described in detail
below, and the inferred parameters for GW150914 are
summarized in Table
I
and Figs.
1
6
.
II. METHOD
Full information about the properties of the source is
provided by the probability density function (PDF)
p
ð
~
θ
j
~
d
Þ
of the unknown parameters
~
θ
, given the two data streams
from the instruments
~
d
.
The posterior PDF is computed through a straightfor-
ward application of Bayes
theorem
[40,41]
. It is propor-
tional to the product of the likelihood of the data given the
parameters
L
ð
~
d
j
~
θ
Þ
, and the prior PDF on the parameters
p
ð
~
θ
Þ
before we consider the data. From the (marginalized)
posterior PDF, shown in Figs.
1
4
for selected parameters,
we then construct credible intervals for the parameters,
reported in Table
I
.
In addition, we can compute the evidence
Z
for the
model under consideration. The evidence (also known as
marginal likelihood) is the average of the likelihood under
the prior on the unknown parameters for a specific model
choice.
The computation of marginalized PDFs and the
model evidence require the evaluation of multidimensional
integrals. This is addressed by using a suite of Bayesian
parameter-estimation and model-selection algorithms tail-
ored to this problem
[42]
. We verify the results by using
PRL
116,
241102 (2016)
PHYSICAL REVIEW LETTERS
week ending
17 JUNE 2016
241102-3
two
independent
stochastic sampling engines based on
Markov-chain Monte Carlo
[43,44]
and nested sampling
[5,45]
techniques. (The marginalized PDFs and model
evidence are computed using the
LALInference
package of the LIGO Algorithm Library (LAL) software
suite
[46]
).
At the detector output we record the data
d
k
ð
t
Þ¼
n
k
ð
t
Þþ
h
M
k
ð
t
;
~
θ
Þ
, where
n
k
is the noise, and
h
M
k
is the
measured strain, which differs from the physical strain
h
k
from Eq.
(1)
as a result of the detectors
calibration
[47]
.In
the frequency domain, we model the effect of calibration
uncertainty by considering
~
h
M
k
ð
f
;
~
θ
Þ¼
~
h
k
ð
f
;
~
θ
Þ½
1
þ
δ
A
k
ð
f
;
~
θ
Þ
× exp
½
i
δφ
k
ð
f
;
~
θ
Þ
;
ð
4
Þ
where
~
h
M
k
and
~
h
k
are the Fourier representation of the time-
domain functions
h
M
k
and
h
k
, respectively.
δ
A
k
ð
f
;
~
θ
Þ
and
δφ
k
ð
f
;
~
θ
Þ
are the frequency-dependent amplitude and
phase calibration-error functions, respectively. These
calibration-error functions are modeled using a cubic spline
polynomial, with five nodes per spline model placed
uniformly in ln
f
[48]
.
We have analyzed the data at the time of this event using
a
coherent
analysis. Under the assumption of stationary,
Gaussian noise uncorrelated in each detector
[49]
, the
likelihood function for the LIGO network is
[20,42]
L
ð
~
d
j
~
θ
Þ
exp

1
2
X
k
¼
1
;
2
h
h
M
k
ð
~
θ
Þ
d
k
j
h
M
k
ð
~
θ
Þ
d
k
i

;
ð
5
Þ
where
h
·
j
·
i
is the noise-weighted inner product
[20]
.We
model the noise as a stationary Gaussian process of zero
mean and known variance, which is estimated from the
power spectrum computed using up to 1024 s of data
adjacent to, but not containing, the GW signal
[42]
.
The source properties are encoded into the two polar-
izations
h
þ
and
h
×
that enter the analysis through Eqs.
(1)
and
(4)
. Here we focus on the case in which they originate
from a compact binary coalescence; we use model
waveforms (described below) that are based on solving
Einstein
s equations for the inspiral and merger of two BHs.
A. BBH waveform models
For the modeled analysis of binary coalescences, an
accurate waveform prediction for the gravitational radiation
h
þ
;
×
is essential. As a consequence of the complexity of
solving the two body problem in general relativity, several
techniques have to be combined to describe all stages
of the binary coalescence. While the early inspiral is well
described by the analytical PN expansion
[27]
, which relies
on small velocities and weak gravitational fields, the
strong-field merger stage can only be solved in full
generality by large-scale NR simulations
[32
34]
. Since
these pioneering works, numerous improvements have
enabled numerical simulations of BBHs with sufficient
accuracy for the applications considered here and for the
region of parameter space of relevance to GW150914
(see, e.g., Refs.
[50
52]
). Tremendous progress has also
been made in the past decade to combine analytical and
numerical approaches, and now several accurate waveform
models are available, and they are able to describe the entire
coalescence for a large variety of possible configurations
[51,53
59]
. Extending and improving such models is an
active area of research, and none of the current models can
capture all possible physical effects (eccentricity, higher
order gravitational modes in the presence of spins, etc.) for
all conceivable binary systems. We discuss the current state
of the art below.
In the Introduction, we outlined how the binary param-
eters affect the observable GW signal, and now we discuss
the BH spins in greater detail. There are two main effects
that the BH spins
S
1
and
S
2
have on the phase and
amplitude evolution of the GW signal. The spin projections
along the direction of the orbital angular momentum affect
the inspiral rate of the binary. In particular, spin compo-
nents aligned (antialigned) with
L
increase (decrease) the
number of orbits from any given separation to merger with
respect to the nonspinning case
[27,60]
. Given the limited
SNR of the observed signal, it is difficult to untangle the
full degrees of freedom of the individual BHs
spins, see,
e.g., Refs.
[61,62]
. However, some spin information is
encoded in a dominant spin effect. Several possible one-
dimensional parametrizations of this effect can be found in
the literature
[21,63,64]
; here, we use a simple mass-
weighted linear combination of the spins
[63,65
67]
χ
eff
¼
c
GM

S
1
m
1
þ
S
2
m
2

·
L
j
L
j
;
ð
6
Þ
which takes values between
1
(both BHs have maximal
spins antialigned with respect to the orbital angular
momentum) and
þ
1
(maximal aligned spins).
Having described the effect of the two spin components
aligned with the orbital angular momentum, four in-plane
spin components remain. These lead to precession of the
spins and the orbital plane, which in turn introduces
modulations in the strain amplitude and phase as measured
at the detectors. At leading order in the PN expansion, the
equations that describe precession in a BBH are
[30]
_
L
¼
G
c
2
r
3
ð
B
1
S
1
þ
B
2
S
2
Þ
×
L
;
ð
7
Þ
_
S
i
¼
G
c
2
r
3
B
i
L
×
S
i
;
ð
8
Þ
where
S
i
is the component of the spin perpendicular to
L
;
overdots denote time derivatives;
r
is the orbital separation;
PRL
116,
241102 (2016)
PHYSICAL REVIEW LETTERS
week ending
17 JUNE 2016
241102-4
B
1
¼
2
þ
3
q=
2
and
B
2
¼
2
þ
3
=
ð
2
q
Þ
, and
i
¼f
1
;
2
g
.It
follows from Eqs.
(7)
and
(8)
that
L
and
S
i
precess around
the almost-constant direction of the total angular momen-
tum
J
. For a nearly equal-mass binary (as we find is the
case for GW150914), the precession angular frequency
can be approximated by
Ω
p
7
GJ=
ð
c
2
r
3
Þ
, and the total
angular momentum is dominated during the inspiral by
the orbital contribution,
J
L
. Additional, higher
order spin-spin interactions can also contribute signifi-
cantly to precession effects for some comparable-mass
binaries
[68,69]
.
The in-plane spin components rotate within the orbital
plane at different velocities. Because of nutation of the
orbital plane, the magnitude of the in-plane spin compo-
nents oscillates around a mean value, but those oscillations
are typically small. To first approximation, one can
quantify the level of precession in a binary by averaging
over the relative in-plane spin orientation. This is achieved
by the following effective precession spin parameter
[70]
:
χ
p
¼
c
B
1
Gm
2
1
max
ð
B
1
S
1
;B
2
S
2
Þ
>
0
;
ð
9
Þ
where
χ
p
¼
0
corresponds to an aligned-spin (nonprecess-
ing) system, and
χ
p
¼
1
to a binary with the maximum
level of precession. Although the definition of
χ
p
is
motivated by configurations that undergo many precession
cycles during inspiral, we find that it is also a good
approximate indicator of the in-plane spin contribution
for the late inspiral and merger of GW150914.
For the analysis of GW150914, we consider waveform
models constructed within two frameworks capable of
accurately describing the GW signal in the parameter space
of interest. The effective-one-body (EOB) formalism
[71
75]
combines perturbative results from the weak-field
PN approximation with strong-field effects from the test-
particle limit. The PN results are resummed to provide the
EOB Hamiltonian, radiation-reaction force, and GW polar-
izations. The Hamiltonian is then improved by calibrating
(unknown) higher order PN terms to NR simulations.
Henceforth, we use
EOBNR
to indicate waveforms
within this formalism. Waveform models that include spin
effects have been developed, for both double nonprecessing
spins
[54,59]
(11 independent parameters) and double
precessing spins
[55]
(15 independent parameters). Here,
we report results using the nonprecessing model
[54]
tuned to NR simulations
[76]
. This model is formulated
as a set of differential equations that are computationally
too expensive to solve for the millions of likelihood
evaluations required for the analysis. Therefore, a fre-
quency-domain reduced-order model
[77,78]
was imple-
mented that faithfully represents the original model with an
accuracy that is better than the statistical uncertainty caused
by the instruments
noise. (In LAL, as well as in technical
publications, the aligned, precessing and reduced-order
EOBNR models are called
SEOBNRv2
,
SEOBNRv3
and
SEOBNRv2_ROM_DoubleSpin
, respectively). Bayesian
analyses that use the double precessing spin model
[55]
are
more time consuming and are not yet finalized. The results
will be fully reported in a future publication.
An alternative inspiral-merger-ringdown phenomeno-
logical formalism
[79
81]
is based on extending frequency-
domain PN expressions and hybridizing PN and EOB with
NR waveforms. Henceforth, we use
IMRPhenom
to
indicate waveforms within this formalism. Several wave-
form models that include aligned-spin effects have been
constructed within this approach
[58,66,67]
, and here we
employ the most recent model based on a fitting of untuned
EOB waveforms
[54]
and NR hybrids
[57,58]
of non-
precessing systems. To include leading-order precession
effects, aligned-spin waveforms are rotated into precessing
waveforms
[82]
. Although this model captures some two-
spin effects, it was principally designed to accurately model
the waveforms with respect to an effective spin parameter
similar to
χ
eff
above. The dominant precession effects are
introduced through a lower-dimensional effective spin
description
[70,83]
, motivated by the same physical argu-
ments as the definition of
χ
p
. This provides us with an
effective-precessing-spin model
[83]
with 13 independent
parameters. (In LAL, as well as in some technical pub-
lications, the model used is called
IMRPhenomPv2
).
All models we use are restricted to circular inspirals.
They also include only the dominant spherical harmonic
modes in the nonprecessing limit.
B. Choice of priors
We analyze coherently 8 s of data with a uniform prior
on
t
c
of width of

0
.
1
s, centered on the time reported by
the online analysis
[1,84]
, and a uniform prior in
½
0
;
2
π

for
φ
c
. We consider the frequency region between 20 Hz,
below which the sensitivity of the instruments significantly
degrades (see panel (b) of Fig. 3 in Ref.
[1]
), and 1024 Hz,
a safe value for the highest frequency contribution to
radiation from binaries in the mass range considered here.
Given the lack of any additional astrophysical constraints
on the source at hand, our prior choices on the parameters
are uninformative. We assume sources uniformly distrib-
uted in volume and isotropically oriented. We use uniform
priors in
m
1
;
2
½
10
;
80

M
, with the constraint that
m
2
m
1
. We use a uniform prior in the spin magnitudes
a
1
;
2
½
0
;
1

. For angles subject to change due to precession
effects we give values at a reference GW frequency
f
ref
¼
20
Hz. We use isotropic priors on the spin orienta-
tion for the precessing model. For the nonprecessing
model, the prior on the spin magnitudes may be interpreted
as the dimensionless spin projection onto
L
having a
uniform distribution
½
1
;
1

. This range includes binaries
where the two spins are strongly antialigned relative to
one another. Many such antialigned-spin comparable-mass
systems are unstable to large-angle precession well before
PRL
116,
241102 (2016)
PHYSICAL REVIEW LETTERS
week ending
17 JUNE 2016
241102-5
entering our sensitive band
[85,86]
and could not have
formed from an asymptotically spin antialigned binary.
We could exclude those systems if we believe the binary is
not precessing. However, we do not make this assumption
here and instead accept that the models can only extract
limited spin information about a more general, precessing
binary.
We also need to specify the prior ranges for the ampli-
tude and phase error functions
δ
A
k
ð
f
;
~
θ
Þ
and
δφ
k
ð
f
;
~
θ
Þ
, see
Eq.
(5)
. The calibration during the time of observation of
GW150914 is characterized by a
1
-
σ
statistical uncertainty
of no more than 10% in amplitude and 10° in phase
[1,47]
.
We use zero-mean Gaussian priors on the values of the
spline at each node with widths corresponding to the
uncertainties quoted above
[48]
. Calibration uncertainties
therefore add 10 parameters per instrument to the model
used in the analysis. For validation purposes we also
considered an independent method that assumes frequency-
independent calibration errors
[87]
, and obtained consistent
results.
III. RESULTS
The results of the analysis using binary coalescence
waveforms are posterior PDFs for the parameters describ-
ing the GW signal and the model evidence. A summary is
provided in Table
I
. For the model evidence, we quote
(the logarithm of) the Bayes factor
B
s=n
¼
Z
=
Z
n
, which
is the evidence for a coherent signal hypothesis divided
by that for (Gaussian) noise
[5]
. At the leading order, the
Bayes factor and the optimal SNR
ρ
¼½
P
k
h
h
M
k
j
h
M
k
i
1
=
2
are
related by ln
B
s=n
ρ
2
=
2
[88]
.
Before discussing parameter estimates in detail, we
consider how the inference is affected by the choice of
the compact-binary waveform model. From Table
I
, we see
that the posterior estimates for each parameter are broadly
consistent across the two models, despite the fact that
they are based on different analytical approaches and that
they include different aspects of BBH spin dynamics. The
models
logarithms of the Bayes factors,
288
.
7

0
.
2
and
290
.
3

0
.
1
, are also comparable for both models: the data
do not allow us to conclusively prefer one model over the
other
[89]
. Therefore, we use both for the Overall column
in Table
I
. We combine the posterior samples of both
distributions with equal weight, in effect marginalizing
over our choice of waveform model. These averaged results
give our best estimate for the parameters describing
GW150914.
In Table
I
, we also indicate how sensitive our results are
to our choice of waveform. For each parameter, we give
systematic errors on the boundaries of the 90% credible
intervals due to the uncertainty in the waveform models
considered in the analysis; the quoted values are the 90%
range of a normal distribution estimated from the variance
of results from the different models. (If
X
were an edge of a
credible interval, we quote systematic uncertainty

1
.
64
σ
sys
using the estimate
σ
2
sys
¼½ð
X
EOBNR
X
Overall
Þ
2
þ
ð
X
IMRPhenom
X
Overall
Þ
2

=
2
. For parameters with bounded
ranges, like the spins, the normal distributions should
be truncated. However, for transparency, we still quote
the 90% range of the uncut distributions. These numbers
provide estimates of the order of magnitude of the potential
systematic error). Assuming a normally distributed error is
the least constraining choice
[90]
and gives a conservative
estimate. The uncertainty from waveform modeling is less
significant than the statistical uncertainty; therefore, we are
confident that the results are robust against this potential
systematic error. We consider this point in detail later in the
Letter.
The analysis presented here yields an optimal coherent
SNR of
ρ
¼
25
.
1
þ
1
.
7
1
.
7
. This value is higher than the one
reported by the search
[1,3]
because it is obtained using a
finer sampling of (a larger) parameter space.
GW150914
s source corresponds to a stellar-mass BBH
with individual source-frame masses
m
source
1
¼
36
þ
5
4
M
and
m
source
2
¼
29
þ
4
4
M
, as shown in Table
I
and Fig.
1
.
The two BHs are nearly equal mass. We bound the mass
ratio to the range
0
.
66
q
1
with 90% probability. For
comparison, the highest observed neutron star mass is
2
.
01

0
.
04
M
[91]
, and the conservative upper-limit for
FIG. 1. Posterior PDFs for the source-frame component masses
m
source
1
and
m
source
2
. We use the convention that
m
source
2
m
source
1
,
which produces the sharp cut in the two-dimensional distribution.
In the one-dimensional marginalized distributions we show the
Overall (solid black), IMRPhenom (blue), and EOBNR (red)
PDFs; the dashed vertical lines mark the 90% credible interval
for the Overall PDF. The two-dimensional plot shows the
contours of the 50% and 90% credible regions plotted over a
color-coded PDF.
PRL
116,
241102 (2016)
PHYSICAL REVIEW LETTERS
week ending
17 JUNE 2016
241102-6
the mass of a stable neutron star is
3
M
[92,93]
. The
masses inferred from GW150914 are an order of magnitude
larger than these values, which implies that these two
compact objects of GW150914 are BHs, unless exotic
alternatives, e.g., boson stars
[94]
, do exist. If the compact
objects were not BHs, this would leave an imprint on the
waveform, e.g., Ref.
[95]
; however, in Ref.
[96]
we verify
that the observed signal is consistent with that predicted
assuming BHs in general relativity. These results establish
the presence of stellar-mass BBHs in the Universe. It also
proves that BBHs formed in nature can merge within a
Hubble time
[97]
.
To convert the masses measured in the detector frame to
physical source-frame masses, we require the redshift of the
source. As discussed in the Introduction, GW observations
are directly sensitive to the luminosity distance to a source,
but not the redshift
[98]
. We find that GW150914 is at
D
L
¼
410
þ
160
180
Mpc. Assuming a flat
Λ
CDM cosmology
with Hubble parameter
H
0
¼
67
.
9
kms
1
Mpc
1
and matter
density parameter
Ω
m
¼
0
.
306
[6]
, the inferred luminosity
distance corresponds to a redshift of
z
¼
0
.
09
þ
0
.
03
0
.
04
.
The luminosity distance is strongly correlated to the
inclination of the orbital plane with respect to the line of
sight
[4,20,99]
. For precessing systems, the orientation of
the orbital plane is time dependent. We therefore describe
the source inclination by
θ
JN
, the angle between the total
angular momentum (which typically is approximately
constant throughout the inspiral) and the line of sight
[30,100]
, and we quote its value at a reference GW
frequency
f
ref
¼
20
Hz. The posterior PDF shows that
an orientation of the total orbital angular momentum of the
BBH strongly misaligned to the line of sight is disfavored;
the probability that
45
°
<
θ
JN
<
135
° is 0.35.
The masses and spins of the BHs in a (circular) binary
are the only parameters needed to determine the final mass
and spin of the BH that is produced at the end of the merger.
Appropriate relations are embedded intrinsically in the
waveform models used in the analysis, but they do not
give direct access to the parameters of the remnant BH.
However, applying the fitting formula calibrated to
nonprecessing NR simulations provided in Ref.
[101]
to
the posterior for the component masses and spins
[102]
,
we infer the mass and spin of the remnant BH to be
M
source
f
¼
62
þ
4
4
M
, and
a
f
¼
0
.
67
þ
0
.
05
0
.
07
, as shown in Fig.
3
and Table
I
. These results are fully consistent with those
obtained using an independent nonprecessing fit
[57]
. The
systematic uncertainties of the fit are much smaller than
the statistical uncertainties. The value of the final spin is a
consequence of conservation of angular momentum in
which the total angular momentum of the system (which
for a nearly equal mass binary, such as GW150914
s
source, is dominated by the orbital angular momentum)
is converted partially into the spin of the remnant black hole
and partially radiated away in GWs during the merger.
Therefore, the final spin is more precisely determined than
either of the spins of the binary
sBHs.
FIG. 2. Posterior PDFs for the source luminosity distance
D
L
and the binary inclination
θ
JN
. In the one-dimensional margin-
alized distributions we show the Overall (solid black), IMRPhe-
nom (blue), and EOBNR (red) PDFs; the dashed vertical lines
mark the 90% credible interval for the Overall PDF. The
two-dimensional plot shows the contours of the 50% and 90%
credible regions plotted over a color-coded PDF.
FIG. 3. PDFs for the source-frame mass and spin of the remnant
BH produced by the coalescence of the binary. In the one-
dimensional marginalized distributions we show the Overall
(solid black), IMRPhenom (blue), and EOBNR (red) PDFs;
the dashed vertical lines mark the 90% credible interval for the
Overall PDF. The two-dimensional plot shows the contours of the
50% and 90% credible regions plotted over a color-coded PDF.
PRL
116,
241102 (2016)
PHYSICAL REVIEW LETTERS
week ending
17 JUNE 2016
241102-7