Black hole puncture initial data with realistic gravitational wave content
B. J. Kelly,1, 2 W. Tichy,3 M. Campanelli,4, 2 and B. F. Whiting5, 2
1Gravitational Astrophysics Laboratory, NASA Goddard Space Flight Center, 8800 Greenbelt Rd., Greenbelt, MD 20771, USA
2Center for Gravitational Wave Astronomy, Department of Physics and Astronomy,
The University of Texas at Brownsville, Brownsville, Texas 78520
3Department of Physics, Florida Atlantic University, Boca Raton Florida 33431-0991
4Center for Computational Relativity and Gravitation,
School of Mathematical Sciences, Rochester Institute of Technology,
78 Lomb Memorial Drive, Rochester, New York 14623
5Department of Physics, University of Florida, Gainsville Florida 32611-8440
(Dated: October 26, 2018)
We present improved post-Newtonian-inspired initial data for non-spinning black-hole binaries,
suitable for numerical evolution with punctures. We revisit the work of Tichy et al. [W. Tichy, B.
Brügmann, M. Campanelli, and P. Diener, Phys. Rev. D 67, 064008 (2003)], explicitly calculating
the remaining integral terms. These terms improve accuracy in the far zone and, for the first
time, include realistic gravitational waves in the initial data. We investigate the behavior of these
data both at the center of mass and in the far zone, demonstrating agreement of the transverse-
traceless parts of the new metric with quadrupole-approximation waveforms. These data can be
used for numerical evolutions, enabling a direct connection between the merger waveforms and the
post-Newtonian inspiral waveforms.
PACS numbers: 04.25.Dm, 04.25.Nx, 04.30.Db, 04.70.Bw
I. INTRODUCTION
Post-Newtonian (PN) methods have played a funda-
mental role in our understanding of the astrophysical im-
plications of Einstein’s theory of general relativity. Most
importantly, they have been used to confirm that the ra-
diation of gravitational waves accounts for energy loss in
known binary pulsar configurations. They have also been
used to create templates for the gravitational waves emit-
ted from compact binaries which might be detected by
ground-based gravitational wave observatories, such as
LIGO [1, 2], and the NASA/ESA planned space-based
mission, LISA [3, 4]. However, PN methods have not
been extensively used to provide initial data for binary
evolution in numerical relativity, nor, until recently (see
[5, 6]), have they been extensively studied so that their
limitations could be well identified and the results of nu-
merical relativity independently confirmed.
Until the end of 2004, the field of numerical relativ-
ity had been struggling to compute even a single or-
bit for a black-hole binary (BHB). Although debate oc-
curred on the advantages of one type of initial data over
another, the primary focus within the numerical rela-
tivity community was on code refinement which would
lead to more stable evolution. Astrophysical realism was
very much a secondary issue. However, this situation
has radically changed in the last few years with the in-
troduction of two essentially independent, but equally
successful techniques: the generalized harmonic gauge
(GHG) method developed by Pretorius [7] and the “mov-
ing puncture” approach, independently developed by the
UTB and NASA Goddard groups [8, 9]. Originally in-
troduced by Brandt & Brügmann [10] in the context of
initial data, the puncture method explicitly factored out
the singular part of the metric. When used in numerical
evolution in which the punctures remained fixed on the
numerical grid, it resulted in distortions of the coordi-
nate system and instabilities in the Baumgarte-Shapiro-
Shibata-Nakamura (BSSN) [11, 12] evolution scheme.
The revolutionary idea behind the moving puncture ap-
proach was precisely, not to factor out the singular part
of the metric, but rather evolve it together with the reg-
ular part, allowing the punctures to move freely across
the grid with a suitable choice of the gauge.
A golden age for numerical relativity is now emerging,
in which multiple groups are using different computer
codes to evolve BHBs for several orbits before plunge and
merger [13, 14, 15, 16, 17, 18, 19, 20, 21]. Comparison of
the numerical results obtained from these various codes
has taken place [22, 23, 24], and comparison with PN
inspiral waveforms has also been carried out with encour-
aging success [5, 6, 25, 26]. The application of successful
numerical relativity tools to study some important as-
trophysical properties (e.g. precession, recoil, spin-orbit
coupling, elliptical orbits, etc) of spinning and/or un-
equal mass-black hole systems is currently producing ex-
tremely interesting new results [27, 28, 29, 30, 31, 32,
33, 34, 35, 36, 37, 38, 39, 40, 41, 42]. It now seems that
the primary obstacle to further progress is simply one of
computing power. In this new situation, it is perhaps
time to return to the question of what initial data will
best describe an astrophysical BHB.
To date, the best-motivated description of pre-merger
BHBs has been supplied by PN methods. We might ex-
pect, then, that a PN-based approach would give us the
most astrophysically correct initial data from which to
run full numerical simulations. In practice, PN results
are frequently obtained in a form ill-adapted to numeri-
http://arxiv.org/abs/0704.0628v2
cal evolution. PN analysis often deals with the full four-
metric, in harmonic coordinates; numerical evolutions
frequently use ADM-type coordinates, with a canonical
decomposition of the four-metric into a spatial metric
and extrinsic curvature.
Fortunately, many PN results have been translated
into the language of ADM by Ohta, Damour, Schäfer and
collaborators. Explicit results for 2.5PN BHB data in the
near zone were given by Schäfer [43] and Jaranowski &
Schäfer (JS) [44], and these were implemented numer-
ically by Tichy et al. [45]. Their insight was that the
ADM-transverse-traceless (TT) gauge used by Schäfer
was well-adapted to a puncture approach. To facilitate
comparison with this earlier work [45], we continue to
use the results of Schäfer and co-workers, anticipating
that higher-order PN results should eventually become
available in a useful form.
The initial data provided previously by Tichy et
al. already include PN information. They are accurate
up to order (v/c)5 in the near zone (r ≪ λ), but the
accuracy drops to order (v/c)3 in the far zone (r ≫ λ)
[here λ ∼ π
r312/G(m1 +m2) is the gravitational wave-
length]. These data were incomplete in the sense that
they did not include the correct TT radiative piece in
the metric, and thus did not contain realistic gravita-
tional waves.
In this paper, we revisit the PN data problem in
ADM-TT coordinates, with the aim of supplying Numer-
ical Relativity with initial BHB data that extend as far
as necessary, and contain realistic gravitational waves.
To do this, we have evaluated the “missing pieces” of
Schäfer’s TT metric for the case of two non-spinning par-
ticles. We have analyzed the near- and far-zone behavior
of these data, and incorporated them numerically in the
Cactus [46] framework. In principle, the most accurate
PN metric available could be used at this step, but it is
not currently available in ADM-TT form.
The remainder of this paper is laid out as follows. In
Section II, we summarize the results of Schäfer (1985)
[43], and Jaranowski & Schäfer (1997) [44] and their ap-
plication by Tichy et al. (2003) [45], to the production
of puncture data for numerical evolution. In Section
III, we describe briefly the additional terms necessary
to complete hTT to order (v/c)4, deferring details to the
Appendix. In Section IV, we study the full data both
analytically and numerically. Section V summarizes our
results, and lays the groundwork for numerical evolution
of these data, to be presented in a subsequent article.
II. ADM-TT GAUGE IN POST-NEWTONIAN
The “ADM-TT” gauge [43, 47] is a 3+1 split of data
where the three-metric differs from conformal flatness
precisely by a TT radiative part:
gij =
ηij + h
ij , (1)
πii = 0. (2)
The fields φ, πij and hTTij can all be expanded in a post-
Newtonian series. Solving the constraint equations of
3+1 general relativity in this gauge, [43, 44] obtained ex-
plicit expressions valid up to O(v/c)5 in the near zone,
incorporating an arbitrary number of spinless point par-
ticles, with arbitrary masses mA. For N particles, the
lowest-order contribution to the conformal factor is1:
φ(2) = 4G
, (3)
where rA =
~x− ~xA is the distance from the field point
to the location of particle A.
In principle hTTij is computed from
hTTij = −δTT klij ✷−1retskl, (4)
where ✷−1ret is the (flat space) inverse d’Alembertian (with
a “no-incoming-radiation” condition [48]), skl is a non-
local source term and δTT klij is the TT-projection oper-
ator. In order to compute hTTij we first rewrite Eq. (4)
hTTij = −δTT klij
∆−1 + (✷−1ret −∆−1)
TT (NZ)
ij − δ
TT kl
ij (✷
ret −∆−1)skl. (5)
Note that the near-zone approximation h
TT (NZ)
ij of h
has already been computed in [43] up to order O(v/c)4
(see also Eq. 12 below). The last term in Eq. (5) is diffi-
cult to compute because
skl = 16πG
pAk pAl
δ(x− xA) +
,l (6)
is a non-local source. However, we can approximate skl
s̄kl =
pAk pAl
B 6=A
nABk nABl
×16πGδ(x− xA). (7)
and show that
hTTij,(div) = −δ
TT kl
ij (✷
ret−∆−1)(skl− s̄kl) ∼ O(v/c)5 (8)
1 We explicitly include the gravitational constant G in all expres-
sions here, as the standard convention G = 1 used in Numerical
Relativity differs from the convention 16πG = 1 employed by
[43, 44].
in the near zone. Furthermore, outside the near zone
ij,(div)
∼ 1/r2, so that hTT
ij,(div)
falls off much faster than
rest of hTTij , which falls off like 1/r. Hence
hTTij = h
TT (NZ)
ij − δ
TT kl
ij (✷
ret −∆−1)s̄kl + hTTij,(div), (9)
where hTT
ij,(div)
can be neglected if we only keep terms up
to O(v/c)4 generally, and O(1/r) at infinity.
The full expression for hTTij for N interacting point
particles from Eq. (4.3) of [43] is:
hTTij = h
TT (NZ)
ij + h
ij,(div) + 16πG
d3~k dω dτ
(2 π)4
pAi pAj
B 6=A
nABi nABj
× (ω/k)
~k·(~x−~xA)−i ω (t−τ)
k2 − (ω + i ǫ)2
. (10)
The first term in (10), h
TT (NZ)
ij can be expanded in v/c
TT (NZ)
ij = h
TT (4)
ij + h
TT (5)
ij +O(v/c)
6. (11)
The leading order term at O(v/c)4, is given explicitly by
Eq. (A20) of [44]:
hTT (4)ij =
mA rA
‖ ~pA ‖2 −5 (n̂A · ~pA)2
δij + 2 piA p
3(n̂A · ~pA)2 − 5 ‖ ~pA ‖2
niA n
A + 12(n̂A · ~pA)n
B 6=A
niABn
AB + 2
rA + rB
niA n
rABrA
+ 3rA
rABrA
, (12)
where sAB ≡ rA + rB + rAB. The other two terms
in Eq. (10) can be shown to be small in the near zone
(r ≪ λ, where the characteristic wavelength λ ∼ 100M
for rAB ∼ 10M). However, hTT (NZ)ij is only a valid ap-
proximation to hTTij in the near zone, and becomes highly
inaccurate when used further afield.
Setting aside these far-field issues, Tichy et al. [45] ap-
plied Schäfer’s formulation, in the context of a black-hole
binary system, to construct initial data that are accurate
up to O(v/c)5 in the near zone. They noted that the
ADM-TT decomposition was well-adapted to the use of
a puncture approach to handle black-hole singularities.
This approach is essentially an extension of the method
introduced in [10]. It allows a simple numerical treatment
of the black holes without the need for excision.
The PN-based puncture data of Tichy et al. have
not been used for numerical evolutions. This is in part
because these data, just like standard puncture data [10,
49, 50, 51], do not contain realistic gravitational waves
in the far zone: h
TT (NZ)
ij does not even vaguely agree
with the 2PN approximation to the waveform amplitude
nor with the quadrupole approximation to the waveform
phase for realistic inspiral.
To illustrate this, we restrict to the case of two point
sources, and compute the “plus” and “cross” polariza-
tions of the near-zone approximation for hTTij :
+ = h
TT (NZ)
θ, (13)
× = h
TT (NZ)
φ. (14)
For comparison, the corresponding polarizations of the
quadrupole approximation for the gravitational-wave
strain are given by (paraphrasing Eq. (3.4) of [52]):
(1+cos2 θ)(πGMfGW)2/3cos(ΦGW), (15)
cos θ(πGMfGW)2/3sin(ΦGW), (16)
where M ≡ ν3/5M is the “chirp mass” of the binary,
given in terms of the total PN mass of the system M =
m1 +m2, and the symmetric mass ratio ν = m1m2/M
The angle θ is the “inclination angle of orbital angular
momentum to the line of sight toward the detector”; that
is, just the polar angle to the field point, when the binary
moves in the x-y plane. ΦGW and fGW are the phase and
frequency of the radiation at time t, exactly twice the
orbital phase Φ(t− r) and orbital frequency Ω(t− r)/2π.
The lowest-order PN prediction for radiation-reaction
effects yields a simple inspiral of the binary over time,
with orbital phasing given by
Φ(τ) = Φ(tc)−
Θ5/8, (17)
Ω(τ) =
Θ−3/8, (18)
where Θ ≡ ν (tc − τ)/5GM , M and ν are given below
(16), and tc is a nominal “coalescence time”. To evaluate
(13-14), we need the transverse momentum p correspond-
ing to the desired separation r12. The simplest expression
for this is the classical Keplerian relation, which we give
parameterized by Ω(τ):
r12 = G
1/3M(MΩ)−2/3, (19)
p = Mν(GMΩ)1/3. (20)
In Fig. 1 we compare the plus polarization of the two
waveforms (13) and (15) at a field point r = 100M ,
θ = π/4, φ = 0, for a binary in the x-y plane, with ini-
tial separation r12 = 10M . The orbital frequency of the
binary is related to the separation r12 and momenta p en-
tering (13) by (19-20). To this level of approximation, the
binary has a nominal PN coalescence time tc ≈ 780M .
As might have been anticipated, both phase and ampli-
tude of h
TT (4)
ij are wrong outside the near zone. This
means that the data constructed from h
TT (4)
ij have the
wrong wave content, but nevertheless these data are still
accurate up to order (v/c)3 in the far zone.
It is evident from the present-time dependence of (12)
that it cannot actually contain any of the past history
of an inspiralling binary. We would expect that a cor-
rect “wave-like” contribution should depend rather on
the retarded time of each contributing point source. It
seems evident that the correct behavior must, in fact, be
contained in the as-yet unevaluated parts of (10). The
requisite evaluation is what we undertake in the next sec-
tion.
III. COMPLETING THE EVALUATION OF hTTij
To move forward, we will simplify (10) and (12) to the
case of only two particles. Then (10) reduces to:
hTTij = h
TT (NZ)
ij + 16πG
p1 i p1 j
~k ·(~x−~x1) +
p2 i p2 j
~k ·(~x−~x2) − G
n12i n12j
~k ·(~x−~x1)
n21i n21j
~k ·(~x−~x2)
· (ω/k)
2 e−i ω (t−τ)
k2 − (ω + i ǫ)2
d3~k dω dτ
(2 π)4
+ hTTij,(div) (21)
TT (NZ)
ij +H
+HTT2ij
−HTT1ij
Gm1m2
2 r12
−HTT2ij
Gm1m2
2 r12
+hTTij,(div), (22)
where
HTTAij [~u] := 16πG
d3~k dω
(2 π)4
[ui uj ]
(ω/k)2
k2 − (ω + i ǫ)2
~k·(~x−~xA(τ)) e−i ω (t−τ). (23)
Here, the “TT projection” is effected using the operator
i := δ
i − ki kj/k2. For an arbitrary spatial vector ~u,
[ui uj ]
TT = uc ud (P
Pij P
= ui uj +
ki kj
u(i kj)
. (24)
Details on the evaluation of these terms are presented
in Appendix A. After calculation, we write the result as
a sum of terms evaluated at the present field-point time
t, the retarded time trA defined by
t− trA − rA(trA) = 0, (25)
and integrals between trA and t,
TTA[~u] = H
TTA[~u; t] +H
TTA[~u; t
TTA[~u; t
A → t], (26)
where the three parts are given by:
0 100 200 300 400 500 600 700 800
-0.0015
-0.001
-0.0005
0.0005
0.001
 (quadrupole)
 (NZ)
FIG. 1: Plus polarization of the quadrupole (black/solid) and
near-zone (red/dashed) strains observed at field point r =
100M , θ = π/4, φ = 0. The binary orbits in the x-y plane,
with initial separation r12 = 10M , and a nominal coalescence
time tc ≈ 780M . Both phase and amplitude of h
TT (4)
very wrong outside the near zone.
TTA[~u; t] = −
rA(t)
u2 − 5 (~u · n̂A)2
δi j + 2 ui uj +
3 (~u · n̂A)2 − 5 u2
niA n
+12 (~u · n̂A)u(i nj)A
, (27)
TTA[~u; t
−2 u2 + 2 (~u · n̂A)2
δi j + 4 ui uj +
2 u2 + 2 (~u · n̂A)2
niA n
−8 (~u · n̂A)u(i nj)A
, (28)
TTA[~u; t
A → t] = −G
(t− τ)
rA(τ)3
−5 u2 + 9 (~u · n̂A)2
δi j + 6 ui uj − 12 (~u · n̂A)u(i nj)A
9 u2 − 15 (~u · n̂A)2
niA n
(t− τ)3
rA(τ)5
u2 − 5 (~u · n̂A)2
δi j + 2 ui uj − 20 (~u · n̂A)u(i nj)A
−5 u2 + 35 (~u · n̂A)2
niA n
. (29)
In Fig. 2, we show the retarded times calculated for
each particle, as measured at points along the x axis, for
the same orbit as in Fig. 1. We also show the corre-
sponding retarded times for a binary in an exactly cir-
cular orbit. Since the small-scale oscillatory effect of the
finite orbital radius would be lost by the overall linear
trend, we have multiplied by the orbital radius.
A. Reconciling with Jaranowski & Schäfer’s h
TT (4)
From the derivation above it is clear that hTTij includes
retardation effects, so it will not depend solely on the
present time. We might even expect that all “present-
time” contributions should vanish individually, or should
cancel out. It can be seen easily from (27) that the “t”
part of the second and third terms of Eq. (22) exactly
cancel out the “kinetic” part (first line) of Eq. (12). Thus,
we can simply remove that line in Eq. (12), and use the
500 1000 1500
particle 1 (circular)
particle 1 (inspiral)
particle 2 (circular)
particle 2 (inspiral)
FIG. 2: Retarded times for particles 1 and 2, as measured by
observers along the x axis at the initial time t = 0, for the
binary of Fig. 1. To highlight the oscillatory effect of the
finite-radius orbit on tr, we first divide by the average field
distance r.
“tr” part instead. One may similarly inquire whether the
“t” parts of the fourth and fifth terms of Eq. (22) above,
TT (pot,now)
ij ≡ −H
Gm1m2
2 r12
n̂12; t
−HTT2ij
Gm1m2
2 r12
n̂12; t
, (30)
also cancel the remaining, “potential” parts of Eq. (12).
The answer is “not completely”; expanding in powers of
1/r, we find:
TT (pot,4)
ij + h
TT (pot,now)
G2m1m2 r12
16 r3
(3 + 14W 2 − 25W 4) δi j − 4 (1 + 5W 2)n12i n12j
−5 (1 + 6W 2 − 7W 4)n1i n1j + 2W (7 + 9W 2) (n12i n1j + n12j n1i)
+O(1/r4),(31)
where W ≡ sin θ cos(φ − Φ(t)), and Φ(t) is the orbital
phase of particle 1 at the present time t. That is, the
“new” contribution cancels the 1/r and 1/r2 pieces of
TT (4)
ij entirely. In the far zone the result is thus smaller
than the hTT
ij,(div)
term which we are ignoring everywhere,
since it is small both in the near and the far zone [43].
We note here two general properties of the contribu-
tions to the full hTTij .
1. In the near zone h
TT (4)
ij is the dominant term since
all other terms arise from (✷−1ret − ∆−1)skl. Thus
all other terms must cancel within the accuracy of
the near-zone approximation.
TT (4)
ij is wrong far from the sources; thus, the
new corrections should “cancel” h
TT (4)
ij entirely, far
from sources. Note, however, that while hij =
−✷−1retskl depends only on retarded time, its TT-
projection hTTij = δ
TT kl
ij hkl has a more complicated
causal structure; E.g. the finite time integral comes
from applying the TT-projection. [Proof: Even if
we had a source given exactly by s̄kl, h
TT (4)
ij would
depend only the present time, hij would depend
only on retarded time, and hTTij would (as we have
computed) contain a finite time integral term.]
Additionally, the full hTTij agrees well with quadrupole
predictions, which we demonstrate in Section IV.
IV. NUMERICAL RESULTS AND INVARIANTS
A. Phasing and Post-Keplerian Relations
It has been known for some time (see for example [53])
that gravitational wave phase plays an even more impor-
tant part in source identification than does wave ampli-
tude. In PN work, phase and amplitude are estimated
somewhat separately; the amplitude requires knowledge
of the time-dependent multipoles, used in developing the
the full metric, while the phase can be relatively simply
approximated from the orbital equations of motion, tak-
ing into account the gravitational wave flux at infinity to
evolve the orbital parameters [54].
The quadrupole waveform introduced for the compar-
ison in Fig. 1 had an amplitude accurate to O(v/c)4
and the simplest available time evolution for the phase.
Waveform phase is a direct consequence of orbital phase.
To lowest order, we could have assumed a binary mov-
ing in a circular orbit (of zero eccentricity) since, up to
2PN order, we can have circular orbits, where the linear
momentum, p, of each particle is related to the separa-
tion r12 by, say, Eq. (24) of [45]. Nevertheless, circular
orbits are physically unrealistic – since radiation reac-
tion will lead to inspiral and merger of the particles –
and Eqs. (17-18) already include leading-order radiation-
reaction effects. Moreover, the phase errors that would
accrue from using purely circular orbits would be larger,
the further from the sources we tried to compute them.
The calculations of section III lead to waveform am-
plitudes that are accurate at O(v/c)4 everywhere. How-
ever, we desire that our initial-data wave content already
encode the phase as accurately as possible. Highly ac-
curate phase for our initial data (via hTT), and hence in
the leading edge of the waveforms we would extract from
numerical evolution, is critical for parameter estimation
following a detection.
For demonstrative purposes, in this section, we will
restrict ourselves to the simplest phasing relations con-
sistent with radiation-reaction inspiral as given by Eqs.
(17-18), while using higher-order PN expressions than
Eqs. (19 -20) for relating the orbit to the phase. For ex-
ample, from [55], we have found to second PN (beyond
leading) order:
r12(Ω)
= (GMΩ)−2/3 − (3− ν)
− (18− 81ν − 8ν
(GMΩ)2/3, (32)
= (GMΩ)1/3 +
(15− ν)
(GMΩ)
(441− 324ν − ν2)
(GMΩ)5/3, (33)
and we note that higher-order equivalents of these can
be computed from [56].
In the numerical construction of initial data, the pri-
mary input is the coordinate separation of the holes. In
placing the punctures on the numerical grid, the separa-
tion must be maintained exactly. To ensure this, we in-
vert Eq. (32) to obtain the exact Ωr corresponding to our
desired r12. Then we use Eq. (18) with t = 0 to find the
coalescence time tc that yields this Ωr. Once we have ob-
tained tc, we then find the orbital phase Φ and frequency
Ω at any source time τ directly from Eqs. (17-18), and
the corresponding separation r12 and momentum p from
Eqs. (32-33), or their higher-order equivalents.
In Fig. 3, we show a representative component of the
retarded-time part of hTTij for both circular and leading-
order inspiral orbits. For both orbits, we use the ex-
tended Keplerian relations (32) and (33); otherwise the
orbital configuration is that of Fig. 1. The coalescence
time is now tc ∼ 1100M . We can see that the cumu-
lative wavelength error of the circular-orbit assumption
becomes very large at large distances from the sources.
This demonstrates that using inspiral orbits instead of
circular orbits will significantly enhance the phase accu-
racy of the initial data, even though circular orbits are
in principle sufficient when we include terms only up to
O(v/c)4 as done in this work. From now on we use only
inspiral orbits.
0 500 1000 1500
circular
inspiral
FIG. 3: The xx component of the full hTTij for a binary with
initial separation r12 = 10M in a circular (black/solid) or
inspiralling (red/dashed) orbit. Both fields have been rescaled
by the observer radius r = z to compensate for the leading
1/r fall-off. The orbital configuration is the same as for Fig.
1, apart from the Keplerian relations, where we have used the
higher-order relations (32-33), yielding tc ∼ 1100M . Note the
frequency broadening at more distant field points.
0 100 200 300 400 500 600 700 800
-0.0015
-0.001
-0.0005
0.0005
0.001
 (quadrupole)
 (full)
 (quadrupole)
 (full)
FIG. 4: Plus and cross polarizations of the strain observed at
field point r = 100M , θ = π/4, φ = 0. Both the quadrupole-
approximation waveform (black/solid and green/dot-dashed)
and the full (red/dashed and blue/dotted) waveforms coming
from hTTij are shown. The orbital configuration is the same as
for Fig. 1.
Next, we compare our full waveform hTTij (expressed
as the combinations h+ and h×) at an intermediate-field
position (r = 100M , θ = π/4, φ = 0) to the lowest-order
quadrupole result. In Fig. 4, the orbital configuration
is the same as for Fig. 1. As one can see, both the
+ and × polarizations of our hTTij agree very well with
quadrupole results, as they should. We demonstrate the
near- and intermediate-zone behavior of the new data on
50 100 150 200 250 300 350 400 450 500
-0.05
-0.05
quadrupole
FIG. 5: Plus and cross polarizations of the strain observed at
t = 0 along the z axis. We show the near-zone (solid/black),
the quadrupole (dashed/red) and full (dot-dashed/green)
waveforms. All waveforms have been rescaled by the observer
radius r = z to compensate for the leading 1/r fall-off. The
orbital configuration is the same as for Fig. 1.
the initial spatial slice in Fig. 5. The quadrupole and
full solutions agree very well outside ∼ 100M . However,
the full solution’s phase and amplitude approach the NZ
solution closer to the sources.
B. Numerical Implementation
After having confirmed that we have a PN three-metric
gij that is accurate up to errors of order O(v/c)
5, and
that correctly approaches the quadrupole limit outside
the near zone, we are now ready to construct initial data
for numerical evolutions. In order to do so, we need the
intrinsic curvature Kij , which can be computed as in
Tichy et al. [45] from the conjugate momentum. The
difference is that here we use the full ḣTTij instead of the
near-zone approximation ḣ
TT (4)
ij to obtain the conjugate
momentum [43]. The result is
Kij = −ψ−10PN
ḣTTij + (φ(2)π̃
+O(v/c)6, (34)
where the error term comes from neglecting terms like
ij,(div)
at O(v/c)5 in hTTij , and where ψPN , π̃
φ(2) can be found in Tichy et al. [45]. An additional
difference is that the time derivative of hTTij is evaluated
numerically in this work. Note that the results for gij
are accurate up to O(v/c)4, while the results for Kij are
0 0.5 1 1.5
 = 10M
 = 20M
 = 50M
 = 100M
0 0.5 1 1.5
FIG. 6: Upper panel: Hamiltonian constraint violation along
the y axis of our new data in the near zone, as a function of
binary separation r12. Lower panel: Momentum constraint
(y-component) violation of the same data along the x axis.
The orbital configuration is that of Fig. 3. Distances have
been scaled relative to r12, so that the punctures are initially
at y/r12 = ±0.5.
accurate up O(v/c)5, because Kij contains an additional
time derivative [45, 57, 58].
Next we show the violations of the Hamiltonian and
momentum constraints computed from gij and Kij , as
functions of the binary separation r12. As we can see in
both panels of Fig. 6, the constraints become smaller for
larger separations, because the post-Newtonian approxi-
mation gets better. Note that, as in [45], the constraint
violation remains finite everywhere, and is largest near
each black hole.
C. Curvature Invariants and Asymptotic Flatness
In analysis of both initial and evolved data, it is often
instructive to investigate the behavior of scalar curva-
ture invariants, as these give some idea of the far-field
properties of our solution. We expect, for an asymptoti-
cally flat space-time, that in the far field, the speciality
index S ≡ 27J 2/I3 will be close to unity. This can
be seen from the following arguments. Let us choose a
tetrad such that the Weyl tensor components ψ1 and ψ3
are both zero. Further, we assume that in the far field
ψ0 and ψ4 are both perturbations of order ǫ off a Kerr
background. Then
S ≈ 1− 3ψ0ψ4
+O(ǫ3), (35)
which is indeed close to one. Note however, that this
argument only works if the components of the Weyl ten-
sor obey the peeling theorem, such that ψ2 ∼ O(r−3),
ψ0 ∼ O(r−5) and ψ4 ∼ O(r−1). In particular, if ψ0 falls
off more slowly than O(r−5), S will grow for large r.
Now observe that ψ0 ∼ O(r−5) ∼ M3/r5 is formally of
O(v/c)6. Thus, in order to see the expected behavior of
S ≈ 1 in the far-field we need to go to O(v/c)6. If we only
go to O(v/c)4 (as done in this work) ψ0 consists of un-
controlled remainders only, which should in principle be
dropped. When we numerically compute S we find that
for our data, S deviates further and further from unity
for large distances from the binary. This reflects the fact
that the so-called “incoming” Weyl scalar ψ0 only falls
off as 1/r3, due to uncontrolled remainders at O(v/c)6,
which arise from a mixing of the background with the
TT waveform.
V. DISCUSSION AND FUTURE WORK
Exploring and validating PN inspiral waveforms is cru-
cially important for gravitational-wave detection and for
our theoretical understanding of black-hole binaries. Our
goal has been to provide a step forward in this under-
standing by building a direct interface between the PN
approach and numerical evolution, along the lines ini-
tially outlined in Ref. [45]. In this paper we have
essentially completed the calculation of the transverse-
traceless part of the ADM-TT metric to O(v/c)4 pro-
vided in [45], yielding data that, on the initial Cauchy
slice, will describe the space-time into the far-field.
We have incorporated this formulation into a numerical
initial-data routine adapted to the “puncture” topology
that has been so successful recently, and have explored
these data’s numerical properties on the initial slice.
Our next step is to evolve these data with moving
punctures, and investigate how the explicit incorporation
of post-Newtonian waveforms in the initial data affects
both the ensuing slow binary inspiral of the sources and
the release of radiation from the system. We note es-
pecially that our data are non-conformally flat beyond
O(v/c)3. We expect our data to incorporate smaller
unphysical initial distortions in the black holes than is
possible with conformal flatness, and hence less spurious
gravitational radiation during the numerical evolution.
We see this as a very positive step toward providing fur-
ther validation of numerical relativity results for multiple
orbit simulations, since it permits comparison with PN
results where they are expected to be reliable. Our initial
data will also allow us to fully evaluate the validity of PN
results for merging binaries by enabling comparison with
the most accurate numerical relativity results.
We expect that further development of these data will
certainly involve the use of more accurate orbital phas-
ing information than the leading order given by Eqs. (17-
18). This information is available in radiative coordinates
(see, e.g. Eq. (6.29) of [59]) appropriate for far-field eval-
uation of the gravitational radiative modes; it may be
possible to produce them in ADM-TT coordinates via
a contact transformation, or by direct calculation (see,
e.g. [60]). For initial separations similar to the fiducial
test case of this paper, r12 = 10M , the order necessary
for clean matching of the initial wave content with the
new radiation generated in evolution should not be par-
ticularly high [26]. As noted, the Keplerian relations
Eqs. (32-33) can easily be extended to higher PN order.
The data presented already allow for arbitrary initial
mass ratios ν; this introduces the possibility of significant
gravitational radiation in odd-l multipoles, together with
associated phenomena, such as in-plane recoil “kicks”.
An interesting future development of these data will be
the inclusion of spin angular momenta on the pre-merger
holes. This will open our initial-data prescription to de-
scribing an even richer spectrum of binary radiation.
Acknowledgments
We would like to thank L. Blanchet and G. Schäfer for
generous assistance and helpful discussion.
M.C., B.K. and B.W. gratefully acknowledge the sup-
port of the NASA Center for Gravitational Wave Astron-
omy (NAG5-13396). M.C. and B.K. also acknowledge the
NSF for financial support under grants PHY-0354867 and
PHY-0722315. B.K. also acknowledges support from the
NASA Postdoctoral Program at the Oak Ridge Associ-
ated Universities. The work of W.T. was supported by
NSF grant PHY-0555644. W.T. also acknowledges par-
tial support from the NCSA under Grant PHY-060040T.
The work of B.W. was also supported by NSF grants
PHY-0245024 and PHY-0555484.
APPENDIX A: DETAILS OF INTEGRAL
CALCULATION
Here we present some more details of the calculations
that lead to the three contributions to Eq. (23): Eqs.
(27-29). Inserting Eq. (24) in the general integral (23),
we can write H
TTA[~u] as a combination of scalar and
tensor terms:
HTTAij [~u] = 16πG
ui uj −
Iij A +
[uc ud
IcdA δij
2 uc u(i
Icj)A +
[uc ud
I cdij A
,(A1)
where the “I” integrals are defined as:
d3~k dω
(2 π)4
(ω/k)2 ei k rA cos θ−i ω T
k2 − (ω + i ǫ)2
, (A2)
d3~k dω
(2 π)4
ki kj
× (ω/k)
2 ei k rA cos θ−i ω T
k2 − (ω + i ǫ)2
, (A3)
i j c d
d3~k dω
(2 π)4
ki kj kc kd
× (ω/k)
2 ei k rA cos θ−i ω T
k2 − (ω + i ǫ)2
. (A4)
Here T ≡ t − τ , and ~rA ≡ ~x − ~xA. We have also taken
our integration coordinates such that ~rA lies in the z di-
rection, so that the dummy momentum vector ~k satisfies
~k · ~rA = k rA cos θ, (A5)
d3~k = k2 dk sin θ dθ dφ. (A6)
Define the unit orthogonal vectors n̂A ≡ (0, 0, 1) , ℓ̂ ≡
(cosφ, sinφ, 0). Then we can write
~k = k cos θ n̂A + k sin θ ℓ̂ ⇒ ~k · ~rA = rA ~k · n̂A.
We can also define a projector tensor onto ℓ̂:
Qa b ≡ δa b − na nb ⇒ Qab = δab − na nb
⇒ QacQcb = Qab , Qab nb = 0 , Qab ℓb = ℓa.
1. Angular integration
We will neglect the A subscript for now, until it be-
comes relevant again. To calculate the integrals (A2-A4),
we begin with the φ integration. The only φ dependence
comes from the ~ℓ parts of the ~k terms. It can be seen
from elementary trigonometric integrals that:
dφ ℓa =
dφ ℓa ℓb ℓc = 0,
dφ ℓa ℓb = π Qa b,
dφ ℓa ℓb ℓc ℓd =
Qa bQc d +Qa cQb d +Qa dQb c
We use these to calculate the φ integrals for Ia bA and
Ia b c dA . Define w ≡ cos θ. Then
dφ 1 = 2 π,
ka kb
= 2 π w2 na nb + π (1− w2)Qa b,
ka kb kc kd
= 2 π w4 na nb nc nd
+6 πw2 (1− w2)Q(a b nc nd)
(1− w2)2Q(a bQc d).
So the next integrals will differ in their θ dependence,
contained in the powers of w above. The θ integrals will
contain the following basic types:
g0(a) ≡
dw eaw = 2
sinh a
, (A7)
g2(a) ≡
dw w2 eaw = 2
sinh a
− 4 cosha
sinh a
g4(a) ≡
dw w4 eaw = 2
sinh a
− 8 cosha
sinh a
− 48 cosha
sinh a
. (A9)
Now Ia b and Ia b c d can be written as the linear combi-
nations:
Ia b =
Qa b I
(na nb − 1
Qa b)K
, (A10)
Ia b c d =
na nb nc nd − 3Q(a b nc nd) + 3
Q(a bQc d)
3Q(a b nc nd) −
Q(a bQc d)
Q(a bQc d) I
. (A11)
I here can be expressed in terms of g0(a) above:
(2 π)3
(ω/k)2
k2 − (ω + i ǫ)2 e
i k r cos θ−i ω T
(2 π)3
ω2 e−i ω T
k2 − (ω + i ǫ)2
g0(i k r)
(2 π)3
ω2 e−i ω T J0. (A12)
The 1/2 factor is because we moved to integrating k over
the whole real line instead of the positive half-line (this
is permissible as gn(a) is an even function of a). K and L
are defined analogously to I, but with extra even powers
of cos θ = w:
(2 π)3
(ω/k)2
k2 − (ω + i ǫ)2
ei k r cos θ−i ω T cos2 θ
(2 π)3
ω2 e−i ω T
k2 − (ω + i ǫ)2 g2(i k r)
(2 π)3
ω2 e−i ω T J2, (A13)
(2 π)3
(ω/k)2
k2 − (ω + i ǫ)2
ei k r cos θ−i ω T cos4 θ
(2 π)3
ω2 e−i ω T
k2 − (ω + i ǫ)2
g4(i k r)
(2 π)3
ω2 e−i ω T J4. (A14)
2. Momentum integration
Now we address the k integrals, defined as:
dk fn(k) =
dk f+n (k) +
dk f−n (k),
where we collect the positive exponents in the gn in the
integrand of f+n (k), and the negative exponents in f
n (k):
f+n (k) ≡
g+n (i k r)/2
k2 − (ω + i ǫ)2
, f−n (k) ≡
g−n (i k r)/2
k2 − (ω + i ǫ)2
We calculate this as the sum of contour integrals of the
“plus” and “minus” integrands (necessary, as the oppo-
site signs require different contours). Each of these has
poles at k = 0, k = k+ ≡ ω + i ǫ, and k = k− ≡ −ω − i ǫ
(the first of these is from the gn). We integrate the “plus”
integrands anticlockwise around the contour C1, and the
“minus” integrands anticlockwise around the contour C2
(see Fig. 7); taking the limit |k| → ∞, the contribu-
tion from the curved segments vanishes, and the residue
theorem gives us:
Jn = 2 π iRes[f
n , k+]− 2 π iRes[f−n , k−]
+π iRes[f+n , 0]− π iRes[f−n , 0]. (A15)
Calculating the residues, we find the values of each of
the Jn:
π ei r (ω+i ǫ)
r (ω + i ǫ)2
r (ω + i ǫ)2
, (A16)
π ei r (ω+i ǫ)
r (ω + i ǫ)2
π ei r(ω+i ǫ) [−2 + 2 i r (ω + i ǫ)]
r3 (ω + i ǫ)4
r3 (ω + i ǫ)4
, (A17)
π ei r (ω+i ǫ)
r (ω + i ǫ)2
4π ei r(ω+i ǫ)
r5 (ω + i ǫ)6
[6− 6 i r (ω + i ǫ)
−3 r2 (ω + i ǫ)2 + i r3 (ω + i ǫ)3
− 24 π
r5 (ω + i ǫ)6
. (A18)
3. Frequency integration
Now we perform the ω integration. Inserting the re-
sults (A16-A18) into (A12-A14) respectively, we see that
each of I, K and L contains a delta function, which we
can extract:
4 π r
[δ(T − r) − δ(T )],
4 π r
δ(T − r) + e−r ǫ
(2 π)3
e−iω (T−r) F2a(ω)
(2 π)3
e−i ω T F2b(ω),
4 π r
δ(T − r) + e−r ǫ
(2 π)3
e−iω (T−r) F4a(ω)
(2 π)3
e−i ω T F4b(ω),
where the new terms on the right-hand side come from
the Jn above, grouped by exponential, as that is what
determines the contours chosen during integration (see
Fig. 7):
F2a(ω) =
π ω2 [−2 + 2 i r (ω + i ǫ)]
r3 (ω + i ǫ)4
F2b(ω) =
2 π ω2
r3 (ω + i ǫ)4
F4a(ω) =
r5 (ω + i ǫ)6
[24− 24 i r (ω + i ǫ)
−12 r2 (ω + i ǫ)2 + 4 i r3 (ω + i ǫ)3
F4b(ω) = −
24 π ω2
r5 (ω + i ǫ)6
Now the residues are as follows (taking the ǫ→ 0 limit):
e−iω (T−r) F2a(ω),−i ǫ
2 π i T
e−iω T F2b(ω),−i ǫ
= −2 π i T
e−iω (T−r) F4a(ω),−i ǫ
4 π i T 3
e−iω T F4b(ω),−i ǫ
= −4 π i T
The only pole is at ω = −i ǫ, so if we can close the contour
in the upper half-plane, we will get zero.
• For T < 0, both the “a” and “b” integrals can be
closed in C1. Result: zero contribution.
• For 0 < T < r, the “a” integrals can be closed
in C1, but the “b” integrals must be closed in C2.
Result: “b” contribution.
• For T > r, both the “a” and “b” integrals must be
closed in C2. But then the “a” and “b” residues
cancel out. Result: zero contribution.
Thus the only interesting contribution happens in the
interval 0 < T < r ⇔ t− r(τ) < τ < t. In this case, the
final integrals yield
(2 π)3
e−iω (T−r) F2b(ω) = −
2 π r3
(2 π)3
e−iω (T−r) F4b(ω) = −
−ω − iǫ
ω + iǫ
FIG. 7: Contours needed to complete integration over k (left) and ω (right).
leading to the final result for K and L:
4 π r
δ(T − r)− 1
4 π r
δ(T ),
4 π r
δ(T − r)−Θ(T )Θ(r − T ) T
2 π r3
4 π r
δ(T − r)−Θ(T )Θ(r − T ) T
We use these to calculate the Ii j and Ii j k l:
Ii j =
ni nj
4 π r
δ(T − r)−Θ(T )Θ(r − T ) T
2 π r3
4 π r
δ(T ) + Θ(T )Θ(r − T ) T
2 π r3
,(A19)
Ii j k l =
ni nj nk nl
4 π r
δ(T − r) −Θ(T )Θ(r − T ) T
− 3Q(i j nk nl) Θ(T )Θ(r − T )
2 π r3
Q(i j Qk l)
4 π r
δ(T ) + Θ(T )Θ(r − T )
. (A20)
4. Time integration
The final integrations will be over the source time τ . The “crossing times” for the two Θ functions are τ = t and
τ = tr, where t is the present field time, and tr the corresponding retarded time defined by (25). Now taking a general
function y(τ), we find that
dτ IA y(τ) =
y(trA)
4 π rA(t
− y(t)
4 π rA(t)
A y(τ) =
niA n
4 π rA
4 π rA
3niA n
A − δ
) (t− τ) y(τ)
4 π rA(τ)3
i j k l
A y(τ) =
niA n
4 π rA
4 π rA
−3Q(i jA n
(t− τ)
2 π rA(τ)3
−niA n
A + 3Q
(t− τ)3
π rA(τ)5
y(τ).
These can now be substituted into the general integral (A1). We write the result as a sum of terms at the present
field-point time t, the retarded time trA, and interval terms between them,
TTA[~u] = H
TTA[~u; t] +H
TTA[~u; t
A] +H
TTA[~u; t
A → t],
TTA[~u; t] = −
rA(t)
ui uj − u
[uk ul
Qk lA δ
2 uk u
[uk ul
TTA[~u; t
ui uj − u
niA n
[uk ul
nkA n
2 uk u
[uk ul
niA n
TTA[~u; t
A → t] = −4G
(t− τ)
rA(τ)3
3niA n
A − δ
[uk ul
3nkA n
A − δk l
2 uk u
A − δj) k
[uk ul
(t− τ)3
rA(τ)5
[uk ul
niA n
A − 3Q
[1] R. Vogt, in Sixth Marcel Grossman Meeting on General
Relativity (Proceedings, Kyoto, Japan, 1991), edited by
H. Sato and T. Nakamura (World Scientific, Singapore,
1992), pp. 244–266.
[2] B. Abbott et al. (LIGO Scientific), Nucl. Instrum. Meth.
A517, 154 (2004), gr-qc/0308043.
[3] P. Bender et al., Tech. Rep. MPQ 233, Max-
Planck-Institut für Quantenoptik (1998), URL:
http://www.lisa-science.org/resources/talks-articles/mission/prephasea.pdf.
[4] K. Danzmann and A. Rudiger, Class. Quantum Grav.
20, S1 (2003).
[5] A. Buonanno, G. B. Cook, and F. Pretorius, Phys. Rev.
D 75, 124018 (2007), gr-qc/0610122.
[6] E. Berti, V. Cardoso, J. A. Gonzalez, U. Sperhake,
M. Hannam, S. Husa, and B. Brügmann (2007), arXiv:gr-
qc/0703053.
[7] F. Pretorius, Phys. Rev. Lett. 95, 121101 (2005), gr-
qc/0507014.
[8] M. Campanelli, C. O. Lousto, P. Marronetti, and
Y. Zlochower, Phys. Rev. Lett. 96, 111101 (2006), gr-
qc/0511048.
[9] J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz, and
J. van Meter, Phys. Rev. Lett. 96, 111102 (2006), gr-
qc/0511103.
[10] S. Brandt and B. Brügmann, Phys. Rev. Lett. 78, 3606
(1997), gr-qc/9703066.
[11] M. Shibata and T. Nakamura, Phys. Rev. D 52, 5428
(1995).
[12] T. Baumgarte and S. Shapiro, Phys. Rev. D 59, 024007
(1999), gr-qc/9810065.
[13] B. Brügmann, W. Tichy, and N. Jansen, Phys. Rev. Lett.
92, 211101 (2004), gr-qc/0312112.
[14] M. Campanelli, C. O. Lousto, and Y. Zlochower, Phys.
Rev. D 73, 061501(R) (2006), gr-qc/0601091.
[15] F. Pretorius, Class. Quantum Grav. 23, S529 (2006), gr-
qc/0602115.
[16] J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz, and
J. van Meter, Phys. Rev. D 73, 104002 (2006), gr-
qc/0602026.
[17] B. Brügmann et al. (2006), gr-qc/0610128.
[18] M. A. Scheel et al., Phys. Rev. D 74, 104006 (2006),
gr-qc/0607056.
[19] P. Marronetti, W. Tichy, B. Brügmann, J. Gonzalez,
M. Hannam, S. Husa, and U. Sperhake, Class. Quant.
Grav. 24, S43 (2007), gr-qc/0701123.
[20] W. Tichy, Phys. Rev. D 74, 084005 (2006), gr-
qc/0609087.
[21] H. P. Pfeiffer, D. A. Brown, L. E. Kidder, L. Lindblom,
G. Lovelace, and M. A. Scheel, Class. Quant. Grav. pp.
S59–S81 (2007), gr-qc/0702106.
[22] J. G. Baker, M. Campanelli, F. Pretorius, and Y. Zlo-
chower, Class. Quant. Grav. 24, S25 (2007), gr-
qc/0701016.
[23] J. Thornburg, P. Diener, D. Pollney, L. Rezzolla,
http://www.lisa-science.org/resources/talks-articles/mission/prephasea.pdf
E. Schnetter, E. Seidel, and R. Takahashi, Class. Quant.
Grav. 24, 3911 (2007), gr-qc/0701038.
[24] NRwaves home page:
https://gravity.psu.edu/wiki NRwaves.
[25] J. G. Baker, S. T. McWilliams, J. R. van Meter, J. Cen-
trella, D.-I. Choi, B. J. Kelly, and M. Koppitz, Phys.
Rev. D 75, 124024 (2007), gr-qc/0612117.
[26] J. G. Baker, J. R. van Meter, S. T. McWilliams, J. Cen-
trella, and B. J. Kelly (2006), gr-qc/0612024.
[27] M. Campanelli, Class. Quant. Grav. 22, S387 (2005),
astro-ph/0411744.
[28] F. Herrmann, D. Shoemaker, and P. Laguna (2006), gr-
qc/0601026.
[29] J. G. Baker et al., Astrophys. J. 653, L93 (2006), astro-
ph/0603204.
[30] M. Campanelli, C. O. Lousto, and Y. Zlochower, Phys.
Rev. D 74, 041501(R) (2006), gr-qc/0604012.
[31] M. Campanelli, C. O. Lousto, and Y. Zlochower, Phys.
Rev. D 74, 084023 (2006), astro-ph/0608275.
[32] M. Campanelli, C. O. Lousto, Y. Zlochower, B. Krishnan,
and D. Merritt, Phys. Rev. D 75, 064030 (2007), gr-
qc/0612076.
[33] J. A. Gonzalez, U. Sperhake, B. Bruegmann, M. Han-
nam, and S. Husa, Phys. Rev. Lett. 98, 091101 (2007),
gr-qc/0610154.
[34] F. Herrmann, I. Hinder, D. Shoemaker, P. Laguna, and
R. A. Matzner (2007), gr-qc/0701143.
[35] M. Campanelli, C. O. Lousto, Y. Zlochower, and D. Mer-
ritt, 659, L5 (2007), revised version has very different
numbers/formulae, gr-qc/0701164.
[36] M. Koppitz, D. Pollney, C. Reisswig, L. Rezzolla,
J. thornburg, P. Diener, and E. Schnetter, Phys. Rev.
Lett. 99, 041102 (2007), gr-qc/0701163.
[37] J. A. Gonzalez, M. D. Hannam, U. Sperhake,
B. Brügmann, and S. Husa, Phys. Rev. Lett. 98, 231101
(2007), gr-qc/0702052.
[38] D.-I. Choi et al. (2007), gr-qc/0702016.
[39] J. G. Baker et al. (2007), astro-ph/0702390.
[40] F. Pretorius and D. Khurana, Class. Quant. Grav. 24,
S83 (2007), gr-qc/0702084.
[41] M. Campanelli, C. O. Lousto, Y. Zlochower, and D. Mer-
ritt, Phys. Rev. Lett. 98, 231102 (2007), arXiv:gr-
qc/0702133.
[42] W. Tichy and P. Marronetti (2007), gr-qc/0703075.
[43] G. Schäfer, Ann. Phys. 161, 81 (1985).
[44] P. Jaranowski and G. Schäfer, Phys. Rev. D 57, 7274
(1998), errata: Phys. Rev. D 63, 029902(E) (2000), gr-
qc/9712075.
[45] W. Tichy, B. Brügmann, M. Campanelli, and P. Diener,
Phys. Rev. D 67, 064008 (2003), gr-qc/0207011.
[46] Cactus Compuational Toolkit,
http://www.cactuscode.org.
[47] T. Ohta, H. Okamura, T. Kimura, and K. Hiida, Prog.
Theor. Phys. 51, 1598 (1974).
[48] V. A. Fock, The Theory of Space, Time and Gravitation,
2nd ed. (Pergamon Press, 1964).
[49] W. Tichy, B. Brügmann, and P. Laguna, Phys. Rev. D
68, 064008 (2003), gr-qc/0306020.
[50] W. Tichy and B. Brügmann, Phys. Rev. D 69, 024006
(2004), gr-qc/0307027.
[51] M. Ansorg, B. Brügmann, and W. Tichy, Phys. Rev. D
70, 064011 (2004), gr-qc/0404056.
[52] L. S. Finn and D. F. Chernoff, Phys. Rev. D 47, 2198
(1993), gr-qc/9301003.
[53] C. Cutler et al., Phys. Rev. Lett. 70, 2984 (1993), astro-
ph/9208005.
[54] W. Tichy, E. E. Flanagan, and E. Poisson, Phys. Rev. D
61, 104015 (2000), gr-qc/9912075.
[55] G. Schäfer and N. Wex, Phys. Lett. A 174, 196 (1993).
[56] R.-M. Memmesheimer, A. Gopakumar, and G. Schäfer,
Phys. Rev. D 70, 104011 (2004), gr-qc/0407049.
[57] N. Yunes, W. Tichy, B. J. Owen, and B. Brügmann,
Phys. Rev. D 74, 104011 (2006), gr-qc/0503011.
[58] N. Yunes and W. Tichy, Phys. Rev. D 74, 064013 (2006),
gr-qc/0601046.
[59] L. Blanchet, Phys. Rev. D 54, 1417 (1996), gr-
qc/9603048.
[60] T. Damour, A. Gopakumar, and B. R. Iyer, Phys. Rev.
D 70, 064028 (2005), gr-qc/0404128.
http://www.cactuscode.org

  We present improved post-Newtonian-inspired initial data for non-spinning
black-hole binaries, suitable for numerical evolution with punctures. We
revisit the work of Tichy et al. [W. Tichy, B. Bruegmann, M. Campanelli, and P.
Diener, Phys. Rev. D 67, 064008 (2003)], explicitly calculating the remaining
integral terms. These terms improve accuracy in the far zone and, for the first
time, include realistic gravitational waves in the initial data. We investigate
the behavior of these data both at the center of mass and in the far zone,
demonstrating agreement of the transverse-traceless parts of the new metric
with quadrupole-approximation waveforms. These data can be used for numerical
evolutions, enabling a direct connection between the merger waveforms and the
post-Newtonian inspiral waveforms.

Black hole puncture initial data with realistic gravitational wave content
B. J. Kelly,1, 2 W. Tichy,3 M. Campanelli,4, 2 and B. F. Whiting5, 2
1Gravitational Astrophysics Laboratory, NASA Goddard Space Flight Center, 8800 Greenbelt Rd., Greenbelt, MD 20771, USA
2Center for Gravitational Wave Astronomy, Department of Physics and Astronomy,
The University of Texas at Brownsville, Brownsville, Texas 78520
3Department of Physics, Florida Atlantic University, Boca Raton Florida 33431-0991
4Center for Computational Relativity and Gravitation,
School of Mathematical Sciences, Rochester Institute of Technology,
78 Lomb Memorial Drive, Rochester, New York 14623
5Department of Physics, University of Florida, Gainsville Florida 32611-8440
(Dated: October 26, 2018)
We present improved post-Newtonian-inspired initial data for non-spinning black-hole binaries,
suitable for numerical evolution with punctures. We revisit the work of Tichy et al. [W. Tichy, B.
Brügmann, M. Campanelli, and P. Diener, Phys. Rev. D 67, 064008 (2003)], explicitly calculating
the remaining integral terms. These terms improve accuracy in the far zone and, for the first
time, include realistic gravitational waves in the initial data. We investigate the behavior of these
data both at the center of mass and in the far zone, demonstrating agreement of the transverse-
traceless parts of the new metric with quadrupole-approximation waveforms. These data can be
used for numerical evolutions, enabling a direct connection between the merger waveforms and the
post-Newtonian inspiral waveforms.
PACS numbers: 04.25.Dm, 04.25.Nx, 04.30.Db, 04.70.Bw
I. INTRODUCTION
Post-Newtonian (PN) methods have played a funda-
mental role in our understanding of the astrophysical im-
plications of Einstein’s theory of general relativity. Most
importantly, they have been used to confirm that the ra-
diation of gravitational waves accounts for energy loss in
known binary pulsar configurations. They have also been
used to create templates for the gravitational waves emit-
ted from compact binaries which might be detected by
ground-based gravitational wave observatories, such as
LIGO [1, 2], and the NASA/ESA planned space-based
mission, LISA [3, 4]. However, PN methods have not
been extensively used to provide initial data for binary
evolution in numerical relativity, nor, until recently (see
[5, 6]), have they been extensively studied so that their
limitations could be well identified and the results of nu-
merical relativity independently confirmed.
Until the end of 2004, the field of numerical relativ-
ity had been struggling to compute even a single or-
bit for a black-hole binary (BHB). Although debate oc-
curred on the advantages of one type of initial data over
another, the primary focus within the numerical rela-
tivity community was on code refinement which would
lead to more stable evolution. Astrophysical realism was
very much a secondary issue. However, this situation
has radically changed in the last few years with the in-
troduction of two essentially independent, but equally
successful techniques: the generalized harmonic gauge
(GHG) method developed by Pretorius [7] and the “mov-
ing puncture” approach, independently developed by the
UTB and NASA Goddard groups [8, 9]. Originally in-
troduced by Brandt & Brügmann [10] in the context of
initial data, the puncture method explicitly factored out
the singular part of the metric. When used in numerical
evolution in which the punctures remained fixed on the
numerical grid, it resulted in distortions of the coordi-
nate system and instabilities in the Baumgarte-Shapiro-
Shibata-Nakamura (BSSN) [11, 12] evolution scheme.
The revolutionary idea behind the moving puncture ap-
proach was precisely, not to factor out the singular part
of the metric, but rather evolve it together with the reg-
ular part, allowing the punctures to move freely across
the grid with a suitable choice of the gauge.
A golden age for numerical relativity is now emerging,
in which multiple groups are using different computer
codes to evolve BHBs for several orbits before plunge and
merger [13, 14, 15, 16, 17, 18, 19, 20, 21]. Comparison of
the numerical results obtained from these various codes
has taken place [22, 23, 24], and comparison with PN
inspiral waveforms has also been carried out with encour-
aging success [5, 6, 25, 26]. The application of successful
numerical relativity tools to study some important as-
trophysical properties (e.g. precession, recoil, spin-orbit
coupling, elliptical orbits, etc) of spinning and/or un-
equal mass-black hole systems is currently producing ex-
tremely interesting new results [27, 28, 29, 30, 31, 32,
33, 34, 35, 36, 37, 38, 39, 40, 41, 42]. It now seems that
the primary obstacle to further progress is simply one of
computing power. In this new situation, it is perhaps
time to return to the question of what initial data will
best describe an astrophysical BHB.
To date, the best-motivated description of pre-merger
BHBs has been supplied by PN methods. We might ex-
pect, then, that a PN-based approach would give us the
most astrophysically correct initial data from which to
run full numerical simulations. In practice, PN results
are frequently obtained in a form ill-adapted to numeri-
http://arxiv.org/abs/0704.0628v2
cal evolution. PN analysis often deals with the full four-
metric, in harmonic coordinates; numerical evolutions
frequently use ADM-type coordinates, with a canonical
decomposition of the four-metric into a spatial metric
and extrinsic curvature.
Fortunately, many PN results have been translated
into the language of ADM by Ohta, Damour, Schäfer and
collaborators. Explicit results for 2.5PN BHB data in the
near zone were given by Schäfer [43] and Jaranowski &
Schäfer (JS) [44], and these were implemented numer-
ically by Tichy et al. [45]. Their insight was that the
ADM-transverse-traceless (TT) gauge used by Schäfer
was well-adapted to a puncture approach. To facilitate
comparison with this earlier work [45], we continue to
use the results of Schäfer and co-workers, anticipating
that higher-order PN results should eventually become
available in a useful form.
The initial data provided previously by Tichy et
al. already include PN information. They are accurate
up to order (v/c)5 in the near zone (r ≪ λ), but the
accuracy drops to order (v/c)3 in the far zone (r ≫ λ)
[here λ ∼ π
r312/G(m1 +m2) is the gravitational wave-
length]. These data were incomplete in the sense that
they did not include the correct TT radiative piece in
the metric, and thus did not contain realistic gravita-
tional waves.
In this paper, we revisit the PN data problem in
ADM-TT coordinates, with the aim of supplying Numer-
ical Relativity with initial BHB data that extend as far
as necessary, and contain realistic gravitational waves.
To do this, we have evaluated the “missing pieces” of
Schäfer’s TT metric for the case of two non-spinning par-
ticles. We have analyzed the near- and far-zone behavior
of these data, and incorporated them numerically in the
Cactus [46] framework. In principle, the most accurate
PN metric available could be used at this step, but it is
not currently available in ADM-TT form.
The remainder of this paper is laid out as follows. In
Section II, we summarize the results of Schäfer (1985)
[43], and Jaranowski & Schäfer (1997) [44] and their ap-
plication by Tichy et al. (2003) [45], to the production
of puncture data for numerical evolution. In Section
III, we describe briefly the additional terms necessary
to complete hTT to order (v/c)4, deferring details to the
Appendix. In Section IV, we study the full data both
analytically and numerically. Section V summarizes our
results, and lays the groundwork for numerical evolution
of these data, to be presented in a subsequent article.
II. ADM-TT GAUGE IN POST-NEWTONIAN
The “ADM-TT” gauge [43, 47] is a 3+1 split of data
where the three-metric differs from conformal flatness
precisely by a TT radiative part:
gij =
ηij + h
ij , (1)
πii = 0. (2)
The fields φ, πij and hTTij can all be expanded in a post-
Newtonian series. Solving the constraint equations of
3+1 general relativity in this gauge, [43, 44] obtained ex-
plicit expressions valid up to O(v/c)5 in the near zone,
incorporating an arbitrary number of spinless point par-
ticles, with arbitrary masses mA. For N particles, the
lowest-order contribution to the conformal factor is1:
φ(2) = 4G
, (3)
where rA =
~x− ~xA is the distance from the field point
to the location of particle A.
In principle hTTij is computed from
hTTij = −δTT klij ✷−1retskl, (4)
where ✷−1ret is the (flat space) inverse d’Alembertian (with
a “no-incoming-radiation” condition [48]), skl is a non-
local source term and δTT klij is the TT-projection oper-
ator. In order to compute hTTij we first rewrite Eq. (4)
hTTij = −δTT klij
∆−1 + (✷−1ret −∆−1)
TT (NZ)
ij − δ
TT kl
ij (✷
ret −∆−1)skl. (5)
Note that the near-zone approximation h
TT (NZ)
ij of h
has already been computed in [43] up to order O(v/c)4
(see also Eq. 12 below). The last term in Eq. (5) is diffi-
cult to compute because
skl = 16πG
pAk pAl
δ(x− xA) +
,l (6)
is a non-local source. However, we can approximate skl
s̄kl =
pAk pAl
B 6=A
nABk nABl
×16πGδ(x− xA). (7)
and show that
hTTij,(div) = −δ
TT kl
ij (✷
ret−∆−1)(skl− s̄kl) ∼ O(v/c)5 (8)
1 We explicitly include the gravitational constant G in all expres-
sions here, as the standard convention G = 1 used in Numerical
Relativity differs from the convention 16πG = 1 employed by
[43, 44].
in the near zone. Furthermore, outside the near zone
ij,(div)
∼ 1/r2, so that hTT
ij,(div)
falls off much faster than
rest of hTTij , which falls off like 1/r. Hence
hTTij = h
TT (NZ)
ij − δ
TT kl
ij (✷
ret −∆−1)s̄kl + hTTij,(div), (9)
where hTT
ij,(div)
can be neglected if we only keep terms up
to O(v/c)4 generally, and O(1/r) at infinity.
The full expression for hTTij for N interacting point
particles from Eq. (4.3) of [43] is:
hTTij = h
TT (NZ)
ij + h
ij,(div) + 16πG
d3~k dω dτ
(2 π)4
pAi pAj
B 6=A
nABi nABj
× (ω/k)
~k·(~x−~xA)−i ω (t−τ)
k2 − (ω + i ǫ)2
. (10)
The first term in (10), h
TT (NZ)
ij can be expanded in v/c
TT (NZ)
ij = h
TT (4)
ij + h
TT (5)
ij +O(v/c)
6. (11)
The leading order term at O(v/c)4, is given explicitly by
Eq. (A20) of [44]:
hTT (4)ij =
mA rA
‖ ~pA ‖2 −5 (n̂A · ~pA)2
δij + 2 piA p
3(n̂A · ~pA)2 − 5 ‖ ~pA ‖2
niA n
A + 12(n̂A · ~pA)n
B 6=A
niABn
AB + 2
rA + rB
niA n
rABrA
+ 3rA
rABrA
, (12)
where sAB ≡ rA + rB + rAB. The other two terms
in Eq. (10) can be shown to be small in the near zone
(r ≪ λ, where the characteristic wavelength λ ∼ 100M
for rAB ∼ 10M). However, hTT (NZ)ij is only a valid ap-
proximation to hTTij in the near zone, and becomes highly
inaccurate when used further afield.
Setting aside these far-field issues, Tichy et al. [45] ap-
plied Schäfer’s formulation, in the context of a black-hole
binary system, to construct initial data that are accurate
up to O(v/c)5 in the near zone. They noted that the
ADM-TT decomposition was well-adapted to the use of
a puncture approach to handle black-hole singularities.
This approach is essentially an extension of the method
introduced in [10]. It allows a simple numerical treatment
of the black holes without the need for excision.
The PN-based puncture data of Tichy et al. have
not been used for numerical evolutions. This is in part
because these data, just like standard puncture data [10,
49, 50, 51], do not contain realistic gravitational waves
in the far zone: h
TT (NZ)
ij does not even vaguely agree
with the 2PN approximation to the waveform amplitude
nor with the quadrupole approximation to the waveform
phase for realistic inspiral.
To illustrate this, we restrict to the case of two point
sources, and compute the “plus” and “cross” polariza-
tions of the near-zone approximation for hTTij :
+ = h
TT (NZ)
θ, (13)
× = h
TT (NZ)
φ. (14)
For comparison, the corresponding polarizations of the
quadrupole approximation for the gravitational-wave
strain are given by (paraphrasing Eq. (3.4) of [52]):
(1+cos2 θ)(πGMfGW)2/3cos(ΦGW), (15)
cos θ(πGMfGW)2/3sin(ΦGW), (16)
where M ≡ ν3/5M is the “chirp mass” of the binary,
given in terms of the total PN mass of the system M =
m1 +m2, and the symmetric mass ratio ν = m1m2/M
The angle θ is the “inclination angle of orbital angular
momentum to the line of sight toward the detector”; that
is, just the polar angle to the field point, when the binary
moves in the x-y plane. ΦGW and fGW are the phase and
frequency of the radiation at time t, exactly twice the
orbital phase Φ(t− r) and orbital frequency Ω(t− r)/2π.
The lowest-order PN prediction for radiation-reaction
effects yields a simple inspiral of the binary over time,
with orbital phasing given by
Φ(τ) = Φ(tc)−
Θ5/8, (17)
Ω(τ) =
Θ−3/8, (18)
where Θ ≡ ν (tc − τ)/5GM , M and ν are given below
(16), and tc is a nominal “coalescence time”. To evaluate
(13-14), we need the transverse momentum p correspond-
ing to the desired separation r12. The simplest expression
for this is the classical Keplerian relation, which we give
parameterized by Ω(τ):
r12 = G
1/3M(MΩ)−2/3, (19)
p = Mν(GMΩ)1/3. (20)
In Fig. 1 we compare the plus polarization of the two
waveforms (13) and (15) at a field point r = 100M ,
θ = π/4, φ = 0, for a binary in the x-y plane, with ini-
tial separation r12 = 10M . The orbital frequency of the
binary is related to the separation r12 and momenta p en-
tering (13) by (19-20). To this level of approximation, the
binary has a nominal PN coalescence time tc ≈ 780M .
As might have been anticipated, both phase and ampli-
tude of h
TT (4)
ij are wrong outside the near zone. This
means that the data constructed from h
TT (4)
ij have the
wrong wave content, but nevertheless these data are still
accurate up to order (v/c)3 in the far zone.
It is evident from the present-time dependence of (12)
that it cannot actually contain any of the past history
of an inspiralling binary. We would expect that a cor-
rect “wave-like” contribution should depend rather on
the retarded time of each contributing point source. It
seems evident that the correct behavior must, in fact, be
contained in the as-yet unevaluated parts of (10). The
requisite evaluation is what we undertake in the next sec-
tion.
III. COMPLETING THE EVALUATION OF hTTij
To move forward, we will simplify (10) and (12) to the
case of only two particles. Then (10) reduces to:
hTTij = h
TT (NZ)
ij + 16πG
p1 i p1 j
~k ·(~x−~x1) +
p2 i p2 j
~k ·(~x−~x2) − G
n12i n12j
~k ·(~x−~x1)
n21i n21j
~k ·(~x−~x2)
· (ω/k)
2 e−i ω (t−τ)
k2 − (ω + i ǫ)2
d3~k dω dτ
(2 π)4
+ hTTij,(div) (21)
TT (NZ)
ij +H
+HTT2ij
−HTT1ij
Gm1m2
2 r12
−HTT2ij
Gm1m2
2 r12
+hTTij,(div), (22)
where
HTTAij [~u] := 16πG
d3~k dω
(2 π)4
[ui uj ]
(ω/k)2
k2 − (ω + i ǫ)2
~k·(~x−~xA(τ)) e−i ω (t−τ). (23)
Here, the “TT projection” is effected using the operator
i := δ
i − ki kj/k2. For an arbitrary spatial vector ~u,
[ui uj ]
TT = uc ud (P
Pij P
= ui uj +
ki kj
u(i kj)
. (24)
Details on the evaluation of these terms are presented
in Appendix A. After calculation, we write the result as
a sum of terms evaluated at the present field-point time
t, the retarded time trA defined by
t− trA − rA(trA) = 0, (25)
and integrals between trA and t,
TTA[~u] = H
TTA[~u; t] +H
TTA[~u; t
TTA[~u; t
A → t], (26)
where the three parts are given by:
0 100 200 300 400 500 600 700 800
-0.0015
-0.001
-0.0005
0.0005
0.001
 (quadrupole)
 (NZ)
FIG. 1: Plus polarization of the quadrupole (black/solid) and
near-zone (red/dashed) strains observed at field point r =
100M , θ = π/4, φ = 0. The binary orbits in the x-y plane,
with initial separation r12 = 10M , and a nominal coalescence
time tc ≈ 780M . Both phase and amplitude of h
TT (4)
very wrong outside the near zone.
TTA[~u; t] = −
rA(t)
u2 − 5 (~u · n̂A)2
δi j + 2 ui uj +
3 (~u · n̂A)2 − 5 u2
niA n
+12 (~u · n̂A)u(i nj)A
, (27)
TTA[~u; t
−2 u2 + 2 (~u · n̂A)2
δi j + 4 ui uj +
2 u2 + 2 (~u · n̂A)2
niA n
−8 (~u · n̂A)u(i nj)A
, (28)
TTA[~u; t
A → t] = −G
(t− τ)
rA(τ)3
−5 u2 + 9 (~u · n̂A)2
δi j + 6 ui uj − 12 (~u · n̂A)u(i nj)A
9 u2 − 15 (~u · n̂A)2
niA n
(t− τ)3
rA(τ)5
u2 − 5 (~u · n̂A)2
δi j + 2 ui uj − 20 (~u · n̂A)u(i nj)A
−5 u2 + 35 (~u · n̂A)2
niA n
. (29)
In Fig. 2, we show the retarded times calculated for
each particle, as measured at points along the x axis, for
the same orbit as in Fig. 1. We also show the corre-
sponding retarded times for a binary in an exactly cir-
cular orbit. Since the small-scale oscillatory effect of the
finite orbital radius would be lost by the overall linear
trend, we have multiplied by the orbital radius.
A. Reconciling with Jaranowski & Schäfer’s h
TT (4)
From the derivation above it is clear that hTTij includes
retardation effects, so it will not depend solely on the
present time. We might even expect that all “present-
time” contributions should vanish individually, or should
cancel out. It can be seen easily from (27) that the “t”
part of the second and third terms of Eq. (22) exactly
cancel out the “kinetic” part (first line) of Eq. (12). Thus,
we can simply remove that line in Eq. (12), and use the
500 1000 1500
particle 1 (circular)
particle 1 (inspiral)
particle 2 (circular)
particle 2 (inspiral)
FIG. 2: Retarded times for particles 1 and 2, as measured by
observers along the x axis at the initial time t = 0, for the
binary of Fig. 1. To highlight the oscillatory effect of the
finite-radius orbit on tr, we first divide by the average field
distance r.
“tr” part instead. One may similarly inquire whether the
“t” parts of the fourth and fifth terms of Eq. (22) above,
TT (pot,now)
ij ≡ −H
Gm1m2
2 r12
n̂12; t
−HTT2ij
Gm1m2
2 r12
n̂12; t
, (30)
also cancel the remaining, “potential” parts of Eq. (12).
The answer is “not completely”; expanding in powers of
1/r, we find:
TT (pot,4)
ij + h
TT (pot,now)
G2m1m2 r12
16 r3
(3 + 14W 2 − 25W 4) δi j − 4 (1 + 5W 2)n12i n12j
−5 (1 + 6W 2 − 7W 4)n1i n1j + 2W (7 + 9W 2) (n12i n1j + n12j n1i)
+O(1/r4),(31)
where W ≡ sin θ cos(φ − Φ(t)), and Φ(t) is the orbital
phase of particle 1 at the present time t. That is, the
“new” contribution cancels the 1/r and 1/r2 pieces of
TT (4)
ij entirely. In the far zone the result is thus smaller
than the hTT
ij,(div)
term which we are ignoring everywhere,
since it is small both in the near and the far zone [43].
We note here two general properties of the contribu-
tions to the full hTTij .
1. In the near zone h
TT (4)
ij is the dominant term since
all other terms arise from (✷−1ret − ∆−1)skl. Thus
all other terms must cancel within the accuracy of
the near-zone approximation.
TT (4)
ij is wrong far from the sources; thus, the
new corrections should “cancel” h
TT (4)
ij entirely, far
from sources. Note, however, that while hij =
−✷−1retskl depends only on retarded time, its TT-
projection hTTij = δ
TT kl
ij hkl has a more complicated
causal structure; E.g. the finite time integral comes
from applying the TT-projection. [Proof: Even if
we had a source given exactly by s̄kl, h
TT (4)
ij would
depend only the present time, hij would depend
only on retarded time, and hTTij would (as we have
computed) contain a finite time integral term.]
Additionally, the full hTTij agrees well with quadrupole
predictions, which we demonstrate in Section IV.
IV. NUMERICAL RESULTS AND INVARIANTS
A. Phasing and Post-Keplerian Relations
It has been known for some time (see for example [53])
that gravitational wave phase plays an even more impor-
tant part in source identification than does wave ampli-
tude. In PN work, phase and amplitude are estimated
somewhat separately; the amplitude requires knowledge
of the time-dependent multipoles, used in developing the
the full metric, while the phase can be relatively simply
approximated from the orbital equations of motion, tak-
ing into account the gravitational wave flux at infinity to
evolve the orbital parameters [54].
The quadrupole waveform introduced for the compar-
ison in Fig. 1 had an amplitude accurate to O(v/c)4
and the simplest available time evolution for the phase.
Waveform phase is a direct consequence of orbital phase.
To lowest order, we could have assumed a binary mov-
ing in a circular orbit (of zero eccentricity) since, up to
2PN order, we can have circular orbits, where the linear
momentum, p, of each particle is related to the separa-
tion r12 by, say, Eq. (24) of [45]. Nevertheless, circular
orbits are physically unrealistic – since radiation reac-
tion will lead to inspiral and merger of the particles –
and Eqs. (17-18) already include leading-order radiation-
reaction effects. Moreover, the phase errors that would
accrue from using purely circular orbits would be larger,
the further from the sources we tried to compute them.
The calculations of section III lead to waveform am-
plitudes that are accurate at O(v/c)4 everywhere. How-
ever, we desire that our initial-data wave content already
encode the phase as accurately as possible. Highly ac-
curate phase for our initial data (via hTT), and hence in
the leading edge of the waveforms we would extract from
numerical evolution, is critical for parameter estimation
following a detection.
For demonstrative purposes, in this section, we will
restrict ourselves to the simplest phasing relations con-
sistent with radiation-reaction inspiral as given by Eqs.
(17-18), while using higher-order PN expressions than
Eqs. (19 -20) for relating the orbit to the phase. For ex-
ample, from [55], we have found to second PN (beyond
leading) order:
r12(Ω)
= (GMΩ)−2/3 − (3− ν)
− (18− 81ν − 8ν
(GMΩ)2/3, (32)
= (GMΩ)1/3 +
(15− ν)
(GMΩ)
(441− 324ν − ν2)
(GMΩ)5/3, (33)
and we note that higher-order equivalents of these can
be computed from [56].
In the numerical construction of initial data, the pri-
mary input is the coordinate separation of the holes. In
placing the punctures on the numerical grid, the separa-
tion must be maintained exactly. To ensure this, we in-
vert Eq. (32) to obtain the exact Ωr corresponding to our
desired r12. Then we use Eq. (18) with t = 0 to find the
coalescence time tc that yields this Ωr. Once we have ob-
tained tc, we then find the orbital phase Φ and frequency
Ω at any source time τ directly from Eqs. (17-18), and
the corresponding separation r12 and momentum p from
Eqs. (32-33), or their higher-order equivalents.
In Fig. 3, we show a representative component of the
retarded-time part of hTTij for both circular and leading-
order inspiral orbits. For both orbits, we use the ex-
tended Keplerian relations (32) and (33); otherwise the
orbital configuration is that of Fig. 1. The coalescence
time is now tc ∼ 1100M . We can see that the cumu-
lative wavelength error of the circular-orbit assumption
becomes very large at large distances from the sources.
This demonstrates that using inspiral orbits instead of
circular orbits will significantly enhance the phase accu-
racy of the initial data, even though circular orbits are
in principle sufficient when we include terms only up to
O(v/c)4 as done in this work. From now on we use only
inspiral orbits.
0 500 1000 1500
circular
inspiral
FIG. 3: The xx component of the full hTTij for a binary with
initial separation r12 = 10M in a circular (black/solid) or
inspiralling (red/dashed) orbit. Both fields have been rescaled
by the observer radius r = z to compensate for the leading
1/r fall-off. The orbital configuration is the same as for Fig.
1, apart from the Keplerian relations, where we have used the
higher-order relations (32-33), yielding tc ∼ 1100M . Note the
frequency broadening at more distant field points.
0 100 200 300 400 500 600 700 800
-0.0015
-0.001
-0.0005
0.0005
0.001
 (quadrupole)
 (full)
 (quadrupole)
 (full)
FIG. 4: Plus and cross polarizations of the strain observed at
field point r = 100M , θ = π/4, φ = 0. Both the quadrupole-
approximation waveform (black/solid and green/dot-dashed)
and the full (red/dashed and blue/dotted) waveforms coming
from hTTij are shown. The orbital configuration is the same as
for Fig. 1.
Next, we compare our full waveform hTTij (expressed
as the combinations h+ and h×) at an intermediate-field
position (r = 100M , θ = π/4, φ = 0) to the lowest-order
quadrupole result. In Fig. 4, the orbital configuration
is the same as for Fig. 1. As one can see, both the
+ and × polarizations of our hTTij agree very well with
quadrupole results, as they should. We demonstrate the
near- and intermediate-zone behavior of the new data on
50 100 150 200 250 300 350 400 450 500
-0.05
-0.05
quadrupole
FIG. 5: Plus and cross polarizations of the strain observed at
t = 0 along the z axis. We show the near-zone (solid/black),
the quadrupole (dashed/red) and full (dot-dashed/green)
waveforms. All waveforms have been rescaled by the observer
radius r = z to compensate for the leading 1/r fall-off. The
orbital configuration is the same as for Fig. 1.
the initial spatial slice in Fig. 5. The quadrupole and
full solutions agree very well outside ∼ 100M . However,
the full solution’s phase and amplitude approach the NZ
solution closer to the sources.
B. Numerical Implementation
After having confirmed that we have a PN three-metric
gij that is accurate up to errors of order O(v/c)
5, and
that correctly approaches the quadrupole limit outside
the near zone, we are now ready to construct initial data
for numerical evolutions. In order to do so, we need the
intrinsic curvature Kij , which can be computed as in
Tichy et al. [45] from the conjugate momentum. The
difference is that here we use the full ḣTTij instead of the
near-zone approximation ḣ
TT (4)
ij to obtain the conjugate
momentum [43]. The result is
Kij = −ψ−10PN
ḣTTij + (φ(2)π̃
+O(v/c)6, (34)
where the error term comes from neglecting terms like
ij,(div)
at O(v/c)5 in hTTij , and where ψPN , π̃
φ(2) can be found in Tichy et al. [45]. An additional
difference is that the time derivative of hTTij is evaluated
numerically in this work. Note that the results for gij
are accurate up to O(v/c)4, while the results for Kij are
0 0.5 1 1.5
 = 10M
 = 20M
 = 50M
 = 100M
0 0.5 1 1.5
FIG. 6: Upper panel: Hamiltonian constraint violation along
the y axis of our new data in the near zone, as a function of
binary separation r12. Lower panel: Momentum constraint
(y-component) violation of the same data along the x axis.
The orbital configuration is that of Fig. 3. Distances have
been scaled relative to r12, so that the punctures are initially
at y/r12 = ±0.5.
accurate up O(v/c)5, because Kij contains an additional
time derivative [45, 57, 58].
Next we show the violations of the Hamiltonian and
momentum constraints computed from gij and Kij , as
functions of the binary separation r12. As we can see in
both panels of Fig. 6, the constraints become smaller for
larger separations, because the post-Newtonian approxi-
mation gets better. Note that, as in [45], the constraint
violation remains finite everywhere, and is largest near
each black hole.
C. Curvature Invariants and Asymptotic Flatness
In analysis of both initial and evolved data, it is often
instructive to investigate the behavior of scalar curva-
ture invariants, as these give some idea of the far-field
properties of our solution. We expect, for an asymptoti-
cally flat space-time, that in the far field, the speciality
index S ≡ 27J 2/I3 will be close to unity. This can
be seen from the following arguments. Let us choose a
tetrad such that the Weyl tensor components ψ1 and ψ3
are both zero. Further, we assume that in the far field
ψ0 and ψ4 are both perturbations of order ǫ off a Kerr
background. Then
S ≈ 1− 3ψ0ψ4
+O(ǫ3), (35)
which is indeed close to one. Note however, that this
argument only works if the components of the Weyl ten-
sor obey the peeling theorem, such that ψ2 ∼ O(r−3),
ψ0 ∼ O(r−5) and ψ4 ∼ O(r−1). In particular, if ψ0 falls
off more slowly than O(r−5), S will grow for large r.
Now observe that ψ0 ∼ O(r−5) ∼ M3/r5 is formally of
O(v/c)6. Thus, in order to see the expected behavior of
S ≈ 1 in the far-field we need to go to O(v/c)6. If we only
go to O(v/c)4 (as done in this work) ψ0 consists of un-
controlled remainders only, which should in principle be
dropped. When we numerically compute S we find that
for our data, S deviates further and further from unity
for large distances from the binary. This reflects the fact
that the so-called “incoming” Weyl scalar ψ0 only falls
off as 1/r3, due to uncontrolled remainders at O(v/c)6,
which arise from a mixing of the background with the
TT waveform.
V. DISCUSSION AND FUTURE WORK
Exploring and validating PN inspiral waveforms is cru-
cially important for gravitational-wave detection and for
our theoretical understanding of black-hole binaries. Our
goal has been to provide a step forward in this under-
standing by building a direct interface between the PN
approach and numerical evolution, along the lines ini-
tially outlined in Ref. [45]. In this paper we have
essentially completed the calculation of the transverse-
traceless part of the ADM-TT metric to O(v/c)4 pro-
vided in [45], yielding data that, on the initial Cauchy
slice, will describe the space-time into the far-field.
We have incorporated this formulation into a numerical
initial-data routine adapted to the “puncture” topology
that has been so successful recently, and have explored
these data’s numerical properties on the initial slice.
Our next step is to evolve these data with moving
punctures, and investigate how the explicit incorporation
of post-Newtonian waveforms in the initial data affects
both the ensuing slow binary inspiral of the sources and
the release of radiation from the system. We note es-
pecially that our data are non-conformally flat beyond
O(v/c)3. We expect our data to incorporate smaller
unphysical initial distortions in the black holes than is
possible with conformal flatness, and hence less spurious
gravitational radiation during the numerical evolution.
We see this as a very positive step toward providing fur-
ther validation of numerical relativity results for multiple
orbit simulations, since it permits comparison with PN
results where they are expected to be reliable. Our initial
data will also allow us to fully evaluate the validity of PN
results for merging binaries by enabling comparison with
the most accurate numerical relativity results.
We expect that further development of these data will
certainly involve the use of more accurate orbital phas-
ing information than the leading order given by Eqs. (17-
18). This information is available in radiative coordinates
(see, e.g. Eq. (6.29) of [59]) appropriate for far-field eval-
uation of the gravitational radiative modes; it may be
possible to produce them in ADM-TT coordinates via
a contact transformation, or by direct calculation (see,
e.g. [60]). For initial separations similar to the fiducial
test case of this paper, r12 = 10M , the order necessary
for clean matching of the initial wave content with the
new radiation generated in evolution should not be par-
ticularly high [26]. As noted, the Keplerian relations
Eqs. (32-33) can easily be extended to higher PN order.
The data presented already allow for arbitrary initial
mass ratios ν; this introduces the possibility of significant
gravitational radiation in odd-l multipoles, together with
associated phenomena, such as in-plane recoil “kicks”.
An interesting future development of these data will be
the inclusion of spin angular momenta on the pre-merger
holes. This will open our initial-data prescription to de-
scribing an even richer spectrum of binary radiation.
Acknowledgments
We would like to thank L. Blanchet and G. Schäfer for
generous assistance and helpful discussion.
M.C., B.K. and B.W. gratefully acknowledge the sup-
port of the NASA Center for Gravitational Wave Astron-
omy (NAG5-13396). M.C. and B.K. also acknowledge the
NSF for financial support under grants PHY-0354867 and
PHY-0722315. B.K. also acknowledges support from the
NASA Postdoctoral Program at the Oak Ridge Associ-
ated Universities. The work of W.T. was supported by
NSF grant PHY-0555644. W.T. also acknowledges par-
tial support from the NCSA under Grant PHY-060040T.
The work of B.W. was also supported by NSF grants
PHY-0245024 and PHY-0555484.
APPENDIX A: DETAILS OF INTEGRAL
CALCULATION
Here we present some more details of the calculations
that lead to the three contributions to Eq. (23): Eqs.
(27-29). Inserting Eq. (24) in the general integral (23),
we can write H
TTA[~u] as a combination of scalar and
tensor terms:
HTTAij [~u] = 16πG
ui uj −
Iij A +
[uc ud
IcdA δij
2 uc u(i
Icj)A +
[uc ud
I cdij A
,(A1)
where the “I” integrals are defined as:
d3~k dω
(2 π)4
(ω/k)2 ei k rA cos θ−i ω T
k2 − (ω + i ǫ)2
, (A2)
d3~k dω
(2 π)4
ki kj
× (ω/k)
2 ei k rA cos θ−i ω T
k2 − (ω + i ǫ)2
, (A3)
i j c d
d3~k dω
(2 π)4
ki kj kc kd
× (ω/k)
2 ei k rA cos θ−i ω T
k2 − (ω + i ǫ)2
. (A4)
Here T ≡ t − τ , and ~rA ≡ ~x − ~xA. We have also taken
our integration coordinates such that ~rA lies in the z di-
rection, so that the dummy momentum vector ~k satisfies
~k · ~rA = k rA cos θ, (A5)
d3~k = k2 dk sin θ dθ dφ. (A6)
Define the unit orthogonal vectors n̂A ≡ (0, 0, 1) , ℓ̂ ≡
(cosφ, sinφ, 0). Then we can write
~k = k cos θ n̂A + k sin θ ℓ̂ ⇒ ~k · ~rA = rA ~k · n̂A.
We can also define a projector tensor onto ℓ̂:
Qa b ≡ δa b − na nb ⇒ Qab = δab − na nb
⇒ QacQcb = Qab , Qab nb = 0 , Qab ℓb = ℓa.
1. Angular integration
We will neglect the A subscript for now, until it be-
comes relevant again. To calculate the integrals (A2-A4),
we begin with the φ integration. The only φ dependence
comes from the ~ℓ parts of the ~k terms. It can be seen
from elementary trigonometric integrals that:
dφ ℓa =
dφ ℓa ℓb ℓc = 0,
dφ ℓa ℓb = π Qa b,
dφ ℓa ℓb ℓc ℓd =
Qa bQc d +Qa cQb d +Qa dQb c
We use these to calculate the φ integrals for Ia bA and
Ia b c dA . Define w ≡ cos θ. Then
dφ 1 = 2 π,
ka kb
= 2 π w2 na nb + π (1− w2)Qa b,
ka kb kc kd
= 2 π w4 na nb nc nd
+6 πw2 (1− w2)Q(a b nc nd)
(1− w2)2Q(a bQc d).
So the next integrals will differ in their θ dependence,
contained in the powers of w above. The θ integrals will
contain the following basic types:
g0(a) ≡
dw eaw = 2
sinh a
, (A7)
g2(a) ≡
dw w2 eaw = 2
sinh a
− 4 cosha
sinh a
g4(a) ≡
dw w4 eaw = 2
sinh a
− 8 cosha
sinh a
− 48 cosha
sinh a
. (A9)
Now Ia b and Ia b c d can be written as the linear combi-
nations:
Ia b =
Qa b I
(na nb − 1
Qa b)K
, (A10)
Ia b c d =
na nb nc nd − 3Q(a b nc nd) + 3
Q(a bQc d)
3Q(a b nc nd) −
Q(a bQc d)
Q(a bQc d) I
. (A11)
I here can be expressed in terms of g0(a) above:
(2 π)3
(ω/k)2
k2 − (ω + i ǫ)2 e
i k r cos θ−i ω T
(2 π)3
ω2 e−i ω T
k2 − (ω + i ǫ)2
g0(i k r)
(2 π)3
ω2 e−i ω T J0. (A12)
The 1/2 factor is because we moved to integrating k over
the whole real line instead of the positive half-line (this
is permissible as gn(a) is an even function of a). K and L
are defined analogously to I, but with extra even powers
of cos θ = w:
(2 π)3
(ω/k)2
k2 − (ω + i ǫ)2
ei k r cos θ−i ω T cos2 θ
(2 π)3
ω2 e−i ω T
k2 − (ω + i ǫ)2 g2(i k r)
(2 π)3
ω2 e−i ω T J2, (A13)
(2 π)3
(ω/k)2
k2 − (ω + i ǫ)2
ei k r cos θ−i ω T cos4 θ
(2 π)3
ω2 e−i ω T
k2 − (ω + i ǫ)2
g4(i k r)
(2 π)3
ω2 e−i ω T J4. (A14)
2. Momentum integration
Now we address the k integrals, defined as:
dk fn(k) =
dk f+n (k) +
dk f−n (k),
where we collect the positive exponents in the gn in the
integrand of f+n (k), and the negative exponents in f
n (k):
f+n (k) ≡
g+n (i k r)/2
k2 − (ω + i ǫ)2
, f−n (k) ≡
g−n (i k r)/2
k2 − (ω + i ǫ)2
We calculate this as the sum of contour integrals of the
“plus” and “minus” integrands (necessary, as the oppo-
site signs require different contours). Each of these has
poles at k = 0, k = k+ ≡ ω + i ǫ, and k = k− ≡ −ω − i ǫ
(the first of these is from the gn). We integrate the “plus”
integrands anticlockwise around the contour C1, and the
“minus” integrands anticlockwise around the contour C2
(see Fig. 7); taking the limit |k| → ∞, the contribu-
tion from the curved segments vanishes, and the residue
theorem gives us:
Jn = 2 π iRes[f
n , k+]− 2 π iRes[f−n , k−]
+π iRes[f+n , 0]− π iRes[f−n , 0]. (A15)
Calculating the residues, we find the values of each of
the Jn:
π ei r (ω+i ǫ)
r (ω + i ǫ)2
r (ω + i ǫ)2
, (A16)
π ei r (ω+i ǫ)
r (ω + i ǫ)2
π ei r(ω+i ǫ) [−2 + 2 i r (ω + i ǫ)]
r3 (ω + i ǫ)4
r3 (ω + i ǫ)4
, (A17)
π ei r (ω+i ǫ)
r (ω + i ǫ)2
4π ei r(ω+i ǫ)
r5 (ω + i ǫ)6
[6− 6 i r (ω + i ǫ)
−3 r2 (ω + i ǫ)2 + i r3 (ω + i ǫ)3
− 24 π
r5 (ω + i ǫ)6
. (A18)
3. Frequency integration
Now we perform the ω integration. Inserting the re-
sults (A16-A18) into (A12-A14) respectively, we see that
each of I, K and L contains a delta function, which we
can extract:
4 π r
[δ(T − r) − δ(T )],
4 π r
δ(T − r) + e−r ǫ
(2 π)3
e−iω (T−r) F2a(ω)
(2 π)3
e−i ω T F2b(ω),
4 π r
δ(T − r) + e−r ǫ
(2 π)3
e−iω (T−r) F4a(ω)
(2 π)3
e−i ω T F4b(ω),
where the new terms on the right-hand side come from
the Jn above, grouped by exponential, as that is what
determines the contours chosen during integration (see
Fig. 7):
F2a(ω) =
π ω2 [−2 + 2 i r (ω + i ǫ)]
r3 (ω + i ǫ)4
F2b(ω) =
2 π ω2
r3 (ω + i ǫ)4
F4a(ω) =
r5 (ω + i ǫ)6
[24− 24 i r (ω + i ǫ)
−12 r2 (ω + i ǫ)2 + 4 i r3 (ω + i ǫ)3
F4b(ω) = −
24 π ω2
r5 (ω + i ǫ)6
Now the residues are as follows (taking the ǫ→ 0 limit):
e−iω (T−r) F2a(ω),−i ǫ
2 π i T
e−iω T F2b(ω),−i ǫ
= −2 π i T
e−iω (T−r) F4a(ω),−i ǫ
4 π i T 3
e−iω T F4b(ω),−i ǫ
= −4 π i T
The only pole is at ω = −i ǫ, so if we can close the contour
in the upper half-plane, we will get zero.
• For T < 0, both the “a” and “b” integrals can be
closed in C1. Result: zero contribution.
• For 0 < T < r, the “a” integrals can be closed
in C1, but the “b” integrals must be closed in C2.
Result: “b” contribution.
• For T > r, both the “a” and “b” integrals must be
closed in C2. But then the “a” and “b” residues
cancel out. Result: zero contribution.
Thus the only interesting contribution happens in the
interval 0 < T < r ⇔ t− r(τ) < τ < t. In this case, the
final integrals yield
(2 π)3
e−iω (T−r) F2b(ω) = −
2 π r3
(2 π)3
e−iω (T−r) F4b(ω) = −
−ω − iǫ
ω + iǫ
FIG. 7: Contours needed to complete integration over k (left) and ω (right).
leading to the final result for K and L:
4 π r
δ(T − r)− 1
4 π r
δ(T ),
4 π r
δ(T − r)−Θ(T )Θ(r − T ) T
2 π r3
4 π r
δ(T − r)−Θ(T )Θ(r − T ) T
We use these to calculate the Ii j and Ii j k l:
Ii j =
ni nj
4 π r
δ(T − r)−Θ(T )Θ(r − T ) T
2 π r3
4 π r
δ(T ) + Θ(T )Θ(r − T ) T
2 π r3
,(A19)
Ii j k l =
ni nj nk nl
4 π r
δ(T − r) −Θ(T )Θ(r − T ) T
− 3Q(i j nk nl) Θ(T )Θ(r − T )
2 π r3
Q(i j Qk l)
4 π r
δ(T ) + Θ(T )Θ(r − T )
. (A20)
4. Time integration
The final integrations will be over the source time τ . The “crossing times” for the two Θ functions are τ = t and
τ = tr, where t is the present field time, and tr the corresponding retarded time defined by (25). Now taking a general
function y(τ), we find that
dτ IA y(τ) =
y(trA)
4 π rA(t
− y(t)
4 π rA(t)
A y(τ) =
niA n
4 π rA
4 π rA
3niA n
A − δ
) (t− τ) y(τ)
4 π rA(τ)3
i j k l
A y(τ) =
niA n
4 π rA
4 π rA
−3Q(i jA n
(t− τ)
2 π rA(τ)3
−niA n
A + 3Q
(t− τ)3
π rA(τ)5
y(τ).
These can now be substituted into the general integral (A1). We write the result as a sum of terms at the present
field-point time t, the retarded time trA, and interval terms between them,
TTA[~u] = H
TTA[~u; t] +H
TTA[~u; t
A] +H
TTA[~u; t
A → t],
TTA[~u; t] = −
rA(t)
ui uj − u
[uk ul
Qk lA δ
2 uk u
[uk ul
TTA[~u; t
ui uj − u
niA n
[uk ul
nkA n
2 uk u
[uk ul
niA n
TTA[~u; t
A → t] = −4G
(t− τ)
rA(τ)3
3niA n
A − δ
[uk ul
3nkA n
A − δk l
2 uk u
A − δj) k
[uk ul
(t− τ)3
rA(τ)5
[uk ul
niA n
A − 3Q
[1] R. Vogt, in Sixth Marcel Grossman Meeting on General
Relativity (Proceedings, Kyoto, Japan, 1991), edited by
H. Sato and T. Nakamura (World Scientific, Singapore,
1992), pp. 244–266.
[2] B. Abbott et al. (LIGO Scientific), Nucl. Instrum. Meth.
A517, 154 (2004), gr-qc/0308043.
[3] P. Bender et al., Tech. Rep. MPQ 233, Max-
Planck-Institut für Quantenoptik (1998), URL:
http://www.lisa-science.org/resources/talks-articles/mission/prephasea.pdf.
[4] K. Danzmann and A. Rudiger, Class. Quantum Grav.
20, S1 (2003).
[5] A. Buonanno, G. B. Cook, and F. Pretorius, Phys. Rev.
D 75, 124018 (2007), gr-qc/0610122.
[6] E. Berti, V. Cardoso, J. A. Gonzalez, U. Sperhake,
M. Hannam, S. Husa, and B. Brügmann (2007), arXiv:gr-
qc/0703053.
[7] F. Pretorius, Phys. Rev. Lett. 95, 121101 (2005), gr-
qc/0507014.
[8] M. Campanelli, C. O. Lousto, P. Marronetti, and
Y. Zlochower, Phys. Rev. Lett. 96, 111101 (2006), gr-
qc/0511048.
[9] J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz, and
J. van Meter, Phys. Rev. Lett. 96, 111102 (2006), gr-
qc/0511103.
[10] S. Brandt and B. Brügmann, Phys. Rev. Lett. 78, 3606
(1997), gr-qc/9703066.
[11] M. Shibata and T. Nakamura, Phys. Rev. D 52, 5428
(1995).
[12] T. Baumgarte and S. Shapiro, Phys. Rev. D 59, 024007
(1999), gr-qc/9810065.
[13] B. Brügmann, W. Tichy, and N. Jansen, Phys. Rev. Lett.
92, 211101 (2004), gr-qc/0312112.
[14] M. Campanelli, C. O. Lousto, and Y. Zlochower, Phys.
Rev. D 73, 061501(R) (2006), gr-qc/0601091.
[15] F. Pretorius, Class. Quantum Grav. 23, S529 (2006), gr-
qc/0602115.
[16] J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz, and
J. van Meter, Phys. Rev. D 73, 104002 (2006), gr-
qc/0602026.
[17] B. Brügmann et al. (2006), gr-qc/0610128.
[18] M. A. Scheel et al., Phys. Rev. D 74, 104006 (2006),
gr-qc/0607056.
[19] P. Marronetti, W. Tichy, B. Brügmann, J. Gonzalez,
M. Hannam, S. Husa, and U. Sperhake, Class. Quant.
Grav. 24, S43 (2007), gr-qc/0701123.
[20] W. Tichy, Phys. Rev. D 74, 084005 (2006), gr-
qc/0609087.
[21] H. P. Pfeiffer, D. A. Brown, L. E. Kidder, L. Lindblom,
G. Lovelace, and M. A. Scheel, Class. Quant. Grav. pp.
S59–S81 (2007), gr-qc/0702106.
[22] J. G. Baker, M. Campanelli, F. Pretorius, and Y. Zlo-
chower, Class. Quant. Grav. 24, S25 (2007), gr-
qc/0701016.
[23] J. Thornburg, P. Diener, D. Pollney, L. Rezzolla,
http://www.lisa-science.org/resources/talks-articles/mission/prephasea.pdf
E. Schnetter, E. Seidel, and R. Takahashi, Class. Quant.
Grav. 24, 3911 (2007), gr-qc/0701038.
[24] NRwaves home page:
https://gravity.psu.edu/wiki NRwaves.
[25] J. G. Baker, S. T. McWilliams, J. R. van Meter, J. Cen-
trella, D.-I. Choi, B. J. Kelly, and M. Koppitz, Phys.
Rev. D 75, 124024 (2007), gr-qc/0612117.
[26] J. G. Baker, J. R. van Meter, S. T. McWilliams, J. Cen-
trella, and B. J. Kelly (2006), gr-qc/0612024.
[27] M. Campanelli, Class. Quant. Grav. 22, S387 (2005),
astro-ph/0411744.
[28] F. Herrmann, D. Shoemaker, and P. Laguna (2006), gr-
qc/0601026.
[29] J. G. Baker et al., Astrophys. J. 653, L93 (2006), astro-
ph/0603204.
[30] M. Campanelli, C. O. Lousto, and Y. Zlochower, Phys.
Rev. D 74, 041501(R) (2006), gr-qc/0604012.
[31] M. Campanelli, C. O. Lousto, and Y. Zlochower, Phys.
Rev. D 74, 084023 (2006), astro-ph/0608275.
[32] M. Campanelli, C. O. Lousto, Y. Zlochower, B. Krishnan,
and D. Merritt, Phys. Rev. D 75, 064030 (2007), gr-
qc/0612076.
[33] J. A. Gonzalez, U. Sperhake, B. Bruegmann, M. Han-
nam, and S. Husa, Phys. Rev. Lett. 98, 091101 (2007),
gr-qc/0610154.
[34] F. Herrmann, I. Hinder, D. Shoemaker, P. Laguna, and
R. A. Matzner (2007), gr-qc/0701143.
[35] M. Campanelli, C. O. Lousto, Y. Zlochower, and D. Mer-
ritt, 659, L5 (2007), revised version has very different
numbers/formulae, gr-qc/0701164.
[36] M. Koppitz, D. Pollney, C. Reisswig, L. Rezzolla,
J. thornburg, P. Diener, and E. Schnetter, Phys. Rev.
Lett. 99, 041102 (2007), gr-qc/0701163.
[37] J. A. Gonzalez, M. D. Hannam, U. Sperhake,
B. Brügmann, and S. Husa, Phys. Rev. Lett. 98, 231101
(2007), gr-qc/0702052.
[38] D.-I. Choi et al. (2007), gr-qc/0702016.
[39] J. G. Baker et al. (2007), astro-ph/0702390.
[40] F. Pretorius and D. Khurana, Class. Quant. Grav. 24,
S83 (2007), gr-qc/0702084.
[41] M. Campanelli, C. O. Lousto, Y. Zlochower, and D. Mer-
ritt, Phys. Rev. Lett. 98, 231102 (2007), arXiv:gr-
qc/0702133.
[42] W. Tichy and P. Marronetti (2007), gr-qc/0703075.
[43] G. Schäfer, Ann. Phys. 161, 81 (1985).
[44] P. Jaranowski and G. Schäfer, Phys. Rev. D 57, 7274
(1998), errata: Phys. Rev. D 63, 029902(E) (2000), gr-
qc/9712075.
[45] W. Tichy, B. Brügmann, M. Campanelli, and P. Diener,
Phys. Rev. D 67, 064008 (2003), gr-qc/0207011.
[46] Cactus Compuational Toolkit,
http://www.cactuscode.org.
[47] T. Ohta, H. Okamura, T. Kimura, and K. Hiida, Prog.
Theor. Phys. 51, 1598 (1974).
[48] V. A. Fock, The Theory of Space, Time and Gravitation,
2nd ed. (Pergamon Press, 1964).
[49] W. Tichy, B. Brügmann, and P. Laguna, Phys. Rev. D
68, 064008 (2003), gr-qc/0306020.
[50] W. Tichy and B. Brügmann, Phys. Rev. D 69, 024006
(2004), gr-qc/0307027.
[51] M. Ansorg, B. Brügmann, and W. Tichy, Phys. Rev. D
70, 064011 (2004), gr-qc/0404056.
[52] L. S. Finn and D. F. Chernoff, Phys. Rev. D 47, 2198
(1993), gr-qc/9301003.
[53] C. Cutler et al., Phys. Rev. Lett. 70, 2984 (1993), astro-
ph/9208005.
[54] W. Tichy, E. E. Flanagan, and E. Poisson, Phys. Rev. D
61, 104015 (2000), gr-qc/9912075.
[55] G. Schäfer and N. Wex, Phys. Lett. A 174, 196 (1993).
[56] R.-M. Memmesheimer, A. Gopakumar, and G. Schäfer,
Phys. Rev. D 70, 104011 (2004), gr-qc/0407049.
[57] N. Yunes, W. Tichy, B. J. Owen, and B. Brügmann,
Phys. Rev. D 74, 104011 (2006), gr-qc/0503011.
[58] N. Yunes and W. Tichy, Phys. Rev. D 74, 064013 (2006),
gr-qc/0601046.
[59] L. Blanchet, Phys. Rev. D 54, 1417 (1996), gr-
qc/9603048.
[60] T. Damour, A. Gopakumar, and B. R. Iyer, Phys. Rev.
D 70, 064028 (2005), gr-qc/0404128.
http://www.cactuscode.org
