of 30
Directly comparing GW150914 with numerical solutions
of Einstein
s equations for binary black hole coalescence
B. P. Abbott
etal.
*
(LIGO Scientific Collaboration and Virgo Collaboration)
(Received 10 June 2016; published 14 September 2016)
We compare GW150914 directly to simulations of coalescing binary black holes in full general relativity,
including several performed specifically to reproduce this event. Our calculations go beyond existing
semianalytic models, because for all simulations
including sources with two independent, precessing spins
we perform comparisons which account for all the spin-weighted quadrupolar modes, and separately
which account for all the quadrupolar and octopolar modes. Consistent with the posterior distributions
reported by Abbott
et al.
[Phys. Rev. Lett.
116
, 241102 (2016)] (at the 90% credible level), we find the data
are compatible with a wide range of nonprecessing
and precessing
simulations. Follow-up simulations
performed using previously estimated binary parameters most resemble the data, even when all quadrupolar
and octopolar modes are included. Comparisons including only the quadrupolar modes constrain the total
redshifted mass
M
z
½
64
M
82
M

, mass ratio
1
=q
¼
m
2
=m
1
½
0
.
6
;
1

, and effective aligned spin
χ
eff
½
0
.
3
;
0
.
2

, where
χ
eff
¼ð
S
1
=m
1
þ
S
2
=m
2
Þ
·
ˆ
L
=M
. Including both quadrupolar and octopolar
modes, we find the mass ratio is even more tightly constrained. Even accounting for precession, simulations
with extreme mass ratios and effective spins are highly inconsistent with the data, at any mass. Several
nonprecessing and precessing simulations with similar mass ratio and
χ
eff
are consistent with the data.
Though correlated, the components
spins (both in magnitude and directions) are not significantly
constrained by the data: the data is consistent with simulations with component spin magnitudes
a
1
;
2
up
to at least 0.8, with random orientations. Further detailed follow-up calculations are needed to determine if
the data contain a weak imprint from transverse (precessing) spins. For nonprecessing binaries, interpolating
between simulations, we reconstruct a posterior distribution consistent with previous results. The final black
hole
s redshifted mass is consistent with
M
f;z
in the range
64
.
0
M
73
.
5
M
and the final black hole
s
dimensionless spin parameter is consistent with
a
f
¼
0
.
62
0
.
73
. As our approach invokes no intermediate
approximations to general relativity and can strongly reject binaries whose radiation is inconsistent with the
data, our analysis provides a valuable complement to Abbott
et al.
[Phys. Rev. Lett.
116
, 241102 (2016)].
DOI:
10.1103/PhysRevD.94.064035
I. INTRODUCTION
On September 14, 2015
09
50
45
UTC, gravitational
waves were observed in coincidence by the twin instruments
of the Laser Interferometer Gravitational-wave Observatory
(LIGO) located at Hanford, Washington, and Livingston,
Louisiana, in the USA, an event known as GW150914
[1]
.
LVC-detect
[1]
,
LVC-PE
[2]
,
LVC-TestGR
[3]
, and
LVC-
Burst
[4]
demonstrated consistency between GW150914
and selected individual predictions for a binary black hole
coalescence, derived using numerical solutions of Einstein
s
equations for general relativity.
LVC-PE
[2]
described a
systematic, Bayesian method to reconstruct the properties of
the coalescing binary, by comparing the data with the
expected gravitational wave signature from binary black
hole coalescence
[5]
, evaluated using state-of-the-art semi-
analytic approximations to its dynamics and radiation
[6
8]
.
In this paper, we present an alternative method of recon-
structing the binary parameters of GW150914, without using
the semianalytic waveform models employed in
LVC-PE
[2]
.
Instead, we compare the data directly with the most physically
complete and generic predictions of general relativity: com-
puter simulations of binary black hole coalescence in full
nonlineargeneralrelativity (henceforthreferred toasnumerical
relativity, or NR). Although the semianalytic models are
calibrated to NR simulations, even the best available models
only imperfectly reproduce the predictions of numerical
relativity, on a mode-by-mode basis
[9,10]
. Furthermore,
typical implementations of these models, such as those used
in
LVC-PE
[2]
, consider only the dominant spherical-harmonic
mode of the waveform (in a corotating frame). For all NR
simulations considered here
including sources with two
independent, precessing spins
we perform comparisons that
account for all the quadrupolar s
pherical-harmonic waveform
modes, and separately comparis
ons that account for all the
quadrupolar and octopolar s
pherical-harmonic modes.
The principal approach introduced in this paper is
different from
LVC-PE
[2]
, which inferred the properties
of GW150914 by adopting analytic waveform models.
Qualitatively speaking, these models interpolate the
outgoing gravitational wave strain
(waveforms) between
the well-characterized results of numerical relativity, as
*
Full author list given at the end of the article.
PHYSICAL REVIEW D
94,
064035 (2016)
2470-0010
=
2016
=
94(6)
=
064035(30)
064035-1
© 2016 American Physical Society
provided by a sparse grid of simulations. These interpolated
or analytic waveforms are used to generate a continuous
posterior distribution over the binary
s parameters. By
contrast, in this study, we compare numerical relativity to
the data first, evaluating a single scalar quantity (the
marginalized likelihood) on the grid of binary parameters
prescribed and provided by all available NR simulations. We
then construct an approximation to the marginalized like-
lihood that interpolates between NR simulations with differ-
ent parameters. To the extent that the likelihood is a simpler
function of parameters than thewaveforms, this method may
require fewer NR simulations and fewer modeling assump-
tions. Moreover, the interpolant for the likelihood needs to
be accurate only near its peak value, and not everywhere in
parameter space. A similar study was conducted on
GW150914 using a subset of numerical relativity wave-
forms directly against reconstructed waveforms
[4]
; the
results reported here are consistent but more thorough.
Despite using an analysis that has few features or code in
common with the methods employed in
LVC-PE
[2]
,we
arrive at similar conclusions regarding the parameters of the
progenitor black holes and the final remnant, although we
extract slightly more information about the binary mass
ratio by using higher-order modes. Thus, we provide
independent corroboration of the results of
LVC-PE
[2]
,
strengthening our confidence in both the employed stat-
istical methods and waveform models.
This paper is organized as follows. Section
II
provides an
extended motivation for and summary of this investigation.
Section
III A
reviews the history of numerical relativity and
introducesnotation to character
ize simulated binaries and their
radiation. Section
III B
describes the simulations used in this
work. Section
III C
briefly describes the method used to
compare simulations to the data; see
PE+NR-Methods
[11]
for
further details. Section
III D
describes the implications of
using NR simulations that in
clude only a small number of
gravitational-wave cycles. Section
III F
relatesthisinvestiga-
tion to prior work. Section
IV
describes our results on the pre-
coalescence parameters. We provide a ranking of simulations
as measured by a simple measure of fit (peak marginalized log
likelihood). When possible, we provide an approximate
posterior distribution over all intrinsic parameters. Using both
our simple ranking and approximate posterior distributions,
we draw conclusions about the range of source parameters that
are consistent with the data. Section
V
describes our results on
the post-coalescence state. Our statements rely on the final
black hole (BH) masses and spins derived from the full NR
simulations used. We summarize our results in Sec.
VI
.In
Appendix
A
, we summarize the simulations used in this work
and their accuracy, referring to the original literature for
complete details.
II. MOTIVATION FOR THIS STUDY
This paper presents an alternative analysis of GW150914
and an alternative determination of its intrinsic parameters.
The methods used here differ from those in
LVC-PE
[2]
in
two important ways. First, the statistical analysis here is
performed in a manner different than and independent of
the one in
LVC-PE
[2]
. Second, the gravitational waveform
models used in
LVC-PE
[2]
are analytic approximations of
particular functional forms, with coefficients calibrated to
match selected NR simulations; in contrast, here we
directly use waveforms from NR simulations. Despite these
differences, our conclusions largely corroborate the quan-
titative results found in
LVC-PE
[2]
.
Our study also addresses key challenges associated with
gravitational wave parameter estimation for black hole
binaries with total mass
M>
50
M
.Inthismassregime,
LIGO is sensitive to the last few dozens of cycles of
coalescence, a strongly nonlinear epoch that is the most
difficult to approximate with analytic (or semianalytic) wave-
form models
[6,12
14]
. Figure
1
illustrates the dynamics and
expected detector response in this regime, for a source like
GW150914. For these last dozen cycles, existing analytic
waveform models have only incomplete descriptions of
precession, lack higher-order spherical-harmonic modes,
and do not fully account for strong nonlinearities.
Preliminary investigations have shown that inferences about
thesourcedrawnusingthese existinganalytic approximations
can be slightly or significantly biased
[9,14
16]
. Systematic
studies are underway to assess how these approximations
influence our best estimates of a candidate binary
sparam-
eters. At present, we can only summarize the rationale for
these investigations, not their results. To provide three
concrete examples of omitted physics, first and foremost,
even the most sophisticated models for binary black hole
coalescence
[6]
do not yet account for the asymmetries
[16]
responsible for the largest gravitational-wave recoil kicks
-2
-1
0
1
2
-0.2
-0.15
-0.1
-0.05
0
0.05
n
w
o
d
g
n
i
R
r
e
g
r
e
M
l
a
r
i
p
s
n
I
Strain (10
-21
)
Time (seconds)
H1 estimated strain, incident
H1 estimated strain, band-passed
FIG. 1.
Simulated waveform
: Predicted strain in H1 for a source
with parameters
q
¼
1
.
22
,
χ
1
;z
¼
0
.
33
,
χ
2
;z
¼
0
.
44
, simulated in
full general relativity; compare to Fig.
2
in
LVC-detect
[1]
. The
gray line shows the idealized strain response
h
ð
t
Þ¼
F
þ
h
þ
ð
t
Þþ
F
×
h
×
ð
t
Þ
, while the solid black line shows the
whitened strain response, using the same noise power spectrum
as
LVC-detect
[1]
.
B. P. ABBOTT
et al.
PHYSICAL REVIEW D
94,
064035 (2016)
064035-2
[17,18]
. Second, the analytic waveform models adopted in
LVC-PE
[2]
adopted simple spin treatments (e.g., a binary
with aligned spins, or a binary with single effective precessing
spin) that cannot capture the full spin dynamics
[19
21]
.A
single precessing spin is often a good approximation,
particularly for unequal masses where one spin dominates
the angular momentum budget
[8,19,22
24]
.However,for
appropriate comparable-mass sources, two-spin effects are
known to be observationally accessible
[16,25]
and cannot be
fully captured by a single spin. Finally,
LVC-PE
[2]
and
LVC-
TestGR
[3]
made an additional approximation to connect the
inferred properties of the binary black holewith thefinal black
hole state
[26,27]
. The analysis presented in this work does
not rely on these approximations: observational data are
directly compared against awiderange ofNR simulations and
thefinal blackholeproperties areextracted directlyfromthese
simulations, without recourse to estimated relationships
between the initial and final state. By circumventing these
approximations, our analysis can corroborate conclusions
about selected physical properties of GW150914 presented in
those papers.
Despite the apparent simplicity of GW150914, we find
that a range of binary black hole masses and spins,
including strongly precessing systems with significant
misaligned black hole spins
[28]
, are reasonably consistent
with the data. The reason the data cannot distinguish
between sources with qualitatively different dynamics is
a consequence of both the orientation of the source relative
to the line of sight and the time scale of GW150914. First, if
the line of sight is along or opposite the total angular
momentum vector of the source, even the most strongly
precessing black hole binary emits a weakly modulated
inspiral signal, which lacks unambiguous signatures of
precession and is easily mistaken for a nonprecessing
binary
[25,29]
. Second, because GW150914 has a large
total mass, very little of the inspiral lies in LIGO
s
frequency band, so the signal is short, with few orbital
cycles and even fewer precession cycles prior to or during
coalescence. The short duration of GW150914 provides
few opportunities for the dynamics of two precessing spins
to introduce distinctive amplitude and phase modulations
into its gravitational wave inspiral signal
[19]
.
Although the orientation of the binary and the short
duration of the signal make it difficult to extract spin
information from the
inspiral
, comparable-mass binaries
with large spins can have exceptionally rich dynamics with
waveform signatures that extend into the late inspiral and
the strong-field merger phase
[8,30]
. By utilizing full NR
waveforms instead of the single-spin (precessing) and
double-spin (nonprecessing) models applied in
LVC-PE
[2]
, the approach described here provides an independent
opportunity to extract additional insight from the data, or to
independently corroborate the results of
LVC-PE
[2]
.
Our study employs a simple method to carry out our
Bayesian calculations: for each NR simulation, we evaluate
the marginalized likelihood of the simulation parameters
given the data. The likelihood is evaluated via an adaptive
Monte Carlo integrator. This method provides a quantita-
tive ranking of simulations; with judicious interpolation in
parameter space, the method also allows us to identify
candidate parameters for follow-up numerical relativity
simulations. To estimate parameters of GW150914, we
can simply select the subset of simulations and masses that
have a marginalized likelihood greater than an observatio-
nally motivated threshold (i.e., large enough to contribute
significantly to the posterior). Even better, with a modest
approximation to fill the gaps between NR simulations, we
can reproduce and corroborate the results in
LVC-PE
[2]
with a completely independent method. We explicitly
construct an approximation to the likelihood that interpo-
lates between simulations of precessing binaries, and
demonstrate its validity and utility.
It is well known that the choice of prior may influence
conclusions of Bayesian studies when the data do not
strongly constrain the relevant parameters. For example, the
results of
LVC-PE
[2]
suggest that GW150914 had low to
moderate spins, but this is due partly to the conventional
prior used in
LVC-PE
[2]
and earlier studies
[5]
. This prior
is uniform in spin
magnitude
, and therefore unfavorable to
the most dynamically interesting possibilities: comparable-
mass binaries with two large, dynamically significant
precessing spins
[25]
. In contrast, by directly considering
the (marginalized) likelihood, the results of our study are
independent of the choice of prior. For example, we find
here that GW150914 is consistent with two large, dynami-
cally significant spins.
Finally, our efforts to identify even subtle hints of spin
precession are motivated by the astrophysical opportunities
afforded by spin measurements with GW150914; see
LVC-
Astro
[31]
. Using the geometric spin prior adopted in
LVC-
PE
[2]
, the data from GW150914 are just as consistent with
an origin from a nonprecessing or precessing binary, as
long as the sum of the components of the spins parallel to
the orbital angular momentum
L
is nearly zero. If the
binary black hole formed from isolated stellar evolution,
one could reasonably expect all angular momenta to be
strictly and positively aligned at coalescence; see
LVC-
Astro
[31]
and
[32]
. Hence, if we believe GW150914
formed from an isolated binary, our data would suggest
black holes are born with small spins:
a
1
¼j
S
1
j
=m
2
1
0
.
22
and
a
2
0
.
25
, where
S
i
and
m
i
are the black hole spins and
masses (
LVC-PE
[2]
). If these strictly aligned isolated
evolution formation scenarios are true, then a low black
hole spin constrains the relevant accretion, angular momen-
tum transport, and tidal interaction processes in the
progenitor binary; cf. Refs.
[32
34]
. On the other hand,
the data are equally consistent with a strongly precessing
black hole binary with large component spins, formed in a
densely interacting stellar cluster (
LVC-Astro
[31]
).
Measurements of the binary black holes
transverse spins
DIRECTLY COMPARING GW150914 WITH NUMERICAL
...
PHYSICAL REVIEW D
94,
064035 (2016)
064035-3
will therefore provide vital clues as to the processes that
formed GW150914. In this work we use numerical
relativity to check for any evidence for or against spin
precession that might otherwise be obscured by model
systematics. Like
LVC-PE
[2]
, we find results consistent
with but with no strong support for precessing spins.
III. METHODS
A. Numerical relativity simulations of binary black
hole coalescence
The first attempts to solve the field equations of general
relativity numerically began in the 1960s, by Hahn and
Lindquist
[35]
, followed with some success by Smarr
[36,37]
. In the 1990s, a large collaboration of US universities
worked together to solve the
Grand Challenge
of evolving
binary black holes
[38
40]
. In 2005, three groups
[41
43]
developed two completely independent techniques that pro-
duced the first collisions of orbiting black holes. The first
technique
[41]
involved the use of generalized harmonic
coordinates and the excision of the black hole horizons, while
the second technique
[42,43]
, dubbed the
moving punctures
approach,
used singularity-avoiding slices of the black hole
spacetimes.
Since then, the field has seen an explosion of activity and
improvements in methods and capabilities; see, e.g.,
Refs.
[44
47]
. Multiple approaches have been pursued and
validated against one another
[10,48,49]
. Binaries can now be
evolved in wide orbits
[50]
,athighmassratiosupto
100
1
[51,52]
, with near-maximal black hole spin
[53
55]
,andfor
many orbits before coalescence
[28,56]
. At sufficiently large
separations, despite small ga
uge and frame ambiguities, the
orbital and spin dynamics evalua
ted using numerical relativity
agrees with post-Newtonian calculations
[6,57
60]
.The
stringent phase and amplitude ne
eds of gravitational wave
detectionandparameterestimationprompted thedevelopment
of revised standards for waveform accuracy
[61,62]
. Several
projects have employed numerical relativity-generated wave-
forms to assess gravitational-
wave detection and parameter
estimation strategies
[12,63
67]
. These results have been used
to calibrate models for the leading-order radiation emitted
from binary black hole coalescence
[6,8,9,13,58,68,69]
.Our
study builds on this past decade
s experience with modeling
the observationally relevant dynamics and radiation from
comparable-mass coalescing black holes.
Inthis and most NR work, theinitial data for a simulationof
coalescing binaries are characterized by the properties and
initial orbit of its two component black holes: by initial black
holes masses
m
1
,
m
2
and spins
S
1
,
S
2
, specified in a
quasicircular orbit such that the (coordinate) orbital angular
momentum is aligned with the
ˆ
z
axis and the initial separation
is along the
ˆ
x
axis. In this work, we characterize these
simulations by the dimensionless mass ratio
q
¼
m
1
=m
2
(adopting the convention
m
1
m
2
), the dimensionless spin
parameters
χ
i
¼
S
i
=m
2
i
, and an initial dimensionless orbital
frequency
M
ω
0
. For each simulation, the orientation-depen-
dent gravitational wave strain
h
ð
t; r;
ˆ
n
Þ
at large distances can
be efficiently characterized by a (spin-weighted) spherical
harmonic decomposition of
h
ð
t; r;
ˆ
n
Þ
as
h
ð
t; r;
ˆ
n
Þ¼
P
l
2
P
l
m
¼
l
h
lm
ð
t; r
Þ
2
Y
lm
ð
ˆ
n
Þ
. To a good first appro-
ximation, only a few terms in this series are necessary to
characterize the observationally accessible radiation in any
direction
[15,70
73]
. For example, when a binary is widely
separated, only two terms dominate this sum:
ð
l;
j
m
jÞ ¼
ð
2
;
2
Þ
. Conversely, however, several terms (modes) are
required for even nonprecessing binaries, viewed along a
generic line of sight; more are needed to capture the radiation
from precessing binaries.
For nonprecessing sources with a well-defined axis of
symmetry, individual modes
ð
l; m
Þ
have distinctive char-
acters, and can be easily isolated numerically and compared
with analytic predictions. For precessing sources, however,
rotation mixes modes with the same
l
. To apply our
procedure self-consistently to both nonprecessing and
precessing sources, we include all modes
ð
l; m
Þ
with the
same
l
. However, at the start of each simulation, the
ð
l; m
Þ
mode oscillates at
m
times the orbital frequency. For
m>
3
,
scaling our simulations to the inferred mass of the source,
this initial mode frequency is often well above 30 Hz, the
minimum frequency we adopt in this work for parameter
estimation. We therefore cannot safely and self-consistently
compare all modes with
l>
3
to the data using numerical
relativity alone: an approximation would be required to go
to higher order (i.e., hybridizing each NR simulation with
an analytic approximation at early times).
Therefore, in this paper, we use all five of the
l
¼
2
modes to draw conclusions about GW150914, which is
necessary and sufficient to capture the leading-order
influence of any orbital precession. To incorporate the
effect of higher harmonics, we repeat our calculations,
using all of the
l
3
modes. We defer a careful treatment of
higher-order modes and the
m
¼
0
modes to
PE+NR-
Methods
[11]
and subsequent work.
B. Numerical relativity simulations used
Our study makes use of 1139 distinct simulations of
binary black hole quasicircular inspiral and coalescence.
Table
II
summarizes the salient features of this set, namely,
the mass ratio and initial spins for the simulations used
here, all initially in a quasicircular orbit with orbital
separation along the
ˆ
x
axis and velocities along

ˆ
y
.
The RIT group provided 394 simulations. The simulations
include binaries with a wide range of mass ratios, as well as a
wide range of black hole angular momentum (spin) magni-
tudes and directions
[18,28,74
76]
, including a simulation
with large transverse spins and several spin precession cycles
which fits GW150914 well
[28]
, as described below. The
SXS group has provided both a publicly available catalog of
coalescing black hole binary mergers
[59]
, a new catalog of
nonprecessing simulations
[77]
, and selected supplementary
B. P. ABBOTT
et al.
PHYSICAL REVIEW D
94,
064035 (2016)
064035-4