June 7, 2016
Directly comparing GW150914 with numerical solutions of Einstein’s equations for binary black
hole coalescence
B. P. Abbott
et al.
∗
We compare GW150914 directly to simulations of coalescing binary black holes in full general relativity, in-
cluding 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 in LVC-PE
[1] (at
the 90% credible level), we find the data are compatible with a wide range of nonprecessing
and precessing
simulations. Followup simulations performed using previously-estimated binary parameters most resemble the
data, even when all quadrupolar and octopolar modes are included. Comparisons including only the quadrupo-
lar modes constrain the total redshifted mass
M
z
∈
[64
M
−
82
M
], mass ratio 1
/
q
=
m
2
/
m
1
∈
[0
.
6
,
1], and
e
ff
ective aligned spin
χ
e
ff
∈
[
−
0
.
3
,
0
.
2], where
χ
e
ff
=
(
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 e
ff
ective spins are highly inconsistent with the data, at any mass.
Several nonprecessing and precessing simulations with similar mass ratio and
χ
e
ff
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 followup 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 LVC-PE
[1].
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[2]. LVC-detect
[2],
LVC-PE
[1], LVC-TestGR
[3], and LVC-Burst
[4] demonstrate
consistency between GW150914 and selected individual pre-
dictions for a binary black hole coalescence, derived using
numerical solutions of Einstein’s equations for general rela-
tivity. LVC-PE
[1] described a systematic, Bayesian method
to reconstruct the properties of the coalescing binary, by com-
paring the data with the expected gravitational wave signa-
ture from binary black hole coalescence [5], evaluated using
state-of-the-art semianalytic 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
[1].
Instead, we compare the data directly with the most physi-
cally complete and generic predictions of general relativity:
computer simulations of binary black hole coalescence in full
nonlinear general relativity (henceforth referred to as numeri-
cal relativity, or NR). Although the semianalytic models are
calibrated to NR simulations, even the best available mod-
∗
Full author list given at the end of the article
els only imperfectly reproduce the predictions of numerical
relativity, on a mode-by-mode basis [9]. Furthermore, typi-
cal implementations of these models, such as those used in
LVC-PE
[1], consider only the dominant spherical-harmonic
mode of the waveform (in a corotating frame). For all NR
simulations considered here—including sources with two in-
dependent, precessing spins—we perform comparisons that
account for all the quadrupolar spherical-harmonic waveform
modes, and separately comparisons that account for all the
quadrupolar and octopolar spherical harmonic modes.
The principal approach introduced in this paper is di
ff
erent
from LVC-PE
[1], which inferred the properties of GW150914
by adopting analytic waveform models. Qualitatively speak-
ing, these models interpolate the
outgoing gravitational wave
strain
(waveforms) between the well-characterized results of
numerical relativity, as provided by a sparse grid of simula-
tions. 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 parame-
ters prescribed and provided by all available NR simulations.
We then construct an approximation to the marginalized like-
lihood that interpolates between NR simulations with di
ff
erent
parameters. To the extent that the likelihood is a simpler func-
tion of parameters than the waveforms, this method may re-
quire fewer NR simulations and fewer modeling assumptions.
Moreover, the interpolant for the likelihood needs to be accu-
rate only near its peak value, and not everywhere in parameter
space. A similar study was conducted on GW150914 using a
arXiv:1606.01262v1 [gr-qc] 3 Jun 2016
2
subset of numerical relativity waveforms directly against re-
constructed waveforms [4]; the results reported here are con-
sistent but more thorough.
Despite using an analysis that has few features or code in
common with the methods employed in LVC-PE
[1], we arrive
at similar conclusions regarding the parameters of the progen-
itor 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 corrobo-
ration of the results of LVC-PE
[1], strengthening our confi-
dence in both the employed statistical methods and waveform
models.
This paper is organized as follows. Section II provides
an extended motivation for and summary of this investiga-
tion. Section III A reviews the history of numerical relativity
and introduces notation to characterize 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
[10]
for further details. Section III D describes the implications
of using NR simulations that include only a small number of
gravitational-wave cycles. Section III F relates this investiga-
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 re-
sults on the post-coalescence state. Our statements rely on the
final black hole masses and spins derived from the full NR
simulations used. We summarize our results in Section VI. In
Appendix A, we summarize the simulations used in this work
and their accuracy, referring to the original literature for com-
plete 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 di
ff
er from those in LVC-PE
[1] in
two important ways. First, the statistical analysis here is per-
formed in a manner di
ff
erent than and independent of the one
in LVC-PE
[1]. Second, the gravitational waveform models
used in LVC-PE
[1] are analytic approximations of particu-
lar functional forms, with coe
ffi
cients calibrated to match se-
lected NR simulations; in contrast, here we directly use wave-
forms from NR simulations. Despite these di
ff
erences, our
conclusions largely corroborate the quantitative results found
in LVC-PE
[1].
Our study also addresses key challenges associated with
gravitational wave parameter estimation for black hole bi-
naries with total mass
M
>
50
M
. In this mass regime,
LIGO is sensitive to the last few dozens of cycles of coa-
lescence, a strongly nonlinear epoch that is the most di
ffi
-
cult to approximate with analytic (or semi-analytic) wave-
form models [6, 11–13]. 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 pre-
cession, lack higher-order spherical-harmonic modes, and do
not fully account for strong nonlinearities. Preliminary inves-
tigations have shown that inferences about the source drawn
using these existing analytic approximations can be slightly
or significantly biased [9, 13–15]. Systematic studies are un-
derway to assess how these approximations influence our best
estimates of a candidate binary’s parameters. 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 mod-
els for binary black hole coalescence [6] do not yet account for
the asymmetries [15] responsible for the largest gravitational-
wave recoil kicks [16, 17]. Second, the analytic waveform
models adopted in LVC-PE
[1] adopted simple spin treatments
(e.g., a binary with aligned spins, or a binary with single e
ff
ec-
tive precessing spin) that cannot capture the full spin dynam-
ics [18–20]. A single precessing spin is often a good approxi-
mation, particularly for unequal masses where one spin domi-
nates the angular momentum budget [8, 18, 21–23]. However,
for appropriate comparable-mass sources, two-spin e
ff
ects are
known to be observationally accessible [15, 24] and cannot
be fully captured by a single spin. Finally, LVC-PE
[1] and
LVC-TestGR
[3] made an additional approximation to connect
the inferred properties of the binary black hole with the final
black hole state [25, 26]. The analysis presented in this work
does not rely on these approximations: observational data are
directly compared against a wide range of NR simulations and
the final black hole properties are extracted directly from these
simulations, without recourse to estimated relationships be-
tween the initial and final state. By circumventing these ap-
proximations, 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, includ-
ing strongly precessing systems with significant misaligned
black hole spins [27], are reasonably consistent with the data.
The reason the data cannot distinguish between sources with
qualitatively di
ff
erent dynamics is a consequence of both the
orientation of the source relative to the line of sight and the
timescale 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, lacking unambiguous sig-
natures of precession and easily mistaken for a nonprecessing
binary [24, 28]. Second, because GW150914 has a large to-
tal 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 [18].
Although the orientation of the binary and the short du-
3
ration of the signal make it di
ffi
cult to extract spin informa-
tion 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, 29]. By utilizing full NR waveforms instead
of the single-spin (precessing) and double-spin (nonprecess-
ing) models applied in LVC-PE
[1], the approach described
here provides an independent opportunity to extract additional
insight from the data, or to independently corroborate the re-
sults of LVC-PE
[1].
Our study employs a simple method to carry out our
Bayesian calculations: for each NR simulation, we evalu-
ate the marginalized likelihood of the simulation parameters
given the data. The likelihood is evaluated via an adaptive
Monte Carlo integrator. This method provides a quantitative
ranking of simulations; with judicious interpolation in param-
eter space, the method also allows us to identify candidate
parameters for followup 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 observationally-motivated threshold
(i.e., large enough to contribute significantly to the posterior).
Even better, with a modest approximation to fill the gaps be-
tween NR simulations, we can reproduce and corroborate the
results in LVC-PE
[1] with a completely independent method.
We explicitly construct an approximation to the likelihood that
interpolates between simulations of precessing binaries, and
demonstrate its validity and utility.
It is well-known that the choice of prior may influence con-
clusions of Bayesian studies when the data do not strongly
constrain the relevant parameters. For example, the results
of LVC-PE
[1] suggest that GW150914 had low to moderate
spins, but this is due partly to the conventional prior used
in LVC-PE
[1] and earlier studies [5]. This prior is uniform
in spin
magnitude
, and therefore unfavorable to the most dy-
namically interesting possibilities: comparable-mass binaries
with two large, dynamically-significant precessing spins [24].
In contrast, by directly considering the (marginalized) likeli-
hood, the results of our study are independent of the choice of
prior. For example, we find here that GW150914 is consistent
with two large, dynamically significant spins.
Finally, our e
ff
orts to identify even subtle hints of spin
precession are motivated by the astrophysical opportuni-
ties a
ff
orded by spin measurements with GW150914; see
LVC-Astro
[30]. Using the geometric spin prior adopted in
LVC-PE
[1], the data from GW150914 are just as consistent
with an origin from a nonprecessing or precessing binary, as
long as 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 pos-
itively aligned at coalescence; see LVC-Astro
[30] and [31].
Hence, if we believe GW150914 formed from an isolated bi-
nary, our data would suggest black holes are born with small
spins:
a
1
=
|
S
1
|
/
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
[1]]. If these strictly-
aligned isolated evolution formation scenarios are true, then a
low black hole spin constrains the relevant accretion, angu-
lar momentum transport, and tidal interaction processes in the
progenitor binary; cf. [31–33]. On the other hand, the data
are equally consistent with a strongly-precessing black hole
binary with large component spins, formed in a densely inter-
acting stellar cluster (LVC-Astro
[30]). Measurements of the
binary black holes’ transverse spins will therefore provide vi-
tal 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
[1], we find results consis-
tent 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 rel-
ativity numerically began in the 1960s, by Hahn and Lindquist
[34], followed with some success by Smarr [35, 36]. In the
1990s, a large collaboration of US universities worked to-
gether to solve the “Grand Challenge” of evolving binary
black holes [37–39]. In 2005, three groups [40–42] developed
two completely independent techniques that produced the first
collisions of orbiting black holes. The first technique [40] in-
volved the use of generalized harmonic coordinates and exci-
sion of the black hole horizons, while the second technique
[41, 42], dubbed “moving punctures approach”, used singu-
larity 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., [43–
46]. Multiple approaches have been pursued and validated
against one another [47, 48]. Binaries can now be evolved in
wide orbits [49]; at high mass ratios up to 100:1 [50, 51];
with near-maximal black hole spin [52–54]; and for many or-
bits before coalescence [27, 55]. At su
ffi
ciently large separa-
tions, despite small gauge and frame ambiguities, the orbital
and spin dynamics evaluated using numerical relativity agrees
with post-Newtonian calculations [6, 56–59]. The stringent
phase and amplitude needs of gravitational wave detection and
parameter estimation prompted the development of revised
standards for waveform accuracy [60, 61]. Several projects
have employed numerical relativity-generated waveforms to
assess gravitational-wave detection and parameter estimation
strategies [11, 62–66]. These results have been used to cali-
brate models for the leading-order radiation emitted from bi-
nary black hole coalescence [6, 8, 9, 12, 57, 67, 68]. Our
study builds on this past decade’s experience with model-
ing the observationally-relevant dynamics and radiation from
comparable-mass coalescing black holes.
In this and most NR work, the initial data for a simulation
of 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 quasicir-
cular orbit such that the (coordinate) orbital angular momen-
tum is aligned with the ˆ
z
axis and the initial separation is along
the ˆ
x
axis. In this work, we characterize these simulations by
4
the dimensionless mass ratio
q
=
m
1
/
m
2
(adopting the con-
vention
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-dependent gravitational
wave strain
h
(
t
,
r
,
ˆ
n
) at large distances can be e
ffi
ciently char-
acterized by a (spin-weighted) spherical harmonic decomposi-
tion of
h
(
t
,
r
,
ˆ
n
) as
h
(
t
,
r
,
ˆ
n
)
=
∑
l
≥
2
∑
l
m
=
−
l
h
lm
(
t
,
r
)
−
2
Y
lm
( ˆ
n
). To
a good first approximation, only a few terms in this series are
necessary to characterize the observationally-accessible radi-
ation in any direction [14, 69–72]. For example, when a bi-
nary is widely separated, only two terms dominate this sum:
(
l
,
|
m
|
)
=
(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 sym-
metry, individual modes (
l
,
m
) have distinctive characters, and
can be easily isolated numerically and compared with ana-
lytic 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 in-
ferred 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 sim-
ulation 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, necessary and su
ffi
cient
to capture the leading-order influence of any orbital preces-
sion. To incorporate the e
ff
ect 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
[10] and subsequent work.
B. Numerical relativity simulations used
Our study makes use of 1139 distinct simulations of bi-
nary black hole quasicircular inspiral and coalescence. Ta-
ble II summarizes the salient features of this set: 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 [17, 27, 73]. The
simulations include binaries with a wide range of mass ra-
tios, as well as a wide range of black hole angular momentum
(spin) magnitudes and directions [17, 27, 73, 74], including
a simulation with large transverse spins and several spin pre-
cession cycles which fits GW150914 well [27], as described
below. The SXS group has provided both a publicly-available
catalog of coalescing black hole binary mergers [75], a new
catalog of nonprecessing simulations [76], and selected sup-
plementary simulations described below. Currently extended
to 310 members in the form used here, this catalog includes
many high-precision zero- and aligned-spin sources; selected
precessing systems; and simulations including extremely high
black hole spin. The Georgia Tech group (GT) provided
406 simulations; see [15] and [77] for further details. This
extensive archive covers a wide range of spin magnitudes
and orientations, including several systematic one- and two-
parameter families. The Cardi
ff
-UIB group provided 29 sim-
ulations, all specifically produced to follow up GW150914 via
a high-dimensional grid stencil, performed via the BAM code
[78, 79]. These four sets of simulations explore the model
space near the event in a well-controlled fashion. In addition
to previously-reported simulations, several groups performed
new simulations (108 in total) designed to reproduce the pa-
rameters of the event, some of which were applied to our anal-
ysis. These simulations are denoted in Table II and our other
reports by an asterisk (*). These followup simulations include
three independent simulations of the same parameters drawn
from the distributions in LVC-PE
[1], from RIT, SXS, and GT,
allowing us to assess our systematic error. These simulations
were reported in LVC-detect
[2] and LVC-Burst
[4], and are in-
dicated by (
+
) in our tables.
The simulations used here have either been published pre-
viously, or were produced using one of three well-tested
procedures operating in familiar circumstances. For refer-
ence, in Appendix A, we outline the three groups’ previously-
established methods and results. For this application, we trust
these simulations’ accuracy, based on their past track record
of good performance. By incorporating simulations of identi-
cal physics provided by di
ff
erent groups, our methods provide
limited direct corroboration: simulations with similar physics
produce similar results.
C. Directly comparing NR with data
For each simulation, each choice of seven extrinsic param-
eters
θ
(4 spacetime coordinates for the coalescence event;
three Euler angles for the binary’s orientation relative to the
Earth), and each choice for the redshifted total binary mass
M
z
=
(1
+
z
)
M
, we can predict the response
h
k
of both of the
k
=
1
,
2 LIGO instruments to the implied gravitational wave
signal. Using
λ
to denote the combination of redshifted mass
M
z
and the numerical relativity simulation parameters needed
to uniquely specify the binary’s dynamics, we can therefore
evaluate the likelihood of the data given the noise:
ln
L
(
λ
;
θ
)
=
−
1
2
∑
k
〈
h
k
(
λ,θ
)
−
d
k
|
h
k
(
λ,θ
)
−
d
k
〉
k
−
〈
d
k
|
d
k
〉
k
,
(1)
where
h
k
are the predicted response of the
k
th detector due to a
source with parameters
λ,θ
;
d
k
are the detector data in instru-
ment
k
; and
〈
a
|
b
〉
k
≡
∫
∞
−∞
2
d f
̃
a
(
f
)
∗
̃
b
(
f
)
/
S
h
,
k
(
|
f
|
) is an inner
product implied by the
k
th’s detector’s noise power spectrum
S
h
,
k
(
f
); see, e.g., [80] for more details. In practice, as dis-
cussed in the next section, we adopt a low-frequency cuto
ff